texts: Can we limit the search range of geocode() function in PostGIS Tiger Geocoder?

texts:
I found the server with only 2 states data loaded is much faster than the server with all states loaded. My theory is bad formatted address that don't have a exact hit at first will cost much more time when the geocoder checked all states. With only 2 states this search is limited and stopped much early.

There is a restrict_region parameter in geocode function looks promising if it can limit the search range, since I have enough information or reason to believe the state information in my addresses input are correct. I wrote a query trying to use one state's geometry as the limiting parameter:

SELECT geocode('501 Fairmount DR , Annapolis, MD 20137', 1, the_geom)
    FROM tiger.state WHERE statefp = '24';

and compared the performance with the simple version

SELECT geocode('501 Fairmount DR , Annapolis, MD 20137',1);

EDIT Note the city and zipcode in this input is wrong which is intentional. If the address is perfect there is no point to limit the search range since the geocoder can go to the right table at first try. Only this attempt is meaningful only when the zipcode or city is wrong.

I didn't find performance gain with the parameter. Instead it lost the performance gain from caching, which usually came from running same query immediately again because all the needed data have been cached in RAM.

Maybe my usage is not proper, or this parameter is not intended to work as I expected. However if the search range can be limited, the performance gain could be substantial, since it's the bad formatted addresses took the most time to geocode, and they also often mess up the already cached data because the geocoder need to search for states, even all my input are in one state and all data can be cached in RAM.


texts: CartoDB - Method to show/hide layers for hatching overlap for > 4 layers?

texts:
Our goal is to show/hide individual hatched maps.
Unfortunately we have more than the 4 max layers allowed in a viz.
We successfully made and example where we show/hide the 4 layers in a viz:
http://maps.massgis.state.ma.us/cartodb_test/cartodb_demo_1_viz_show_hide_animals_chkbox_legend.html

Next, we tried to have an example where we have 2 viz, each with 4 layers and try to show/hide all 8 items on the map.  They all draw on the map at first and clear, but then the show/hide doesn't properly work. We tried two methods:

1) Try to put the second viz into a different objects called layers2
http://maps.massgis.state.ma.us/cartodb_test/cartodb_2_viz_show_hide_animals_chkbox_method1.html
2) Try to put the second viz into the same layers object
http://maps.massgis.state.ma.us/cartodb_test/cartodb_2_viz_show_hide_animals_chkbox_method2.html

Is there a way to combine > 1 viz with > 1 layer each and do show/hide.  Alternatively, we could have 1 layer per viz and show/hide each viz somehow?


texts: How to load vector files from file path in QGIS?

texts: We have around 500 files in different folders. For QC purpose, we've selected 25 files(we have the list of files/path in excel) out of those 500. Is it possible to open all those 25 files(it might be a shp or tab) in one shot using the file path in excel sheet ?


texts: How can I overlay two raster without turning them in vector format in Python to get zone wise statistics?

texts:
I have two raster re-sampled and same extent
one is elevation class raster from DEM reclassification
classified in 10 zone with 500m interval
and another is snow cover fraction raster has snow fraction value 0-100 and some reassigned classes like 201 for cloud 205 for water etc.

I want to have a table of elevation class-wise snow cover fraction. How to calculate that in Python? How to sort the array and make query?

I do not want to convert the input raster format to vector format because in the later stage I need the output file in raster format actually the array to read in to build contingency matrix.

I have tried with from rasterstats import zonal_stats but it needs to have one file in vector format however I have two raster data.


texts: Get all residential addresses within bounding box

texts:
How can I get all the residential addresses within a bounding box?  I tried the query from this question (reproduced here):

way
  ["building"~"residential|house|terrace|detached|apartments"]
  ({{bbox}});
(._;>;);
out;

The bounding box is whatever that is visible in the map area.

Whatever map I try this out on in Overpass Turbo, nothing is returned.


texts: What might cause OpenLayers plugin to hang on windows 8.1?

texts:
I'm trying to add OpenStreetMap layer to the map through OpenLayers plugin and then application hangs.
The problem occures on Windows 8.1 QGIS 2.0.1. On Windows 7 everything works correct.
Does anyone know what might be problem and how to fix it or how to add OSM another way?


texts:
I am trying to build a model in arcmap that (1) converts a raster to polygon (2) dissolves the polygons by landcover class (3) adds a field that calculates area of the landcover polygons. I want to combine all the polygons of different landcover classes into a single record for each class so I have the total area for the class instead of multiple polygons for each class with their own area.
I know how to manually merge records in Editor, but how do I build it into my model?


texts: Learning about ArcMap's Tracking Analyst tool?

texts:
I was examining the ArcMap toolboxes to learn more about their functions and found the Tracking Analyst tool.

Does anybody have examples of how can this function is used and when best to use it?


texts: osm2po is missing edges while importing to pgRouting

texts:
I am using osm2po for importing pbf file to PostgreSQL but when I find route using built-in algorithms I found that some of the edges are missing so route does not match with road.

Does anybody know why this is happening?


texts: Seeking shapefile for Australia and New Zealand?

texts:
I was trying to find a shapefile which has Australia and New Zealand together (states) so that I could convert to jVectorMap for ANZ.

Where I can find it?


texts:
I have been learning about the TIGER and the National Map datasets for mapping roads in the United States.

One limitation is that as far as I can see there is no data on the number of lanes. The main available attribute seems to be "road class" which is things like "Freeway", "Collector", "Local Road" etc, which does not really indicate the number of lanes the road has. For my purpose I am mostly interested in knowing the number of lanes and whether it is a divided highway or not (ie has a central barrier), or if it has a median, like most interstates do. Is there any way to get this information?


texts: Efficient preallocation of simple features in R (sf r)

texts:
I have a use case that requires looping over a large number of API calls to populate a simple features dataset. Read/write is very slow when I preallocate a large empty collection initially, which goes against what I though the whole point of preallocating was for.
Simple features are likely different from a memory standpoint, as ex ante one doesn't know if the geometry for each feature will be "small" or "large". Still, it would be great if there was a way to speed this up.
Minimal reprex:
rm(list = ls())
library(sf)
library(tidyverse)

## Initialize size to 100 rows, populate with 100 features ##
start_time100 <- Sys.time()
rsize <- 100

route <- st_sf(id = 1:rsize, geometry = st_sfc(lapply(1:rsize, function(x) st_linestring() )), crs = 4326)

for (i in 1:100) {
  sln <- rbind(c(runif(1),runif(1)), c(runif(1),runif(1)), c(runif(1),runif(1)))
  route_temp <- st_linestring(sln, dim = "XY") %>%
    st_sfc(crs = 4326)
  (route[i,] <- st_sf(id = i, geometry = route_temp))
}

route <- route[!st_is_empty(route),]
end_time100 <- Sys.time()

## Initialize size to 10000 rows, populate with 100 features ##
start_time10k <- Sys.time()
rsize <- 10000

route <- st_sf(id = 1:rsize, geometry = st_sfc(lapply(1:rsize, function(x) st_linestring() )), crs = 4326)

for (i in 1:100) {
  sln <- rbind(c(runif(1),runif(1)), c(runif(1),runif(1)), c(runif(1),runif(1)))

  route_temp <- st_linestring(sln, dim = "XY") %>%
    st_sfc(crs = 4326)
  (route[i,] <- st_sf(id = i, geometry = route_temp))
}

route <- route[!st_is_empty(route),]
end_time10k <- Sys.time()

end_time100 - start_time100
end_time10k - start_time10k

Running on my machine gives:
> end_time100 - start_time100
Time difference of 0.2343311 secs
> end_time10k - start_time10k
Time difference of 9.813453 secs

Modifying the reprex based on the suggestion from @mdsumner is much faster (prob fast enough for my use case), but still sees substantial slowdown for editing within the "larger" object. Replacing the main portion of the code:
route_list <- vector("list", rsize)
route_id <- vector("numeric", rsize)

for (i in 1:100) {
  sln <- rbind(c(runif(1),runif(1)), c(runif(1),runif(1)), c(runif(1),runif(1)))
  route_list[[i]] <- st_linestring(sln, dim = "XY")
  route_id[i] <- i
}

route_list <- route_list %>%
  st_sfc(crs = 4326)
route <- st_sf(id = route_id, geometry = route_list)

