texts: Is there a measure that account for the the location diversity and trajectory distance of agents?

texts:
I am currently working on project to assess the success of musicians. For each musician I know when and where (coordinates,city,county) they performed in a gig. I assume a successful musician is the one who travel get frequent bookings in diverse countries. I like to account for the number of gigs, and number of cities and countries. For example, the musician A only performs only in US might have a longer travel distance, compare to the musician B who performed in multiple countries in Europe.

So far I thought about computing the total travel distance by summing up the distance between the cities for consecutive gigs. This doesn't account for the location diversity as one artists can move frequently between few cities.

I also thought to create a polygon of a musician trajectory and use the "area" of the polygon as the success measure. However, it seems duplicated cities create a holes in the polygon and reduce the area. For example in the following example I was hoping geom1 (3.05402) have a higher area than geom2 (8.695). When I remove the the second visit to Cologne in the geom1 the area changes to 11.87. If I remove the duplicated from each trajectory I can not control for the number of gigs any more.

geom1 = Polygon([city_points_obj['cologne'],city_points_obj['berlin'],city_points_obj['paris'],city_points_obj['cologne'],city_points_obj['milan']]) #area: 3.05402

geom2 = Polygon([city_points_obj['cologne'],city_points_obj['berlin'],city_points_obj['stuttgart'],city_points_obj['cologne']]) #area: 8.6952

Since I have a very limited (almost none) education on in GIS analysis, I though there might be some other measures that I'm not aware of, or get some suggestions from experts.


texts: Using internal QGIS icons / banner in own plugin

texts:
I have a problem using the internal symbols/icons in my plugin. For example I want to use the "save" Symbol in a QToolButton or the "cancel" symbol in a QPushButton.

I designed the GUI in QT 5 Designer. I also added the symbols to the button. After I compiled the ui file with pyuic5 I changed the absolute path to the symbols to the relative path in my plugin folder after putting the downloaded svg files (like its discribed here and there) in the plugin directory.

Running the file QgsAlkOsmCompare_dialog.py in the terminal, produces this:
All icons show up in their correct position (I cutted the screenshot because of its size).

Opening my plugin in QGIS is functional, but the icons disappear. Even if I run the same py file in the  ./QGIS/QGIS3/profiles/default/python/plugins folder the plugin doesn't show the icons. I tried it with png file instead but nothing changed.

Is it impossible to show buttons with icons in QGIS? I guess no, because I have seen icons in some of the GDAL processing functions too.

Can I import these internal QGIS icons in my GUI file or is it compiled to the qgis.exe/qgis or some other kind of resource file?


texts: Switching between python versions for Arcpy

texts:
I’m developing a multi-module project in Python. I’m new to GIS and my boss doesn’t have a strong technical background. I’ve been instructed to use the Arcpy module to do functions such as intersect. The project requirements state I need to use Python 3.3.
Is there a contradiction between using the Arcpy module and Python 3.3?
I’m using Liclipse (Eclipse styled for Python) and in one module use the library lxml and in another use Arcpy. In properties, if I put the python C:\Python33\python.exe on top then lxml works but Arcpy doesn’t and if I put pythonArcGis C:\Python27\...  on top then Arcpy works and lxml doesn’t. I don’t really understand what’s going on…

Questions:
Can Arcpy work with Python 3.3?
If yes, how do you get both Arcpy and lxml to work at the same time?
If no, then is it possible to launch only the module that uses Arcpy as Python 2.7?


texts: What is the proper way to calculate distance in km from degrees

texts:
For a genetic analysis of Australian humpback dolphins I am currently trying to assess whether or not genetic differentiation correlates with geographic distance in water. For this I have done the following:

Get bathymetric data using the R package marmap

getNOAA.bathy(lon1 = 145, lon2=154, lat1 = -28, lat2 = -17, resolution = 1)

Convert it to a raster and subsequently to a transition matrix with low cost for moving in water and prohibitively high cost for traversing land

Get the shortest distance in water between each pair of points using gdistance::shortestPath, and converting this into a distance matrix

I am now faced with the following: The lengths I am getting from shortestPath are in the same dimension as the coordinate system I am using, which is decimal degrees. So my question is: What is the proper way to convert this distance matrix into km?
I understand these degrees must be projected somehow but I am currently failing at figuring out how exactly I can do this.
my centroids are:
data.frame(lon = c(146.1333, 146.8698, 148.7083, 150.8709, 151.3147, 152.9791, 153.0279, 153.1816), lat = c(-18.33419, -19.22425, -20.55407, -23.43784, -23.84044, -25.26728, -25.81851, -27.35664))

and the corresponding distance matrix is:
1.3629804
3.9769987  2.7628686
7.7416532  6.5480830  4.0348631
8.1950482  6.9809181  4.4882582  0.8048330
10.4037986  9.1896685  6.6970086  3.0371065  2.3537371
10.9785875  9.7644574  7.2650158  3.6118972  2.9285279  0.6491764
12.7441285 11.5573399  9.0441220  5.3774364  4.7008506  2.4420589  1.7487997


texts: Creating a subset of XY data by most recent date

texts:
I have a large set of well test data that I would like to pull just the most recent test data from. The point data has several years of tests in it, at co-incident points. All of the overlapping points have identical x & y coordinates.I was wondering how best to delete the older test data for each point, leaving me with just one point for each xy.

The date of test is based off of 2 fields [SAMP_DATE] (which is data type "DATE") and [SAMP_TIME] (which is data type "String"). There are multiple samples taken on the same date (at different times), so I will need to use both fields.

Hoping someone can help me figure this out via a repeatable definition query, as I have to pull out the same data for different sample tests (eg. Arsenic, Nitrate, etc.).

Here is what my attribute table looks like, if it helps:

I attempted Using the sort & delete identical method described in the comments below by @Bera & @ FelixIP, but did not achieve the result I wanted. Here is a photo of the results using that method (for one point) The original table is on right & results of Sort/Delete identical processing on the left.


texts: Popup in Sidebar With Hosted Feature Layer, not ArcGIS Web Map

texts:
I see this example of how to put the contents of a popup into a sidebar, but I can't really follow it. I'm not sure how to implement it for my site, which uses feature layers hosted on ArcGIS Online, rather than a web map like the example.

It seems like the general idea is to turn off the popup, listen for selection changes, and then grab the infoWindow content and put it into a specified DIV. Am I right here?

How can I do this without using a web map?


texts: How to delete a map from the GRASS interface?

texts: I want to delete a map that I generated from GRASS. I tried removing a map from the "Map Layers for each display" tab in the GRASS GIS Layer Manager, but it only removed the map from the list and not deleted. I want to delete the map totally that it won't appear as an option in when I load the maps in the workspace (Workspace>Load Map Layers). How is it done?


texts: Getting values from multiple rasters that overlap polygon using ArcGIS Desktop

texts:
I have a series of rasters from consecutive years that represent different levels of fire severity. Some of these rasters have pixels that overlap each other in the different years. I have a set of polygons that represent specific species' habitat and I need to see what percent of each polygon is affected by different levels of fire severity.

For example, my rasters have values ranging from around -11000 to 3000 and I only want to see how much of 115-308, 308-643, and 643+ are in each polygon. I would like to just get rid of everything below 115 but the Reclassify tool is not updating the attribute field like I though it was supposed to, but instead creates a new TIFF with only 4 values and rows.

How can I do this?

I'm using ArcGIS Desktop 10.7


texts: MySQL Query to fetch all cities and POI's near the route

texts:
I'm implementing alghorithm to search POI near the route (array of Lat, Long points).
My problem is that I don't know how to calculate distance beetween line segment and POI.
Look at picture that I attached. Can somebody give me some formula for fast calculation this distance?

Notes:
Using PHP and MySQL with Spatial Extension.

Look at picture

Assumptions:
- I have database of POI with lat, long -
- I have array of points (lat, long) in correct (drive) order - route
- The calculations will be made on max 20 km area so flat surface assumption is ok.
- The POI must be returned in order that the driver will be passing the route.

Any ideas?


texts: Accelerating the code for random selection of polygons

texts:
I have a shapefile with thousands of features (polygons). I need to randomly sample one polygon at every 1 sq. km. area if at least one polygon is available. I wrote the following code which works. However, it is very slow (took few days to complete).

import numpy as np, pandas as pd, geopandas as gpd
from shapely.geometry import Polygon

fi = 'test.shp'
gdf = gpd.read_file(fi)[['geometry']]
print (gdf.shape)
crs = gdf.crs

ULX, LRY, LRX, ULY = gdf.total_bounds
ts = 0.008983157 #1km
lats, lons = np.mgrid[ULY:LRY-ts:-ts, ULX:LRX+ts:ts]

ulx_s = lons[:-1,:-1].ravel()
uly_s = lats[:-1,:-1].ravel()
lrx_s = lons[1:,1:].ravel()
lry_s = lats[1:,1:].ravel()
tiles = [[(ulx,uly),(ulx,lry),(lrx,lry),(lrx,uly)] for ulx, uly, lrx, lry in zip(ulx_s, uly_s, lrx_s, lry_s)]
print(len(tiles))

#### this is the slowest part ###########
outputs = []
for tile in tiles:
    poly = Polygon(tile)
    ok = gdf[gdf.geometry.intersects(poly)]
    if ok.shape[0] >= 1:
        out = ok.sample(1)
        outputs.append(out)
###########

out_df = pd.concat(outputs, axis=0)
out_gdf = gpd.GeoDataFrame(out_df, geometry='geometry')
out_gdf.crs = crs

out_gdf.to_file(fo, layer='test', driver='ESRI Shapefile')

How can I do it faster?


texts: Which cell size should I use in the LecoS plugin?

texts: LecoS asks for the cellsize when calculating landscape statistics. However, there is no such information in the metadata of my input raster layer. Pixel size is 0.00126545,-0.00126497, but that can't be the required information? What am I doing wrong?


texts: Could not find Easy custom labeling tool in plugin manager

texts:
I could not find the tool "Easy custom labeling" in Plugin Manager.

My colleague can install the tool for QGIS in her PC but she uses Windows 7 while I use windows 10. Could anyone advise me how to get that tool installed for windows 10?


texts: How do I select shapefiles to be converted to GeoJSON in folder with multiple shapefiles?

texts:
How do I only select .shp's to read and then convert to .geojson. The homefolder has 3 different shapefiles and this code currently converts every single file in the folder to a .json. Just trying to speed up the running process. I have some other issues too with naming conventions that I will deal with next. Among other things, this is a work in progress.

homefolder = "D:/UW/Winter_2016/501/pythonProject/publicLandsData/newShapefiles/"

# read the shapefiles
filenames = os.listdir(homefolder)

for filename in filenames:
  reader = shapefile.Reader(homefolder + os.path.splitext(os.path.basename(filename))[0])
  print filename
  fields = reader.fields[1:]
  field_names = [field[0] for field in fields]
  buffer = []
  print filename + "buffer dictionary done"
  for sr in reader.shapeRecords():
      atr = dict(zip(field_names, sr.record))
      geom = sr.shape.__geo_interface__
      buffer.append(dict(type="Feature", \
      geometry=geom, properties=atr))

# write the GeoJSON file
#currently it saves into the wrong directory..and all one shapefile, might be nice to be in multiple geojson's
#list directory "homefolder" in variable filenames
  filenames = os.listdir(homefolder)

  for filename in filenames:
#for loop to rename local government to loca
    geojson = open(((os.path.splitext(os.path.basename(filename))[0]) + ".json"), "w")
    #os.rename(os.path.join(homefolder, filename), os.path.join(homefolder, (filename.replace('.shp', 'shp'))))
    geojson.write(dumps({"type": "FeatureCollection",\
      " features": buffer}, indent=2) + "\n")
    geojson.close()