Gives run times as below (and I double checked that it's not the route_id component). sfheaders + template sounds useful, though I haven't tried that yet.
> end_time100 - start_time100
Time difference of 0.01904988 secs
> end_time10k - start_time10k
Time difference of 0.09023905 secs


texts: Splitting geometry collection into its items using Python and QGIS API

texts:
I'm getting some geometries via intersecting lines by polygons (somewhat blindly - I don't really know if intersection exists and actually I don't need to know) using python and QGIS API: some_line.intersection(some_polygon).

To my mind I may get several types of geometry:

None (if there is no intersection).
Point
Linestring
Multilinestring
Geometry collections (when there are several points and/or lines)

but the geometry I need in the end is a Linestring (Polyline) or several Linestrings (not MultilineStrings). I can identify geometry type using wkbType(). It returns integer, and this integer corresponds to this order.

How can I split MultilineStrings and Geometry collections and get the parts they are consist of (preferably using QGIS API)?

I wasn't able to find such function in QGIS API, maybe just missed it? There is also GetGeometryRef() function in OGR that I believe could help me, but I can't figure out how to implement it in my case.


texts:
I have measures many dozens of lines and every now and then there is a point wrongly numbered. By sorting on Orthographic Height in the points tab I can easily find these, for example:

4009 : -0.2039 (1 measurement)
3009 : -0.2033 (Many measurements)

Obviously 4009 should be assigned to 3009. In the level-proctab you can reassign a measured tippled by right clicking the point you want to reassign. Right now I have to search within all my measured lines for point 4009 to reassign it in level-proc. I should be able to do this in points or at the very least find out what line it was measured in  automated way.

Same goes for doubles that deviate to much. If the average limit is exceeded Geo-Office automatically discards it that specific point. But in most cases that specific point is incorrectly numbered and need to be assigned to another point ID.

So my question in a nutshell:

How can I assign a point to another ID from the Pointstab?
Or how can I quickly find the line where this point was measured to modify it in the line in Level-Proc?

My current method is rather devious. I am writing down each point that needs to be changed, it are many dozens ID's so I cannot remember them all. Then I'm looking through all my lines in level-Proc to find either the incorrect ID or a ID with a incorrect height (also an incorrect ID essentially). In level proc I can assign ID's to other ID's so I fix it here. I'm really hoping to do this in Points, changing ID's is possible but only if they do not exist and 99 out of 100 times I need to assign it to a existent ID.


texts:
I want to perform IDW interpolation using R, and more specifically, using the idw command from the gstat package. Also, I want to use the optimal number of neighbors (point locations) and distance power, which will be determined by what combination of "idp" and "nmax" produces the lowest RMSE in leave-one-out-cross validation. For this reason, after each run of IDW with different "nmax" and "idp" parameters, I perform the leave-one-out cross validation. In order to do this, I have written a double for loop in R that tests each "idp" parameter for each number of "nmax" neighbors, but my script can't produce the result I want.

power = seq(from = 1, to = 1.6, by = 0.1)
neigh = seq(from = 4, to = 5, by = 1)
# January
for(i in power) {
  for(j in 1:length(neigh)
    jan_idw = idw(Jan ~ 1, st, grd, nmax = power, idp = neigh)
    list_jan_idw = c(list_jan_idw,jan_idw )
    krige.cv(Jan ~ 1, st, nfold = 52)
  }
}

For each run, I want the results to be stored in a data frame or matrix.


texts: Distance to Multiple Neighboring Sectors

texts:
I have two shapefiles - one of neighborhood polygons and another of centroids within those neighborhoods:

For each centroid I would like to calculate:

Distance to the next closest neighborhood (and the neighborhood's identifying attributes)
Distance to the 2nd closest neighborhood (and the neighborhood's identifying attributes)
Distance to the 3rd closest neighborhood (and the neighborhood's identifying attributes)

Is there a simple tool or plugin for QGIS that would allow me to calculate these values for each centroid? The desired output is a table with 3 rows for each centroid and columns of the centroid ID, neighborhood ID, and distance to that neighborhood.


texts: Heat Maps only for feature layer (why not for dynamic layer ?) in ArcGIS Java script API 3.0/3.2

texts:
I am referring this link to generate the heat map. In this example they used feature layer and I tried using dynamic layer but getting an error (that no map is displaying).

I researched on Google regarding this but not able to find such information .

So is there any specific requirement/purpose to use feature layer instead of dynamic layer??

If anyone know this please let me know….

Thanks…

Software platform : ArcGIS 10.1 (Server & Desktop), ArcGIS JSAPI 3.0/3.2


texts: Obtaining the distance from the QGIS Measure Tool?

texts:
I have built a plugin for feeding a Spatialite database with an input mask.
There are two fields which represent the min/max dimension of an object.

In the moment I call the Measure Tool (iface.actionMeasure().trigger()) with the two Push Buttons "measure", click the dimension and type the value with my keyboard into the fields (QLineEdit).

Is there the possibility that after I do the right click (finish measuring) with the Measure Tool the value will be automatically committed to the field (min or max, QLineEdit)?

Call Measure Tool (it works)
Left Click: begin measuring
and 4. at once: Right click and commit Total value of Measure Tool to the min field (QLineEdit).

UPDATE 1:

One day after asking this question I've got the idea to have a linestring layer which I could use as a measure tool (ST_Length(geom)). The advantage is that you could have the full control of manipulating and obtaining data from it. A disadvantage is a additional layer in the project.

A workflow could be:

click measure: toggle editing, activate "Add Feature"
Create linestring that represents the min dimension of an object
Confirm changes, calculate length of linestring    (ST_Length(geom)), commit length to QLineEdit
Close Dialog and remove linestring (DELETE FROM )

But my actual question remains. Is there a way to get the distance value from the Measure Tool directly?

Update 2:

With the help if Jakob I can update my QLineEdit with the Total value of the Measure Tool until I close the Measure Tool. Then the value is set to '0,000 m'. I solved this problem by using this:

def text_changed(str_measure):
    if str_measure == '0,000 m':
        pass
    else:
        self.min.setText(str_measure)

It works, but I don't know if this is a proper solution.


texts: Trim the grid in boundaries of a shapefile

texts:
I have extracted the London Shapefile  from https://data.gov.uk/dataset/6cdebf5d-c69b-4480-8c9c-53ab8a816b9d/statistical-gis-boundary-files-for-london and uploaded London_Borough_Excluding_MHW.shp in QGIS to draw a grid . I am able to draw a grid using the extent layer of the shapefile uploaded . I get the grid mapped to layer as below

How can I go from the first visualisation to the below one , where the grid is fitted into the shapefile below and not lying on top of it

I tried clip but it changes nothing , Any pointers will be helpful as you can find i am a novice user of QGIS . The grid is 100 m * 100m in the 2nd visualisation so its finer than 1st Grid Visualisation


texts: Classifying powerlines in LiDAR using ArcGIS Pro?

texts:
I have a LAS dataset that includes powerlines. The dataset is unclassified, but I have access to a feature class that includes the powerlines as lines.

Is there a tool or method in ArcGIS Pro that lets me classify the points as powerlines and still lets me keep the other points?

I have access to 3D analyst.


texts: Creating trade flow maps in ArcGIS for Desktop (or any free software)?

texts:
I am trying to represent trade flows on a world map where the size and direction of arrows indicate the exports/imports of a given country.

This type of map is commonly found in magazines and newspapers. However, I find it extremely difficult to make.

Is there any function in ArcGIS for Desktop or any free software that would help me in making it?


texts: How to compare 2000 Census block group data with 2010 block group data?

texts:
I want to calculate the percent change for median household income in Los Angeles at the block group level from 2000 to 2010. Someone told me I have to reaggregate 2000 block data to 2010 block group geography, do a spatial join based on centroid, and that I have to download "crosswalk" or relationship data from the U.S. Census comparing 2000 to 2010. I downloaded the relationship .txt file with 2000 California Census tabulation data and the 2010 tabulation data. Right now I have the unjoined data files. I don't know what to do next. What do I join together?

1) Can someone give me step-by-step instructions on how to calculate the percent change at the block group level? :/

2) If calculating it at the block group level is too complicated, I may just change my level of analysis to census tract level. How do I calculate the percent change at the census tract level?


texts: Unable to connect to "Database Connections/.sde" from python

texts:
I just upgraded to 10.3.1 and now I am having issues running my administrative tasks set with Python scripts.

I am unable to use the string "Database Connections/some.sde", i.e:

arcpy.ListUsers("Database Connections/some.sde) Returns the following error:

Traceback (most recent call last):
  File "<pyshell#0>", line 1, in <module>
    userList = arcpy.ListUsers("Database Connections/some.sde")
  File "C:\Program Files\ArcGIS\Server\ArcPy\arcpy\__init__.py", line 1058, in ListUsers
    return utils.listofnamedtuples(gp.listUsers(sde_workspace), 'user')
  File "C:\Program Files\ArcGIS\Server\ArcPy\arcpy\geoprocessing\_base.py", line 504, in <lambda>
    return lambda *args: val(*gp_fixargs(args, True))
TypeError: Could not open SDE workspace.

I can use the full path string for some tasks like ListUsers, AcceptConnections,etc., but not others like reconcile and post, e.g.:

"C:\Users\%username%\AppData\Roaming\ESRI\Desktop10.3\ArcCatalog\some.sde"

Fails without throwing error at: arcpy.ReconcileVersions_management("C:\Users\%username%\AppData\Roaming\ESRI\Desktop10.3\ArcCatalog\some.sde","ALL_VERSIONS","DBO.QC","DBO.Chris;DBO.Joe","LOCK_ACQUIRED","NO_ABORT","BY_OBJECT","FAVOR_EDIT_VERSION","POST","KEEP_VERSION","#")

The full script is here, which worked great prior to 10.3 upgrade: reconcile and post script

In ArcCatalog I have created these database connections, they work fine, I can browse in Catolog, and the .sde files exist in my AppData/Roaming/ESRI subdirectory. I have tried copying all my old connection strings from the 10.2 into my roaming profile. Is there some environmental variable that needs to be repaired? Do I need to re-install Server and Desktop on this machine?

Question: How can I repair my "Database Connections" to run Python scripts


texts: Recently I have updated my postgresql version. I have also upgraded pgAdmin from 3 to 4. But I was unable to load shapefile. I didn't find any plugin called shapefile loader. How can I get the shapefile loader plugin in pgAdmin 4?


texts: Setting opacity info in ArcGis Api for Javascript overwrites featureLayer.renderer.colorInfo

texts:
I'm having trouble setting the color of a feature layer in ArcGIS api for javascript. I have a feature layer, which, on the load event I set its renderer's colorInfo and opacityInfo. Both work individually, but when I try to set a specific color, and set the opacity of the color together, the color reverts to a default color from the feature layer. I've searched the documentation and tried everything--not a lot of examples pertinent to simple stuff like this. It's like colorInfo and opacityInfo conflict with each other.

<script>
require(["stuff"],function(stuff){

map = new Map("div",center, zoom, etc);

var featureLayer = new FeatureLayer("url",{
    id: "id",
    infoTemplate: esri.InfoTemplate({title: "${title}",content:     "${content}"}),
    outFields:["*"],
    visible: true
));

featureLayer.on("load", function(){
featureLayer.renderer.setColorInfo({field: 'my_field', colors: ['#mycolor']});
featureLayer.renderer.setOpacityInfo({opacityValues: [0.25]}); //removes mycolor and sets it to some weird default color
});

});

All I'm really trying to do is define a specific color and set its opacity. I've also tried using a SimpleRenderer, but this seems to conflict with a toggle menu I'm using, so trying to stick with the feature layers own renderer methods.


texts: Python error message using QGIS "load them all" plugin

texts:
I have many thousands of shapefiles organized into folders for residential neighborhoods (state/city/neighborhood).  I installed "load them all" plugin but get a python error message.

  An error has occured while executing Python code:

  Traceback (most recent call last):   File
  "C:/Users/realmarkable1/.qgis2/python/plugins\loadthemall\LoadThemAllDialog.py",
  line 251, in load
      self.apply()   File "C:/Users/realmarkable1/.qgis2/python/plugins\loadthemall\LoadThemAllDialog.py",
  line 193, in apply
      filterList )   File "C:/Users/realmarkable1/.qgis2/python/plugins\loadthemall\LoadFiles.py",
  line 197, in init
      bLayersOff, bDoNotEmpty, bSort, bReverseSort, numLayersToConfirm )   File
  "C:/Users/realmarkable1/.qgis2/python/plugins\loadthemall\LoadFiles.py",
  line 29, in init
      self.tree = Tree( baseDir, False ) # Don't create a parent group   File
  "C:/Users/realmarkable1/.qgis2/python/plugins\loadthemall\LoadFiles.py",
  line 246, in init
      self.root = QgsProject.instance().layerTreeRoot() AttributeError: 'QgsProject' object has no attribute 'layerTreeRoot'

Python version:
2.7.4 (default, Apr  6 2013, 19:54:46) [MSC v.1500 32 bit (Intel)]

QGIS version:
2.2.0-Valmiera Valmiera, c3a2817

Python path:

['C:/PROGRA~2/QGISVA~1/apps/qgis/./python/plugins\\processing', 'C:/PROGRA~2/QGISVA~1/apps/qgis/./python', u'C:/Users/realmarkable1/.qgis2/python', u'C:/Users/realmarkable1/.qgis2/python/plugins', 'C:/PROGRA~2/QGISVA~1/apps/qgis/./python/plugins', 'C:\\PROGRA~2\\QGISVA~1\\bin\\python27.zip', 'C:\\PROGRA~2\\QGISVA~1\\apps\\Python27\\DLLs', 'C:\\PROGRA~2\\QGISVA~1\\apps\\Python27\\lib', 'C:\\PROGRA~2\\QGISVA~1\\apps\\Python27\\lib\\plat-win', 'C:\\PROGRA~2\\QGISVA~1\\apps\\Python27\\lib\\lib-tk', 'C:\\PROGRA~2\\QGISVA~1\\bin', 'C:\\PROGRA~2\\QGISVA~1\\apps\\Python27', 'C:\\PROGRA~2\\QGISVA~1\\apps\\Python27\\lib\\site-packages', 'C:\\PROGRA~2\\QGISVA~1\\apps\\Python27\\lib\\site-packages\\PIL', 'C:\\PROGRA~2\\QGISVA~1\\apps\\Python27\\lib\\site-packages\\win32', 'C:\\PROGRA~2\\QGISVA~1\\apps\\Python27\\lib\\site-packages\\win32\\lib', 'C:\\PROGRA~2\\QGISVA~1\\apps\\Python27\\lib\\site-packages\\Pythonwin', 'C:\\PROGRA~2\\QGISVA~1\\apps\\Python27\\lib\\site-packages\\Shapely-1.2.18-py2.7-win32.egg', 'C:\\PROGRA~2\\QGISVA~1\\apps\\Python27\\lib\\site-packages\\wx-2.8-msw-unicode', 'C:\\PROGRA~2\\QGISVA~1\\apps\\qgis\\python\\plugins\\fTools\\tools', 'C:/Users/realmarkable1/Desktop']


texts: Loading shapefile using ArcPy in Pycharm?

texts:
I want to use ArcPy and pycharm to do zonal statistics. When I loaded shapefile just as I did in the python console in arcgis, the pycharm identifies the shapefile as string type and I got the ERROR 000865: Input value raster: D:/target/New File Geodatabase.gdb\elev does not exist. So how to fix this? Below is my code:

import arcpy
from arcpy import env
from arcpy.sa import *
arcpy.env.workspace = "D:/target/New File Geodatabase.gdb"
featureClass = "D:/target/New File Geodatabase.gdb/target"
inZoneData = "target.shp"   # shapefile
zoneField = "FID"           # the FID field of target shapefile
inValueRaster = arcpy.Raster("elev.img")
arcpy.CheckOutExtension("spatial")
outZonalStatis1 = arcpy.sa.ZonalStatistics(inZoneData,"Id",inValueRaster,"MAXIMUM")


texts: Handle add new feature event and/or access feature before commit?

texts: Is it possible with PyQGIS to catch "add new feature" event before commiting changes on the layer? or at least to access new features before commit.


texts: PostGIS Clipping a Raster with a Geotiff

texts:
I believe I am doing a common operation where I am trying to clip a a Raster with a Geometry.  They both share the same srid and I can successfully perform st_clip and st_intersection with the crop option set to true.

The result is what I expected which a raster where all the values that are outside of the geometry have been given a null value.

What I can't figure out is if I have a rectangular geometry, how I can align the cropping such that the resulting raster is the same shape without any null values?


texts: I am Getting higher stream number for lower order streams in my analysis on ArcGIS,  order 6 has more streams than order 5. is this common in morphometric analysis?


texts: App using latitude and longitude, should decimal degrees be separate a selection from metric?

texts: I am coding a mobile app that displays coordinates in parts of the interface. We currently have an option to toggle between Imperial and Metric. But I am unsure if decimal degrees (vs D/M/S) should be its own option or if it is something that should be tied to the Imperial/Metric display option as well?


texts: Publishing ArcGIS Server service without cache does not show map on mobile?

texts:
We're publishing a whole country data set using ArcGIS Server and we need to display it on web and on mobile applications.

The issue we're facing is that the data is too big so we're not creating cache while publishing the service and here we're facing a strange issue i.e. the service is not showing map on mobile when we try to publish it using mobile SDK's while the service do show map on web when we use ArcGIS JavaScript API.

The question is, how we can make it possible to show maps on mobile if we don't want to create cache while publishing the service?


texts: Saving calculated fields to reuse in other QGIS projects

texts:
I am using QGIS 3.4.

I want to make attribute tables with many calculated fields that link to styles, and reuse the same structure for many projects.

For example, when I draw a pipe (line vector), and enter attributes like diameter and flow, these should calculate a range of derived values that affect how the style is presented on the map.

In a new project when I make a new geopackage layer and apply the "Irrigation Pipe" style, how can I make sure the field calculations are in the table and auto-generated?

I could make a function for each type of field calculation, but would I have to manually load the functions to each field every time I made a new attribute table?

QGIS novice here, it's been 20 years since I did lots of GIS!


texts: PostGIS, ST_Intersects & LIMIT performance question

texts:
I have a pretty standard table with one geometry column and a datetime column. Geom is indexed using gist and table is clustered + VACUUM ANALYZE. The Docker image is kartoza/postgis:12.0.
I consider the following query (bbox intersection):
SELECT product_id
FROM products
WHERE
    ST_Intersects(
        geom,
        ST_MakeEnvelope(%s, %s, %s, %s, 4326)
    )
LIMIT 1000

I have noticed the following behavior: query runs < 10ms when the bbox is below 1 deg² (returning 1000 rows), but the performance degrades as the AOI surface increases (for example, 250ms if the bbox is 20 deg x 20 deg).
I would have expected the opposite as PostGIS is more likely to find quickly 1000 products intersecting a bbox as the bbox is larger.
So I guess the LIMIT clause is not quite effective in my query.
How can I correct this behavior, and implement a robust & fast pagination?


texts: difference in behavior between merge from SAGA and merge from QGIS geoalgorithms

texts:
I'd like to know what I'm doing wrong with the merge algorithms in Python.

I have several shapefiles (polygons) ending in 'LNA', 'LNC', 'LDA', and 'LDC', and I'd like to merge them depending on their names (as part of a bigger pyQgis script). The expected output is 4 files, one for each suffix.
Here is the code:

import processing
from PyQt4 import QtGui
layer = qgis.utils.iface.mapCanvas().layers()
destpath = str((QtGui.QInputDialog.getText(None, "File path", "without shp suffix"))[0])

for typeCarte in ('LNC', 'LNA', 'LDC', 'LDA'):
    layerList = [lyr for lyr in QgsMapLayerRegistry.instance().mapLayers().values() if (lyr.name().endswith(typeCarte) and (lyr.name().startswith('0')))]

    processing.runalg("qgis:mergevectorlayers", layerList, destpath)

    iface.addVectorLayer(
        destpath + ".shp",
        typeCarte + '_temp_Fusion',
        "ogr")

The processing works but only with the files ending in 'LNC', that is the first item in my list. The LDA_temp_Fusion, LNA_temp_Fusion and so forth are just identical.
If I replace now this processing bit with the following line

   processing.runalg('saga:mergelayers', layerList,True,True,destpath)

I get the expected result. What is different in those algorithms that make them behave differently? Or is it just me?


texts: How to enable a Leaflet Draw tool programatically

texts:
I'm using the Leaflet Draw control, and would like to enable/activate a particular tool (Draw Polygon) using JavaScript, so that this tool is active when the page loads. Is this possible?

https://stackoverflow.com/a/17144318/906814 shows how to bypass the control and start drawing directly, but I'd like instead to programatically enable the existing tool since it gives easy access to additional functionality.


texts: VectorSource() as import module - is it possible?

texts:
How to achieve use VectorSource() as import module and pass layerUrl as parameter inside object?
e.g.
new VectorSource({format:  new GeoJSON(), loader: function(), my_layer_name: layerName})

I have many wfs layers and trying simplicity the code.
But have no idea how to pass my own parameter into OpenLayers object VectorSource({}).
Any tips are welcome.
It's a code snipped
    const layerName = "v_alarm_co";

    const vSource = new VectorSource({
          format: new GeoJSON(),
          loader: function(extent) {
            const layerUrl =
              "/wfs?service=WFS&version=1.1.0&request=GetFeature&typename=" +
              layerName +
              "&outputFormat=application/json&srsname=EPSG%3A3857&bbox=" +
              extent.join(",") +
              ",EPSG:3857";
            xhr.open("GET", layerUrl, true);
            const onError = () => {
              vSource.removeLoadedExtent(extent);
            };
            xhr.onerror = onError;
            xhr.onload = () => {
              if (xhr.status == 200) {
                vSource.addFeatures(vSource.getFormat().readFeatures(xhr.responseText));
              } else {
                onError();
              }
            };
            xhr.send(null);
          },
          strategy: bbox
        });


texts: Layers with same CRS not projecting correctly in QGIS

texts:
In QGIS, I have multiple layers that are not lining up, despite the fact that I have enabled "on the fly" projection and clicked "Save As" to reproject specific layers to CRS EPSG: 26915 (NAD83/UTM zone 15N). I've also tried "Define Current Projection" from "Geometry Tools." No matter what I do, I can't seem to get layers to line up. The oddest part about it is that about half of my shapefiles line up in one part of the map, and the other half line up in a completely different part. All of the previous questions I have found on this topic have not helped me solve this issue, or I'm too dense to understand.

The layers I'm working with are a variety (polygon, point and line), as well as .shp and .kml.


texts: Is there a function or tool that will add vertices to polylines where endpoints are snapped to edges? (ArcMap)

texts: My friend and I are creating a walking network, and we have been having issues getting the routing to work. We just realized that the reason why our routes are messed up is because we snapped our routes to edges as well as endpoints. The routing works fine where the lines are snapped at endpoints because they share a common vertex, but when snapped to edges only one line has a vertex at that location, and the network doesn't see them as connected. Is there a tool or function that will automatically create a matching vertex on the other line?


texts: What is the function of GDA2020 dynamic formats?

texts:
Comments against another question on GDA2020 identified that while EPSG 9 includes static GDA2020 support, there are apparently some dynamic (temporal dependency) variant(s).

How do the GDA2020 dynamic format(s) work (in practice)?


texts: Esri leaflet plugin only showing North America tile?

texts:
I am using the esri-leaflet plugin to connect to an internally hosted ArcGIS Server.  I have tried using a basemap or tilemap implementation but get the same results with both. The problem is that I am only seeing North America and South America with the same tile starting over again at the prime meridian.  I can find the tiles for europe and Asia by looking at them on the arcgis server so I know they are there.

My ArcGIS Server (set up by someone long gone) is saying we are using spatial reference of 4326.

Leaflet Esri plugin states this -

  Your map service must be published using the Web Mercator Auxiliary Sphere tiling scheme (WKID 102100/3857) and the default scale option used by Google Maps, Bing Maps and ArcGIS Online. Esri Leaflet will not support any other spatial reference for tile layers.

Am I able to use this Esri leaflet plugin or do I need the proj4Leaflet with a new CRS projection setting?

I think I am getting tripped up on the tiling scheme of this pre-existing ArcGis server and don't know how to get a leaflet 2D map to work with it.

Some example code -

var crs = L.CRS.EPSG4326;
var map = L.map("map").setView([1,1],1); //just picked a random spot.
L.esri.tiledMapLayer({
  url: myUrl,
  crs:crs,
  minZoom:0,
  maxZoom:12
}).addTo(map);


texts: Leaflet+D3 how to show multiple layers as plane with markers on it

texts: I'm trying to draw map which looks similar to 3 Dimensional. It has 3 layers bottom most is base layer which is OSM map layer. 2 and 3 is svg layer. as shown in below diagram, I'm able to tilt/rotate the map to using css transform functionality. However I'm not able to give background color to svg element using d3 and there by fill the layer with custom color to represent a plane and even not able to transform the polygon lat/long to x,y and z(offset) on svg using d3 such that polygon looks like placed on corresponding layer. Below is the digram I'm trying to achieve.


texts: QGIS qt.qpa.plugin: Could not load the Qt platform plugin "windows" in "" even

texts:
I run often Metshape software from QGIS app, however I recently encounter a problem.
I use this script:
import sys
import os
script_path = "D:/script.py"
commandline = '"C:/Program Files/Agisoft/Metashape Pro/metashape.exe"' + " -r " + script_path
timer_value = 10
os.system(commandline + " & timeout " + str(timer_value))

This script run perfectly with python 3.6 but in QGIS app, I have this error :

qt.qpa.plugin: Could not load the Qt platform plugin "windows" in "" even though it was found.
This application failed to start because no Qt platform plugin could be initialized. Reinstalling the application may fix this problem.
Available platform plugins are: direct2d, minimal, offscreen, windows.

Can anyone guess ?


texts: Displaying a line-chart along with some text fields inside a popUpTemplate.

texts:
I am building a sample web-mapping application using the ESRI 4.5 API. I have a single marker with a popUp on it. When an user clicks on the pop-up I want some text fields along with a line chart inside the popUp template. Right now, I am able to do both individually but I am not finding a way to display both of them together.

<!DOCTYPE html>

<style>
    html,
    body,
    #mapDiv {
        height: 100%;
        width: 100%;
        margin: 0;
        padding: 0;
    }

</style>

<script>
    require([
        "esri/Map",
        "esri/views/MapView",
        "esri/layers/FeatureLayer",
        "esri/views/SceneView",
        "esri/PopupTemplate",
        "esri/geometry/Point",
        "esri/Graphic",
        "esri/symbols/SimpleMarkerSymbol",
        "dojox/charting/Chart2D",
        "dojo/domReady!"
    ], function(Map, MapView, FeatureLayer, SceneView, Point, SimpleMarkerSymbol, Graphic, Chart2D) {

        var map = new Map({
            basemap: "streets",
            ground: "world-elevation"
        });

        var view = new SceneView({
            container: "mapDiv",
            map: map,
            scale: 50000000,
            center: [77.580643, 12.972442]
        });

        var point = {
            type: "point",
            longitude: 77.580643,
            latitude: 12.972442
        };

        var markerSymbol = {
            type: "simple-marker", // autocasts as new SimpleMarkerSymbol()
            color: [226, 119, 40],
            outline: { // autocasts as new SimpleLineSymbol()
                color: [255, 255, 255],
                width: 2
            }
        };

        tagValues = [8, 12, 23, 543, 1, 324, 10]

        var c = dojo.create("div", {
            id: "demoChart"
        }, dojo.create('div'));

        var chart = new dojox.charting.Chart2D(c);
        chart.addPlot("default");
        chart.addSeries("Series A", tagValues);
        chart.addAxis("y", {
            vertical: true,
            fixLower: "major",
            fixUpper: "major",
            title: "Temperature",
            min: 0
        })
        chart.addAxis("x", {
            fixLower: "major",
            fixUpper: "major"

        });

         var lineAtt = {
            Name: "Keystone Pipeline",
            Owner: "TransCanada",
            Length: "3,456 km"
        };

        // Create a graphic and add the geometry and symbol to it
        var pointGraphic = new Graphic({
            geometry: point,
            symbol: markerSymbol,
            attributes: lineAtt,
            popupTemplate: { // autocasts as new PopupTemplate()
                title: "{Name}",
                //How do I change the line below to display the chart along with some text fields.
                content: c
            }
        });

        view.graphics.add(pointGraphic);
        chart.render();

    });