I have tried using basename, selecting only ".shp" files with .endwith() etc to no avail. Suggestions?

I found this shp to geojson conversion code at the link below:
https://gist.github.com/frankrowe/6071443


texts: How to use POST method with Overpass QL and OpenLayers?

texts: I am using OpenLayers and Overpass QL with GET method to retrieve data from the Overpass server. Some URLs tend to be very long and hence get rejected by the Overpass server. To get around this, I am thinking of using POST method to query the Overpass server using OpenLayers. Is there a way to use POST method with Overpass QL (not XML) in OpenLayers? And if yes, what (syntax) goes into the url property of OpenLayers.Layer.Vector?


texts: How to convert gps receiver native format to google lat long or XY format

texts:
I have received some GPS co-ordinates in the following format like - "17267690N 078279340E"..

How to convert the co-ordinates to show it on Google map.


texts: ArcGIS REST service is returning error when distinct param is supplied

texts:
ArcGIS REST service (10.4.1) is returning the error

Unable to complete operation.
Unable to perform query operation.

When ReturnDistinct=true and InputGeometry is supplied.
This operation works correctly if ReturnDistinct=false

Full query params as below

query?where=
&objectIds=
&time=
&geometry=%7B"xmin"%3A16712664.599773334%2C"ymin"%3A-2674043.1889597364%2C"xmax"%3A16772285.481835691%2C"ymax"%3A-2639761.1817738814%2C"spatialReference"%3A%7B"wkid"%3A102100%7D%2C"cache"%3A%7B"_parts"%3A%5B%7B"extent"%3A%7B"xmin"%3A16712664.599773334%2C"ymin"%3A-2674043.1889597364%2C"xmax"%3A16772285.481835691%2C"ymax"%3A-2639761.1817738814%2C"spatialReference"%3A%7B"wkid"%3A102100%7D%7D%2C"frameIds"%3A%5B0%5D%7D%5D%7D%7D
&geometryType=esriGeometryEnvelope
&inSR=
&spatialRel=esriSpatialRelIntersects
&distance=
&units=esriSRUnit_Foot
&relationParam=
&outFields=
&returnGeometry=false
&maxAllowableOffset=
&geometryPrecision=
&outSR=
&gdbVersion=
&returnDistinctValues=true
&returnIdsOnly=false
&returnCountOnly=false
&returnExtentOnly=false
&orderByFields=
&groupByFieldsForStatistics=
&outStatistics=
&returnZ=false
&returnM=false
&multipatchOption=
&resultOffset=
&resultRecordCount=
&f=html


texts: extending existing line in QGIS by adding nodes

texts: I am trying to extend an existing line by adding nodes to it.  Is this possible to do in QGIS without creating a whole new line?


texts: Autocompletion for VS and Eclipse when debugging QGIS plugins

texts:
I have been trying to get autocompletion to work in first Eclipse and afterwards tried with Visual Studio. I have been following these guidelines / Stack Questions

IDE settings for writing and debugging plugins

import qgis.core problem with QgsFeature and QgsGeometry ImportError: DLL load failed: The specified procedure could not be found

Setting up a QGIS Plugin Development Environment in Visual Studio and PTVS

In both environments I managed to be able to debug the python code step by step, but iI was not able to get any code completion.

I suspect that its the same issue in both cases. I have been trying to add all possible paths to the PYTHONPATH. Both in Eclipse and in windows environment variables.

According to the guideline I did open Eclipse using a bat script to init windows specific things for QGIS.

call "C:\Program Files\QGIS 2.18\bin\o4w_env.bat"
set PATH=%PATH%;C:\Program Files\QGIS
2.18\apps\qgis\bin;C:\ProgramData\Oracle\Java\javapath
C:\software\eclipse\java-oxygen\eclipse\eclipse.exe

I made a similar script for VS without success.

The configuration I have been setting up for the python env. looks like the following.

What could I be missing that causes Eclipse not to recognize libs such as qgis.core?

Could it be the same reason that causes VS to not recognize it?


texts: Which projection to use to map three Canadian provinces?

texts: I was making a map of Alberta, Saskatchewan, and Alberta, thinking I would use a UTM projection, however, these provinces cover zones 12-16, so this likely isn't ideal. What other options are there?


texts: Is there a way on Atlas to put a Margin around feature but also use Predefined scales?

texts:
In QGIS, how can I put a margin around a feature but also use predefined scales in Atlas? Like an in-between of the two options. I have a feature that fits in the map but as you can see in the image it almost reach the edges leaving some annotation out of the viewmap, so I would like QGIS to use a smaller scale to better accommodate the drawing.

I need it to be the same coverage layer with the features - not another layer with a buffer - due to the fact that I have a symbology that depends on it, like the expression that I use to show only the features that touch the atlas feature in more than two vertices.
Tried to create a buffer with "Geometry Generator" around the feature thinking the program would understand that I also want it to appear and so shrink the drawing, but it didn't work. Is there a way to do that?


texts: Getting coordinates in Projected Coordinate System using ArcGIS API for JavaScript?

texts:
I have map built using the ArcGIS API for Javascript 4.8. The basemap is the ESRI "streets" basemap in wkid 102100 (latestWkid 3857).

How do I get the coordinates in the following to be output as wkid 2810 (NAD_1983_HARN_StatePlane_Minnesota_North)?

Here's the code:

mapView.on("pointer-move", function(event){
    console.log(mapView.toMap(event));
});


texts: Cokriging: zero distance cross-semivariance in undersampled cases [gstat error]

texts:
Trying to cokrige two variables that are not perfectly colocated (one has sparser measurements than the other), I faced an issue that I'll illustrate with the following MRE.

With the meuse dataset, consider we krige lead as principal variable assisted by copper measurements as auxiliary variable. We subset the lead, to simulate sparser measurements and hence, a need for cokriging with the denser copper measurements:

library("gstat")
library("sp")
data("meuse")
coordinates(meuse) = ~x+y

g <- gstat(NULL, data=meuse[1:80,], formula=lead ~ 1) # subset
g <- gstat(g, data=meuse, formula=copper ~ 1)
v <- variogram(g)
plot(v) # zero distance semivariance dot in panel var1.var2
g <- fit.lmc(v=v, g, vgm("Sph")) # error

  Error in fit.variogram(x, m, fit.ranges = fit.ranges, ...) :
    fit.method 7 will not work with zero distance semivariances; use another fit.method value

You can see in the bottom-left panel of the plot a zero distance dot, which subsequently causes fit.lmc() to fail.

Now, if no subset is done (somewhat limiting the interest of cokriging, no?), everything works fine. This because, the zero distance dot in the cross-semivariogram does not appear in this case:

g <- gstat(NULL, data=meuse, formula=lead ~ 1)
g <- gstat(g, data=meuse, formula=copper ~ 1)
v <- variogram(g)
plot(v) # no zero dist dot on var1.var2 panel
g <- fit.lmc(v=v, g, vgm("Sph")) # fit fine with default fit.method

There is no reason for this dot to appear only in the subsetted cases, right?

Another example of this can be seen in this 2009 exercise document by Edzer Pebesma. In section 8.13, he refers to this as the "undersampled case". But the provided code is not working any more, likely for the reason mentioned above.

Is there is a simple way around this (hopefully) temporary bug?

PS: I moved this from Cross-validated because it is software-related.


texts: Merging multiple rasters (Very large dataset)

texts:
I am working on project to create a DSM raster for whole Kentucky. I will have over 46000 tiff tiles and each one is about 10 Mb and one band. Normally, I use Mosaic to New Raster Tool but this time data is very large. I am planning to divide the project into smaller pieces but I have no idea what size should I set the each piece.
I will use python code with ArcGis Desktop 10.6 to complete project. Considering my computer has 2.7 Ghz (i7) processor with 16 Gb RAM, How many tiles I need to use at one time? like 1000? 10000? or what number?


texts:
So I'm in the middle of learning PostGIS/PostgreSQL.

I've a polygon table ans I've a point table. The polygon layer contains a foreign key that points to another table. For the example the Foreign key points to the type of polygon. So TypeFK in the polygon table  references TypeID in the type_tbl;

I can do a:

select type
from polygon as a, type_tbl as b
where a.typefk = b.typeid;

Now I know the polygon is a type of tree.

I can do a:

select typefk
from polygon
join points on ST_Contains(polygon.geom, points.geom);

I know what typefk they fall in.

I'm ultimately trying to return the type of polygon the points are located within.

So - gurus of PostGIS/SQL - how do I combine the two?

I've been sitting here too long without a break unfortunately and am just not getting it.


texts: Is the PostGIS st_intersection function viable for very large data-sets?

texts:
I am running an area intersection query on two large multipolygon tables (a few million rows).

SELECT r.geom
FROM public.my_table r
WHERE EXISTS (
    SELECT gid
    FROM public.my_table_intersection
    WHERE
        r.class = '2' AND
        st_intersects(r.geom, geom) AND
        ST_IsValid(r.geom) AND
        ST_IsValid(geom) AND
        (st_area(st_intersection(r.geom, geom))/st_area(r.geom)) > 0.1
);

The code appears to work with small test files but I have left it for ages with my main files with no result yet. I know that such processing on very large datasets can be very slow, is this sort of operation within PostGis's capabilities or am I pushing it too far? Is there perhaps an alternative way of doing such processing on very large data-sets? would ArcGIS be quicker?

EDIT:
The table row size is 3.7 million for my_table and then 2.1 million for my_table_intersection with indexes on geom for both.

The query plan is:

"Nested Loop Semi Join  (cost=0.41..50013102.40 rows=33 width=275)"

"  ->  Seq Scan on my_table r  (cost=0.00..1153671.74 rows=295608 width=275)"

"        Filter: ((modal_clas = 3::numeric) AND st_isvalid(geom))"

"  ->  Index Scan using my_table_intersection_geom_gist on my_table_intersection  (cost=0.41..165.22 rows=8 width=277)"

"        Index Cond: (r.geom && geom)"

"        Filter: (st_isvalid(geom) AND _st_intersects(r.geom, geom) AND ((st_area(st_intersection(r.geom, geom)) / st_area(r.geom)) > 0.1::double precision))"


texts:
After launching the GeoServer, it prepares for map display for 60 seconds and processes data(layers, image pyramids).  Further, the processing of requests is fast.If no requests are received from the client to geoserver for a long time, the situation with long processing is repeated.

What does the "org.geotools.gce.imagemosaic] - Checking" and for what?