</script>


texts: Exporting image from QGIS precisely on layer extent?

texts:
I have a map of the word, the bottom layer is a world raster layer from Natural Earth and on top of it several vector layer.

Using QGIS 3.4.4, I need to export image without a background, it needs to be precisely on the raster extent.

When I use the Layout manager (Print composer) I cannot get a perfect result, there is always a background.

So, I am looking for a way to perfect extent choosing.


texts: How to Create a Catchment Area of a Polygon (Wetland) using DEM in QGIS or Whitebox GAT

texts:
I am looking for a way to calculate the upstream catchment area where water would flow into a given wetland.

I have delineated the boundary of a wetland (creating a polygon) and would like to determine the upslope catchment area and downslope area of the wetland polygon. My wetland polygon is not on a stream network.

In QGIS, I used the SAGA Rasterize Tool to rasterize my polygon - making sure that my new rasterized wetland has the same projection, resolution, and extent as my DEM. Then I proceeded to use the SAGA Upslope Area tool using my rasterized wetland as the "target area," and my DEM for the Elevation. However, I keep getting the error "Input layers do not have the same grid extent."

As an alternative, I tried using Whitebox GAT using the process described in here:
How to create a catchment from a point using a DTM raster in QGIS?

However, this process is for a point vector, whereas, I am looking for something that can calculate the catchment for a specific polygon or area.

I'm fairly new to using QGIS and very new to using Whitebox. I have some experience with ArcGIS as well, but don't have access to a licence.


texts: Load One Layer from ArcGIS Server MapServer

texts:
I am pretty new to ArcGIS Server and am building a map applications using the Javascript API.  I can easily load the map layer below by adding it as a ArcGISDynamicMapServiceLayer

http://www.parkcounty.org/ArcGIS/rest/services/Roads_Map/MapServer/

As you can see there are many layers as part of the above url.  What if I just want to load the roads layer as a ArcGISDynamicMapServiceLayer?  I've tried replacing the above url with the road specific URL below and that does not seem to work

http://www.parkcounty.org/ArcGIS/rest/services/Roads_Map/MapServer//2

Is there something else I can do as the front end developer or is it up to those that are hosting and managing the server to create a unique url for the roads?


texts: Creating proportional symbols in QGIS that are then graduated by another variable

texts:
I have a table that features two variables that I want to visualise - total population and percentage of total population affected. I need to find a way of including both of these variables in the visualisation and can't see an easy way to do it.
I can create a proportional symbol map easily by graduating the total population by size, but I then want to colour the proportional symbols based on the percentage of affected population. I can't find an obvious way to do that.
I can't manually change the colourings, as the percentages  don't relate directly to the total populations so I need to find some way of writing a query that says create proportional symbols based on the total population and then graduate those based on the percentage affected. I'm using QGIS 3.4.1 if it matters.


texts:
I would like to project a raster using a shapefile read in via the "vect" command. I have made the following work:
r <- rast(xmin=-110, xmax=-90, ymin=40, ymax=60, ncols=40, nrows=40)
r3<-terra::project(r,crs(west))

(where "west" is a shapefile read in with vect)
But what I want to verify is that the crs command used to extract the coordinate reference system from the shapefile is passing on the WKT format? In other words, given the changes to GDAL, if this ref system happened to be in a datum other than WGS84, would this be an appropriate way to pass on the shape's coordinate reference system to the project command? And similarly, passing on the PROJ4 string would NOT be the correct method...


texts: How to check that a point is outside of a route?

texts:
I am newbie in PostGIS and also GIS particularly.

I need to check that a vehicle is outside of a route. This route has been defined in advance, and of course it stores in a table: RROUTEDETAIL. This table contains: ROUTE_ID, LONGITUDE and LATITUDE.

So, when vehicle leaves the route they should be on, the system will give a warning message and then record it into other table.

I imagine to check the point I get with every point store in the RROUTEDETAIL, one by one with ST_DWITHIN. If all check get FALSE value, then the system give a warning message.

Is this the best approach or is there a better idea?

For example, this is the point where vehicle leaves: -7.539079, 109.322891 with value tolerance: 0.5. Then I write:

    select st_dwithin((select st_makeline(ST_Point(-7.539079, 109.322891))),
        (select st_makeline(ST_Point(foo."TRACKERDATA_LATITUDE",foo."TRACKERDATA_LONGITUDE"))
            from (select ht."TRACKERDATA_LATITUDE",
                        ht."TRACKERDATA_LONGITUDE"
                        from "HTRACKERDATA" ht
                        left join "RVEHICLE" rv on ht."VEHICLE_ID" = rv."VEHICLE_ID"
                        where ht."VEHICLE_ID" = 106 and
                        ht."TRACKERDATA_GPSTIME" BETWEEN '2012-01-01 08:00'
                        AND '2012-01-01 17:00'
                        order by ht."TRACKERDATA_GPSTIME", ht."TRACKERDATA_REPORTTIME") as foo),
        0.5);

It always return FALSE. Then I check again with ST_DISTANCE, it gives 99999999999.9 (double precision).

How to fix that?


texts: Creating multipoint from point using PostGIS?

texts:
I am quite new to PostGIS.

I need take some points from an existing table and create a multipoint and save this.

How would i go about this?

The only way i can think without creating a new table is:

Alter table
alter column geom TYPE geometry (MULTIPOINT,3857) USING ST_Multi(geom);

if i had to create a table what would the statement look like when using the x,y of the points in the point table?

Ive been asked to then find the centre of the mutlipoint, I know i need to use  ST_Centroid but what would the code look like?


texts: Is possible to add a web map service like streetdirectory to ArcMap/QGIS?

texts: I would like to know if there is some chance to get web map like streedirectory.sg to arcmap or qgis. I need to use this map as a basemap that's why i'm looking for a way how to get this map to ArcMap. Not sure whether arcgis server is available or not. According my research, i didn't find solution so far.


texts: OSMBuildings building:levels are not working

texts:
I am integrating Leaflet with OSMBuildings and i am using Classic 2.5D. Everything is working and it is showing 3D building on 2D Map, But The Building levels are not working (I am gone through this docs http://wiki.openstreetmap.org/wiki/Simple_3D_buildings).

Findout in HTML code as mentioned in **"building:levels": 12**, I have added this line to my geojson property.

Here is my HTML code :

<!DOCTYPE html>

    <meta charset="UTF-8">
    <title>GeoJSON</title>
    <link rel='stylesheet prefetch' href='http://cdn.leafletjs.com/leaflet-0.7.3/leaflet.css'>
        <link rel="stylesheet" href="css/style.css">

Objects on the map are build from custom GeoJSON.

        Edit the data or paste your own examples here.

<textarea>
{
  "type": "FeatureCollection",
  "features": [{
    "type": "Feature",
    "geometry": {
      "type": "Polygon",
      "coordinates": [[[72.8999412059784,19.0702453117883],[72.8998178243637,19.0699157606889],[72.9002845287323,19.0700019510397],[72.9003381729126,19.0703416419859],[72.8999197483063,19.0702351717642],[72.8999412059784,19.0702453117883]]
      ]
    },
    "properties": {
      "color": "rgb(255,200,150)",
      "height": 150,
      **"building:levels": 12**
    }
  }]
}
</textarea>

<button onclick="setGeoJSON()">Apply changes</button>

and my JS file is :

var map = new L.Map('map');
map.setView([52.52111, 13.40988], 16, false);
new L.TileLayer('http://{s}.tiles.mapbox.com/v3/osmbuildings.kbpalbpk/{z}/{x}/{y}.png', {
  attribution: 'Map tiles &copy; MapBox',
  maxZoom: 18,
  maxNativeZoom: 20
}).addTo(map);
var osmb = new OSMBuildings(map);
//************************************************************
var geoJsonBlock;
function getCenter(geoJson) {
  var geometry, len = 0, lat = 0, lon = 0;

  for (var i = 0, il = geoJson.features.length; i < il; i++) {
    geometry = geoJson.features[i].geometry.coordinates[0];
    if (geometry[0][0][0]) {
      geometry = geometry[0];
    }

    len += geometry.length-1;
    for (var j = 0, jl = geometry.length-1; j < jl; j++) {
      lat += geometry[j][1];
      lon += geometry[j][0];
    }
  }

  return { lat:lat/len, lon:lon/len };
}

function setGeoJSON() {
  try {
    var geoJson = JSON.parse(geoJsonBlock.value);
  } catch(e) {}
  osmb.set(geoJson);
  var center = getCenter(geoJson);
  map.setView([center.lat, center.lon], 17)
}

geoJsonBlock = document.querySelector('textarea');
setGeoJSON();

Can anybody help how i am displaying my 3D building to levels, roofs, sphere?


texts:
I find that when things are not working as they should repairing the geometry on everything is a good first point of troubleshooting so I run a few lines of python:

mxd = arcpy.mapping.MapDocument("CURRENT")
for lyr in arcpy.mapping.ListLayers(mxd):
    arcpy.RepairGeometry_management(lyr)

Of course this only works with certain layers, so to make this applicable when there are rasters I've added:

if lyr.isFeatureLayer == True:

But I still have the problem when there are SDE layers. How can I either test for SDE or test for layers that are a Shapefile or a feature class in a Personal or File Geodatabase?


texts: Selecting raster pixels based on dem range

texts:
I have two raster layers that are equal extent and resolution: one that is qualified as nonforest tree cover that is made up of pixels with either the value of 1 or NA, and one DEM.

What I am hoping to do is to essentially create a varying buffer around each pixel of value 1 in the nonforest tree cover raster layer that includes +200 and -200 meters of elevation.

I imagine this would have to be accomplished by going through each pixel where the value is 1 (non-NA) and progressively testing each adjacent pixel to see if it falls within 200 meters of the original pixel elevation and continuing from each within range pixel until the desired elevation range is exceeded, however I am unsure of how best to do this. My desired outcome is a raster layer where all of the relevant pixels share one value.


texts: Choice approaches for passing a selected feature to a method in ArcObjects? (C#)

texts:
It pains me to ask such a basic question, but after several hours (not exaggerating) of trying to figure this out and failing, I come to GIS.se

What I am trying to do is pass a selected feature and/or layer to a given method. In this example, it is to add a field to the selected layer in the ToC. The snippet calls for a feature class that is a member of the IFeatureClass interface (if my vocabulary is off, please correct me):

public void AddFieldToFeatureClass(IFeatureClass featureClass, IField field)

I create an IField easily enough:

IField statusField = new FieldClass();

But what is the preferred, cleanest, or most effective method to pass an IFeatureClass that is the currently selected layer in the TOC?


texts: Running QGIS Python code from the R rPython package

texts:
The following script works fine from python or ipython at the command line:

from qgis.core import *
from qgis.gui import *
import sys
sys.path.append("/usr/share/qgis/python/plugins")
from PyQt4.QtCore import *
from PyQt4.QtGui import *
app = QgsApplication([], True)
QgsApplication.initQgis()
from processing.core.Processing import Processing
Processing.initialize()
import processing

layer = QgsVectorLayer('ger.shp','ger','ogr')

Where ger.shp is a shapefile in the current folder. After this, layer.isValid() is True. No warnings or messages.

But if I try and run it from the rPython package, I get an error message and the layer is not valid:

library(rPython)

qtest <- function(){
    code = "from qgis.core import *
from qgis.gui import *
import sys
sys.path.append('/usr/share/qgis/python/plugins')
from PyQt4.QtCore import *
from PyQt4.QtGui import *
app = QgsApplication([], True)
QgsApplication.initQgis()
from processing.core.Processing import Processing
Processing.initialize()
import processing

layer = QgsVectorLayer('ger.shp','ger','ogr')
"
    rPython::python.exec(code)
}

qtest()
QObject::connect: Cannot connect (null)::raiseError( QString ) to QgsVectorLayer::raiseError( QString )

and then:

python.get("layer.isValid()")
[1] FALSE

It is identical code and identical python version, so I think it can only be in the way the code is run - perhaps the python interpreter is initialised differently from the way embedded python is run in rPython? I'm stuck.


texts:
I have a DEM in ArcGIS 10.1 that I would like to get into STL format for a 3D printer. Does anyone know some software that will allow me to make this conversion?

I have to use the STL format because that is the format required by my university's 3D printer.


texts: Shapefile incompatibility between ArcGIS Online and QGIS 3.4 for a point feature layer

texts:
I created a hosted feature point layer ('Threatened Plant Locations') in an ArcGIS Online webmap, collected some data using Collector for GIS app on my iPad, synced data back to my ArcGIS Online webmap, exported (as .shp and GeoJSON files) and downloaded data to my PC and then imported .shp and GeoJSON files to my QGIS 3.4 project map. I noticed that the points in my QGIS map are up to 10 metres away from their collected 'true/actual' positions. The points are correctly placed in the ArcGIS Online webmap and when I uploaded the shapefile to ESRI support, they informed me that the points loaded correctly in ArcMap as well.
Note the point feature layer source data is in EPSG: 3857 WGS84 (ArcGIS Online/Collector data is always WGS84) as is the source CRS for the public imagery layer I am using as a basemap in QGIS (NSW LPI). The CRS I am using for my QGIS Project map is GDA94 MGA56 (EPSG:28356) as I want all layers rendered/displayed to this CRS. I thought this might be a projection issue but when I exported the shapefile in order to create a new reprojected shapefile (with WGS84 CRS) there was still no change in the point positions on my QGIS map.
There was a similar question asked a number of years ago but the only answer related to reprojecting shapefile which I have already tried. Not sure it is a projection issue?
Any there any ArcGIS Online/QGIS users out there?


texts: Intersect shapefile with raster, then determine proportion of each raster value?

texts:
Using the arcpy module in Python, I want to find the proportion of each pixel value within the boundary of an overlaying polygon.  (The polygon is in a shapefile.)

For example, there are 10 possible values in the raster, and intersecting the polygon 17% of the pixels are class1, 8% are class2, etc.

This may involve use of the Zonal Statistics tool, but apparently I'm having trouble working out the syntax programmatically.


texts: Huge point cloud laser data in PostGIS - Storing and Processing it

texts:
I wonder, how it's possible to store huge sets of laser scanned point cloud data in PostGIS, with the time-aspect for processing it in mind. I know, there exists a geometry-object Point in PostGIS. But as far as I know it saves each point in a new tupel, which can make searching for any certain point a very slow process, if a few millions or more of them are stored.

I found a paper from HSR Universtiy of Applied Sciences Rapperswill, discussing this topic. It suggests three ways to store such data: Whole data in one tupel, Each point in one tupel or Splitting Data into Blocks which are referenced by info-tables, holding the extends of each block. As the third way seems the most useful for locating stored points, I wonder if anyone already has made some experiences with it?

The paper can be found here: http://wiki.hsr.ch/Datenbanken/files/pgsql_point_cloud.pdf

Last but not least, I stumpled across a project on github, which seems to deal with point cloud manners in PostgeSQL. Unfortunately not much information about it around the net. So the same question here: Has someone already made some experiences with it? Is it usable for such purposes?

Project can be found here: https://github.com/pramsey/pointcloud

I would also be glad to hear about other suggestions, ideas or experiences, if there are any. But i must admit, that non-commercial solutions are prefered.


texts: Convert a ee.Geometry() into a Leaflet polygon

texts:
I'm using the Leaflet and gee API in a Python notebook.
In this notebook I'd like to display a featureCollection as a layer on a Leaflet map.
#retreive the asset
asset_descripsion = 'Glad_' + country_selection + '_boundaries'
folder = ee.data.getAssetRoots()[0]['id'] + '/'

asset = ee.FeaturCollection(folder+asset_descripsion)

#create the map
m = Map(center=center, zoom=zoom)
dark_matter_layer = basemap_to_tiles(basemaps.CartoDB.DarkMatter)
m.add_layer(dark_matter_layer)

#add the layer
asset_geom = ee.FeatureCollection(asset).geometry().getInfo()['coordinates']#polygon = Polygon(
    locations=asset_geom,
    color="green",
    fill_color="green"
)
m.add_layer(polygon);
m

Problem the way leaflet and gee construct polygon is not compatible so my Map end up like that

instead of the Argentina boundaries that I have in GEE

Does anyone know how to perform this trick ?


texts: Loop clipping Imagecollection by features of a vector in Earth Engine and exporting all tifs

texts:
I am looking for a method to clip all the images in an ImageCollection through all the features of a vector (world borders), and then save the resulting tifs in my Google Drive.

So far I wrote this script, which does not work but it gives an idea of what I would like to do. A part for looping over the save with the corresponding filename is missing, as I am not sure how to write it.

var coll = ee.ImageCollection('NOAA/VIIRS/DNB/MONTHLY_V1/VCMSLCFG').select('avg_rad');
var geometry = ee.FeatureCollection("users/gg/gadm");

var loop = function(image, feature){
    var collection_roi = image.clip(feature);

  // missing part to save each output to drive

    }

var stacked = loop(coll, geometry);

Edit:

Using the following code I did manage to create two-argument loop. However, I end up with a list of imagecollections. I am not sure how to save on my drive all the images of each imagecollections in the list.

var coll = ee.ImageCollection('NOAA/VIIRS/DNB/MONTHLY_V1/VCMSLCFG').select('avg_rad');

var geometry = ee.FeatureCollection("users/gg/gadm");

// loop over features and each image of
var loop = function(image){
 var output = geometry.map(function(geom){
 // Clip the image to the polygon geometry
    var ic = image.clip(geom);
    return ic;
})
  return output
};

var outputd = coll.map(loop)

var listOfImages = outputd.toList(outputd.size());


texts: When to replace Google Maps with MapServer?

texts:
I have a contract handed to me to create a google map interface for a real estate company where the area they are operating is moving sales around in extremely high volume.  To date, I have created 12 maps of polygons of zip code boundaries.  The 12 maps are attached to a single, multi-column spreadsheet in Google fusion tables, and each map uses a different column to create a value-based fill color.

Now I'm being asked to create heatmaps for each zip code.  Unfortunately we are looking at thousands of markers per zip code, and potentially we will need to heatmap zipcode against zipcode.  The Proof Of Concept zipcode has over 16,000 markers (lots of homes).

My question is this:  At what point do you tell your client that google maps will not work this far for you, and you MUST think about switching to a GIS server (ie. mapserver) to manage your interface?

Of note, he is renting these maps out as service for real estate managers around the area, so the traffic on this website will be high.


texts: Update records with ArcPy update cursor? Geodatabase format

texts:
I try to update thousands of rows when all records get the same value.
My code has a loop that runs thousands of times.
Is it possible for the program to make the change to all the columns in one command and save run time ?

For example, php+mySQL can do this:

mysql_query("UPDATE all_records SET Xcenter = $val1, Ycenter = $val2 WHERE ID > 100");

This is my code:

    for row in rows:
    row.setValue("Xcenter", val1)
    row.setValue("Ycenter", val2)
    rows.updateRow(row)


texts:
I have this map with two controls

<!DOCTYPE html>

    <title>Simple Map</title>
    <link rel="stylesheet" href="http://openlayers.org/en/v3.17.1/css/ol.css" type="text/css">

I want to select + or -, and when it is selected, I want to zoom in or zoom out when I click on my map.

How can I do this?


texts: Can ArcGIS Pro use Python toolboxes written for ArcMap?

texts:
Can the new ArcGIS Pro use existing ArcToolbox scripts?

Or to put it another way, will I have to rewrite my existing ArcMap toolbox & python scripts if I want them to run in ArcGIS Pro?


texts:
Is there a solution to fill a popup with content getting from a AJAX call, without JQuery in OpenLayers 5?
The content is in *.txt files, or I can handle it to send the content as a JSON or HTML via PHP script.  I searched for a solution to send all the data for all the popups not by the first web page call, but only when a popup is opened.
My code:
var feature = new ol.Feature({
    type: 'click',
    desc: popupcontentstring,
    geometry: new ol.geom.Point(ol.proj.fromLonLat(coordarray_coord))
});
feature.setStyle(style);
sourceVector.addFeature(feature);
}

popupcontentstring is the variable with the content for the popup. That I will replace with a AJAX request.


texts: I am trying to change the default point layer symbol to a set one.  For example I have a fire hydrant point layer that is represented by a default point symbol when I import it into my map.  Is there a way I can change the default symbol to one I select and it will continue to be the default point symbol?


texts: Create certain number of Clusters, equal sum of Z value