2020-02-17 09:36:05,056 DEBUG [org.geotools.gce.imagemosaic] - Checking file:/usr/share/geoserver/bin/../data_dir/data/yamap/z18/O-36/O-36_2-5/0/
isHidden:false
exists:false
isFilefalse
canRead:false
canWritefalse
canExecute:false
isAbsolute:true
lastModified:0
length:0
2020-02-17 09:36:05,056 DEBUG [org.geotools.gce.imagemosaic] - Checking file:/usr/share/geoserver/bin/../data_dir/data/yamap/z18/O-36/O-36_2-5/0/
isHidden:false
exists:false
isFilefalse
canRead:false
canWritefalse
canExecute:false
isAbsolute:true
lastModified:0
length:0
2020-02-17 09:36:05,057 DEBUG [org.geotools.gce.imagemosaic] - Checking file:/usr/share/geoserver/bin/../data_dir/data/yamap/z18/O-36/O-36_2-5/0/
isHidden:false
exists:false
isFilefalse
canRead:false
canWritefalse
canExecute:false
isAbsolute:true
lastModified:0
length:0
2020-02-17 09:36:05,057 DEBUG [org.geotools.gce.imagemosaic] - Checking file:/usr/share/geoserver/bin/../data_dir/data/yamap/z18/O-36/O-36_2-5/0/
isHidden:false
exists:false
isFilefalse
canRead:false
canWritefalse
canExecute:false
isAbsolute:true
lastModified:0
length:0
2020-02-17 09:36:05,057 DEBUG [org.geotools.gce.imagemosaic] - Checking file:/usr/share/geoserver/bin/../data_dir/data/yamap/z18/O-36/O-36_2-5/0/
isHidden:false
exists:true
isFiletrue
canRead:true
canWritetrue
canExecute:false
isAbsolute:true
lastModified:1580307269992
length:367708
2020-02-17 09:36:05,057 DEBUG [org.geotools.gce.imagemosaic] - Checking file:/usr/share/geoserver/bin/../data_dir/data/yamap/z18/O-36/O-36_2-5/0/
isHidden:false
exists:true
isFiletrue
canRead:true
canWritetrue
canExecute:false
isAbsolute:true
lastModified:1580307269992
length:367708
2020-02-17 09:36:05,057 DEBUG [org.geotools.gce.imagemosaic] - Trying to load properties file from URL:file:/usr/share/geoserver/bin/../data_dir/data/yamap/z18/O-36/O-36_2-5/0/0.shp
2020-02-17 09:36:05,057 DEBUG [org.geotools.gce.imagemosaic] - Trying to load properties file from URL:file:/usr/share/geoserver/bin/../data_dir/data/yamap/z18/O-36/O-36_2-5/0/0.shp
2020-02-17 09:36:05,057 DEBUG [org.geotools.gce.imagemosaic] - Checking file:/usr/share/geoserver/bin/../data_dir/data/yamap/z18/O-36/O-36_2-5/0/
isHidden:false
exists:false
isFilefalse
canRead:false
canWritefalse
canExecute:false
isAbsolute:true
lastModified:0
length:0
2020-02-17 09:36:05,057 DEBUG [org.geotools.gce.imagemosaic] - Checking file:/usr/share/geoserver/bin/../data_dir/data/yamap/z18/O-36/O-36_2-5/0/
isHidden:false
exists:false
isFilefalse
canRead:false
canWritefalse
canExecute:false
isAbsolute:true
lastModified:0
length:0
2020-02-17 09:36:05,057 DEBUG [org.geotools.gce.imagemosaic] - Checking file:/usr/share/geoserver/bin/../data_dir/data/yamap/z18/O-36/O-36_2-5/0/
isHidden:false
exists:false
isFilefalse
canRead:false
canWritefalse
canExecute:false
isAbsolute:true
lastModified:0
length:0
2020-02-17 09:36:05,058 DEBUG [org.geotools.gce.imagemosaic] - Checking file:/usr/share/geoserver/bin/../data_dir/data/yamap/z18/O-36/O-36_2-5/0/
isHidden:false
exists:false
isFilefalse
canRead:false
canWritefalse
canExecute:false
isAbsolute:true
lastModified:0
length:0


texts: Setting CRS for comma delimited layer using QGIS

texts:
I am using QGIS 3x.
I would like to overlay my .csv comma delimited layers in QGIS 3x properly. I have got the Ordnance Survey Maps as the raster players, which CRS is EPSG 27700 1936 British National Grid.
I overlaid these maps with my .csv data, where I set EPSG 4326: WGS 84. As a result, the data is visible, although there are some mismatches, as you can see below:
<img src="https://i.stack.imgur.com/UA8w9.jpg" alt="enter i
mage description here" />
In the following example, you can see, that the difference is quite big, as the pole appears on the other side of the road effectively. I theoretically know, that the mismatch is driven by different CRS between the comma-delimited layers (WGS 84) and the map layer (British National Grid).
When I change the CRS for my comma-delimited layers, then they don't appear in my location anymore. Their new place is somewhat at the beginning of the coordinate system, as you can see below:

I would like to make sure, that all these layers match to each other as much as possible.
What should I do then?


texts: How to specify shapefile variables length in R?

texts:
I am using R to edit some shapefiles with readOGR and writeOGR functions. I can edit the name and type of my shapefile variables, but not the length.

Example: i want that my variable of type string takes a length of 20. For now it takes the default length of 80 (i visualize the shapefile under QGIS).

Any hint ?


texts: How to calculate if the location of geomorphological locations are likely to be in the proximity of geological faults

texts:
I have two shapefile data. One is point data with 800+ observations of certain geomorphologic features in karst and other is a line shapefile which represents faults (geology). I want to know if there is any relationship between those two, or more specifically I would like to know if there are significantly more points near the faults.

Do you know of any functions that would help me with this?


texts: Counting points inside polygons using PostGIS?

texts:
I have 2 vector data in my Postgis DB: one is a point layer (punti) with about 200000 elements and the other is a polygon (province) with 4 elements.
They cover the same area and I need to find the number of points for each polygon.

This is the query:

province.gid, count(*) AS totale FROM punti, province WHERE
st_intersects(province.the_geom,punti.the_geom) GROUP BY province.gid;

This is the EXPLAIN ANALYZE

  HashAggregate  (cost=35.35..35.36 rows=1 width=4) (actual time=23660.787..23660.789 rows=4 loops=1)
     ->  Nested Loop  (cost=0.00..35.34 rows=1 width=4) (actual time=0.307..23532.479 rows=199977 loops=1)
           Join Filter: _st_intersects(province.the_geom, punti.the_geom)
           ->  Seq Scan on province  (cost=0.00..1.04 rows=4 width=36) (actual time=0.003..0.009 rows=4 loops=1)
           ->  Index Scan using punti_the_geom_idx on punti  (cost=0.00..8.31 rows=1 width=100) (actual time=0.038..76.254 rows=69541 loops=4)
                 Index Cond: (province.the_geom && punti.the_geom)
   Total runtime: 23660.850 ms

What do you think about this total runtime?


texts: Using GeoJSON in Openlayers 4 and Geoserver

texts:
First time trying to draw an Openlayers map with GeoJSON instead of WMS, and not having much luck. My layer is hosted in Geoserver and is named 'sections', in a workspace called 'basic_map'. This is my code so far:

var vectorLayer = new ol.layer.Vector({
    source: new ol.source.Vector({
    url: 'http://localhost:8082/geoserver/basic_map/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=basic_map:sections&maxFeatures=50&outputFormat=application%2Fjson',
    format: new ol.format.GeoJSON(),
    crossOrigin: 'anonymous',
    })
    });

 var mapcentre= ol.proj.transform([157.83430,-33.07209],
'EPSG:4326', 'EPSG:3857');

  var map = new ol.Map({
    layers: [

      vectorLayer
    ],
    target: 'map',
    view: new ol.View({
      center: mapcentre,
      zoom: 14
    })
  });

The URL link is just copied from the layer preview, which I'm not sure on as some examples use wfs instead.

The developer console also shows an error saying 'Cross-Origin Request Blocked: The Same Origin Policy disallows reading the remote resource at http://localhost:8082/geoserver/basic_map/ows?.........'


texts:
In my code, I have displayed the statistics using hardcode i.e information is fed manually in a list and displayed.
I would like to display all statistics correctly using for loop by iself automatically.

//Create a circle by drawing a 2000 m buffer around a point and saving it to variable roi.
var roi = ee.Geometry.Point([-76.147, 43.046]).buffer(20000);

//Create a variable that points to a particular Landsat image in the Earth Engine
var myimage = ee.Image('LANDSAT/LT05/C01/T1_TOA/LT05_015030_20100531');

// Compute the mean .
var meanDict = myimage.reduceRegion({
  reducer: ee.Reducer.mean(),
  geometry: polygon,
  scale: 90
});

// Compute the median.
var medianDict = myimage.reduceRegion({
  reducer: ee.Reducer.median(),
  geometry: polygon,
  scale: 90
});

// Compute the mode.
var modeDict = myimage.reduceRegion({
  reducer: ee.Reducer.mode(),
  geometry: polygon,
  scale: 90
});

// Reduce the image to get a one-band maximum value image.
var maxValue = myimage.reduce(ee.Reducer.max());

//min value
var minValue = myimage.reduce(ee.Reducer.min());

// Extract the bands from collection and make a list
var bands= ee.List([
    [ 'B1','0.45','0.092377477','0.086616925','0.086616925'],
    ['B2','0.51','0.085137095','0.076246974','0.076246974'],
    ['B3','0.59','0.062736553','0.053236987','0.044920224'],
    ['B4','0.67','0.324248892','0.334150672','0.357299032'],
    ['B5','0.88','0.161820343','0.157588621','0.153408903'],
    ['B6','1.65','0.3812067','0.295290527','293.7689209'],
    ['B7','2.29','0.076728936','0.064302692','0.052958585'],
    ['B8','0.5','672.0124447','671','672']
]);

var featList=bands.map(function(list){
  return ee.Feature(null).set({bands:ee.List(list).get(0),wavelength:ee.List(list).get(1),
      mean:ee.List(list).get(2),median:ee.List(list).get(3),mode:ee.List(list).get(4)})
})

var features = ee.FeatureCollection(featList)

//Export table to Drive;
Export.table.toDrive({
  collection: features, //Select FeatureCollection
  folder: 'GEE', //Set folder name
  fileFormat: 'CSV' //Set file format
});


texts: FUSION Microsoft visual runtime library C++ blank error

texts:
I am trying to run the DTM script on 73 laz files (each file >500 MB) and a big single ground filter.laz file > 1 GB that has 303239327 points. The laz files ran just fine with the ground filter script, but as soon as the DTM script begins to read the ground filter file, I get this error. I have a feeling that this is happening due to a memory/CPU processing issue, although upon checking the load,the task manager only shows around 2 GB of RAM in use by FUSION and 54% single core use on an i7. The DTM script otherwise runs fine on smaller size ground filter file. Is there a way to allocate more resources to FUSION? I usually run the batch file through Notepad++. Should I try LTK.exe for this particular script on a big ground filter laz file? In case the script is needed, it is as follows:
:: FUSION/LDV
set FUSION_DIR=%FUSION_DIR%;C:\FUSION\;

:: LAZ file
set LAZ_DIR=D:\LAS_Files\laz

:: Processing
set PROC_DIR=C:\03_processing
cd %PROC_DIR%

:: Project Name
set PROJ=ABCD

:: #--------------------------------------------------------------------------------------------------#
:: # B) ############ Processing #############
:: #--------------------------------------------------------------------------------------------------#

:: 3) # Digital terrain model (DTM): Gridded surface model using collections of random points.
mkDir 03_dtm
C:\FUSION\GridSurfaceCreate /smooth:5 %PROC_DIR%\03_dtm\%PROJ%_dtm.dtm 1 M M 1 11 2 2 %PROC_DIR%\02_groundfilter\%PROJ%_ground.laz
C:\FUSION\dtm2ascii %PROC_DIR%\03_dtm\%PROJ%_dtm.dtm

CMD /k

UPDATE 1
It seems to be a Microsoft visual runtime error


texts: if-else logic expression in field calculator

texts:
I want to calculate the value of a new field using an if-else condition with field calculator in modelbuilder. Want I want to do is as follows: if my new field is empty, then use the expression entered by the user to calculate the value of the field; else do nothing and keep it as it is.

What I have done so far is the following:


texts: what is the rightclick handler in openlayers3

texts: I used rightclick handler with openlayers2 and now I am upgrading my website to openlayers3 and I cannot find this handler in ol3


texts: Resize and relocate Alaska and Hawaii on US county map with R

texts:
How can I resize, relocate and frame Alaska and Hawaii on a US county map, to look more like a traditional map of the fifty states (and DC)?

I've loaded Esri's US county shapefile into R:

library(rgeos)
library(maptools)
# Download from www.baruch.cuny.edu/confluence/display/geoportal/ESRI+USA+Data and unzip
counties.shp <- readShapeSpatial("counties.shp")

I've processed and merged the shapes and data:

library(ggplot2)
counties.shp.f <- fortify(counties.shp, region = "FIPS")
cty_map <- merge(counties.shp.f, cty_pop_msa[,c(2:3,5:6)], by.x="id", by.y="cty_geoid", all.x=TRUE)
cty_map <- cty_map[order(cty_map$order),]

I've plotted a choropleth map of MSA counties:

library(ggmap)
ggplot() +
  geom_polygon(data = cty_map,
               aes(x = long, y = lat, group = group, fill = metro),
               color = "black", size = 0.25) +
  coord_map(xlim=c(-200,-50), ylim=c(25,70)) + theme_nothing(legend = TRUE)

However, Alaska and Hawaii are represented to scale and in their proper locations (I guess), which makes the contiguous US barely legible in A4 format, Hawai invisible, and Alaska too emphasized :

What I need is something more like this, in terms of Alaska's and Hawaii's relative sizes and positions, e.g. sized to fit south-west of Texas and framed:


texts:
I have a map (raster, tif file) with 7 land use classifications and 6 colors representing each classification.

The 7 classifications are

Land use 1
Land use 2
Land use 3
Land use 4
Land use 5
Residential
Urban

I would like to produce a binary exclusion layer (raster) that excludes land use 1.

To do this, I am trying to use the reclassify tool. However, the categories in the table are strange.

Under old values

0-1
1-2
3-4
4-6

There is no 'Residential' or 'Urban category' in this table.

So, how can I properly reclassify the data so that 'land use 1' is reclassified as 0 and everything else 1 (complies with the binary exclusion criteria).

I'm using Arcmap 10.6


texts:
I have two Feature collections: One that includes the location of dams and one that provides the basin area of a river within a country. I want to use gee to obtain a list of dams constructed in a basin area (I don't need dams outside of basin areas). In the end I want to have a table that includes all properties of the dam feature collection and add the name of the country ("CCODE") and the name of the river basin ("BCODE") from the basin dataset.
I used a spatial join. But when I export the data I loose the information on the basin dataset (I guess it should be in 'scenes'). This is my code:

//Load dams
var GrandDams = ee.FeatureCollection('users/basins/GRanD_dams_v1_3');

//load basins
var BCU = ee.FeatureCollection('users/basins/TFDD_BCU');

var spatialFilter = ee.Filter.intersects({
  leftField: '.geo',
  rightField: '.geo',
  maxError: 10
});

var saveAllJoin = ee.Join.saveAll({
  matchesKey: 'scenes',
});

var BCU_dams = saveAllJoin.apply(GrandDams, BCU, spatialFilter);

Export.table.toDrive({
  collection: BCU_dams,
  description: 'Dams',
  fileFormat: 'CSV',
});

How can I add the country and basin name from the basin dataset to the exported table?


texts: Viewing GeoServer tile layers in QGIS through WMS

texts:
I have an issue with viewing GeoServer 2.13-generated tile cache layers in QGIS 2.18.18 on Windows 10. The tile cache layers in GeoServer are being viewed in QGIS via an Add WMS/WMTS Connection, with a WMS connection URL of http://localhost:8080/geoserver/wms?tiled=true. GeoServer is configured with the Direct WMS Integration option selected, which allows WMS requests served through GeoServer to be cached as if they were received and processed by GeoWebCache. When viewing tiled layers from GeoServer with a URL like http://localhost:8080/geoserver/geospace/wms?tiled=true&service=WMS&version=1.1.0&request=GetMap&layers=geospace:Boulder_3857&styles=&bbox=-1.1722528E7,4858919.5,-1.1704111E7,4880837.0&width=1024&height=900&srs=EPSG:3857&format=application/openlayers, the performance is near-instantaneous, as you would expect. However, in QGIS, after adding a GeoServer tiled layer per the WMS connection URL of http://localhost:8080/geoserver/wms?tiled=true, a zoom to layer command takes way too long (~10 seconds for ~230K polygons) to bring up the initial view. After that zoomin/out commands are speedy, but whenever reverting back to the overall layer view the same delay is encountered.

Is this a function of how QGIS handles WMS tile layers? Is there a workaround to speed up the initial view?


texts: Newbie here trying to figure out how to work with Excel data in ArcMap. I am currently using ArcMap ArcGIS for Personal Use. Before my online ArcGIS Online subscription ran out I managed to add a layer of Excel real estate data and the coordinates were as I intended. This was due to the address geocoding feature I suppose. After displaying x y data, I would geocode address and it would work fine. However, the online free trial ran out I have not been able to geocode addresses inside of ArcMap. I'm not sure what I'm doing wrong. I have tried everything and the data is showing on the map halfway across the world. Any suggestions?


texts: How to deal Leaflet map control checkboxes?

texts:
I add some geojson layers to the map control once a layer is loaded and the control can show all the layer names. But I notice that all the checkboxes are unchecked by default. So that once I click the control to uncollapse it, all the layers disappear. And I have to check them one by one to display them again.

The code is like:

ctrl = L.control.layers([aBaseLayer]); ctrl.addOverlay(aGeoJsonLayer, 'ageojsonlayer');

Then I can control the aGeoJsonLayer, but unlike the aBaseLayer, the aGeoJsonLayer is not checked by default in the control panel. So every time after the layer is loaded, when I click the control, the layer will be gone until I check it in the control panel.  Any way to fix this?


texts:
I'm getting crazy with grass gis and python. I want to call a grass function from an external python script (using spyder for instance).

I've installed QGIS + GRASS with the OSGeo4W64 installation suite (Desktop Express installation mode).

My code is the following:

I've specivied which line is working and which is not in the code:

import os
import sys

GISBASE = os.environ['GISBASE'] = 'C:\\OSGeo4W64\\apps\\grass\\grass-7.2.2'
#GRASS_PYTHON = os.environ['GRASS_PYTHON'] = os.path.join(GISBASE,'extrabin','python.exe')
#PYTHONHOME = os.environ['PYTHONHOME'] = os.path.join(GISBASE,'Python27')
#GRASS_PROJSHARE = os.environ['GRASS_PROJSHARE'] = os.path.join(GISBASE,'share','proj')

#os.environ['PATH']+=os.pathsep +os.path.join(GISBASE,'extrabin')
os.environ['PATH']+=os.pathsep +os.path.join(GISBASE,'bin')
os.environ['PATH']+=os.pathsep +os.path.join(GISBASE,'lib')
LD_LIBRARY_PATH= os.environ['LD_LIBRARY_PATH'] = os.path.join(GISBASE,'lib')

os.environ['PATH']+=os.pathsep +"C:\\OSGeo4W64\\bin"

gisdb = os.path.join(os.path.expanduser("~"), "grassdata")

location = "lacLeman"
mapset = "PERMANENT"

grass_pydir = os.path.join(GISBASE, "etc", "python")
sys.path.append(grass_pydir)

import grass.script as gscript
import grass.script.setup as gsetup
#this line works
rcfile = gsetup.init(GISBASE, gisdb, location, mapset)

## this line works
print(gscript.read_command('r.info',map='elevation',flags='g'))
## this line doesn't (that's what I want ultimately to do)
#gscript.run_command('r.horizon', elevation='elevation',step=1,coordinate=(512102.844311,148320.733533))
## this line doesn't work
gscript.read_command('g.region', flags='p')

My Output Errors:

I got two successive windows error (see screenshot below):

Additionally I code an error in my Ipython shell:

  UnicodeDecodeError: 'ascii' codec can't decode byte 0xe9 in position 45: ordinal not in range(128)

But I think this error is just due to the errors mentioned above.

What did I miss?

I'm running on windows 10.


texts: Coordinate Transformation in Gdal/Ogr (EPSG:4326 to EPSG:3068)

texts:
I think gdal/ogr(/proj4) doesn't recognize epsg:3068 right?

Here is what I tried:

public static void transformSRS(){
        // creating EPSG:4326
        SpatialReference sourceSRS = new SpatialReference();
        sourceSRS.ImportFromEPSG(4326);
        // creating EPSG:3068
        SpatialReference targetSRS = new SpatialReference();
        sourceSRS.ImportFromEPSG(3068);
        // printing both SRS's
        System.out.println(sourceSRS); // does get printed
        System.out.println(targetSRS); // doesn't get printed!?
        // creating a point in EPSG:3068 coordinates 850,-150
        Geometry point = Geometry.CreateFromWkt("POINT(850 -150)");
        // assigning EPSG:3068 to the point so the system knows that it is in EPSG:3068
        point.AssignSpatialReference(sourceSRS);
        // trying to transform it to EPSG:4326 but Exception occurs.
        point.TransformTo(targetSRS);
}

and this is the output (first System.out.println() is working but the second one prints null and than the exception at point.TransformTo(targetSRS) ):

GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]

ERROR 6: No translation for an empty SRS to PROJ.4 format is known.
Exception in thread "main" java.lang.RuntimeException: OGR Error: General Error
    at org.gdal.ogr.ogrJNI.Geometry_TransformTo(Native Method)
    at org.gdal.ogr.Geometry.TransformTo(Geometry.java:409)
    at org.geotools.mavenGdal.Test.transformSRS(Test.java:880)
    at org.geotools.mavenGdal.Test.main(Test.java:109)


texts: Clip boundaries of heatmap according to another feature

texts:
Here is a heatmap in Leaflet (using heatmap.js and leaflet-heatmap.js) with hot spots that extend beyond the basemap boundaries:

How would someone clip this? In QGIS or ArcMap you would follow a certain workflow, but this is not available in Leaflet.

We can set a GeoJSON layer of just the Georgia state, but how would one tell the heatmap not to cross the boundaries of this GeoJSON layer?

Slightly similar question (not using Leaflet), unanswered: Limiting OpenLayers heatmap to vector feature boundaries


texts:
I have a map over Sweden with plenty of data represented as circles of various sizes. The map is saved as PDF from Illustrator, fully editable.
I wonder if there is any way I could get these circles into Tilemill so that I can make a searchable georeferenced map? I don't know the projection but I would guess WGS84 or SWEREF99.

One way that comes to mind is to extrakt the geodata from a geotiff of Sweden (if I find one with the same projection), place PDF-map on top of the geotiff in Illustrator/Photoshop and then use gdal to insert the geodata again.

Any ideas if this is a good way, or if there are better methods?


texts: Measuring distance between points yet within polygons using ArcGIS Desktop?

texts:
I have two data sets. One has points of populations and the other has one merged polygon of all suitable habitat types.

I want to find the distance between all the point to one another, but without exiting the polygon. So the route would have to go around any holes in the polygon. I cannot add network to the polygon so I do not know if this is possible.

I am using ArcGIS Desktop.


texts: Finding polygon self-intersection with PyQGIS on QGIS 3.4

texts:
My script:

from qgis.core import *
import qgis.utils
import processing, os, sys, glob

polys = QgsVectorLayer("/Users/john/shapefile/layer.shp","intersection","ogr")
output = "/Users/john/shapefile/intersection.shp"
processing.run("saga:polygonselfintersection", polys, output)

Generates the return message:

if param not in parameters or parameters[param.name()] is None:
TypeError: argument of type 'QgsVectorLayer' is not iterable


texts: Connecting to non geodatabase MS Access database using ArcObjects

texts:
I am trying to connect to an MS Access database from a program I have written in VB.Net with ArcObjects 10.0.  My MS Access database is not a personal geodatabase so apparently I have to connect to Access using an OLEDBWorkspaceFactory.

Dim propset As IPropertySet2 = New PropertySet
propset.SetProperty("CONNECTIONSTRING", "Provider=Microsoft.Jet.OLEDB.4.0;Data Source=" & txtbox1.Text & ";User Id=admin;Password=;")

Dim oleWorkspace As IWorkspaceFactory2 = New OLEDBWorkspaceFactory
Dim m_workspace As IWorkspace2