texts:
Divide polygons into *n* number of groups of equal counts with ArcGIS 10.2 is very similar to my question. But instead of equal number of polygons, I have multiple points (schools) with a z value (students in a special program). I am in 10.2. I have spatial analyst.

I need a certain number of groups (based number of program aides) that each have about the same sum of all their point's z value and are grouped spatially.  The idea is that each aid goes administers the program at schools in an area and all of the aids have about the same number of students(z value) to help.

I have tried grouping analysis to see of it could get me close and no dice.

I know I can do this manually but I have a feeling this analysis will be required again for a different program so I would like to automate it as much as possible.

Attached is a proportional symbol map showing the range of z value and the general location of schools. I want to group it into 3 groups, same sum of z values.

EDIT/UPDATE

I think I will just visually (and with input from program director) pick three Central locations and create service areas using network analyst. Once I have finalized the appropriate service areas I will export them into a polygon shapefile and sum the number of z values for the school that are in the nearest service areas. Any schools left in a middle zone will be distributed to the group with the lowest z value sum.


texts: How can I merge roads layers that do not align perfectly?

texts:
I have multiple roads layers for a given area and would like to merge them into one layer since they each have some unique roads the other layer does not show. However, the roads that these layers share in common do not align perfectly. Therefore, if I just merge the two layers, all of those roads will be shown twice in slightly different places.

Is there any way to make some sort of tolerance where it won't duplicate roads displayed a few meters apart, for example? Or should I just merge them anyway and deal with a messier and somewhat more confusing layer when they're merged?


texts: Getting local incidence angle from Sentinel-1 GRD image collection in Google Earth Engine

texts:
I am using python API in google earth engine (colab) to retrieve local incidence angle from sentinel-1 GRD data. However, the .select ('') function only allows the VV, VH, HH, HV, and 'angle' band selection. The angle band is viewing angle of the sensor, and I need a local incidence angle (retrieved as a by-product of the terrain correction step).
Does anyone have the answer?
Thanks


texts: How to create map tiles out of 25 images?

texts:
I am an absolute beginner of GIS and have the following situation:

I need to create an offline map for Android,
I have 25 pictures that represent a region that needs to be mapped in my application,
I have the top left and the bottom right coordinates of the map (of the 5x5 grid comprised of the 25 images).
I am thinking of using the Nutiteq SDK manage the map in my application (in order to add markers and user position).

The Nutiteq Maps SDK supports: MBTiles, MGMaps, GeoTIFF and other formats.

How do I create a usable format file out of 25 images and top left, bottom right coordinates?


texts: Exporting 3 layers together as JPEG using Data Driven Pages

texts:
I'm trying to use data driven pages to export JPEGs of features that share a common attribute.
E.g. I have two point layers, "beach" and "light". And a polyline that connects them, called "vector".
Each of these features has a common field "Waypoint". I want to export a JPEG of all 3 layers at once that have the same waypoint value, e.g. "1".
The problem that I am encountering is that data driven pages are focused on one "index" layer. So I need to try to figure out the logic to use in ArcPy to overcome this.
Any suggestions?
https://support.esri.com/en/technical-article/000013432


texts: What is the maximum number of vertices for polygons & polylines in PostGIS?

texts: What is the maximum number of vertices for polygons and polylines (including multipart) in PostGIS?


texts: Fresh Install of ArcMap 10.4 Stalls on "Loading Document"

texts:
I have just installed the trial version of the ArcGIS for Desktop which contains ArcMap on my Windows 7 Pro 64-bit.

I haven't been able to launch ArcMap, ArcScene or ArcGlobe yet. As of ArcMap it loads that initial banner and then stalls on Loading Document, then after a minute closes.

I haven't ever installed or used ArcGIS. I have already tried:

Remove the ESRI folder from %APP_DATA%/Roaming
Remove the ESRI key from the Registry
Run as administrator

Don't know what is the problem and support won't answer if I don't buy the product.


texts: How to determine which stream end point is at the basin outlet?

texts: I am working with streams, stream end points, elevation raster and I have created basins using the watershed tool in ArcGIS.  From that data, I would like to determine for each basin, which single point is at or nearest the outlet as shown by the arrow in the example below.


texts: How to add group layers/TOC to ArcGIS web map/web application?

texts:
I am new to arcgis online and am trying to create a web map that has some kind of grouped layers. I have tried adding multiple shape files that are all zipped together but because together the files are too large to upload onto the map I have to resort to some other way to do this.

I have come across TOC coding for API and I have never used API before. I was wondering if there was a way to code the TOC or group layers into a web map or web application (without using a template) through the google chrome developer tool? I want to create my web app through the Web Application Builder if possible.

Does anyone know any way to use the developer tool kit on chrome to do this?


texts: Performing TIFF to MBTILES conversion on Mac?

texts:
I have a large number of tiff files (map tiles) along with tab and tfw file types as part of an Ordnance Survey download. This is for an offline iOS app. I managed to get the tiffs into QGIS with correct projections and save as a virtual raster. What I need to do now is somehow convert or save it to an mbtiles format.

How can this be done on a Mac?


texts: CartoDB Twitter code bl.ocks non-functional with namedmap type viz

texts:
I'm working with the code example shown here http://bl.ocks.org/javisantana/6029822
When I substitute Javi's viz for a namedmap map type viz it fails with the following error
"TypeError: r is not a constructor cartodb.js:18"
Swapping out the js include from the  dev version to the  production version
the map loads, but the info windows fails populate or get fully styled.
My Question - Does anyone see what might need to be changed in order for this code to function properly with namedmaps?
If not any help in understanding the difference between namedmap and layergroup type vizs would be very helpful.


texts: Do pad functions of QGIS Field Calculator have bug?

texts:
I have this expression in the field calculator:

regexp_replace('1', '^(\\d+)', rpad('\\1', 3, '0') ) gives 01 rather than 001 as expected.

the *pad functions work fine when called outside regexp_replace

Weird!   This has to be a bug? or am I missing something?

If it is a bug I'll log something in the issue tracker.

BTW I was delighted when I found I could use complex expressions including functions in the third parameter of regexp_replace!   kudos to the QGIS developers!


texts: A strange phenomenon on ST_Intersection() in PostGIS


texts:
I am attempting to create a new raster layer based on the values of one raster layer using a formula that is conditional based on the values in another layer.

Two source layers:

> ndvi
class       : RasterLayer
dimensions  : 7741, 7611, 58916751  (nrow, ncol, ncell)
resolution  : 30, 30  (x, y)
extent      : 580185, 808515, 4191885, 4424115  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : C:/Users/xxx/AppData/Local/Temp/RtmpScc57V/raster/r_tmp_2019-04-22_112901_2168_61926.grd
names       : layer
values      : -1, 1  (min, max)

> pv
class       : RasterLayer
dimensions  : 7741, 7611, 58916751  (nrow, ncol, ncell)
resolution  : 30, 30  (x, y)
extent      : 580185, 808515, 4191885, 4424115  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : C:/Users/xxx/AppData/Local/Temp/RtmpScc57V/raster/r_tmp_2019-04-22_112936_2168_54474.grd
names       : layer
values      : 0, 1  (min, max)

Layers each have the same dimensions, projection, and extent.

I want to create a new layer, e, that has the value:

0.97 when ndvi < 0.2
0.004 * pv + 0.986 when ndvi >= 0.2 and ndvi <= 0.5
0.99 when ndvi > 0.5

My statement is:

e <- overlay(ndvi, pv, fun = function(ndvi, pv)
    {
        ifelse(ndvi < 0.2, 0.97, ifelse(ndvi > 0.5, 0.99, 0.004 * pv + 0.986))
    }
)

However, when I execute this, I receive the error:

  Error in ifelse(ndvi > 0.5, 0.99, 0.004 * pv + 0.986) :    argument
  "pv" is missing, with no default

How can I write this function so that it is compatible with overlay()?  Is another approach more appropriate?

I'm somewhat new to R.

Other permutations attempted:

calc_e <- function(ndvi, pv)
    {
        ifelse(ndvi < 0.2, 0.97, ifelse(ndvi <= 0.5, 0.004 * pv + 0.986, 0.99))
    }

e <- overlay(ndvi, pv, fun = calc_e)

Same error was returned.

Then I tried:

ndvi_2 <- ndvi
ndvi_2[ndvi_2 < 0.2] <- 0.92

ndvi_m <- ndvi
ndvi_m[ndvi_m >= 0.2 && ndvi_m <= 0.5] <- 0.004 * pv + 0.986

ndvi_5 <- ndvi
ndvi_5[ndvi_5 > 0.5] <- 0.99

The first and third subset assignments ran.  However, I received an error on the second one "ndvi_m[ndvi_m >= 0.2 && ndvi_m <= 0.5] <- 0.004 * pv + 0.986"

  Error in ndvi_m >= 0.2 && ndvi_m <= 0.5 : invalid 'x' type in 'x && y'


texts:
I am building a map for the northeastern U.S.A. The map background needs to be either an altitude map or a mean annual temperature map. I have two rasters from Worldclim.org which give me these variables but I need to clip them to the extent of the states I am interested in.
How do I do this?
This is what I have so far:
#load libraries
library (sp)
library (rgdal)
library (raster)
library (maps)
library (mapproj)

#load data
state<- data (stateMapEnv)
elevation<-raster("alt.bil")
meantemp<-raster ("bio_1.asc")

#build the raw map
nestates<- c("maine", "vermont", "massachusetts", "new hampshire" ,"connecticut",
  "rhode island","new york","pennsylvania", "new jersey",
  "maryland", "delaware", "virginia", "west virginia")

map(database="state", regions = nestates, interior=T,  lwd=2)
map.axes()

#add site localities
sites<-read.csv("sites.csv", header=T)
lat<-sites$Latitude
lon<-sites$Longitude

map(database="state", regions = nestates, interior=T, lwd=2)
points (x=lon, y=lat, pch=17, cex=1.5, col="black")
map.axes()
library(maps)                                                                  #Add axes to  main map
map.scale(x=-73,y=38, relwidth=0.15, metric=T,  ratio=F)

#create an inset map

 # Next, we create a new graphics space in the lower-right hand corner.  The numbers are proportional distances within the graphics window (xmin,xmax,ymin,ymax) on a scale of 0 to 1.
  # "plt" is the key parameter to adjust
    par(plt = c(0.1, 0.53, 0.57, 0.90), new = TRUE)

  # I think this is the key command from http://www.stat.auckland.ac.nz/~paul/RGraphics/examples-map.R
    plot.window(xlim=c(-127, -66),ylim=c(23,53))

  # fill the box with white
    polygon(c(0,360,360,0),c(0,0,90,90),col="white")

  # draw the map
    map(database="state", interior=T, add=TRUE, fill=FALSE)
    map(database="state", regions=nestates, interior=TRUE, add=TRUE, fill=TRUE, col="grey")

The elevation and meantemp objects are the ones that need to be clipped to the area extent of the nestates object.


texts: Visualising layer properties of a GeoJSON polygon layer in Leaflet

texts:
In a Leaflet GeoJSON polygon layer, I want a popup information box, showing attributes of each of the polygon properties. So the final output should show the Country Name, Area, Population etc properties accessed from the GeoJSON layer of the particular polygon; something like the image.

For this reason what should I write inside Popup?

gj = folium.GeoJson(open("country.geojson",encoding = "utf-8-sig").read())
popup = folium.Popup('Country name: ')
popup.add_to(gj)
gj.add_to(map)


texts: Remembering selected items on server round trip

texts:
I am using GeoServer, GeoExt, PostGIS, etc (all the standard stuff).

I have a map displayed in my browser and performing a box select adds relevant items to my selection vector layer.  In addition to this a popup box appears with a summary of selected items and against each layer summary is a hyperlink which when clicked causes a server round trip (this round trip is necessary as it displays the main detail of each selected item).  The server round trip results in the map being re-displayed but now my selected items are no longer selected.

Can I maintain the list of selected items and keep them highlighted after the map re-draw?

I could probably cobble something together using cookies but was wondering if there is an easier solution!


texts: Geoserver - how clear memory after tile-seeding and normal wmts queries

texts:
I am currently doing a lot of seeding of tiles in GeoServer and I noticed that after the tiling process is over the memory that is required during the task is not cleared. The same thing happens after a WMTS query, the memory usage goes up but isn't cleared after the request has been fulfilled. The memory is only cleared after I restart GeoServer.

I am running it in Tomcat7 and I am using the following JAVA memory options:

JAVA_OPTS="-Djava.security.egd=file:/dev/./urandom -Djava.awt.headless=true -Xmx8192m -XX:MaxPermSize=512m -XX:+UseConcMarkSweepGC"

I have memory BlobStores turned off as well as "Use memory mapped buffers" and "Cache and reuse memory maps" in the Data Stores. I did notice that when you edit a layer and go to "Tile caching" there is an option to "Enable In Memory Caching for this Layer", which is selected and I can't turn it off as it's greyed out. Does anyone know what the solution to this problem is?


texts: I have hypsometric map/model of mountain which I get by contour lines with Z value. Can you explain how to visualize 3d model in ArcGIS? I need to get something like nviz command in GRASS.


texts: Exported files from GeoTracker (kml,kmz.gpx) are not visible after certain zoom in QGIS

texts: I have problem with track recorded with GeoTracker on my Android. Exported track, either in any of given formats: gpx, kml or kmz, is not visible after certain zoom level - scale 1:6.141 (I am using QGIS v2.18.5), and I must note that this is only happening with linestring, not points. Option to set scale visibilty for selected layers in QGIS is not enabled, so QGIS is not a problem. I also exported all the files into shapefile just to try if it would maybe work, but it is still not visible, because logically I guess it withdraws attributes from the original file. I ve been checking out visibilty in files itself, however it states 1. I didnt have problem with track I ve recorded after this one, where configuration and options were exactly the same. I need to be able to edit track somehow.


texts: How to find which table has layer in mapinfo?

texts:
Is it possible to find if table has layer Or print list of layers?
 Edit::
I have a table when I selete some data from table like:

Select * from earth where population >100

Then it create new table but if want to hide layer closing table would not work so to hide layer I need to know table that are used to create layer?


texts: Converting from Gauss-Kruger to Bing coordinates / projection

texts:
I am trying to convert a shapefile with ogr2ogr from Gauss-Kruger (zone 4) to Bings coordinate system / projection using the advice given in http://alastaira.wordpress.com/2011/01/23/the-google-maps-bing-maps-spherical-mercator-projection/:

ogr2ogr -overwrite -s_srs EPSG:31468 -t_srs Bing-Projection.txt out.shp in.shp

Where Bing-Projection.txt contains

PROJCS["Popular Visualisation CRS / Mercator",
 GEOGCS["Popular Visualisation CRS",
  DATUM["WGS84",
    SPHEROID["WGS84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7059"]],
  AUTHORITY["EPSG","6055"]],
 PRIMEM["Greenwich", 0, AUTHORITY["EPSG", "8901"]],
 UNIT["degree", 0.0174532925199433, AUTHORITY["EPSG", "9102"]],
 AXIS["E", EAST], AXIS["N", NORTH], AUTHORITY["EPSG","4055"]],
PROJECTION["Mercator"],
PARAMETER["semi_minor",6378137],
PARAMETER["False_Easting", 0],
PARAMETER["False_Northing", 0],
PARAMETER["Central_Meridian", 0],
PARAMETER["Latitude_of_origin", 0],
UNIT["metre", 1, AUTHORITY["EPSG", "9001"]],
AXIS["East", EAST], AXIS["North", NORTH],
AUTHORITY["EPSG","3785"]]

But I get the following output:

ERROR 6: No translation for Mercator to PROJ.4 format is known.
Failed to create coordinate transformation between the
following coordinate systems.  This may be because they
are not transformable, or because projection services
(PROJ.4 DLL/.so) could not be loaded.

See this gist for the full output.


texts: Executing original features of merged layers using QGIS

texts:
I am using QGIS 3.10.
Following this query:
QGIS 3.10 problem with merging layers & multipart promotion
I've managed by merging the layers featuring the same geometry.
Unfortunately, the effect is not enough for me, since the newly merged layer contains its own feature (especially with styling).
Basically I merged two layers with completely different styling, as per below:

but my merged layer has its own stylisation, as it worst only 1:

I want to have the layers merged, but with the original styles kept.
Is it possible?

I tried to categorize these layers as per the image below:

creating the custom styles in the "categorized section". Unfortunately in the effect, my shape disappeared.

I don't know how to do now.
I set "Area" as a value, although I noticed, that my object is gone once I switch on the "Categorized" option.


texts: Round Raster to next higher or lower int

texts:
I am using ArcGIS 10.1.

I am searching for a solution to calculate a raster.
The values are like:
0.2
0.6
5.3
i.e.
How can I now get 0.2 as 0 and 0.6 to 1 and 5.3 to 5? In the end integer numbers.
The "roundup" and "rounddown" functions are worse to do this.

Thanks for any hints.

Luigi


texts: Symbol Level in ArcGIS 9: 'Broken' look of symbols, how to resolve

texts:
I am using symbol levels in ArcMap to style roads so that I have a banded feel. My problem is that at intersections and sharp turns, the cartography gives a broken effect in the road. How can this be resolved. I have tweaked and searched without success.

.

 Unwanted effect pointers in blue.

UPDATE EDIT:

After further investigation I had an effect closer to what I wanted but with some side effect. From left to right; First image shows (i)the settings in ArcMap (Symbol Property Editor) , (ii) Resulting effect solving the 'craked effect' but undesirable capping at the end and (iii) the desired effect at the end of a road.

.


texts:
What is the difference between NAD CORS96 (Illinois East 1201) and the "normal NAD1983?

We have data that was GPS'ed in CORS96, but data was imported into ArcMap using the "normal" NAD1983?  What effect, if any, will it have on data? For example, will it affect the location and the area calculations of the data?


texts:
I am developing a web app in geodjango (django) framework but I need to import raster files.

Is it possible a geodjango app read raster files from a directory?


texts: Converting Shapefiles to Mapbox Vector Tiles

texts:
We are currently hosting OSM vector tiles on our mapping server (courtesy osm2vectortiles and tilesever-gl). Now I'm trying to convert our own shapefiles into the same Mapbox Vector Tile format so we can add them as an overlay to our maps.

I believe one way this can be accomplished is by importing the shapefile to PostGIS and then creating vector tiles from that, but I'd like to avoid having to use PostGIS as a middleman if I can, and instead go directly from Shapefile->Mapbox Vector Tiles (which will be stored in an MBTiles container). These will then be served via tileserver-gl.

Does any such tool exist that can convert shapefiles (1-2GB) directly to vector tiles, and ideally perform simplification while it does so?


texts: PostGIS - Fast area computation on large overlapping (or not) polygon dataset