If oleWorkspace IsNot Nothing Then m_workspace = oleWorkspace.Open(propset, 0)

The above code works except for one minor detail.  When I call oleWorkspace.Open(propset, 0) I get the ArcGIS Database Connection dialogue asking what driver to use and for connection information (which I have already provided in my propertyset).

Is there a way to avoid this dialogue?


texts: How to display a visualization created with CartoDB in a Joomla article?

texts:
The code in the article is:

<div id="map" style="width:600px; height:400px;"> </div>
<script>
    function main() {
        cartodb.createVis('map','http://tsubaki.cartodb.com/api/v2/viz/4a9159ca-467e-11e4-80a0-0e4fddd5de28/viz.json' , {
            shareable: true,
            title: true,
            description: false,
            search: false,
            tiles_loader: false,
            center_lat: 0,
            center_lon: 0,
            zoom: 3
        })
        .done(function(vis) {
            map = vis.getNativeMap();
        })
        .error(function(err) {
            console.log(err);
        });
    }
    window.onload = main;
</script>

but when I load the web page of the article I get these errors with Firebug:

  creating a layer without type
  ..."geo/gmaps/gmaps_cartodb_layer.js","geo/gmaps/gmaps.js","ui/common/dialog.js","u...
  cartodb.js (riga 15)

  TypeError: e is undefined
      ...&")+r.join("&")},template:function(e,t){return e.replace(/{ *([\w_]+) *}/g,fun...
  cartodb.js (linea 7, col 30868)

  TypeError: n.remove is not a function
  ...ed[e]=t;if(!t)n=this.interaction[e],n&&(n.remove(),this.interaction[e]=null);els...
  cartodb.js (linea 18, col 18443)

I don't understand...

Please someone help me!

This is the link of the page. Can you help me with the errors?
http://www.locattivitaecon.altervista.org/index.php/mineral-deposits-map


texts: How can I insert a bounding box or polygon in a PostGIS geometry field?

texts:
I have a field in my table called 'geometry' that was created thus:

CREATE TABLE "public"."sites" (
    "id" int8 DEFAULT nextval('sites_id_seq' :: regclass) NOT NULL,
    #... A few MORE fields
    "geometry" "public"."geometry",
    CONSTRAINT "sites_pkey" PRIMARY KEY ("id")
) WITH (OIDS = FALSE);
CREATE INDEX "index_sites_on_geometry" ON "public"."sites" USING gist (geometry);

I can create a point with the command:

UPDATE sites
SET geometry = ST_SetSRID (
    ST_MakePoint (44.28, - 96.68, 495),
    4326
),
WHERE
    ID = 12;

I would like to create a bounding box, but being new to the syntax have not been able to get the correct syntax.

My best attempt:

UPDATE sites
SET geometry = ST_SetSRID (
    ST_3DMakeBox (
        ST_MakePoint (44.2, - 96.6, 400),
        ST_MakePoint (44.3, - 97.0, 400)
    ),
    4326
)
WHERE
    ID = 1123;

Gives the error:

  Column has Z dimension but geometry does not

I get this error even after trying