texts:
I'm co-developing an PostGIS based API that will work with complex and potentially large sets (>10k) of polygons, extracted from a large table following custom, user defined filters. I have precomputed the area of every polygon and put it on an 'area' field. Sometimes the selected polygon set will have no overlapping at all, most of the times there will be ~10% overlapping. I need to efficiently compute  the total, non-overlapping area of the filtered dataset.
The easy answer is to use st_Area(st_Union(geom)::geography)), but it will be inefficient on large datasets with little overlapping, and completly unnecessary if there's no overlapping (I'd just need to add up the 'area' field).
I wonder if there is way to create a spatial function (let's call it smart_area) that will:

Evaluate if there's any overlapping at all on the input dataset
If there's no overlapping, it will return the sum of the 'area' field
If there's overlapping, it will compute the area of the overlapping polygons and then add it to the summed area of the 'area' field of the rest of the polygons.

To generate a working example, I have modified the routine presented on http://boundlessgeo.com/2014/10/postgis-training-creating-overlays/ :
DROP TABLE IF EXISTS pts;
DROP TABLE IF EXISTS circles;
-- CREATE RANDOM POINTS
CREATE TABLE pts AS
  WITH rands AS (
  SELECT generate_series as id,
         random() AS u1,
         random() AS u2
  FROM generate_series(1,10000)
)
SELECT
  id,
  ST_SetSRID(ST_MakePoint(
    10*u1,
    10*u2), 4326) AS geom
FROM rands;
-- CREATE CIRCLES WITH AREA (in km²)
CREATE TABLE circles AS
  SELECT id, ST_Area(ST_Buffer(geom, .01)::geography)/1000000 as area, ST_Buffer(geom, .01) AS geom
    FROM pts;
-- COMPUTE TOTAL AREA
SELECT SUM(area) FROM circles;
-- COMPUTE TOTAL NON-OVERLAPPING AREA
SELECT ST_Area(ST_Union(geom)::geography)/1000000 FROM circles

This piece of code will generate 10.000 circles, some of them will overlap. The second query, which computes the non-overlapping area (the area that I want to compute) took around 10 seconds on my computer, while the first one is almost instantaneous, but inaccurate.
The 'smart_area' function would be a mix of both queries, as it would add up the 'area' field of the non-overlapping polygons and compute the st_union area of the overlapping ones. I guess (I hope) that this function will be way faster than the second query.
Any insight on this problem?


texts:
I am trying to digitize curved roads in QGIS but cannot find the right tool.

I want somthing like end point arc segment in arcmap as shown in figure below.


texts:
I am trying to assign a value from one field of my attribute table to a new field based on the 3D field...

The code block I am thinking about is as follows:

def conditional (!FarmCateg!):
   if (FarmCateg=='4'):
       return
        '''(and here I want to return values from "Percent Area" Field (float)):'''
   elif( ): return "0"

Expression: conditional(!FarmCateg!)


texts: I'm unable to find and run NVIZ and now I've read several places that it's not functioning in the latest versions of QGIS? Is this a temporarily bug or definite?


texts: Quickly temporarily disabling snapping in QGIS 3?

texts:
In some other applications I use you can use one of the "shift" keys to disable snapping for a single operation.

Is there an equivalent for the new node tool?


texts: Purchasing Points of Interest data in USA?

texts:
I have thousands of POIs with addresses and they are geocoded.  The lat longs are on the road.

I have Tele Atlas POI data, but it also shows the locations closer to the road.

Is there any company which sells POI data (for the USA) and the exact location(lat long) on the ground?


texts: Protecting KML/KMZ files created in ArcMap from being copied?

texts: Is there a way to lock/protect a kml/kmz file that is created in ArcMap from being copied in Google Earth and/or geoprocessed in ArcMap when you want to give to third parties?


texts: ST_ValueCount() returning multiple rows of same VALUE

texts:
I want to calculate the proportion of each pixel value in a raster after clipping a much larger raster to my area of interest. However, after clipping and saving the raster from PostgreSQL database/PostGIS, the table created does not calculate the pixel number or proportion as expected.

Step 1: clip larger raster (Minnesota gridded soils) based on extent of "plotbounds" (polygon geometry table):

CREATE TABLE out.eonr_gaylord_2012_soils AS
SELECT ST_Clip(rast, p.geom)
FROM gssurgo_mn.mapunitraster_mn_10m, temp.eonr_gaylord_2012_plotbounds_data AS p
WHERE ST_Intersects(p.geom, rast);

Step 2: Calculate the proportion of each pixel

SELECT DISTINCT pvc.VALUE, SUM(pvc.COUNT) AS tot_pix, 100 * SUM(pvc.COUNT) / SUM(SUM(pvc.COUNT)) OVER () AS prop
FROM out.eonr_gaylord_2012_soils, ST_ValueCount(rast, 1) AS pvc
GROUP BY pvc.VALUE, pvc.COUNT
ORDER BY pvc.VALUE

Output:

Why are there multiple rows (6) for each value?

I can load this raster into QGIS and it displays fine. I tried exporting to my local hard drive as a geotiff (from QGIS), then reloaded into PostgreSQL using raster2pgsql in the command prompt, and ran the same query to calculate pixel number and proportion (with the imported geotiff raster), and it calculates as expected:

When looking at the Info from QGIS DB Manager, the original raster (the one that doesn't behave as expected) has a few differences from the imported raster (the one that works):

Pages: 3 vs 1
Rows: 128 vs 1
Spatial ref: [none] vs "WGS 84 / UTM zone 15N (32615)"
Constraints: [none] vs many (e.g., enforce_max_extent_rast, enforce_out_db_rast, etc.)

Version info:

SELECT version()

PostgreSQL 9.5.10, compiled by Visual C++ build 1800, 64-bit

SELECT PostGIS_full_version()

POSTGIS="2.3.5 r16110" GEOS="3.6.2-CAPI-1.10.2 4d2925d" PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 2.2.2, released 2017/09/15" LIBXML="2.7.8" LIBJSON="0.12" RASTER


texts: What is the order of parameter (lat/long) in python UTM package for conversion into x, y coordinate?

texts:
Currently, I am trying to convert (lat, long) into (x, y) coordinates in Python. There are some packages like pyproj or utm to do that. For example:

import utm
lat = 41.8836718276551
lon = 12.4877775603346
x, y, zn, zl = utm.from_latlon(float(lat), float(lon))

My question is that what is the order of the parameter? is it (lat, lon) or (lon, lat)? Even If I use pyproj to do the same conversion, then what will be the order?


texts: Prior art on rhumb line and great ellipse areas?

texts:
I'm adding the computation of the area betweeen a rhumb line or a great
ellipse and the equator to GeographicLib.  My formulations are
briefly given in here and here.  These formulations are
similar to the one used by GeographicLib in computation geodesic areas.
With these basic areas known, it will be possible to compute the area of
a polygon whose sides are a mixture of geodesics, rhumb lines, or great
ellipses.

Is there any "prior art" I should reference?  I'm particularly
interested in methods which apply to the ellipsoid and can deal with
long segments without slicing the area into many strips.  How do the
commercial tools, ArcGIS and Blue Marble, solve these problems?


texts: Using Overlay Analysis (e.g., Union) with point data (point shapefile)?

texts:
Is it possible to compare a point shapefile with a DEM raster in which the Aspect tool has been applied?

I am reading about Overlay Analysis but cannot find a definitive answer.

The attached photo is an example of my point dataset (enlarged black dots) shown in front of the study area with Aspect applied. I would like to see how the shapefile is distributed by aspect (e.g., How many points are on a south-facing aspect, how many on a southwest-facing aspect, etc.). With over 200 points, I am trying to avoid doing this manually. I am using ArcMap 10.3.


texts: Rename outputs from QGIS processing saga algorithms

texts:
I have a processing script which chains together several saga processing tools like this:

outputs_SAGARASTERCALCULATOR_3=processing.runalg('saga:rastercalculator', inputvx,[inputvy],'(atan2(a,b)) * (180/pi())',True,7,None)
outputs_SAGAGRIDVALUESTOPOINTS_1=processing.runalg('saga:gridvaluestopoints', [outputs_SAGARASTERCALCULATOR_3['RESULT'],PolygonExtent,None,True,0,OutputName)

The calculator process defaults to calling the raster calculator output'RESULT', which means it ends up using the same interim file when you run a batch process - they are all called result regardless of the inputs.

How do I change the name of the calculation result to be something other than 'result'?

Here's the python console for the calculator tool - I assume I can alter both 'Result:' and 'Name:' but it won't seem to allow me to add any other parameters into the tool.

        Grid system: 10; 2601x 1500y; 383472.948554x 557408.826780y
        Grids: 1 object (Tyne10mBaselineUpdated02MainFlood0100Vx)
        Grids from different Systems: 1 object (Tyne10mBaselineUpdated02MainFlood0100Vy)
        Result: Result
        Formula: (atan2(a,b)) * (180/pi())
        Name: Calculation
        Take Formula: no
        Use NoData: yes
        Data Type: 4 byte floating point number


texts: QGIS update from 2.14 to 2.16 failed on Ubuntu 14.04

texts:
Until yesterday my whole qgis and python-qgis worked fine with version 2.14 under Ubuntu 14.04.
After a simple sudo apt-get update and sudo apt-get upgrade qgis failed to update... Several attempts to reinstall it with different versions failed.

Current error message (for version 2.14.3) goes like this:

qgis-providers (2.14.3+dfsg-2~trusty1) is configured ...
/usr/lib/qgis/crssync: symbol lookup error: /usr/lib/libgdal.so.20: undefined symbol: GEOSDelaunayTriangulation_r

I also tried to reinstall all gdal libraries - not sure if I got all?!

Maybe it is somehow related to this one: QGIS 2.14 fails to install on Ubuntu 14.04 ?

UPDATE:

Trying to install version 2.16.0 it gives different but similar error:

/usr/lib/qgis/crssync: symbol lookup error: /usr/lib/libqgis_core.so.2.16.0: undefined symbol: GEOSNearestPoints_r

UPDATE Sources:

For version 2.16.0 there are these sources in /etc/apt/sources.list:

deb http://qgis.org/debian trusty main
deb-src http://qgis.org/debian trusty main
deb http://qgis.org/debian trusty main
deb-src http://qgis.org/debian trusty main

which are all the same, right?

When I do apt-cache policy qgis it gives me:

qgis:
  Installiert:           1:2.16.0+20trusty
  Installationskandidat: 1:2.16.0+20trusty
  Versionstabelle:
 *** 1:2.16.0+20trusty 0
        500 http://qgis.org/debian/ trusty/main amd64 Packages
        100 /var/lib/dpkg/status
     2.0.1-2build2 0
        500 http://de.archive.ubuntu.com/ubuntu/ trusty/universe amd64 Packages


texts:
I get results using pgRouting I cannot explain.
I created a topology  using osm2pgrouting. I then want to find the shortest way through this topology, for someone walking on foot.
To calculate the costs I wrote a function. This function works, sort of.

CREATE OR REPLACE FUNCTION costs_foot(valueHighway text, costs double
RETURNS double precision AS
$$
SELECT CASE
    WHEN $1 IN ('motorway', 'motorway_link', 'motorway_junction', 'trunk', 'trunk_link', 'primary', 'primary_link', 'secondary', 'tertiary', 'bus_guideway', 'cycleway', 'byway', 'secondary_link', 'tertiary_link', 'lane', 'opposite_lane', 'opposite', 'road', 'unclassified', 'grade1', 'grade2', 'grade3', 'grade4', 'grade5')
    THEN -1 -- lines in question
    ELSE costs -- lines in question
    END
$$
language 'sql';

And then I calculate the route

CREATE TABLE route_foot AS
SELECT route.*, w.the_geom, w.length_m FROM pgr_dijkstra('
    SELECT gid AS id,
         source,
         target,
     costs_foot(str_name, cost_s) AS cost,
         reverse_cost_s AS reverse_cost
    FROM ways',
    pgr_pointToEdgeNode('ways', ST_SetSRID(
                    ST_Point(-73.930397, 40.783351),
                    4326), 0.01
            ),
    pgr_pointToEdgeNode('ways', ST_SetSRID(
                    ST_Point(-73.882022, 40.852214),
                    4326), 0.01),
    directed := false) AS route
LEFT JOIN ways w
ON route.edge = w.gid
ORDER BY seq;

That gives me route which goes over highways, motorways, ... basically everything I am trying to avoid with -1.

However, when I change the function to this:

CREATE OR REPLACE FUNCTION costs_foot(valueHighway text, costs double
RETURNS double precision AS
$$
SELECT CASE
    WHEN $1 IN ('motorway', 'motorway_link', 'motorway_junction', 'trunk', 'trunk_link', 'primary', 'primary_link', 'secondary', 'tertiary', 'bus_guideway', 'cycleway', 'byway', 'secondary_link', 'tertiary_link', 'lane', 'opposite_lane', 'opposite', 'road', 'unclassified', 'grade1', 'grade2', 'grade3', 'grade4', 'grade5')
    THEN costs -- lines in question
    ELSE -1 -- lines in question
    END
$$
language 'sql';

I get a route that only choses footways, pedestrians, ...
Looking at the logic in the SQL Function it should be the other way around.

Why is that?


texts: Entry point not found Error initiating GRASS GIS from Python

texts:
I cannot run any GRASS tools from within Python 2.7.15. While the GRASS environment, grass.script etc. are loading correctly and I can run basic checks like

gscript.read_command('g.gisenv')
gscript.read_command('g.mapset', flags='p')

but receiving this error (screenshot) every time I tried to run tools like

gscript.read_command('g.proj', flags='p')
gscript.read_command('g.region', flags='p')

Everything was running smoothly before. If you want to check the startup script that I'm using, please see here). Not duplicating here as already posted earlier. This problem started after I fresh-installed Anaconda 2 with Python 2.7.15 to fix some broken Python libraries.

Already tried fresh installing GRASS 7.4, restarting the machine after each install and uninstall. But not fixing this issue. Need to run the previously working scripts desperately. Any thoughts?


texts: How does the projection of a dataset affect the results of spatial analysis?

texts:
I am planning on conducting a multi-criteria-analysis with a number of different source datasets. A number of the datasets are found in different coordinate systems and projections. Do I need to project these datasets into the same coordinate system before conducting this analysis? My data currently covers more than one UTM Zone, usually they fall within one. Or should all of the data be in the same coordinate system, but unprojected?

Since there is such a large amount of data I would like it to be as efficient as possible.

Although I am utilizing Esri tools for this problem my question is not specific to those. I am working with a dataset which covers 650,000 square kilometres of the NWT and NU in Canada. I've decided to use spatialreference.org/ref/epsg/nad83-nwt-lambert for all of the data. I was wondering if it would be better to leave the data unprojected or pick a projection for all of the data to be in.


texts: Combining/getting elevations from .asc files or .tiff with x and y locations

texts:
I am new to R and still getting to know the environment.

I am doing a project in animal movement. I have x and y locations (plus IDs, dates, hours) in an Excel table, converted into .csv as well. I downloaded the data from LiDAR and another web page for the area the animals were moving around at, so I got .asc and .tiff files. Now I need to combine the x and y locations with either of these files to get an elevation for each x and y (elevation will be one of the environmental covariates).

I tried "elevatr" package/guidance but it's only for US and it only shows how to download it from the internet.

Has anyone any idea how to merge or combine these files and get elevations for each location?

I believe there's a very simple way of doing it but as I said, I'm new to this coding and wandering in the dark most of the time.


texts: Is there a Swipe tool for ArcGIS JavaScript API?

texts:
I was researching the ESRI JavaScript API, and was wondering if there is a sample which shows layer swiping in the browser.

Is there such a Swipe sample?