wrapping the ST_SetSRID(...) in a select statement as some online examples show (though the select statement returns a 2D result)
using ST_Envelope(ST_Makeline as suggested by @mattmakesmaps

My apologies that I am new to this. I have been looking through the documentation and googling postgis insert bouding box etc, but at my level, this isn't clear.


texts: Georeference an unreferenced image using rasterio in Python

texts:
I have unregistered/unreferenced 6 band TIFF images that I have been trying to export as georeferenced GeoTIFFs. The GPS coordinates of the corners of each image is known and each image is 1024 height and 1280 width. I attempted the GDAL solution in https://stackoverflow.com/questions/55681995/how-to-georeference-an-unreferenced-aerial-image-using-ground-control-points-in but was receiving "item has no get_attr_" errors at the setGCP command. Attempt 1:

import rasterio as rio
from rasterio import Affine
from rasterio.warp import reproject, Resampling, calculate_default_transform
from rasterio.transform import from_origin
from rasterio.control import GroundControlPoint

im = rio.open('path to a tif')
arr = im.read()
dst_fn1 = outdir + 'out.tif'

tl = GroundControlPoint(0, 0, -83.70113410256013, 42.307951446432604)
bl = GroundControlPoint(1024, 0, -83.69940501521428, 42.307603183805234)
br = GroundControlPoint(1024, 1280, -83.698829074736, 42.3091785425499)
tr = GroundControlPoint(0, 1280, -83.70055820297041, 42.309526812647555)
gcps = [tl, bl, br, tr]

src_crs = im._get_crs()  ## This is a blank object. No crs
dst_crs = {'init': 'epsg:4326'}
src_tran = im.transform
dst_transform, dst_width, dst_height = calculate_default_transform(a, crs, 1280, 1024, gcps=gcps)
dest_ds = np.empty(shape=(arr.shape[0], dst_width, dst_height))
reproject(sand, dest_ds, src_crs=dst_crs, dst_crs=dst_crs,
      src_transform=src_tran, dst_transform=dst_transform,
      resampling=Resampling.nearest)

This method does not throw any errors but the resulting array is only zeroes. Also exported array is not rotated to match where the corners should be. Perhaps this is the wrong expectation but I thought that the data would be inside the yellow box and the outside would be nodata values to square it off.

Attempt 2: ( This is essentially the method found here https://www.earthdatascience.org/courses/earth-analytics-python/lidar-raster-data/reproject-raster/)

with rio.open('path to a tif') as src:
    transform, width, height = calculate_default_transform(src.crs, dst_crs, src.width, src.height, gcps=gcps)
    kwargs = src.meta.copy()
    kwargs.update({'crs': dst_crs, 'transform': transform, 'width': width, 'height': height})
        with rio.open('path/output.tif', 'w', **kwargs) as dst:
            for i in range(1, src.count + 1):
                reproject(
                    source=rio.band(src, i),
                    destination=rio.band(dst, i),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=transform, dst_crs=dst_crs,
                    resampling=Resampling.nearest)

Results in this error "rasterio._err.CPLE_AppDefinedError: The transformation is already "north up" or a transformation between pixel/line and georeferenced coordinates cannot be computed for TIFFPATH. There is no affine transformation and no GCPs. Specify transformation option SRC_METHOD=NO_GEOTRANSFORM to bypass this check.


texts: Reduce squarish effect on group of polygon

texts:
I am trying to represent the average speed of some ship in the north of Mediterranean sea.

I am interested in aggregate the speed information by zone.

The raw data I am working with are .shp with the following information:
Lat_first;Long_first;average_speed

the representation is the following:

I have managed to create a vectorial grid and then calculate the average speed in every polygon as shown here and spatial join the two layers:

the result is the following picture:

which has correct data but is too squared, I would like to smooth the edges to create something more similar to a heat map or interpolation.

Do anybody have dome ideas ? I am currently using QGIS but I would change if someone has a simpler procedure on some other software.


texts:
Bottom Line Issue: The distances I get from NNJoin for the shortest distance from a point to the nearest polygon are several hundreds or even thousands of meters shorter than the distance I calculate using the "Measure Line" tool in QGIS. Any idea why these might be different?

More Details: First, I imported a shapefile of polygons (census designated places in US) and a text file of points (long-lat locations). Next, I projected them both into WGS/UTM Zone 11N (the relevant UTM zone where they are located). Finally, I used the NNJoin plugin to find the distance between each point and the nearest shapefile. However, the distances that NNJoin returns seem to be several hundred or even thousands of meters shorter than the distances I calculate using the Measure Line tool in QGIS (or doing the calculations using Google Maps). For example, NNJoin says the distance from a point called Brown's Gultch (lat: 41.834721, long:
-115.97) to Owyhee, Nevada (nearest census designated place) is 2,391 meters. By contrast, the "Measure Line" tool indicates it is around ~2,600 meters.

Any idea what I am doing wrong?


texts: Make sidebar div active on Marker click (and vice versa)?

texts:
What's the best way to change an element to active when you click on a marker.

For example, when you click on a marker on a Leaflet map, how to make it so the sidebar element that's related to it will be active?

I can write the sudo code but I'm missing pieces.

markers on click
$(#sidebar div)
   make (related) li element active

edit:
Basically, I'm trying to mimic what a lot of map directory sites like Yelp or Foursquare do.  You click on the sidebar list item and the marker becomes active, and vice versa.  This makes it so a user and can read info about the marker while seeing where the place is.

Mapbox had a cool (now deprecated) Foursquare example. I went through the code for that but barely understood it.  That said, the example gets the point across.

When you mouseover a listitem in yelp, the marker goes black.


texts: Runaway memory usage with Event Layers in ArcGIS Engine and Desktop?

texts:
I have an ArcEngine application, based off an MXD with an axMapControl.  Several of the layers in the map are event layers built with XY data in a SQL table.  It has routing info applied from a routing table too.

My problem is that when I pan the map, the memory usage grows and grows.  The same behavior is exhibited in ArcMap.  This inevitably leads the target computer to run at a crawl as memory is depleted, and hard disk swapping all the time.

I'm assuming that the map control is caching the spatial side as new data is drawn on the map.

Is there any way to manage this memory usage?  Perhaps flush from time to time, or disable all together?  At this stage slower drawing would be preferable as it's simply unworkable in this scenario.

UPDATE:
Further investigation suggests that this is always the case, though on low spec (CPU) computers, when panning a lot, the CPU simply can't manage/flush the cache as fast as it grows.  On faster machines, you can see the memory climb, but it returns to a lower level pretty quickly too.  So this one's simply chalked up to slow CPU's not being able to keep up with all the map control cleverness.  It would still be see a good explanation of how this all works behind the scenes.


texts: Using ESRI Query Layer to query OSM data in PostGIS database connection?

texts:
I have managed to connect to my OSM data using database connections in ArcGIS Desktop.  I have the OSM line, point, polygon and roads tables.  Now I am trying to query the tables using the Add New Query Layer.

I can load all data into ArcMap, but I would like to query the OSM tags.  The column tags is currently a Unknown type.  Therefore all my queries seems to fail.

Does anyone know how I should query the tags column here?


texts: Leaflet: Capturing points for features inside circle

texts:
I'm using the leaflet draw tool to create a circle and then finding the points inside the circle. My data is from GeoJSON, and the code works as far as figuring out how many points are inside, but I want to then capture the point in a new array, test.push.layer doesn't work, it seems to create a null object, not a layer.

if (type == 'circle') {
    var test = [];

    var theCenterPt = layer.getLatLng();
    var theRadius = layer.getRadius();
    var counter_points_in_circle = 0;

    // Loop through each point in GeoJSON file  //var allPoints = L.geoJson(data);
    allPoints.eachLayer(function (layer) {
        var theLayer = layer;
        // Lat, long of current point
        layer_lat_long = layer.getLatLng();
        // Distance from our circle marker To current point in meters
        distance_from_centerPoint = layer_lat_long.distanceTo(theCenterPt);
        // See if meters is within radius
        if (distance_from_centerPoint <= theRadius) {
            counter_points_in_circle += 1;

            test.push.layer;     //Fails

        }

    });
    alert(counter_points_in_circle);  //Works.
    alert(test.features.length);    //test.features is undefined..
    alert("results "+JSON.stringify(test));  //results []
    alert("results "+JSON.stringify(test.toGeoJSON()));  //test.toGeoJSON is not a function
}


texts: ogr2ogr add geometry from x y when importing to postgis

texts:
Im importing a xlsx with x y coordinates to a postgis database, with the following command:

ogr2ogr -overwrite -progress -nln data -nlt POINT -f PostgreSQL "PG:host=Localhost user=postgres password=Postgres dbname=database" -nln data C:\data\data.XLSX

It successfully add the data to my database. With some sql included in the command i will now like to take the x y columns and update the (empty) geometry column, so that the x y data is converted to geometry.

I tried to add:

-sql "SELECT AddGeometryColumn ('database', 'data', 'geom', 26913, 'POINT', 2)" -dialect sqlite

Just to try to add / define some geometry to my table.

Any suggestion on where to look or what to do?


texts:
I can create a NetCDF file (using the Python netCDF4 library) and load it in QGIS. However, as the size of these files gets sometimes quite big, I try to compress my data. This works very well using

createVariable(h, 'float32', ("y","x"), zlib=True,complevel=3,chunksizes=(1,len(x)),fill_value=np.nan)

However, playing around a bit showed that chunking both dimensions (e.g. chunksizes=(100,100) or chunksizes=(len(y),len(x))) reduced the resulting file size by more than 50%. Unfortunately, in this case I am not able any more to load this grid correctly to QGIS. The only chunksize that works is chunking only the rightmost variable (1,len(x)). I guess, this has something to do with the GDAL. About QGIS says:

QGIS version 2.14.12-Essen
Compiled against GDAL/OGR 1.10.1
Running against GDAL/OGR 1.11.2

I'm working on Ubuntu14.04 and installed QGIS via the debian-ltr PPA. Is there any way to make other chunksizes work too? Would an update to QGIS 2.18 help or is this an issue with GDAL?


texts: XY Coordinates for Points in odd format, not matching any projection

texts:
I am trying to match up some data without success. I have a state shapefile from the TIGER/Line and a series of X/Y Coordinates from a St. Louis, Missouri police department (CSV).

I am using an OpenStreetMap basemap in QGIS 2.14.8 to verify the data is in the right place.

The Tiger/line shapefile is in NAD83.

The XY coordinates are in what I can only guess is U.S. Survey feet? Example: "903948.39, 1033477.17". I do not know what CRS they are in. When I try to project them to..well..anything, they do not match up to the Tiger/line shapefile. I assume this is because the units are so wonky.

Is it possible to display these points in a format that will actually coincide with the Tiger/line shapefile? I'm not sure how to convert them since the units are wonky. Simply changing the projection to NAD83 does not help.

If it is possible, will this solve the problem and how do I do it?


texts:
I have a large .gdb of all blocks in the United States, available from Census. I would like to create a point shape file of all census tract centroids in the US.

I have QGIS with FileGDB enabled, as well as R with rgdal, so reading the data is not a technical problem. My issue is that reading the layer into R has taken 2 hours and used up 56 Gb (I have 64) of memory, and it's not done. I know exactly what I want to do:

Create a tract field by pulling the first 7 characters out of the GEOID field
Dissolve the blocks on this field
Calculate the centroids of the new geographies
Export as a points shape file, or even as another layer in the database.

QGIS shows the data, but the field calculations take just about as long. It seems that I should be able to send these commands into the database rather than load the database and try to do it interactively. I might need some hand holding to do this, or at least good pointers to the right manuals.

On the other hand, is there a way to pull the database out in chunks? I am running Ubuntu 14.04 on a 12-core machine, so parallelization is always a good option...


texts: Changing Layer Order with ArcPy in ArcGIS Pro

texts:
Working on a Python script that does some routine table management of a couple of feature classes. In one step I define a table join and then export that to a GDB. When that runs, it will add it to the Geodatabase and also add the new Feature Class as a layer to the table of contents. It always adds to the top of the Table of Contents, but I need it to be added to the bottom as the last item.
The problem that I am running into is that in my code, I define several local variables using an index from the Table of Contents like old_sites = mapx.listLayers()[3]
If I have this hardcoded and a new layer is added to the top, this layer is no longer indexed at position 3, it's now 4, but the code won't know that.
Is there a way to change the default add position to the bottom?


texts: QGIS python plugin appends values every time it's started

texts: I've written a small python plugin (using plugin builder as base) to simplify the expression process for a given layer (just a few comboboxes that are used to create a query). The combo box populating is done in the def run(self): part with a simple .addItem loop. However when I close the plugin and reopen it again, the values are added again. So if it's "a, b, c" and I close/open the plugin again, it's "a, b, c, a, b, c" in the box.. can I stop/avoid that? Or should I just wrap it in a case that sets a flag on the first run?


texts: OpenLayers 3: Missing Tiles when map is wrapped horizontally

texts:
I am encountering a strange error when my map is being wrapped horizontally. See image below for the problem. I have three layers in the map: the grid, countries, and the colors. As you can see the grid is being generated. Also I am using Geoserver 2.6.1 to generate the tiles. Any insight on the problem would be helpful.

Code for map creation:

this.map = new ol.Map({
    controls: ol.control.defaults().extend([
        new ol.control.ScaleLine()
    ]),
    target: 'map',
    view: new ol.View({
        projection: projection,
        center: [0, 0],
        zoom: 2
    })
});

The layers get added dynamically but here is a general way I create the layers:

        var layer = new ol.layer.Tile({
            id: i,
            title: layerName,
            layerName: layerName,
            layerIndex: i,
            source: new ol.source.TileWMS({
                url: this.geoServerUrl,
                serverType: 'geoserver',
                params: {
                    'LAYERS' : layerName,
                    'TILED': true,
                },
            })
        });


texts: GDAL 3.0.4 invalid coordinate transformation result via ogr2ogr

texts:
I'm using command
ogr2ogr -f GML -s_srs "+proj=latlong +datum=WGS84 +axis=neu +wktext" -t_srs "+proj=latlong +datum=WGS84 +axis=enu +wktext" ./file.gml ./file.shp

to transform file format and change long/lat order to lat/long order.
GDAL 2.2.3 works fine, but GDAL 3.0.4 doesn't change the order.
What am I doing wrong?


texts: Map Matching between two geometry (line string)

texts:
I have two tables in PostGIS, both having line string as the geometry type.  I wanted to extract a similar or exact matching street. Can you please suggest query or algorithm for it?

Table 1

column name(Id Geom1) data(1 0102000020E610000002000000590E3EDEF5AB23409255B6B4BF1B4A4031197DBBDBAB23404AF663EEB51B4A40)

Table 2

column name(Id Geom2) data(1  10102000020E61000000300000057900B73277D2340A1675831011E4A4094C55801197D23406A204C40021E4A40B0CFF7AE9C7C234079AAE8B4131E4A40)

I have tired ST_HausdorffDistance(table1.Geom1, table2.Geom2)

but not able to get the desired result, my dataset is big, table one has 2370 record and table two 23125 records.

Can you please some effective way to do map matching in Python or Using SQL in Postgis.


texts:
This question seems simple but I couldn't get it done for most of the layers I tested:

This is the adress that is predefined. It works with a small subset
ogc.bgs.ac.uk/cgi-bin/BGS_Bedrock_and_Superficial_Geology/wms

I tried several others but keep getting  Tool Execution Failed

http://ows.mundialis.de/services/service?
http://ows.terrestris.de/osm/service?
http://www.webatlasde.de/arcgis/services/Maps4BW/MapServer/WMSServer?
http://sg.geodatenzentrum.de/wms_webatlasde.light

These all work in QGIS but not in SAGA:

I understood it as follows (and so worked the first example):

Open Import a Map Via Web Map Service (WMS)
Insert the address, leave all other fields
Hit Okay
The image pops up again (see below), enter extents and select bands
Hit Okay

Works for me with the given default address but not with the others. Any ideas?


texts: QGIS (qgs file) meaning of annotationform and editform tag

texts:
What is the meaning of tags annotationform and editform in the qgs file? The reason I ask is because there is only file path within. Mostly relative path as dot, sometimes path to something like ../../../My_Project/GIS. And by something I mean I have no idea why this path. Is it supposed to set path for some kind of form or whatever might be saved there or why? There is definitely not saved anything apart of qgs projects in that directory.

If I remove the paths from within these tags completely or leave simple dot is it ok? What exactly is it going to affect?


texts: Display combinations of multiple fields containing same data type for each feature class point

texts:
I am trying to display small mammal trapping data in ArcMap 10.3.1.

I am dealing with a point feature class.  There are 3 relevant fields in the attribute table, one for each night of trapping, containing categorical data (type of species caught).  Each night, the traps can catch 1 of any 4 species.  So, some points will have a capture of the same species across the 3 nights, while others might capture a different species each night.  There might be any combination of species caught during the 3 nights of trapping for any particular point.

I am interested in displaying which species were caught at which trap.  I don't care about the order or if a trap caught the same species more than once.

I am having difficulty finding an efficient way to display these data.  I've tried "Unique values, many fields" in the symbology tab of properties box, but that gives all permutations while I am only interested in relevant combinations.

I think one efficient display would be to assign each species a circle with a unique color and diameter.  This way, multiple species could be seen for each point if one trap catches multiple species, and the legend would only have to show 1 symbol for each species.

Another problem I cannot resolve is how to collapse the information when one point has duplicate data across fields.  In other words, how to collapse if one trap catches the same species for multiple nights.


texts: GeoServer GeoWebCache caching wrong resolution tiles from USGS WMS

texts:
I have a simple OpenLayers 3 map with a topo layer that is cascaded from a USGS WMS through my GeoServer / GeoWebCache system.  The tile zoom level resolution is correct when the GeoWebCache is turned off, but I get the wrong resolution tiles from the cascaded USGS link when I turn on the caching with the GeoServer GeoWebCache.    For example, at zoom level 13, the tiles should show 50-foot contours like this:

When the caching for the layer with GeoWebCache is turn on, the tiles showing 20-foot contours for level 14 are rendered into level 13 like this:

It is the same for all zoom levels.  When caching is turned on, the current zoom level gets the tiles for the next zoom level.

The zoom level tile resolution problem started when the USGS moved this layer and their other cached layers to Amazon Cloud on July 1, 2017.

Why does this happen just when the caching is turned on for the integrated GeoWebCache or when going through the separate stand alone GeoWebCahce?


texts: shapely Accessing Data from Geometry collections

texts:
I am new to shapely and python. I recently read Splitting Polygon by Linestring in GeoDjango?  and below is the resulting answer.

from shapely import wkt
from shapely.ops import linemerge, unary_union, polygonize

POLY = "POLYGON ((34.67491149902344 31.59900710035676, 34.85000610351562 31.59900710035676, 34.85000610351562 31.73867905688433, 34.67491149902344 31.73867905688433, 34.67491149902344 31.59900710035676))"
LINE = "LINESTRING (34.64401245117188 31.63292168314889, 34.80812072753906 31.75911546882192)"

poly = wkt.loads(POLY)
line = wkt.loads(LINE)

merged = linemerge([poly.boundary, line])
borders = unary_union(merged)
polygons = polygonize(borders)
for p in polygons:
    print(p)

The print generates two polygon collections.

How do I access the arrays?


texts: Calculating elevation profile along line from a DEM?

texts:
Given a DEM (Digital Elevation Model) and a line (x1,y1) --> (x2,y2) (paired coordinates (x,y)), how can one calculate the elevation profile or cross-section projected on that line?

I am looking for an open source code that I can use in my project or pseudo-code that I can use as a guide.


texts:
I am trying to add two different rasters using the Con function. I am iterating the raster layer using the for loop to obtain two rasters belonging to the same species and their formats are in grid and tif format. On running the below scripts i am not able to execute the con function. If anyone has a way to solve this issue, let me know your suggestion.

# Import arcpy module

import arcpy

import os

from arcpy.sa import Con

# Check out any necessary licenses

arcpy.CheckOutExtension("spatial")

# Overwrite the ouptut

arcpy.env.overwriteOutput = 1

# set environmental settings # Location of the rasters

output_workspace = "C:\\GIS_Analysis_programming\\GAPanalysisData\\Output"

allMamRows = arcpy.SearchCursor("C:\\GIS_Analysis_programming\\GAPanalysisData\\Flmammal.dbf")

# loop through mammal list and load code to list

     mamList = []

     i=0

     for row in allMamRows:

         mamList.append(row.getValue("CODE"))

       # print mamList[i]

         i = i +1

# Local variables:
# Creating a list of rasters from the workspace.
# Iterations through (only 10) mammals

   j=0

   for species in mamList:

    if j < 10 :

    species_str = str(species)

    rastersHab = arcpy.ListRasters("hab_species_str", "GRID")

    rastersCnty = arcpy.ListRasters("cnty_species_str", ".tif")

    species_richness = os.path.join(output_workspace,  "_richness" + species_str + ".tif")

    (Con(("%rastersHab%" > 0)&("%rastersCnty%" == 1),1,0),species_richness)

    print j

    j = j+1


texts: CORS lets one ArcGIS JavaScript (REST) request through but not the other

texts:
I have a problem with Cross-Origin-Requests invoked by the ArcGIS Javascript API.

I'm able to access my web app without any issues if I use localhost:8080. But a REST call (by the ArcGIS JS API) fails if I use mymachine.local.

The exportImage REST call works properly but the query REST call invoked by featureLayer.selectFeatures fails with a Not allowed by Access-Control-Allow-Origin error.

I've gone through http://blogs.esri.com/esri/arcgis/2011/09/29/cross-origin-resource-sharing-cors-with-the-arcgis-api-for-javascript/ but I don't think there's a problem with the CORS setup as one request goes through and the other one fails.


texts: Trigger event after KML Layer completely loaded

texts:
I want to set the bounds of my map after a Vector (KML) Layer is loaded.

kmllayer.events.register('refresh', kmllayer, function(evt) {
    map.zoomToExtent(kmllayer.getDataExtent());
});

In the documentation there if something like the afterAdd function. But it does not work at all.

I am using version 2 of Openlayers.


texts: create a vegetation mask using Landsat 7

texts:
I'm working on a very large study area. The aim is to create a vegetation mask representing only pines. I'm using Landsat 7 images. Since the recordings were made in the winter it's pretty easy to distinguish the conifers from the rest, because everything else is covered with snow.

Unfortunately I do not get a satisfactory result. What is the best methodology to get good results? I've tried to use a NDVI. Unfortunately, the results are bad.

I have also tried to use a tasseled cap transformation. Unfortunately, I have no experience using TC, so I'm not sure how to interpret/use the results.

I am using Erdas Imagine and ArcGIS.

Do you have any tips on how I can solve the problem?


texts:
I have many CSV file in the form of lon, lat, Z that in would like to import into gdal for further processing.  I need to do this on the fly in a python script so any step involving gdal_grid or anything else from the command line is not an option.

It seemed that creating a vrt xml would solve the problem but I am unable to get gdal to load the data.

Loading data this way seems to be such a common task I can't believe there is not a way to do this.


texts:
For debugging purposes I would like to print the content of some variables in my PYGT-function, like:

print(myGeometrylength)

...

The variable has some content but nothing is printed?

I just get it out by

qgis.utils.iface.messageBar().pushMessage("blabla", u"blba :"+str(myGeometrylength))

Has anybody an idea how to get the print-command working?


texts: Can't make a pop-up appear with thumbnail nor link (ArcGIS Online)

texts:
I'm trying to enrich my web map by giving pop-up to a layer with a thumbnail with link to larger picture.
I've uploaded two files up to my AGOL account. I'm looking to use those files as pop-ups in a service layer.

I've manged to open them via the 'My Content' GUI. Still I could not open them when copy-pasting the url as an attribute and then using it as a pop-up

I think the problem is some how related to the url's suffix: when a file extension is included the pop-up appears OK(*.jpg); when the url don't contain file extension suffix it doesn't appear.

How to make a pop-up appear with a valid url (for the thumbnail and a link)?


texts: Calculate centre point of SpatialGrid object

texts:
I am using R and have a SpatialGrid object similar to the example from this question:

Overlay a spatial polygon with a grid and check in which grid element specific coordinates are located

## read shapefile
library("rgdal")
shp <- readOGR("nybb_13a", "nybb")

proj4string(shp)  # units us-ft
# [1] "+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333
# +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +datum=NAD83
# +units=us-ft +no_defs +ellps=GRS80 +towgs84=0,0,0"

### define coordinates and convert to SpatialPointsDataFrame
poi <- data.frame(x=c(919500, 959500, 1019500, 1049500, 1029500, 989500),
              y=c(130600, 150600, 180600, 198000, 248000, 218000),
              id="A", stringsAsFactors=F)
coordinates(poi) <- ~ x + y
proj4string(poi) <- proj4string(shp)

### define SpatialGrid object
bb <- bbox(shp)
cs <- c(3.28084, 3.28084)*6000  # cell size 6km x 6km (for illustration)
                            # 1 ft = 3.28084 m
cc <- bb[, 1] + (cs/2)  # cell offset
cd <- ceiling(diff(t(bb))/cs)  # number of cells per direction
grd <- GridTopology(cellcentre.offset=cc, cellsize=cs, cells.dim=cd)
grd
# cellcentre.offset 923018 129964
# cellsize           19685  19685
# cells.dim              8      8

sp_grd <- SpatialGridDataFrame(grd,
                           data=data.frame(id=1:prod(cd)),
                           proj4string=CRS(proj4string(shp)))

I have a SpatialGrid rather than SpatialGridDataFrame, so:

new_sp_grd <- SpatialGrid(grd,proj4string=CRS(proj4string(shp)))

I now want to calculate the centre coordinates of each grid cell, so that I can calculate the distance between this point and a variety of SpatialLines and SpatialPolygons objects.

Can anyone suggest how to do this? There is help for finding the centre of a polygon, but not a grid, and functions such as gCentroid from the rgeos package do not work on SpatialGrids.


texts: Getting CSV file of Wikipedia geolocations?

texts:
Where do I get a CSV file containing (a reference to) all Wikipedia articles with lon/lat coordinates?

I am aware of Retrieving Wikipedia's coordinate information? which refers to this page:
https://de.wikipedia.org/wiki/Wikipedia:WikiProjekt_Georeferenzierung/Hauptseite/Wikipedia-World/en

At one time that page offered a CSV download but it is no longer present, only the option to download a PostGIS data dump.

I prefer to get the dewiki geolocations


texts:
I have been reading a lot about PostGIS, but I am struggling to find a good beginner's guide.

Which tutorials/books do you recommend?


texts: How to sum attributes from different vector grids?

texts:
I am using QGIS and have a vector grid covering part of South Europe. The grid is made of rectangles each with equal size 0.01 degree. So there are thousands of features defined by four values : left, right, top, bottom, which resulted from a standard 'Create grid' task. Also there is a field Count which shows some integer value for each grid cell.

Then I clipped my polygons - about 100 which represent species occurrence in small areas from the big grid. Each polygon I received now has some rectangles (1 to several hundred) defined by coordinates and a value in the count field.

How can I unite all the polygons and sum any overlapping values in each cell? I mean if three polygons occupy the upper left grid then it receives the value 3, etc.

This should be probably better with rasterization but I don't know how to make the pixels the same as these grids which were given to me. CRS is the same for every polygon.


texts: Marking zonal statistics areas with a point

texts:
I'm working with Lidar data, and my end goal is to get a collection of points representing cliff bands, with key information attributed to the individual points. So far I've been able to identify my cliffs create an area around the cliff and used this zone to find the range of elevation data around the cliff (approximate max height). I now have one raster that has a bunch of blobs with different pixel values.

I need help marking a point in the middle of each of these blobs that contains the raster value of their parent blob.


texts: Randomly limiting number of features QGIS displays for specific point layer?

texts:
I have a map with various point layers (e. g. shops, bars, hotels). Some of these contain many features and the map is overcrowded – particularly in urban areas.

Hence, I would like to (randomly) reduce the number of displayed features for specific layers (ideally not by changing the underlying data).
This is possible for labels (Limit number of features to be labeled to... and Discourage labels from covering features in the Layer Properties under Labels -> Rendering).

Basically I am looking for an option such as: “Limit number of features to be shown to...”

What I have:

What I want (e. g. limited to 5 displayed shopping-features):


texts: How to specify one layer among many with getWMSFeatureInfo?

texts:
I am using OpenLayers getWMSFeatureInfo with multiple point layers.

When I click on a point on map, how to specify which layer that point belongs to?


texts: Zonal statistics as table fails using Geoprocessor and C#

texts: The above code snippets fails while calculating the Zonal statistics as table using Geo processor in C#, butit works in ArcMap. Any ideas where/what I did wrong?


texts: How it is possible to create a 3D cube with bars showing some quality, located at certain points with geographic coordinates? Is it possible to do this with some qgis plugin or, maybe, any other comparable GIS software (arcGIS or so)? The desired map would look like as follows (arbitrary picture from search results). Instead of (x-y) parameters there should be geographic coordinates, and some GIS map in the X-Y plane, and colored bars with different heights would be located sporadically (not on a grid) in a set of points, given by their coordinates.


texts: Proximity analysis of posted speed limit

texts:
I want to perform a proximity search around my current position (moving vehicle) to find the links on the roads that have unique posted speed limits associated. I am imagining a big circular region with, say 1 mile radius, around my location as the center, and I want to request the static speed limit for each link or a segment of the road.

Right now i am using the rest-api query method to find out the speed limit at my current location on a given link.However, I want to know the speed limit of the connecting road, if I take the next exit.

The closest match to what i am trying to achieve is on this link: (https://developer.here.com/platform-extensions/documentation/platform-data/topics/resource-search-proximity.html). Has anyone tried this method before?

A solution is given at (https://stackoverflow.com/questions/37679799/here-maps-api-speed-limit-with-multiple-waypoints) where it is advised to use tiles and layers concept.

UPDATE: After some research, i figured out that indeed using tile method is the way to go. Additionally, knowing what HERE maps layers contain is fundamental to designing the algorithm. I found this link really helpful (https://developer.here.com/platform-extensions/documentation/platform-data/topics/example-layers.html).

It is indeed a very interesting puzzle, to know what each layer has to offer and then try querying for the right attributes from a specific layer based on the "functional class" of your choice.


texts: Assign Z value to new coordinates based on interpolated response surface in R

texts:
I've created a response surface of sorts based on interpolated z values (bottom of post).

I'm looking for the most efficient way to assign the respective z value to a new point (actually, a dataset of new points), based on the corresponding lat long (which is what the original grid is based on).

Anyone have a recommendation for a command and/or package I may look up?

The grid is not a shapefile so the spatial joining packages I've seen have not been any use.


texts: QGis QSpatialite plugin corrupting dates when importing layers

texts:
I'm using Qspatialite version 7.0.4 in QGis 2.14. We use QSpatialite to save selected polygons from an MSSQL (2008) layer into spatialite databases.

Since Qgis 2.14 we have noticed that any dates held in the SQL table as a date data type are corrupted when the data is imported into a spatialite database through QSpatialite.

Original Data from MSSQL layer within QGis:

Data after importing into spatialite database using QSpatialite:

N.B the field name is the data type used to create the table in MSSQL.

I have tried using the DB manager and this fails stating 'Error 5 - creation of fields failed'

We currently have a workaround that involves using QGis 2.8 to import the layer and then opening it through QGis 2.14 to continue editing, this doesn't corrupt the data.

Does anyone have any experience of this? or any ideas of other things to try?

Please ask if you need any more information!


texts: What is causing syntax error in arcpy NetCDF Raster script

texts:
I'm trying to save the raster layer in the stand-alone python GUI but I keep getting a syntax error.

I'm using the code from Save To Layer File (Data Management).
And the error that shows up is: the LA MESA part of below:

in_features = "spain_1901_1.shp"
out_layer0 = "spain_1901_1"
where_clause = '"NAME" = 'LA MESA''
workspace = "C:/output"

now I've replaced their fill in things from the website with what I need but I keep getting an error saying that LA MESA isn't the correct syntax even when I replace it with my own name.

please help!

EDIT: Here's the entire code:

#MAKE RASTER FILE
# Import system modules
import arcpy

# Set local variables
inNetCDFFile = "C:/Users/ESS/Desktop/temp/tmpy1901m01.nc"
variable = "tmp"
XDimension = "lon"
YDimension = "lat"
outRasterLayer = "C:/Users/ESS/Desktop/out/spain_1901_1.shp"
bandDimmension = ""
dimensionValues = ""
valueSelectionMethod = ""

# Execute MakeNetCDFRasterLayer
arcpy.MakeNetCDFRasterLayer_md(inNetCDFFile, variable, XDimension, YDimension, outRasterLayer, bandDimmension, dimensionValues, valueSelectionMethod)

#SAVE IT!!!!!!
# Name: SaveToLayerFile_Example2.py
# Description: Saves an inMemory layer to a file on disk
# Author: ESRI

# Import system modules
import arcpy
from arcpy import env

# Set workspace
env.workspace = "C:/Users/ESS/Desktop/out"

# Set local variables
in_layer = "spain_1901_1"
out_layer = "spain_1901_1.lyr"

#MakeFeatureLayer variables
in_features = "spain_1901_1.shp"
out_layer0 = "spain_1901_1"
where_clause = '"NAME" = \'LA MESA\''
workspace = "C:/Users/ESS/Desktop/out"

try:
    # Execute MakeFeatureLayer
    arcpy.MakeFeatureLayer_management(in_features, out_layer0, where_clause, workspace)

    # Execute SaveToLayerFile
    arcpy.SaveToLayerFile_management(in_layer, out_layer, "ABSOLUTE")
except:
    print arcpy.GetMessages()

##THE CALCULATIONS
# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *

# Set environment settings
env.workspace = "C:/Users/ESS/Desktop/out"

# Set local variables
inZoneData = "es2001.shp"
zoneField = "FIDS"
inValueRaster = "spain_1901_1"
outTable = "spain_1901_1_data.dbf"

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Execute ZonalStatisticsAsTable
outZSaT = ZonalStatisticsAsTable(inZoneData, zoneField, inValueRaster,
                                 outTable, "NODATA", "ALL")


texts: Passing Arguments Between Widgets with Mixedin Templates - CMV Builder

texts:
I am attempting to manipulate a custom widget (Table) I built. I first initialize the widget in the viewer.js class and attach it to the 'sidebarBottom' pane in the viewer. Since the 'Find' widget is a titlePane and attached to the sidebar, I need a way to pass attributes in between each widget.

if I initialize the custom widget in the viewer.js class I can successfully attach it to the bottom pane, but I do not know how to pass dojo/dijit arguments to it via the Find widget. I have attempted to use dojo.byId() with no success.

If I initialize my custom widget from inside of the Find widget using the second block of code, I can successfully set() and get() attributes. However, I am unable to attach the widget to the "sidebarBottom" pane defined in the viewer.js declaration.

Table Widget Section in Viewer.js:

  table: {
            include: true,
            id: 'table',
            type: 'contentPane',
            path: 'gis/dijit/Table',
            placeAt: 'bottom',
            options: {
                map: true,
            }
        },

Table Widget Init in 'Find' Widget Class

        var fOut = query.url + "/" + query.layerIds;
        var table;
            table = new Table({
                include: true,
                open: true,
                FeatureInput: fOut,
                map: this.map,
                //placeAt: this.bottom,
            });
           alert(table.get('FeatureInput'));
           alert(table.get('open'));


texts: Calculating maximum value in Raster Calculator of ArcGIS for Desktop?

texts:
I have MODIS NDVI data of four months and I have to calculate NDVImax.

I am using the Raster Calculator of ArcGIS and have been using various algorithms like [MAX ("june2000", "july2000" , "august2000", "september2000")] or [(MAX ("june2000", "july2000" , "august2000", "september2000")]

I tried the same in ERDAS Imagine with Global MAXIMUM function, of which i have attache an image, but I am not able to get any results.

Can anyone help?


texts: Applicable OpenLayers settings for WFS layer

texts:
I am trying to display a WFS layer on an OpenLayers map. Here are the appropriate portions of GETCAPABILITIES response.

<ows:Operation name="GetFeature">
  <ows:DCP>
    <ows:HTTP>
      <ows:Get xlink:href="http://localhost/MapServer/Default.aspx"/>
      <ows:Post xlink:href="http://localhost/MapServer/Default.aspx"/>
    </ows:HTTP>
  </ows:DCP>
  <ows:Parameter name="resultType">
    <ows:Value>results</ows:Value>
    <ows:Value>hits</ows:Value>
  </ows:Parameter>
  <ows:Parameter name="outputFormat">
    <ows:Value>text/xml; subtype=gml/3.1.1</ows:Value>
  </ows:Parameter>
</ows:Operation>
<ows:Parameter name="srsName">
  <ows:Value>EPSG:4326</ows:Value>
</ows:Parameter>
<ows:Constraint name="DefaultMaxFeatures">
  <ows:Value>20000</ows:Value>
</ows:Constraint>

The feature I am trying to display:

<wfs:FeatureType>
  <wfs:Name>Data_ACTIVE</wfs:Name>
  <wfs:Title/>
  <wfs:Abstract/>
  <ows:Keywords/>
  <wfs:DefaultSRS>EPSG:4326</wfs:DefaultSRS>
  <wfs:OutputFormats>
    <wfs:Format>text/xml; subtype=gml/3.1.1</wfs:Format>
  </wfs:OutputFormats>
  <ows:WGS84BoundingBox>
    <ows:LowerCorner>-102.01545715332 28.970645904541</ows:LowerCorner>
    <ows:UpperCorner>-71.021598815918 43.1671752929688</ows:UpperCorner>
  </ows:WGS84BoundingBox>
</wfs:FeatureType>

Here is the code I am using:

$(document).ready(
    function () {
        // allow testing of specific renderers via "?renderer=Canvas", etc
        var renderer = OpenLayers.Util.getParameters(window.location.href).renderer;
        renderer = (renderer) ? [renderer] : OpenLayers.Layer.Vector.prototype.renderers;
        var geographic = new OpenLayers.Projection("EPSG:4326");
        var mercator = new OpenLayers.Projection("EPSG:900913");

        map = new OpenLayers.Map({
            div: "map",
            layers: [
                new OpenLayers.Layer.WMS("OpenLayers WMS",
                    "http://vmap0.tiles.osgeo.org/wms/vmap0",
                    { layers: "basic" }
                ),
                new OpenLayers.Layer.Vector("GML", {
                    strategies: [new OpenLayers.Strategy.Fixed()],
                    protocol: new OpenLayers.Protocol.WFS({
                            url: "http://localhost/MapServer/Default.aspx",
                            featureType: "Data_ACTIVE",
                            featureNS: "http://www.tstgis.org/gml",
                            version: "1.1.0",
                            geometryName: "line"
                        }),
                        renderers: renderer
                    })
            ],
            zoom: 15
        });

        var bb = new OpenLayers.Bounds(-179.821327209473, 12.1057098342161, -56.5289154052734, 78.1442901657839);
        map.zoomToExtent(bb);
});

I see a post request being made to the server for the layer and a result coming back but it is not displayed at all. Why would the layer not show? Any help would be appreciated.

NOTE
This test script is part of an HTML shell that is running under localhost/mapserver, so it rules out the familiar cross-domain issue.


texts:
I have a short python script, called test:

import fiona
with fiona.open("C:\\Essai_compil", 'r') as l_Obs_shp:
    for shp in l_Obs_shp:
        print shp['geometry']

I want to compile with py2exeor PyInstaller. OS: Windows7, Python 2.7 (32-bits),
I've got erros. With py2exe, the 'log.txt' said :

Traceback (most recent call last):
File "test.py", line 1, in <module>
File "fiona\__init__.pyc", line 78, in <module>
File "fiona\collection.pyc", line 7, in <module>
File "fiona\ogrext.pyc", line 12, in <module>
File "fiona\ogrext.pyc", line 10, in __load
ImportError: DLL load failed: %1 n’est pas une application Win32 valide.
ImportError: DLL load failed: %1 is not a valid Win32 application.

With PyInstaller, log says:

File "ogrext.pyx", line 12 in init fiona.ogrext (src/fiona/ogrext.c:19107)
ImportError: cannot import name ogrinit

Problem are related to fiona module or some ogr stuff... I supposed this is not a 32/64 bits problem; all the modules I've used are 32-bits. More surely a DLL problem... Any ideas ?


texts:
I have a QgsLayoutAtlas,  I need to capture the geometry of each sheet within a  loop
project = QgsProject.instance()
layout = project.layoutManager().layoutByName('Layout')
atlas=layout.atlas ()
hojas=atlas.count ()
i=0
for i in range(hojas):
    atlas.seekTo(i)
    #HERE I NEED THE GEOMETRY OF @atlas_geometry
    i=+1

Thanks


texts: Converting polygon to lines without duplicate edges?

texts:
This seems like it should be easy, but I'm not figuring it out: How do I convert a polygon layer into a line layer representing the polygon edges, where each edge belongs only to a single feature (or part)? When I convert from polygons to lines, each edge is represented twice (because each edge belongs to two adjacent polygons). Even if I dissolve the features, the duplicates remain as separate parts, so I cannot just dissolve, then split into singleparts. This causes problems with symbolization, for example if I want to use a dashed line style.

My thought was to split each polygon boundary into separate edge segments, as shown in these two questions:
Split polygon into lines
How do I split polygons into line segments?

…and then to remove duplicates, for example using v.clean with rmdupl. But I can't figure out how to do either of those steps.

For splitting segments, either the answer isn't in those two linked threads, or else I'm not understanding the explanation. Using v.split doesn't work, because I can only select to split by length, or by number of vertices, not by intersections with other line segments.

For removing duplicates, rmdupl isn't working; the duplicate edges always remain after I run the algorithm. It seems that the edges aren't recognized as duplicate geometries, even though they have identical points. (Running snap first didn't help; the layer appears to be topologically correct already.)

What am I missing?


texts: How to fix overlapping polygons using ArcGIS for Desktop?

texts:
I'm using ArcMap 10.3.

I have a layer of different land use types. After calculating the percentage of each land use in the catchment using tabulate intersection, I found out that the percentage of total land use is larger than 100% which indicates that polygons are overlapping each other. I managed to find all overlapping polygons using intersect. The result of intersect is shown below.

Could you please let me know what is the best way to fix all these overlap?


texts: Spatio-temporal interpolation in R or ArcGIS?

texts:
I'm trying to calculate average rainfall value from a number of points using Inverse Weighted Distance tool in ArcGIS 9.3.

My problem is that: each point has it own time series, therefore the interpolation process should be able to carry out for all years (kind of iteration so to speak).

Following is a sample attribute table:

ID X Y Name Rain1990 Rain1991 Rain1992 Rain1993 .... Rain2010

1 xx1 yy1 AA 1210 1189 1863 1269 ......
2 xx2 yy2 BB 1492 1502 2187 1923 ......
......

Could anybody show me how to do that?

Edit 1: I finally did this by using C++ code which required ArcGIS mask grid, data files & locations of all the points.

Edit 2: I recently used R to do this interpolation task. You can use either hydroTSM, gstat or spacetime packages. Few example links below:

http://spatial-analyst.net/wiki/index.php?title=Spatial_interpolation_exercises_%28NL%29

http://www.geostat-course.org/Topic_Bivand_2012

Edit 3: Added a working example below for future readers