texts: Producing a high-res jpeg with grid from GEOTIFF and Shapefile in Python

texts:
I have a bunch of 4-band GeoTIFFs that I would like to save as high resolution JPEG images. But, in the output image, I also need a separate point .shp to be overlayed and labeled with one of the fields in the shapefile. Additionally, I would like to have a coordinate grid overlayed on the image.

I have attached a screenshot (the 4-meter JPG is too large to share) of a subset (upper left corner) of an example that was generated using ENVI's QuickMap tool, which gives you the option of adding a coordinate grid. The points were overlayed using Arc, I believe.

For each geoTIFF (input), there will be two output jpgs. One with a 4-3-2 (meaning, 4th band = Red, 3rd band = Blue, 2nd band = Green) color composite and another with a 3-2-1 color composite, each with the coordinate grid and labeled points overlayed. ENVI's QuickMap tool works well to convert a geotiff raster to a color composite JPG map with the coordinate grid overlay, but I need to find a way to automate the entire process (including adding the points/labels) as there will be hundreds to thousands of inputs.

I can use GDAL, ENVI/IDL, QGIS or ArcGIS, but if there's a way to do it all in GDAL that would be great. If necessary, I can split the process up to do different steps with different tools.


texts: Setting up CRS for Mars images taken by a satellite

texts:
I have a PDS image of the Valles Marineris taken by the Mars Orbiter Missions (MOM) Mars Colour Camera (MCC). I have to generate a depth map and I am in the process of creating a 3D model.
In the provided metadata I have latitude and longitudes of 9 reticle points on the given image.
When I try to set the provided CRS I am prompted with the following error:

No transform is available between EPSG:4326 - WGS 84 and Unknown CRS: GEOGCRS["unknown",DATUM["unknown",ELLIPSOID["unkno….
proj_create_operations: Source and target ellipsoid do not belong to the same >celestial body

The CRS that I tried were the following:

Mars_2000_(Sphere) ESRI:104971
GCS_Mars_1979 ESRI:104904.

Are there any suggestions for setting up a CRS or any other way to generate a depth map?


texts: Calculating distance between geographical centroids and points?

texts:
I have one zone system and two layers. The 1st layer is the geographical centroids of the zone system (each zone has 1 geographical centroid), while the 2nd layer displays the rail stations. I would like to calculate and display the distance between the centroids and the rail stations. Ideally, I want each station to be connected to the centroid of the zone that it belongs (see pic attached).

In essence, it is the shortest distance, but in some cases, one rail station might be closer to a geographical centroid of another zone and QGIS will give me this distance, which I don't want. In addition, one zone might have more than 1 rail stations. In that case, I want only the shortest distance.

Any idea on how to perform this task?


texts:
I have to compare geometries anywhere in the world; mainly between 65N and 40S but also in the Arctic latitudes. I need "commercial accuracy" - not military or navigational. For example does a line intersect with a polygon and - if so - so they intersect for 50 km rather than 55 km. I store the lat and lng of  these geometries.

Can anyone please tell me how to select the most suitable cartesian projection (from spatial_ref_sys) to transform the coordinates into before compare any line with any polygon ?

I have tried to use the lng of the centre of the line to select one of the SRIDs between 32601-32660 and 32701-32760. As many of you will know, these are based on the equator every 6 degrees and the 601-660 are for northern hemisphere while the 701-760 are for the southern hemisphere.

To test whether this works, I check whether the ST_Intersects() function detects that a point in Asia at 36.9417, 75.4083 intersects with a polygon of Asia after both are transformed with SRID 32643 whose central_meridian is 75 degrees. Sadly the ST_Intersects() function doesn't detect that the point is right in the middle of Asia, but when I use SRID 32644 whose central_meridian is 81 degrees, ST_Intersects() function detects an intersection between point and polygon.

There are many other examples of my "calculated" SRID being the wrong one. Maybe I'm overcomplicating ? I originally used SRID 4326 anywhere worldwide before testing with ST_Intersects(), and things seem fairly accurate except in high northern latitudes, so I was hoping to minimize inaccuracies by using different cartesian SRIDs at different longitudes.

Perhaps I need to find SRIDs which centre on the latitude of the geometries - not just the longitude ? Maybe I'm overcomplicating ?

Many thanks for your initial welcoming response. I hope that gives good information.


texts: ArcGIS Desktop equivalent to GRASS r.geomorphon algorithm?

texts: I would like to know if in ArcGIS Desktop 10 there is some tool/algorithm which can be used to perform the same analysis that r.geomorphon does in GRASS 7.


texts:
I am an ArcGIS Desktop user and a layman to cloud computing. Our GIS data is currently hosted on a server and recently our IT dept started to look at possibilities to get rid of this server and move all the GIS data (vector & raster) to Amazon cloud. I am not sure how ArcGIS can connect to the data in cloud.

Is it an easy setup or does it require extra stuff?

I googled a bit and most articles I found are related to ArcGIS server. We don't have ArcGIS Server here.


texts:
I am an archaeologist, and I am trying to add a geologic cross-section to my ArcScene map, but I have so far been unable to do so. I have created two (basic) cross-sections using points and lines (see picture), but I have unable to convert this into a polygon so I can display the profile. I have tried to use the 3D editor, but It will not create a vertical polygon and, I have also tried the line to polygon tool, but that creates an empty polygon.
.


texts: Creating line in mean position respect to position of other lines in QGIS

texts:
I have a series of shapefiles (let's say 3 shape files), each of which contains a line representing the coastline during a certain period of prehistory.

red line = coastline when sea level was -16 m than now
green line = coastline when sea level was -14 m than now
blue line = coastline when sea level was -5 m than now

The attributes table of the first shape file (which is referred to the moment in which sea level was 14 m lower than now) is made of 3 columns and 5369 rows and looks like this:

ID    GRIDCODE    Sea level
2247  0           -14
2244  0           -14
2245  0           -14
2241  0           -14
2243  0           -14

The attributes table of the second shape file (which is referred to the moment in which sea level was 5 m lower than now) is made of 3 columns and 6285 rows and looks like this:

ID    GRIDCODE    Sea level
1622  0           -5
1619  0           -5
1620  0           -5
1616  0           -5
1618  0           -5

The attributes table of the third shape file (which is referred to the moment in which sea level was 16 m lower than now) is made of 3 columns and 5382 rows and looks like this:

  ID    GRIDCODE    Sea level
  2352  0           -16
  2347  0           -16
  2348  0           -16
  2344  0           -16
  2345  0           -16

What I'd obtain is a new shapefile of a line representing the mean position of the coastline, calculated on the basis of the three initial shapefiles.


texts: Adding Layer package to Web Map in ArcGIS Online?

texts:
I am using ArcMap 10.4.1, and created a layer package of a raster file, but it does not allow me to add it to a web map in ArcGIS . It shows up in "My Content", but only allows me to download it or view it in ArcMap.

It is 6MB in size.

If there is some kind of size limit, how can I load this into a web map?

It is reclassified data from the Hansen Forest cover loss dataset. I also published a web service of the map, but it still doesn't display the raster data.


texts:
We are seeing .prj files from our users that don't look like the normal WKT formatted ones we've seen up to this point, instead, they match the format illustrated at "Esri Grid format".  When we attempt to use these files in GeoTools, it fails to parse these prj files - it expects them to be WKT.

An example that looks like NZ Transverse Mercator looks like this:

Projection    TRANSVERSE
Datum         AI_NZGD00
Spheroid      GRS80
Units         METERS
Zunits        NO
Xshift        0.0
Yshift        0.0
Parameters
0.9996 /* scale factor at central meridian
 173  0  0.0 /* longitude of central meridian
   0  0  0.0 /* latitude of origin
1600000.0 /* false easting (meters)
10000000.0 /* false northing (meters)

A quick look through some GDAL source says that these might be 'old style .prj files' from Arc/Info version 7 - but maybe they're still being proliferated for some reason?  Is there an option in newer versions of Arc to output the old version?

Or maybe there's some GeoTools code that can read these files, or a way to get gdal to convert them?

UPDATE I started trying to port the gdal code to use with GeoTools, but it looked like a total headache.  I've ended up using GDAL JNI to parse the prj file with a kludgy but workable solution.  The importFromESRI method doesn't look usable as it seems to want a pointer to some text, so instead I create an almost empty ascii grid and create a sidecar prj file containing the esri prj data I want to parse.   I do this like so:

  public CoordinateReferenceSystem parsePrj(String prjContents) {
    try {
      // prepare an asc file with the prj sidecar filled with the prj data to parse
      File asc = File.createTempFile("dummy", ".asc");
      File prj = Paths.get(asc.getParent(), Files.getNameWithoutExtension(asc.getName()) + ".prj").toFile();
      Files.write(dummyAscContents.getBytes(), asc);
      Files.write(prjContents.getBytes(), prj);

      // We can't get at `SpatialReference_ImportFromESRI(long, SpatialReference, Vector)` directly, but it'll be used
      // when opening a bitmap dataset - so we open a tiny dataset with a prj file loaded with the required contents...
      Dataset dummyDataset = gdal.Open(asc.getAbsolutePath(), gdalconstConstants.GF_Read);

      if (dummyDataset == null) {
        // not an expected user facing error, no need for i18n
        throw new RuntimeException("GDAL configuration error - could not open dummy asc");
      }
      // ... so we get get the wkt out from here
      String wkt = dummyDataset.GetProjection();

      if (Strings.isNullOrEmpty(wkt)) {
        throw new RuntimeException("could not parse prj - " + prjContents, e);
      }

      return CRS.parseWKT(wkt);

    } catch (Exception e) {
      ...
    }
  }


texts: Enclosing building models from bottom (Missing GroundFeatures) using FME?

texts:
I have 3D building model of the city and I would like to export it to proper CityGML model.

I am trying to find out, how to make an intersection of building and terrain in order to make "TerrainIntersectionCurve".

So literally, I am looking for some steps how to make intersection of terrain and building models using FME?


texts:
I am trying to write SQL statement in WHERE Clause FME in Navigator but it does not work.

WHERE Clause:

NAME IN SELECT CONCAT(''$(',SELECT REPLACE ('Roman Ann Jack',' ',','),')'')

This sql code suppose to be NAME IN '$(NAME)' but NAME is actually list of names derived from published parameter defined as Multi Choice Value. I think I have to replace EMPTY CHARACTER between names to comma character and that should be recognizable as list for FME SQL?!
I asked this question before, and got answers that this suppose to work in FME2012 or by writing python script.

But at the moment I am curious about this sql statement and if it is syntax and logic ok?!
Can anyone check if whole or just part

SELECT CONCAT(''$(',SELECT REPLACE ('Roman Ann Jack',' ',','),')'')

is working?

Thank you very much


texts: PostGIS Rasters now loading differently with QGIS 3.14+

texts:
In QGIS version 3.10 when loading in a raster from my PostGIS DB it worked fine and when hoovering over for info gave this message

Now with version 3.16(anything over 3.14) I am getting this when on info hoover

For some reason the table column now contains 'schema'.'table' followed by the contents of my table column which says 'rast' / (rast)
When loading in this data it now has large black squares where blank 'no data' values used to be in 3.10 and also when clipping this dataset down to a smaller area I am getting the following error message about a not recognised dataset name.
GDAL command: gdalwarp -dstnodata -99.0 -r near -of GTiff "dbname='Data' host=10.0.0.17 port=5432 user='David' password='password' sslmode=disable srid=0 table=\"1m_LiDAR\".\"ardstraw_02_02_2012\" (rast)" C:/Users/david.laverty/AppData/Local/Temp/processing_kCAqKD/2020d69a712c4b7680bc3e15de2f4b9d/OUTPUT.tif GDAL command output: ERROR 4: dbname='Data' host=10.0.0.17 port=5432 user='David' password=password sslmode=disable srid=0 table="1m_LiDAR"."ardstraw_02_02_2012" (rast)' does not exist in the file system, and is not recognized as a supported dataset name.
Does anyone know how to fix this issue or revert back to the old way of reading the data as to not give these errors when using Version 3.14+. As this is I can no longer use any of my rasters stored in my DB and would like to know if anyone else has experienced this or can offer me a solution to the issues I am facing
System is running on Windows 10 and was working fine the past 2 years until these most recent version updates starting with 3.14


texts: Raster to ASCII tool in ArcMap returning different max and min values

texts:
I am using Raster to ASCII conversion tool in ArcMap 10.7.1, and the issue is reproducible in 10.6.1, and 10.8.1.
I have two raster files in TIFF format.
When I used the tool for one raster the output .asc file returned the exact max and min values as that of the raster, but for another .tif the converted .asc file returned inaccurate max and min values.
The correct input and output converted .asc

Here we can see the difference in the max and min values

The issue seems to be data specific, but is the deciding factor the radiometric resolution or something else?


texts: Labeling classes for small roads on large extent maps in ArcMap with SQL statement?

texts:
I have many small sub divisions that won't label nice, (lot of overlapping) unreadable.

While using the bullet leader symbol.

What would be the best statement that would place the labels in a clear spot near the sub division that would point to the road and also look clear and aligned horizontal?


texts: Direct connection from geotools to arcgis 10.3

texts:
Is it possible to direcly connect to ESRI 10.3 database using geotools ?

I use geotools 13.0 and have ESRI database set up on oracle.
Previously I was able to connect to SDE using gt-arcsde-plugin but since ESRI 10.3 there is no SDE service.
Whats more, there is no jsde_sdk.jar provided which is required by the geotools arcsde plugin.

I was trying to set up direct connection - but did not succeed. I get 'ClassNotFoundException' pointing that there is no: com.esri.sde.sdk.client.SeException.


texts: Postgres: Remove Duplicate Rows returned by ST_Dwithin query

texts:
I am having trouble finding a query that returns one row for every unique vehicle crash that is within 100 feet of a select group of road features. My query currently returns a row for every collision that is within 100 feet of these roads, but because some crashes occurred within 100 feet of multiple road features, there are duplicate rows for some crashes. How can I remove these duplicates? I only one one row for each collision.

So far I have:

SELECT crashes.* , roads.geom, ST_Distance(crashes.geom, roads.geom)
FROM crashes, roads
WHERE ST_Dwithin(crashes.geom, roads.geom, 100)


texts: R and osmar : the doc's example does not work for me, "Opening and ending tag mismatch"

texts:
I can't make osmar work anymore. I tried to reproduce the documentation's example, and even that produces a very long error code. I could not find the same error anywhere. Any idea ?

The code is :

require(osmar)
api <- osmsource_api()
box <- corner_bbox(11.579341, 48.15102, 11.582852, 48.1530)
gschw <- get_osm(box, source = api)

The error is :

  Opening and ending tag mismatch: BR line 6 and FONT
  AttValue: " or ' expected
  attributes construct error
  Couldn't find end of Start Tag TABLE line 9
  AttValue: " or ' expected
  attributes construct error
  Couldn't find end of Start Tag FONT line 11
  Opening and ending tag mismatch: BR line 14 and FONT
  Opening and ending tag mismatch: BR line 13 and TD
  Opening and ending tag mismatch: TD line 10 and TR
  AttValue: " or ' expected
  attributes construct error
  Couldn't find end of Start Tag HR line 19
  AttValue: " or ' expected
  attributes construct error
  Couldn't find end of Start Tag HR line 24
  Opening and ending tag mismatch: BR line 23 and FONT
  Opening and ending tag mismatch: BR line 22 and TD
  Opening and ending tag mismatch: BR line 20 and TR
  AttValue: " or ' expected
  attributes construct error
  Couldn't find end of Start Tag FONT line 34
  EntityRef: expecting ';'
  Opening and ending tag mismatch: br line 36 and FONT
  Opening and ending tag mismatch: br line 36 and i
  Opening and ending tag mismatch: br line 36 and B
  Opening and ending tag mismatch: BR line 35 and TD
  Opening and ending tag mismatch: B line 33 and TR
  Opening and ending tag mismatch: i line 33 and TABLE
  Opening and ending tag mismatch: TD line 32 and blockquote
  Opening and ending tag mismatch: TR line 32 and FONT
  Opening and ending tag mismatch: FONT line 18 and BODY
  Opening and ending tag mismatch: TD line 17 and HTML
  Premature end of data in tag TR line 17
  Premature end of data in tag TR line 10
  Premature end of data in tag blockquote line 8
  Premature end of data in tag FONT line 5
  Premature end of data in tag BODY line 4
  Premature end of data in tag HTML line 1
  Error: 1: Opening and ending tag mismatch: BR line 6 and FONT
  2: AttValue: " or ' expected
  3: attributes construct error
  4: Couldn't find end of Start Tag TABLE line 9
  5: AttValue: " or ' expected
  6: attributes construct error
  7: Couldn't find end of Start Tag FONT line 11
  8: Opening and ending tag mismatch: BR line 14 and FONT
  9: Opening and ending tag mismatch: BR line 13 and TD
  10: Opening and ending tag mismatch: TD line 10 and TR
  11: AttValue: " or ' expected
  12: attributes construct error
  13: Couldn't find end of Start Tag HR line 19
  14: AttValue: " or ' expected
  15: attributes construct error
  16: Couldn't find end of Start Tag HR line 24
  17: Opening and ending tag mismatch: BR line 23 and FONT
  18: Opening and ending tag mismatch: BR line 22 and TD
  19: Opening and ending tag mismatch: BR line 20 and TR
  20: AttValue: " or ' expected
  21: attributes construct error
  22: Couldn't find end of Start Tag FONT line 34
  23: EntityRef: expecting ';'
  24: Opening and ending tag mismatch: br l


texts: Extracting regions from NDVI using range of values in Google Earth Engine

texts:
I want to identify locations in an area of interest where the NDVI values are maximum and to achieve that I need to extract regions in NDVI image based on the range of values i.e 0-0.5 classA and 0.51-1 classB. So the output image will have two classes(calss A and B) and many regions based on the connectivity of the similar class pixels.

Following figure illustrates what I want to achieve.

At first there is the image with raw values
Image is classified into two classes based on the given ranges
Then three regions are detected in image based on the connectivity of the similar class pixels
Finally the pixels with the highest NDVI are extracted.

I need to automate this process for a very large therefore I am using Google Earth Engine.
I have searched about this kind of regioning a lot but could not find anything specific to this.
Is is possible to do this in Google Earth Engine or any function available in GEE or will I have to extract the data for the region of interest and then do all this by own?


texts: Transparency feature on my QGIS 3.10 is Opacity

texts: When adding a raster layer to a vector layer, I am unable to make the top layer transparent. On getting into Layer Panel properties and checking it shows Opacity slider. Adjusting opacity at anything between 0 - 30% does not help. How do I manage this?


texts: Fastest way to join points table to polygons tables in ORACLE

texts:
I have two tables in 11g oracle db, one with points geometry, one with polygons. Currently I need to enrich the points table daily with the name of the polygon (country) that the point is inside, but what I'm using takes way too long. How can I improve performance? Here's what I'm doing:

    BEGIN
    FOR rec IN (
        SELECT country_name, SHAPE
        FROM territory
        ORDER BY country_name
        )
    LOOP
        UPDATE points_table
        SET Country = rec.country_name
        WHERE (Country is null) AND (SDO_INSIDE(Geometry, rec.SHAPE) = 'TRUE';
    END LOOP;
END
;


texts:
I have GIS data in a PostGIS data base and would like to generate a density map and display them on a Leaflet/Google Maps.

Question: What is the recommended/typical method for doing so, without transferring ALL (millions) the data points within the map's viewpoint to the Leaflet/Gmaps map? No user intervention (ie. Import data into ArcGIS to generate the map) should be required unless automation is possible via CLI etc.

Currently using Postgresql/PostGIS as the spatial database and Leaflet on the browser.


texts: QGIS Points in Polygon Query Not Working

texts:
I'm trying to count points within a polygon but it's not working.

I think it has something to do with the co-ordinate system.

The points are from a text file and the line are from a dxf file. I had to set the co-ordinates for both as one was lat long and the other was an old Irish grid co-ordinate system.

When I buffered and dissolve the lines the resulting layer defaulted to the original co-ordinates.

I reset the line file to the correct co-ordinates and ran the query to count points within a polygon. The resulting layer defaulted to its original co-ordinates and all the values were null.


texts: I'm looking at the transparency window of the layer properties of my "Heat map". How does the No Data Value effect a raster layer.  I'm assuming an initial setting of -9999 sets it out of the way, so when should I consider setting it.


texts: I'm adding an id in QGIS with the function @row_number and it is adding the numbers descending, as you see it in the picture. But i need them exactly in the reverse order.


texts: How to save selected features from a postgis database into a new mini database

texts: I have a postgis(v2.1.3)/postgresql(v9.3) database with geometry (geom1, geom2, etc.) and attribute (attr1, attr2, attr3, etc.) tables.  The attribute tables are linked to the geometry tables using foreign keys.  I use QGIS to view all these tables.  I want to select features in QGIS (either using a rectangle, polygone or manual select) from the map and then create a mini-database using only these selected features (from the geometry tables) and the relevant records from the attribute tables.  Is this possible?  Ideally I want the mini-database to be in spatialite/sqlite but I can find a way to convert from postgis to spatialite using ogr2ogr.


texts: What is the difference between Aspect and Flow Direction tools?

texts:
There are two seemingly identical Spatial Analysis tools:

(1) Surface -> Aspect:

  Derives aspect from a raster surface. The aspect identifies the
  downslope direction of the maximum rate of change in value from each
  cell to its neighbors.

  Aspect can be thought of as the slope direction. The values of the
  output raster will be the compass direction of the aspect.

(2) Hydrology -> Flow direction:

  Creates a raster of flow direction from each cell to its steepest
  downslope neighbor.

Manual comparison

To test this I created an aspect layer and a flow direction layer, converted the aspect layer to the format of flow direction using Raster Calculator:

64*(("aspect" >= 0)  &  ("aspect" < 22.5)) +
128*(("aspect" >= 22.5) & ("aspect" <67.5)) +
1*(("aspect" >=67.5) & ("aspect" <112.5)) +
2*(("aspect" >=112.5) & ("aspect" <157.5)) +
4*(("aspect" >=157.5) & ("aspect" <202.5)) +
8*(("aspect" >=202.5) & ("aspect" <247.5)) +
16*(("aspect" >=247.5) & ("aspect" <292.5)) +
32*(("aspect" >=292.5) & ("aspect" <337.5)) +
64*(("aspect" >=337.5) & ("aspect" <=360))

and then compared to the flow direction layer:

("aspect_conv" == "flowdir")

The vast majority of the resulting layer is 1 (true) and I suppose the few areas with 0 is because of flat aspects or nulls.

It therefore appears that the only difference is that "Flow Direction" returns unique values (powers of 2), while aspect returns a continuous range of floats. Is that all?


texts: Lizmap URL highlight and zoom to feature on layer

texts:
I would like to pass an id or searchstring and a layername in the URL when calling Lizmap, so that Lizmap highlights the feature on the specified layer and zooms to it. I read a lot about all the cool features Lizmap aready offers, but I didn't find anything in the well written documentation.
I was thinking about something in that direction:
https://157.xxx.xx.xx/lizmap/index.php/lizmap/service/?repository=lizmap&project=myproject&layer=mylayer&id=feature_id
Does anybody know if its possible or if there is a workaround?


texts:
I have two grids.  Cells in each grid have an attribute.  Size of grid cells differs between the two grids and the two grids do not overlay each other perfectly.  I would like to create a new layer, or group of polygons, or grid, that results from joining the two existing grids.  Each cell in the new grid should retain the attribute from each of the original grids.  Area not occurring in both original grids can be discarded from the new grid.  Such discarded areas would only have one attribute (and one missing observation for the second attribute).

How can I do this is R?

So far I have tried over in the sp package and unionSpatialPolygons.  I have also considered the PBSmapping package, which allows for creating intersections of polygons.  This post describes how to join polygons within a layer, but perhaps not how to join polygons in different layers: Joining polygons in R

I have written code to compute the (x,y) coordinates of the four corners of each cell that results from joining the two original grids.  That code is at the very bottom of this post. Perhaps I can create a new polygon layer from the resulting data frame.  However, I suspect that process has already been implemented in a package.

Below is code to create the two original grids, assign an attribute to each cell in each grid and plot the two grids.

There will be 16 cells in the new grid.  Those cells generally will not be squares, but rather rectangles of various shapes and sizes.

library(rgdal)
library(maptools)
library(raster)
library(rgeos)
library(gridExtra)

setwd('c:/users/mmiller21/simple R programs')

set.seed(1234)

# create grid 1

grd1 <- GridTopology(c(2,2), c(1,1), c(4,4))
polys1 <- as(grd1, "SpatialPolygons")

# assign projection to grid

proj4string(polys1) = CRS("+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0")

# create fake atttribute data for each grid cell

poly.data1 = data.frame(f1=runif(length(row.names(polys1)), 1, 10))
row.names(poly.data1) <- paste0('g', 1:length(row.names(polys1)))

# convert grid to a SpatialPolygonsDataFrame:

poly.df1 = SpatialPolygonsDataFrame(polys1, poly.data1)

# create grid 2

grd2 <- GridTopology(c(3.8,3.8), c(2,2), c(2,2))
polys2 <- as(grd2, "SpatialPolygons")

# assign projection to grid

proj4string(polys2) = CRS("+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0")

# create fake atttribute data for each grid cell

poly.data2 = data.frame(f2=runif(length(row.names(polys2)), 100, 200))
row.names(poly.data2) <- paste0('g', 1:length(row.names(polys2)))

# convert grid to a SpatialPolygonsDataFrame:

poly.df2 = SpatialPolygonsDataFrame(polys2, poly.data2)

jpeg(filename = "joined.grid.jpeg")

plot(poly.df2, col='blue', xlim=c(1,8), ylim=c(1,8))
text(coordinates(poly.df2), labels=row.names(poly.df2))
axis(1, lty = 0, lwd = 0)
axis(2, lty = 0, lwd = 0)
plot(poly.df1, add=TRUE)
text(coordinates(poly.df1), labels=row.names(poly.df1))

dev.off()

Here is code to compute the (x,y) coordinates of each cell in the new grid:

x1 <- seq(1.5,5.5,1)
y1 <- seq(1.5,5.5,1)

x2 <- seq(2.8,6.8,2)
y2 <- seq(2.8,6.8,2)

x <- c(x1,x2)
x <- sort(x)
x
# [1] 1.0 2.0 2.8 3.0 3.8 4.0 4.8 5.0 5.8 6.8

new.x <- x[x>= min(x2) & x <= max(x1)]

y <- c(y1,y2)
y <- sort(y)
# [1] 1.0 2.0 2.8 3.0 3.8 4.0 4.8 5.0 5.8 6.8

new.y <- y[y>= min(y2) & y <= max(y1)]

z <- merge(new.x, new.y)

plot(z$x, z$y)


texts: How to Export NDVI image using Fusion Table as featurecollection boundary in GEE python API

texts:
I am using python api of google earth engine to create NDVI raster. Now I want to export the raster in my google drive. But during exporting it is giving the following error:

  ee.featurecollection.FeatureCollection object at 0x03D29EF0> is not
  JSON serializable

I am using the following piece of code for the entire work

import datetime
import ee
import ee.mapclient
import json
from ee.batch import Export
ee.Initialize()

boundary = ee.FeatureCollection('ft:1rbhNtC1TqDBvY9Rt2BZR-DjhpIPuC3nU5kmz49WW')

image = (ee.ImageCollection('COPERNICUS/S2').
         filterBounds(boundary).
         filterDate(datetime.datetime(2018,01,24),
                     datetime.datetime(2018,02,28)).
         sort("CLOUD_COVERAGE_ASSESMENT").
         limit(20))

mosaic = image.mosaic();
cloudBitMask = ee.Number(2).pow(10).int();
cirrusBitMask = ee.Number(2).pow(11).int();

def maskS2clouds(image): #(
  qa = mosaic.select('QA20');
  #// Both flags should be set to zero, indicating clear conditions.
  mask = qa.bitwiseAnd(cloudBitMask).eq(0)and(qa.bitwiseAnd(cirrusBitMask).eq(0));
  return image.updateMask(mask);

cloudMasked = image.filterBounds(boundary).map(maskS2clouds)
median = cloudMasked.median()

clip = median.clip(boundary)

ndvi = clip.normalizedDifference(['B8', 'B4'])

task_config = {
    #'image:': ndvi,
    'description': 'PythonExample',
    'scale': 30,
    'region': boundary
    }

task = ee.batch.Export.image(ndvi, 'nd', task_config)

Please help me how I can covert FeatureCollection as json object? My per my understanding from the error I might have to convert the fusion table FeatureCollection into JSON object


texts:
I have a 6 acre permaculture farm, and I would like to be able to plant something at a location, log the location, and be able to find it again to within about a foot accuracy.

I do not need the position of these plants globally:  I only need them relative to each other.

Is there any cheap way to accomplish this?  For instance, could I put some kind of transmitter on two fence posts (at fixed locations), and triangulate my position relative to them in a cheap way?  Are there off the shelf solutions for this?

My apologies if this is not the correct place for this question.


texts: Building Python add-in button to calculate elevation change

texts:
I'm building a python add in button that calculates elevation change. I want the user to be able to click the button, click one cell on a DEM, click another cell, and then display a message that has the calculated elevation change between the two cells. Below is my code. When I install it and click on the newly created button in ArcMap nothing happens. My cursor doesn't even change to the one I selected.

Can anyone tell me what my code is missing?

import arcpy
import pythonaddins

class CollectElevations(object):
    """Implementation for ElevationDifference_addin.CollectElevations
(Button)"""
    def __init__(self):
        self.enabled = True
        self.checked = True
        self.shape="None"
        global clickCount
        self.cursor=3
        clickCount=0

    def onMouseDownMap(self,x,y,button,shift):
        global clickCount
        rasters = arcpy.ListRasters()

        if clickCount==0:
            for raster in rasters:
                self.result=arcpy.GetCellValue_management(raster, "{} {}".format(x, y), "1").getOutput(0)
                clickCount=1

        elif clickCount==1:
            for raster in rasters:
                self.result2=arcpy.GetCellValue_management(raster, "{} {}".format(x, y), "1").getOutput(0)
                message = "Elevation difference is"+str(" ")+str(self.result-self.result2)
                pythonaddins.MessageBox(message, "Results:")
                clickCount=0


texts: Python Error with QGis2Web on QGIS 3.6.3 Noosa

texts:
I'm receiving a Python Error when trying to create a web map using QGIS2Web and for the life of me cannot figure out what to do.

Running QGIS 3.6.3 Noosa
Windows 10 Professional

KeyError: 100
Traceback (most recent call last):
  File "C:/Users/Rodney/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\qgis2web\qgis2web.py", line 60, in run
    self.dlg = MainDialog(self.iface)
  File "C:/Users/Rodney/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\qgis2web\maindialog.py", line 156, in __init__
    self.autoUpdatePreview()
  File "C:/Users/Rodney/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\qgis2web\maindialog.py", line 320, in autoUpdatePreview
    self.previewMap()
  File "C:/Users/Rodney/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\qgis2web\maindialog.py", line 323, in previewMap
    preview_file = self.createPreview()
  File "C:/Users/Rodney/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\qgis2web\maindialog.py", line 286, in createPreview
    dest_folder=utils.tempFolder()).index_file
  File "C:/Users/Rodney/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\qgis2web\olwriter.py", line 91, in write
    folder=dest_folder)
  File "C:/Users/Rodney/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\qgis2web\olwriter.py", line 130, in writeOL
    popup, json, restrictToExtent, extent, feedback, matchCRS)
  File "C:/Users/Rodney/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\qgis2web\utils.py", line 220, in exportLayers
    iface, extent, precision, crs, optimize)
  File "C:/Users/Rodney/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\qgis2web\utils.py", line 233, in exportVector
    cleanLayer = writeTmpLayer(layer, restrictToExtent, iface, extent)
  File "C:/Users/Rodney/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\qgis2web\utils.py", line 162, in writeTmpLayer
    uri = TYPE_MAP[layer.wkbType()]
KeyError: 100

Python version: 3.7.0 (v3.7.0:1bf9cc5093, Jun 27 2018, 04:59:51) [MSC v.1914 64 bit (AMD64)]
QGIS version: 3.6.3-Noosa Noosa, 0c5774c068

Python Path:

C:/PROGRA~1/QGIS 3.6/apps/qgis/./python
C:/Users/Rodney/AppData/Roaming/QGIS/QGIS3\profiles\default/python
C:/Users/Rodney/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins
C:/PROGRA~1/QGIS 3.6/apps/qgis/./python/plugins
C:\Program Files\QGIS 3.6\bin\python37.zip
C:\PROGRA~1\QGIS 3.6\apps\Python37\DLLs
C:\PROGRA~1\QGIS 3.6\apps\Python37\lib
C:\Program Files\QGIS 3.6\bin
C:\PROGRA~1\QGIS 3.6\apps\Python37
C:\PROGRA~1\QGIS 3.6\apps\Python37\lib\site-packages
C:\PROGRA~1\QGIS 3.6\apps\Python37\lib\site-packages\win32
C:\PROGRA~1\QGIS 3.6\apps\Python37\lib\site-packages\win32\lib
C:\PROGRA~1\QGIS 3.6\apps\Python37\lib\site-packages\Pythonwin
C:/Users/Rodney/AppData/Roaming/QGIS/QGIS3\profiles\default/python
C:\Users\Rodney\AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugins
C:/Users/Rodney/Google Drive/QGis


texts: Draw postgis geometries in google maps api

texts:
Postgis have a base with a table of data lines and other points.
I'm doing a java applet to query these tables and space to draw on google maps.
My question is: what better way to spend these geometries to google maps? When I give the select the bank which is the best format to return the geomtrias to draw on the map, I had thought of GeoJSON.
Does anyone have any example code or article addressing postgis + google maps api?


texts: How to create a buffer on a linestring with varying or increasing/decreasing width

texts:
I want to buffer rivers/streams for a topographic map. I'm using a PostgreSQL/PostGIS database. Usually, I create buffer with ST_Buffer with a constant width.

For example, the width of a river buffer should be increasing from the spring to the estuary with a distinct value.

Is there a way to create buffer with an increasing or decreasing width starting from a distinct point and ending at a distinct point?

UPDATE

With the help of simplexico I tried to create some geometries that could help to create a thinning buffer around a linestring:

First I create the buffer lines parallel to the river geom and insert them into a table called offsetcurve:

WITH geometries AS (
    SELECT
        gid,
        (ST_Dump(geom)).geom AS geom
    FROM waterways
        WHERE name='river_name')

    SELECT
        gid,
        ST_OffsetCurve(geom, 10, 'quad_segs=4 join=mitre') AS pos_curve,
        ST_Reverse(ST_OffsetCurve(geom,-10, 'quad_segs=4 join=mitre')) AS neg_curve
    FROM geometries;

I have to use the reverse function for a further step.

The next step is to dump the points from the buffer lines and the river geom and insert them into a table called offsetpoints:

SELECT
    offsetcurve.gid,
    (ST_DumpPoints(offsetcurve.pos_curve)).geom AS pos_curve_dump,
    (ST_DumpPoints(offsetcurve.neg_curve)).geom AS neg_curve_dump,
    ST_SetSRID((ST_DumpPoints(geometries.geom)).geom, 31468) AS river_dump
FROM offsetcurve, geometries

After that I create lines between the buffer line points and the river geom points and store them in a table called offsetazimuth:

SELECT
    offsetpoints.gid,
    ST_Azimuth(offsetpoints.pos_curve_dump, river_dump) AS pos_curve_azimuth,
    ST_MakeLine(offsetpoints.pos_curve_dump, river_dump) AS pos_azimuth_line,
    ST_Azimuth(offsetpoints.neg_curve_dump, river_dump) AS neg_curve_azimuth,
    ST_MakeLine(offsetpoints.neg_curve_dump, river_dump) AS neg_azimuth_line
FROM offsetpoints

Visualizing all new geometries, I have the following picture:

The red and blue points are dumps from the buffer line.
The green points are the river geom dumps.

The blue and red points are connected to the green points with the offsetazimuth line.

I can change the position of the buffer points (e.g. for one buffer line) along the azimuth line with:

SELECT
    ST_MakePoint(
        ST_X(st_endpoint(pos_azimuth_line)) - (ST_X(st_endpoint(pos_azimuth_line))-st_x(st_startpoint(pos_azimuth_line)))/6,
        st_y(st_endpoint(pos_azimuth_line)) - (st_y(st_endpoint(pos_azimuth_line))-st_y(st_startpoint(pos_azimuth_line)))/6
                ) as geom
FROM offsetazimuth

The results are new points with a constant distance from the river geom:

The first step of a final solution has to be that the first point has the distance like (buffer_width(B)-buffer_width(A))/count(river_dump_points) (see first picture) from the river geom on the offsetazimuth line (see the last but one query).

For example:
The first new point has a distance of (200m-100m)/100=1m from the river geom, the second one ((200m-100m)/100)*2=2m and so on. 200m is the maximum width,100m is the minimum width of the wanted buffer. 100 is the number of dumped points from the river geom.

It's the question if this (increasing or decreasing distance value for every buffer line point) is possible with PostGIS!?

The next idea could be a function: setting parameters like the gradient or the length where the buffer has to take place.


texts: Displaying generated lines out of map extent using QGIS?

texts:
I have drawn an OD flow map using QGIS geometry generator, using the following expression :

make_line( centroid( $geometry ),
centroid( geometry( get_feature( 'Layer2', 'Attribute1', attribute(
           $currentfeature, 'Attribute2') ) ) ) )

I got everything fine, defining the width of my lines according to attributes.

But when I zoom to the "hub", lines linked to places outside of map canvas doesn't show up. I guess it's because the features are not supposed to appear on the map.

I tried to uncheck the "Clip features to canvas extent" option in the advanced menu, both in the "Geometry generator" and in the "Line" levels, but it still doesn't appear.

I would like to have all my lines appearing, even when zoomed on the hub.


texts: Can someone tell me how, or if it's possible, to send street centerline updates to Bing for inclusion in their base maps.


texts: Looking for a way to create dynamically updating point layers in QGIS (linked to spreadsheet)

texts:
I'm looking for a way to create a dynamic link between a QGIS point layer and a spreadsheet that contains point coordinates.
For example, I have a spreadsheet of location points for a species of plants, and I want points to automatically display in QGIS when I find more locations (and add the XY to an excel sheet). Currently, I have to reimport the CSV to a new point layer each time I want to update a map.
Is there already a plugin that exists to do this in QGIS? I know that this is (or at least, was) possible in ArcMap.


texts: Calculate the common area between polygons using python?

texts:
I am trying to calculate the common "empty" area between the neighbouring polygons in order to calculate the porosity volume (the red area in the image)

My idea was, first to use Convert FeatureToVertices then to calculate the shortest distance between each point to the next polygon and then trying to generate the new points (blue points) and finally to generate the common polygon from these points. I could not find the way to write the code and I have only succeeded to do the first two steps. Any idea how to manage this problem?


texts: Using scale and offset values of MODIS products

texts:
I'm working on NDVI data from MODIS_13Q1 product. As far as I know, NDVI data has scale and offset values. I used these values to get geophysical values of pixels.

My problem is after processing these data (linear interpolation and spline fitting missing values), the value of pixels is float, not integer as original data, and the size of processed datasets when saved is huge (from 100 MB to 700-800 MB).
How can I use scale and offset values to convert value of pixels back to integer, and does it cause losing data?


texts: CompositeSymbol that includes TextSymbols and SimpleFillSymbol does not draw map

texts:
Following is the code snippet that does not draw both TextSymbols.

private function MapClick(evt:MapMouseEvent):void
        {
            var cs:CompositeSymbol = new CompositeSymbol();
            var textSymbol1:TextSymbol = new TextSymbol("hello");
            textSymbol1.textFormat= new TextFormat("Arial",5);
            var textSymbol2:TextSymbol = new TextSymbol("hello 123");
            textSymbol1.textFormat= new TextFormat("Arial",10);

            var myPolygon:Polygon = new Polygon(
                [[
                    evt.mapPoint,
                    new MapPoint(evt.mapPoint.x+5, evt.mapPoint.y+8),
                    new MapPoint(evt.mapPoint.x+10, evt.mapPoint.y+12)
                ]], new SpatialReference(4326));
            cs.symbols= [textSymbol1, textSymbol2, new SimpleFillSymbol(SimpleFillSymbol.STYLE_BACKWARD_DIAGONAL,0xFF00FF,1)];
            var m_graphic:Graphic = new Graphic(myPolygon , cs );
            gLayer.add( m_graphic );
        }

TextSymbols accept point geometry and I am assigning them with multipoint geometry I think that is where problem lies. I need to draw TextSymbols all around the polygon so how can I accomplish it. I am using FlexViewer 2.5.

Any help would be highly appreciated.
Thanks


texts: Summarize values from a raster brick by latitude bands in R

texts:
I would like to sumarize (sum and/or mean) values from layers on a raster brick based on latitude bands (been using raster package). Been looking at the package documentation and didn't find a way to do it. Can anyone give any advice on how to do it?

I guess one option would be to use extract by spatial polygons? Or is there a way to do it with another raster where the values would be classes of the latitudes of cells?


texts: Seeking alternatives to pgRouting for open source routing / network analysis?

texts:
I was wondering if there are any alternatives to pgRouting tools.

Thinking about it for a bit, I think the only alternative is to use R or RPY to do the analysis.

Basically I'm trying to solve routing problems based on nonPostGIS databases. As most of my data is actually stored in a different format; however, it's not to hard to convert. I'm finding that pgRouting to my dumb silly brain is difficult to implement. So perhaps a non-database solution that I can easily tweak would be ideal for me


texts: Weighted heatmap analysis- density of points, vs density of weight variable

texts:
Was just wondering if anyone could help me with a problem. I'm trying to map wellbeing scores (out of 10), using locations from a survey of 15,000 households, using the HEatmap command in QGIS. It's an unclustered sample, but obviously the density of points is dependent on the population density of the area. I'm using wellbeing score as a weight in the heatmap, but unfortunately every map I produce is effectively a map of population density. I have tried multiple ways of weighting the wellbeing score to factor in the differing densities of points (so scores in areas with only a few data points are weighted, up and scores in dense areas are weighted down), but the heatmap looks effectively the same every time.
Is there any sensible way of doing this?
Thanks.


texts: How to open .dbf table exported from ArcMap 10.7 with Access 2013 or later Microsoft Office versions?

texts:
I have exported an attribute table as .dbf file with ArcMap 10.7 but if I try to link it to an Access 2013 database (it works fine with Access 2007) I get a message error  "Could not find installable ISAM".

It happens with any .dbf files and on different computers. Only Access 2007 seems to work.

I am not too happy with the embedded reporting tool on ArcMap and would like to extract the table and use the features on the newer Microsoft Access to then print a report.

Is there any other ways I can extract the attribute table and link it to a newer version (not 2007) of Microsoft Access database?


texts: How to create a (offline) custom map from OSM data?

texts:
I am a beginner with GIS software, but know some basics of the OSM format.

So I want to create a custom map for a game. I already created a heightmap from DEM data, but I need to color the water bodies in black. The way I wanted to do this is to create a custom map from OSM data, which shows only the water bodies and then overlay it with the heightmap using Gimp or something similar. The Region I want to map is approximately 70km * 70km and will probably be scaled down later.

So I thought this would be simple, but after several hours and trying out several tools I am still not anywhere near the desired result. It seems like the available tools are straight out of the 90's and dying when used with the amount of data I want to process or the development has been discontinued and some bugs make it impossible to use.

I tried the following tools:

Merkaartor (0.18.3 on Linux): It takes more than 10 minutes to load the .osm file I have downloaded.
QGIS (2.18.16, Linux): This seemed to be the most promising solution so far. I managed to download the OSM data, create the database and layers. Unfortunately for unknown reasons some water bodies are simply not shown. The current Windows version of QGIS just hangs in the last step of the setup and then is missing some .dll file.
Maperitive: I have read that it should be able to do this, however the rule system is quite complex and it seems a bit of an overkill.
TileMill (v 0.10.1, Windows): This was another promising software, but soon it turned out that there were some bugs and I could not really figure out how to get it to do what I want. Overall the project seems to be discontinued and refers to Mapbox Studio.
Mapbox Studio: At least I could export a png with the data I wanted. However as I cannot export an unlimited amount of data I would rather have an alternative. Plus the export dies when i request the area with 7k x 7k pixels.

Is there an easy way to archive what I want that I have overlooked? As a software developer all of this software seems so odd to me. It seems so hard and time consuming to archive a seemingly simple objective.

EDIT: Added my experience with Maperitive.


texts: Finding format/type of the vector geometry

texts:
Does anybody know what format is the following element?
01000000804AE9B78F8F0D416258AE25F8E45141

It looks like WKB or PostGIS geometry/geography but I didn't manage to convert it.
It has to be point in the centre of the city of Sofia, Bulgaria. The coordinate system should be EPSG:32635. There are some other attributes of the SRS:
FalseX=-5120900; FalseY=-9998100; XYUnits=450445547.391054; FalseZ-1000; ZUnits=2147483.645; FalseM=-100000; MUnits=10000; IsHighPrecision=1; XYTolerance=0.001; ZTolerance=0.001; MTolerance=0.001

But I don't know if I should create user defined CRS with these elements.
The coordinates of this point given in another table are (198; 4733). If I use them with the standard EPSG:32635 the point is located far more South.


texts: Why do I get error "Runtime: b'no arguments in initialization list'" when running to_crs?

texts:
I have a geopandas dataframe that has POINT geometry in NAD83. I want to use to_crs to change them to espg: 4326. When I first bring in the feature class as a gpd dataframe, it shows no crs, so I use a function to find the source crs and then assign it to the dataframe (my understanding is that to use to_crs it has to already be assigned to a crs).

bldg_feat_cent_df = gpd.read_file(bldg_cent_cb_intersects_gdb, layer=state_bldg_centroid)
print(bldg_feat_cent_df.crs)
>>>None
bldg_feat_cent_df.crs = {'init' :'+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs'}
print(bldg_feat_cent_df.crs)
>>>'+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs '

Now that its in a crs, I simply want to apply to_crs. I've tried this a couple of different ways, including inplace=True but I still get the same error:

bldg_reproj = bldg_feat_cent_df.to_crs({'init': 'epsg:4326'})

> RuntimeError: b'no arguments in initialization list'

I had a similar issue trying to use pyproj and it seems like the 2 are related. Conda shows that my pyproj version is 1.96. How can achieve my goal and resolve this issue?

Update: As requested, here's the output when I run the above line:


texts: Cannot get OSM tags when downloading data with osmnx functions which are based on Nominatim requests

texts:
Context
I am trying to download some footway graph from OpenStreetMap (OSM) using osmnx.
Here is a piece of code:
import os, sys
import pandas as pd
import geopandas as gpd
import matplotlib.cm as cm
import matplotlib.colors as colors
import networkx as nx
import osmnx as ox
ox.config(log_console=True, use_cache=True)
ox.__version__

G_test = ox.osm_net_download(polygon=None,
                             north=46.95489, west=7.46006,
                             south=46.95335, east=7.46132,
                             network_type='walk'
                             )

The ox.osm_net_download() method, is based on the Overpass API, so if I check what's in G_test and search for this OSM way 295557051 I can see that the OSM tags are shipped with the data:
{'type': 'way',
    'id': 295557051,
    'nodes': [2992861879,
     2992861880,
     (...) # bunch of vertex coordinates here, skipped for clarity
     2992864399],
    'tags': {'footway': 'sidewalk', 'highway': 'footway', 'surface': 'paved'}
}

The tags are of course the same as the one found by making a direct call to the Overpass API or the OSM API (on which the previous seems to rely) at this address for this way:
https://www.openstreetmap.org/way/295557051
Up to this point, everything is fine.
But the ox.osm_net_download() method doesn't return a networkx graph object, it returns a Python list.
Therefore, to download a graph, I usually query the data through any of the proposed graph_from_*() methods, for example a graph_from_bbox*() as:
G = ox.graph_from_bbox(north=46.95490, west=7.46000,
                       south=46.95330, east=7.46230,
                       network_type='walk',
                       retain_all=True)

Issue
But when downloading the data using the latter method, the OSM tags no longer come with the data!
This is probably due the fact that the graph_from_*() methods are not based on the Overpass API but rather on the Nominatim API.
This may be related to: Cannot get response from nominatim API on some OSM objects
If I explore the result of the G graph, searching for the OSM way 295557051:
G_raw.edges(keys=True, data=True)

One of the may results (as they are many edges composing this OSM feature) look like this:
(2992861879, 2992861880, 0, {'osmid': 295557051, 'highway': 'footway', 'oneway': False, 'length': 4.821}),

For a better comprehension, we can convert it to a GeoDataFrame using the ox.graph_to_gdfs() method:
G_gdf = ox.graph_to_gdfs(G,
                         nodes=True,
                         edges=True,
                         node_geometry=True,
                         fill_edge_geometry=True)

And print its content:
G_gdf[1][G_gdf[1]['osmid']==295557051]

(OSM tags are not there.)
This problem is also propagated when saving the graph to an ESRI Shapefile using ox.save_graph_shapefile():
ox.save_graph_shapefile(G, filename=os.path.join('~/test_graph'))

Indeed, when I open the resulting Shapefile, for example in QGIS, a query on that same way doesn't show the OSM tags in the attribute table.

Question
I wonder where the tags are lost in the process and how to solve for this issue?
Note:
I was first thinking of converting the G graph to a GeoDataFrame, merging it with G_test which I also convert to a GeoDataFrame and then convert back that merged GeoDataFrame to a networkX graph, but I went into a strange error because the retrieved graph doesn't have the right structure:
TypeError: float() argument must be a string or a number, not 'NoneType'

Or to directly write the GeoDataFrame to PostGIS:
AttributeError: 'GeometryDtype' object has no attribute 'base'

But GeoPandas doesn't handle (yet) writing data to PostGIS:
https://github.com/geopandas/geopandas/pull/1248
So I'm stuck.


texts: Why is MapServer adding horizontal artifacts when projecting EPSG:4326 data to EPSG:3857?

texts:
We're running MapServer with a map that has just one layer, political boundaries from a SHP file.  This layer is using a basic EPSG:4326 projection, and the map itself is using EPSG:3857.  The resulting images generated by MapServer CGI queries have horizontal lines at certain latitudes:

Beyond this problem, the map is okay and lines up perfectly with the Open Street Map in OpenLayers.  The MAPFILE we're using:

MAP
  WEB
          METADATA
        "wms_title" "Map Server"
        "wms_enable_request"  "*"
        "wms_srs"             "EPSG:3857"
           END
  END

  PROJECTION
    "init=epsg:3857"
  END

  IMAGETYPE      PNG
  EXTENT -20037508.34 -20037508.34 20037508.34 20037508.34
  SIZE           400 300
  SHAPEPATH      "../data"
  IMAGECOLOR     71 245 242

  LAYER # States polygon layer begins here
    NAME         Borders
    DATA         ne_10m_admin_1_states_provinces_shp
    STATUS       OFF
    TYPE         POLYGON

  PROJECTION
    "init=epsg:4326"
   END

    CLASS
      NAME       "Borders"

      STYLE
        COLOR        43 156 51
        OUTLINECOLOR 32 32 32
      END
    END
  END # States polygon layer ends here

END

I've verified that the proj.4 EPSG file has the proper definition of EPSG:3857, and don't understand what's happening.  Any insight would be greatly appreciated.


texts:
I have a problem with the function buffers in QGIS. I prepared a layer with bus stops and buffers for them (500 m). Then I joined layer with buffers and layer with density of population. I want to know the number of inhabitants living within 500 m from the bus stop. That's why I use join attributes by location. Everything came out well, but the problem is the overlapping buffers. "Connected people" are in two or three buffers at once. And I need to know the number of inhabitants outside of these buffers. Is there any possibility of division between buffers or converted, so as not counted them a few at a time?

Orange and brown spots are inhabitants.


texts:
I have several DEMs as geotiff files which I load into R using raster::raster. Now some of these DEMs contain patches of missing values, represented by NA. I would like to interpolate these missing values from the surrounding cells using Kriging.

One approach I tried was the raster::interpolate function. However, the examples in the documentation use spatial points data frames as predictors (from the "meuse" data set) and I can't figure out how it works for my raster.

I have also seen this, but the approach suffers from the same problem (interpolation of a raster from points). A similar example is this.

Several other questions on https://gis.stackexchange.com address similar issues, however none of them provided a working solution for me.

What I also tried is the R package filling. After converting my raster to a matrix, several functions therein work for me, yet I am not sure if these methods are appropriate for the problem at hand.

I have the impression this should be a pretty straightforward issue, yet I could not solve it so far despite considerable effort. Any help would be appreciated.


texts: ArcObjects: with transformation reprojected tiles are drawn with spaces

texts:
I have used this Code from Kirk to choose a transformation if the spatial reference do not match when adding a layer with tiles. The tiles already got automatically reprojected to the existing spatialreference. I only added the IGeoTransformation and the transformation direction to the raster before being reprojected. As result the tiles get drawn but do not match and leave blank areas like this (it's getting worse with higher zoom levels):

I do unterstand that this has to do with the transformation and is caused by the performed rotation and translation.

How can I fix this?


texts: How to draw circles with a radius based on lat/long data in DB

texts:
I have a . sqlite file with many columns, two of those columns being lat and long.

I've been searching around for a bit now (don't think I want to be using BuildCircleMbr from spatialite?) but I'm not entirely too sure how to go about it.

I'm thinking it has something to do with GeomFromText but can't get the query set up correctly.

Final product would be having these circles and display them on a map using MapServer/OpenLayers.


texts: Converting geoalchemy2.elements.WKBElement to WKT

texts:
I have a geoalchemy2.elements.WKBElement that comes from a query on a Postgres/PostGIS view. I want to convert this to WKT, or otherwise encode it so that it is (Geo)JSON serialisable. How can I accomplish this?

Assuming similarity to shapely, I thought that geoalchemy2.elements.WKBElement would have a wkt property or class method, but it does not. The __str__ method returns the binascii.hexlify representation of the WKB. In the source for the geoalchemy2.elements there appears to be little else to use.

I'm about to build GDAL with python bindings to attempt to read the WKB and then output as WKT, but this feels like more work than I expected.


texts: Geom and ST_AsText(geom) gives different intersection results for same geometries

texts:
I am working to get the count of points intersecting with a convex hull. When I use the query in 1 it gives a count different than count in query 2. Only difference is that in query 1 I am using directly geom whereas in 2, I am converting geom to text and back to geom using ST_GeomFromText()
1:
select count(*) from table
where st_intersects(
        geom,
        (select st_convexhull(st_collect(geom)) a
         from  table)
       );

2:
select count(*)
from table
where st_intersects(
       geom,
       (select st_geomfromtext(
                 st_astext(
                   st_convexhull(st_collect(geom))
                   ),
                4326)
        from table)
      );

Does this conversion using ST_AsText change the geometry somehow? The reason to test this is that I need to send the convex geometry from web frontend as text and convert it to geom in PostgreSQL.


texts: How to extract NMEA outputs from an external bluetooth GPS to an Android device

texts:
I am requested to create an Android app that have to connect to an external Bluetooth GPS and read the NMEA output and visualize the location on ad digital map. I am a new comer in this aspect. Grateful if anyone can suggest how can i connect the GPS and read the NMEA file? Please suggest any good reference materials / source code if available.

I am instructed to simply import the coordinates from the non-professional grade GPS receiver (i.e. GoPass GPT-810 Bluetooth GPS Trip Recorder), and visualize them on top of a geotiff. How can the android "locationlistener" get the coordinates from the GNSS receiver?


texts: What projection is used in Argentina cadastral maps from 1901?

texts:
I am looking to georeference and digitize the cadastral maps for Argentina from 1901 and I need some help in determining the projection which was used in the creation of this maps.
For more detailed maps I can use land marks and roads but in some maps I really need to use the grid from the map I have no idea on the projection and google does not offer any help on this.


texts: Aggregate multiple image class values into a single class (Earth Engine)

texts:
I have created 20 clusters from unsupervised classification of a specific region in Google Earth Engine. Now I want to reclass them. For example, let's say class values 10, 4, and 16 are forest, I want to combine these class values into a class that represents forest. Similarly, I want to do this for water, agricultural land, etc. How can I do that?

link code

var poi=ee.Geometry.Point([78.032188,30.316496]).buffer(35000);

var dataset = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
                  .filterDate('2015-01-01', '2015-12-31');
var trueColor321 = dataset.select(['B3', 'B2', 'B1']);
var tru=trueColor321.reduce(ee.Reducer.mean());
var truclip= tru.clip(poi);
var trueColor321Vis = {
  min: 0.0,
  max: 0.4,
  gamma: 1.2,
};
Map.setCenter(78.032188,30.316496, 9);
Map.addLayer(truclip, trueColor321Vis, 'True Color (321)');

var dataset1=ee.ImageCollection('LANDSAT/LC08/C01/T1_32DAY_NDVI').filterDate('2015-01-01','2015-12-31');
var NDVI= dataset1.select('NDVI');
var ndvimean=NDVI.reduce(ee.Reducer.mean());
var ndviclip=ndvimean.clip(poi);
var visparam={
  min: 0.0,
  max: 1.0,
  palette: [
    'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
    '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
    '012E01', '011D01', '011301'
  ]
};
Map.addLayer(ndviclip,visparam,'NDVI',0);

//clustering
var training=ndviclip.sample({
  region : poi,
  scale : 30,
  numPixels : 5000
});
var clusterer=ee.Clusterer.wekaKMeans(20).train(training);

var result = ndviclip.cluster(clusterer);

Map.addLayer(result.randomVisualizer(),{},"clusters");

var subset= result.select("cluster").eq(4).selfMask();
Map.addLayer(subset,{},'only 19');


texts: How to make a private OSM (Save personal information to) a map using OpenLayers or TileMill or mapbox

texts:
OSM (OpenStreetMap) is great for making a real time public map for general purposes but what if I want to add private info to a map such as the location and date of my first kiss or a note about every location I posted fliers etc.

I'd like to know how to add private information to OSM.

I have heard that TileMill or OpenLayers are best for this but cannot find a video tutorial on adding private map data to a private layer and being able to back that info up and update it or transfer it to others to edit or update etc...

As I research more I'm seeing that TileMill/Mapbox are what I am looking for.

Which makes no sense. It is expected that OSM has this functionality and what does the name TileMill have to do with OSM or maps at all? I think a lot is lost in new possible community members coming to OSM for this functionality and being totally confused and sent to something called TileMill.

Still haven't found any video tutorials for TileMill/Mapbox or any clear info that these are in fact what I am searching for.


texts:
Goal is to group polygon (1,2,4) and polygon (3) based on overlap. All polygons are part of the same shapefile layer.

In ArcMap I can simply do dissolve and uncheck Create multipart features.

However dissolve in geopandas requires you to set an attribute to dissolve on. What would be the easiest alternative in geopandas?


texts: Raster Function for NDVI Band Calculation - Divide by Zero

texts:
I am using rasterio to open image files related to 2 specific bands:

b4 = rio.open(r)
b8 = rio.open(poly)

I then read the bands into a python array, I think it's a numpy array but not sure?:

#Read Red(b4) and NIR(b8) as arrays
red = b4.read()
nir = b8.read()

Then I perform band based calculation against them via the following:

#Calculate NDVI
ndvi = (nir.astype(float) - red.astype(float)) / (nir + red)

Now, this of course can have the denominator (nir + red) be equal to zero, which will cause a divide by zero issue that I believe is messing with the results.

Is there a way to set the instances where the denominator is equal to zero to have a value of -1 or something?  How are you handling those in your own raster functions?


texts: How to paste a feature multiple times at once?

texts: I have a feature class and I want to paste one of the records multiple times. I know this will duplicate the record but that is my intention. I also know that I can simply copy and paste the record and it will duplicate it. I was wondering if there was a way to duplicate the record 10 times at the same time instead of copying and pasting over and over.


texts: How is dateline a problem when calculating distances?

texts:
From http://postgis.net/workshops/postgis-intro/geography.html, I can understand how poles might be an issue but how do datelines create problems?

  The distance between points get larger as problem areas like the poles or the international dateline are approached.


texts: QGIS - Map won't allow me to move around, only shows bottom half of Africa

texts:
I just installed QGIS today, so I'm very new to it. I've installed the Quick Map Services plugin and tried adding Google Maps, Google Satellite, Bing Maps, etc., and none of them seem to let me pan beyond this section shown here.

It's the bottom half of Africa. Whenever I try to pan somewhere else, it shows blank map around this section and snaps back to here.

Any idea why this is happening?


texts:
I've been able to export features as KML using GeoTools. However, I would like to insert a name attribute at the top level (e.g. Document, though a Folder would be fine too). (Example of what I'd like at the end.)

This is what I am getting so far:

<?xml version="1.0" encoding="UTF-8"?>
<kml:kml xmlns:xs="http://www.w3.org/2001/XMLSchema"
         xmlns:xal="urn:oasis:names:tc:ciq:xsdschema:xAL:2.0"
         xmlns:kml="http://www.opengis.net/kml/2.2"
         xmlns:atom="http://www.w3.org/2005/Atom">
    <kml:Document id="area-418e1090-6ca8-4e34-a563-9c36db53cb35">
        <!-- would like "name" here -->
        <kml:Placemark id="shape-1">
            <kml:ExtendedData/>
            <kml:Polygon>
                <kml:outerBoundaryIs>
                    <kml:LinearRing>
                        <kml:coordinates>0.0,0.0 0.0,20.0 20.0,20.0 20.0,0.0 0.0,0.0</kml:coordinates>
                    </kml:LinearRing>
                </kml:outerBoundaryIs>
            </kml:Polygon>
        </kml:Placemark>
        <kml:Placemark id="shape-2">
            ... et cetera

        </kml:Placemark>
    </kml:Document>
</kml:kml>

The code for that I've scraped together from a bunch of examples, it boils down to this:

ArrayList<SimpleFeature> list = new ArrayList<>();
addShapes(entity, list);

DefaultFeatureCollection featureCollection = new DefaultFeatureCollection(id("area", entity.getUuid()), null);
featureCollection.addAll(list);

OutputStream outputStream = httpServletResponse.getOutputStream();
encoder.encode(featureCollection, KML.kml, outputStream);

And:

private void addShapes(AbstractAreaEntity entity, ArrayList<SimpleFeature> list) {
    List<Geometry> shapes = entity.getShapes();
    for (int i = 0; i < shapes.size(); i++) {
        Geometry geometry = shapes.get(i);
        list.add(new SimpleFeatureImpl(singletonList(geometry), shapeType, new FeatureIdImpl(id("shape", i + 1))));
    }
}

This is what I would like:

<?xml version="1.0" encoding="UTF-8"?>
<kml:kml xmlns:xs=...>
    <kml:Document id="area-418e1090-6ca8-4e34-a563-9c36db53cb35">
        <kml:name>My Document With Placemarks</kml:name> <!-- this! -->
        <kml:Placemark id="shape-1">
            ... et cetera

        </kml:Placemark>
    </kml:Document>
</kml:kml>

It appears that Encoder.encode only works with a FeatureCollection implementation, and it does not seem I can set properties or attributes on it (like name).


texts:
I have whole bunch of GPS coordinates collected from a device and sometimes they shift around 10m. Is there anyway to filter/modify points and snap them to sidewalks? So far I've written some javascripts to map those points to the center line of the nearest driving path but I still want them to be on sidewalks.

I noticed ArcGIS can display street curbs. Is there any way to extract the curbs? I cannot finde APIs for curbs at all.

There are some questions on this site regarding making sidewalk maps but no real solutions are available yet. I would prefer utilizing existing data sets like shape file and curbs to map points instead of adding sidewalks on google maps manually.


texts:
I need a few thousands of street distances between cities in germany.
I know that there are a few web APIs where I get the distance thru a REST request, but I wonder if I can do that locally on my Linux PC. I saw the OSRM project and libosrm seems to do what I want.
Are there any other tools doing that?
What do I have to download and what do I have to install?
(I prefer to program in python)


texts: "Polygon With Hole" gives a 'NoneType' error when assigning it to a shapely shape in Python

texts:
I am using Python and shapely to do certain spatial calculations (difference, intersections etc).

Bur before I can do the calculations I first need to get the object into a shapely object.
so I iterate through the shapes in a feature class using arcpy.da.searchcursor:

with arcpy.da.SearchCursor(OurCADLAYER, ['SHAPE@'], whereclause) as onerows:
        for onerow in onerows:

and then assign the shapely object here:

            Shape2 = shape(onerow[0])

and this is where I get the error message:
TypeError: 'NoneType' object has no attribute getitem

PS. THIS ONLY HAPPENS when the shape is a Polygon that has holes (interior rings). no issues with normal polygons..


texts: How to find the height information from LiDAR point cloud data

texts:
I need to find the biophysical parameters(biomass) from the LiDAR point cloud data. For that I need the height of the trees from the data. Currently I am having the general height with respect to sea level(400 range).

I tried with lasheight, and also tried to create dem and dsm and from them chm but nothing is working with the data set. I have attached my las data too below.

Data

I tried with lasheight, and took the output, its supposed to give the height information as 0 for ground and for the trees for example 10m or something. but I'm getting the same 400 range values for ground and trees.

I tried with las2dem and tried to save its output as tif image.(I am programing using R), while reading that tif image to R. Its giving error. Its not taking that tif image


texts: Adding text and legend to image collection animation in Earth Engine

texts:
I have a GIF image generated using the ui.Thumbnail() function

I wanted to embed the label below, inside the GIF image.

The GIF is generated using ui.Thumbnail and the label is generated using the ui.Panel. I Want to combine them together into a single GIF image, so that I could share the GIF link directly.
This is the very first time I am using Google Earth Engine. If this is not the right way to do this, what other way do I have to achieve the same result?


texts: How to rename duplicated layers in QGIS3?

texts:
I duplicated a layer ("Buildings") in QGIS 3.4.4 and with the Query Builder under Source I selected to just display schools. Subsequently I renamed this layer "Schools" by rightclicking on the layer, choosing Rename Layer.

The new name of the layer appears but once I tick the box to toggle visibility it jumps back to its previous name ("Buildings"). I'm also unable to rename the layer under Source (the 'displayed as' field is greyed out), see screenshot below.

How can I rename this layer permanently?

Edit 21-02-2019: I found that in other projects the layer could be simply renamed or saved with a different name just as was suggested in the answers below. I have not been able to do it with the project I had worked on and suspect a bug of some sort with the imported layers.


texts:
I have a folder of GML files that I want to be able to merge to create one GML file that I can then view.
I have tried this using the code described but searching for GML in this website without result.

The code I used was:

mkdir merged
for %f in (*.gml) do (
if not exist merged\merged.gml (
ogr2ogr -f "GML" merged\merged.gml %f)
else (
ogr2ogr -f "GML" -update -append merged\merged.gml %f -nln Merged )
)

For clarity the files are for the OS mastermap topography data.

Forgot to say I get the following error message several times when I run the code.

  C:\OSMM_GMLS>(if not exist merged\merged.gml (ogr2ogr -f "GML"
  merged\merged.gml  M00001177-SX6590-5c296.gml )  else (ogr2ogr -f
  "GML" -update -append merged\mer ged.gml M00001177-SX6590-5c296.gml
  -nln Merged  ) ) FAILURE: Unable to open existing output datasource `merged\merged.gml'.


texts: Where can i find free mapnik map styles with focus on waterways?

texts:
I am using Google maps for my map tiles. These tiles are very road oriented. As i am making an application for waterways i would like to change my map so that its main focus is on waterways.

I have a little bit experience with making custom map tiles with Mapnik, but it still requires me to make a lot of styling decisions and styling is not my strongest point. So it would be great if there are some free stylings available that let me render map tiles with a focus on waterways, like the ANWB (Dutch) water map.

I have found a site where i could change and test styles on the fly, but they had to serve the tiles and i think it also cost money to use their service. I am able to host such a tile service myself once i have the correct tiles.


texts: After creating a Visualization I don't see the Publish button on the top right of the page, instead I see a Share button. Where is the Publish button?


texts: Check if shapefile Polygon has a valid number of points - Python

texts:
I'm trying to create a polygon from ShapefileReader.shapes() object in python using the shapely.wkt library. But I'm having the followeing error:

IllegalArgumentException: Invalid number of points in LinearRing found 3 - must be 0 or >= 4

Is there anyway to create a valid Polygon before parsing it? I know that ShapefileReader.shapes() returns a String object, so how can I check this String to see if it will create a valid polygon, or even better, if there's a way to force it to create a valid polygon?


texts:
I am building a map in leaflet, rendering a variety of layers on top of the map.
One layer contains geotagged images fetched from a folder (with Leaflet.Photo).
Some of the images GPS locations are not accurate enough and I would want to re-position them.
My question is:

Is there a open source tool that allows repositioning of images on a
map and edits the exif data accordingly?
Would this be possible to integrate this into the Leaflet map itself (e.g have the popup contain an 'edit' button that allows write access
to an images metadata)?


texts: What is the second argument in QgsProcessingUtils.mapLayerFromString()?

texts:
I am new to Qgis3. I am trying to get a layer in the Python console. And I have to put to arguments in the function QgsProcessingUtils.mapLayerFromString(). The first argument is self-explanatory, the string. But I don't understand what the second argument is (the QgsProject).
Could someone clarify me what that second argument is?


texts: What has prevented a "better" format than shapefiles from being widely adopted?

texts:
Reading through this answer got me thinking: Why hasn't a better format than the shapefile been widely adopted?

Shapefiles have (as I'm sure most of us have noticed) many annoying quirks, shortcomings and restrictions. I won't go through the trouble of finding every complaint that has ever been made about the format, but common complaints I see are:

Multiple files are needed for each layer
File name limitations
Field-length limitations
Tables are no longer editable in MS Excel

With these (and other) complaints in mind, what has prevented a "better" format from being widely adopted? Sure there are alternatives such as personal and file geodatabases, but these aren't as portable/shareable as shapefiles.


texts: QGIS Joining Parcel Polygons where overlay of Building Footprints crosses Multiple Parcels

texts:
QGIS 3.4.4

Polygon Layer 1 = Land Parcels

Polygon Layer 2 = Building Footprints

Some buildings exist across several parcels. How can I use the information in Layer 2 (red in the image below) to identify and join all parcel polygons in Layer 1 (black in the image below) that are beneath the same building?

And, in the process, assign the data from Layer 2 (area, height, elevation) to the newly joined parcel polygon?

I realize that I could probably assign the Object ID from Layer 2 to Layer 1 and sort -- but I feel like that might not work where multiple buildings are involved.


texts: Garbage Truck Fleet Route Optimization Using PGRouting?

texts:
I am trying to create an optimized route plan for a townships garbage truck fleet for garbage pickup. I have tried using ESRI's Network Analyst within ArcMap but have found the analysis to be to complex and the results returned are not as desired. There are approximately 2000 stops within each zone within the township with 3 trucks per zone.
Recently I have come across the PGRouting SQL extension and would like to know if this type of project would be to complex for PGRouting to run? Can PGRouting handle the large volume of pickups?


texts: Finding path between points on limited road network

texts:
I am struggling to find a path between health facilities in Madagascar.

I think that because the road network (from OpenStreetMap) is limited, the functions in QGIS, such as v.net.distance and shortest path (point to point) can't find a path.

The OSM tools plugin can't find a path either.

When I use the online routing mapper plugin, however, a path is always found.

Unfortunately, you can't import a shapefile, but you have to click on the coordinates to find a route. Since I have multiple health sites between which I would like to get a path, this is barely doable. My programming skills are very limited.

Is there a way to do it more efficiently or another plugin I should try or am I may be wrong to think the problem is due to a limited network?


texts: Why does OpenLayers load OSM-like tiles, but does not show them?

texts:
I created a tile server with OSM-data in the German style. The data is quite similar to the official tile server for openstreetmap.de, provided by Fraunhofer IOSB.

http://tile.iosb.fraunhofer.de/tiles/osmde/${z}/${x}/${y}.png

If I include the standard OSM-server

http://a.tile.openstreetmap.org/${z}/${x}/${y}.png

with the following constructor

new OpenLayers.Layer.OSM("OSM",myOSMURL,{...});

I get a valid map. The OpenCycleMap tile server works as well.

If I want to use my own tile server or the one from Fraunhofer, the browser successfully gets an respond to the GET-request from the tile server and loads the tiles needed, but doesn't show them in the map.

Following the steps here, I generated a map with local data successfully.  Nevertheless, if I replace the link of the local data with a link to Fraunhofer or my tile server, the browser still loads the tiles, but doesn't show them.

// This is the layer that uses the locally stored tiles
//var newLayer = new OpenLayers.Layer.OSM("Local Tiles", "tiles/${z}/${x}/${y}.png", {numZoomLevels: 19, alpha: true, isBaseLayer: false});
var newLayer = new OpenLayers.Layer.OSM("Local Tiles", "http://tile.iosb.fraunhofer.de/tiles/osmde/${z}/${x}/${y}.png", {numZoomLevels: 19, alpha: true, isBaseLayer: false});
map.addLayer(newLayer);
// This is the end of the layer

Image in full resolution: http://i.stack.imgur.com/73jOX.jpg

How do I get the tiles from my tile server or Fraunhofer to my map?


texts: Label or font colors on an imagery background examples

texts:
I have been in searching of a website that was done by somebody from ESRI that wrote about...if I remember right...it was in either on the ArcGIS Mapping Center or in the ESRI Technical support article....

There was an article that an author from ESRI wrote about an examples of colors of Labels or font on over a satellite imagery in the background...

Can't remember what it was about the article but I believe it was meant to give the users an examples of colors to show or visualize you the colors with the image in the background.

I had that site saved it to my bookmark back then before ESRI decided to redesign their website and it may have move or remove from their server.

I would like to get that information restore again but I can't find it on the internet...

Anyone remember it?

P.S. I wanted to modify this and add :  This is mean to help my end-user to distinct the colors of the labels on the image backgound and I know ESRI wrote that article..

And one more thing from the Mapbox : https://blog.mapbox.com/designing-labels-for-satellite-and-aerial-imagery-7036163aa21

But it lacks all the different colors...


texts: How to check rasters with at least on cell falling within a polygon boundaries

texts:
I have some independent rasters (named r1, r2, r3, r4, and so on), congruent each other concerning reference system, extent and resolution. Each raster has cells that score the following values: 1, 2, 3, or NA.

Moreover, I have a layer represented by a polygon (named "boundary").

Is there a method to count and list (i.e. identify) all rasters which have at least one cell with no NA value falling within the polygon boundaries?

I'm working with R: could you provide me the appropriate code, please?


texts: How to identify interior rings and set them to empty in a shapefile to avoid warning

texts:
import geopandas as gpd
from pyproj import CRS

I'm working with a lot of shapefiles.

I've checked their geometry to be valid .is_valid.
I have updated to geopandas 0.81 shapely 1.7.0 python 3.6.10
I have reprojected my coordinated systems to projected one geo.to_crs(CRS.from_epsg(31983))

Still, everytime I check the file for a simple geo.head() query, I get the warning from geopandas which is great, but has very scarce documentation...
UserWarning: Only Polygon objects have interior rings. For other geometry types, None is returned.
"Only Polygon objects have interior rings. For other "

I have searched around quite a lot. Found this #issue here https://github.com/geopandas/geopandas/issues/1331 which is closed and did not help.

In the end, what I think is wrong is that there is at least one interior.ring inside my shape which is not an empty list.

How can I set a specific interiors line to []
And in order to set it automatically, I would also need a test that returns which observations have interior rings.
I'm using PyCharm 2020.2
EDIT You can find the data here: https://github.com/BAFurtado/firms_growth_by_APs/ Download all 5 example.* files (shapefile format)


texts: Improving Trimble GPS Pathfinder Office and TerraSync coordinate systems workflow?

texts:
We have GPS survey projects all over the country (CONUS and OCONUS) and I am wondering what the best practice is for having to deal with multiple coordinate systems. I can go from WGS84 UTM 18N meters on one project to NAD83 TX SP E feet to multiple others in one day.

We're using GeoXH receivers and here's my current flow:

On one site, I start with GIS source data in say NAD83 NC State Plane feet.
In PFO: Options -> Coord Sys -> set to NAD83 NC SP feet. Import -> select the shapefile and set the coord sys to NAD83 NC SP feet. Open the imp file. Data Transfer -> Send -> Select the file.
Then in TerraSync: Setup -> Coord Sys -> NAD83 NC SP feet. Set the imp as a background file, they collect the data, email it back (or toss the GeoXH on my desk).
Back in PFO, I set the coord sys to whatever their data is in (sometimes it's the same - yay!), transfer, open, diff correct (*question on that later) Export -> set coord system, and back into ArcGIS it goes.

How do I simplify or streamline this workflow?

Would it be better just to reproject the source shapefiles to lat/lon, import and leave PFO and TerraSync in lat/lon, and then reproject in ArcGIS (back to SP or UTM, whatever) to avoid changing the coord sys in PFO and TerraSync back and forth? And is the projection info stored inside the SSF? I know the data only come in as lat/lon from a sat.

It it baffling that what are considered "projects" in PFO do not remember something as essential as the set coordinate system or am I missing a setting?


texts:
I merged four rasters using GDAL > Miscellaneous Raster > Merge
The result is a black image (only 0 value everywhere on it) of four images extents.
I've executed this geo-algorithm before without any problems... until today.
What is going on?

Aditional data: OS Linux Ubuntu 20.04
QGIS version 3.16.7-Hannover
QGIS code revision 04717438


texts: When using a Wacom Cintiq tablet and pen, can you pan while creating a feature using the freehand polygon plug-in in QGIS?

texts: Is there a way to pause the free-hand tool in QGIS so that i can pan and move over the image and then continue my polygon?


texts: How can I upload more than 10 images in a category map?

texts:
I did a category map with 30 different categories. For each, I would like to upload an image but the "wizards" permits me to upload only ten images and then one image for "others".

Is it possible to upload more?


texts: GRASS GIS producing different slope than GDAL?

texts:
When producing a slope in GRASS GIS 7.4 the results are different than with QGIS (GDAL) and ArcMap for the same input data and same parameter (Z ratio is 1.0 and degrees). QGIS and ArcMap generate almost the same output. In GRASS, I tried both for DCELL and FCELL and the results are same. The original layer has a resolution on 1 m and the GRASS region resolution is set to 0.5 m.

Code GRASS: r.slope.aspect elevation=DMR@PERMANENT slope=slope

Code GDAL: gdaldem slope .../DMR.tif .../slope.tif -of GTiff -b 1 -s 1.0

Do I do something wrong?
I am attaching the screenshots of the situation.

Left - GRASS calculated slope and right - GDAL calculated one


texts: From Raster to Masked and Filtered Polygon Vector

texts:
I have a set of rasters (geologic risk) from which I want to extract zones with a specific color, clipping them with a polygon (administrative borders), outputting a shapefile.
I'm using QGis (and GDAL from command line).
I know I can use clipper, polygonize and show only features with the value I need, but:

it takes ages to polygonize the raster, and I don't want much of the information in it;
clipper doesn't work because the polygon border is larger than each raster, but fits inside the sum of all rasters
to make things complicated, the zone I would like to extract don't have the same color value in all rasters (but do have the same rgb color)

How can I solve it in a fast way?


texts: How can I change the class of one pixel in a classified image?

texts:
I have a classified image (I did it in ERDAS) with six categories. There are some misclassified pixels (especially high-density urban category). Is there any way that I can choose one pixel and then change its class (for instance change bareland to high-density urban) in ERDAS or Arcmap?
I do appreciate your help. Thanks.


texts: Loading WFS GeoJSON from Geoserver in Openlayers

texts:
I am trying to load a WFS JSON layer into my web map. I am running the web page on a MAMP server. The layers are loaded from Geoserver and from a POSTGIS database.

I have 2 WMS layers being loaded from Geoserver and they are displayed on the map without a problem. However, when I try to load the WFS layer, it is not displayed on the map. I have read several similar issues here on StackExchange, but the solutions offered do not really work.

I have tried the WFS Example from Openlayers and it works fine and is displayed on my map. The key lines are:

var vectorSource = new ol.source.Vector({
    format: new ol.format.GeoJSON(),
    url: function(extent) {
      return 'https://ahocevar.com/geoserver/wfs?service=WFS&' +
          'version=1.1.0&request=GetFeature&typename=osm:water_areas&' +
          'outputFormat=application/json&srsname=EPSG:3857&' +
          'bbox=' + extent.join(',') + ',EPSG:3857';
    },
    strategy: ol.loadingstrategy.bbox
});
var vector = new VectorLayer({
    source: vectorSource
});

When I tried to adjust the example to load my own layer, it was not displayed. Here is my code:

var vectorSource = new ol.source.Vector({
format: new ol.format.GeoJSON(),
url: function(extent) {
  return 'localhost:8080/geoserver/postgis_test/ows?service=WFS&' +
         'version=1.0.0&request=GetFeature&' +
         'typeName=postgis_test:clusterresult&maxFeatures=500&' +
         'outputFormat=application%2Fjson' +
         'bbox=' + extent.join(',') + ',EPSG:31370';
},
strategy: ol.loadingstrategy.bbox
});

var vector = new ol.layer.Vector({
   source: vectorSource
});
map.addLayer(vector);

It is also not:
- A CORS problem (this has been solved)
- a problem of MAMP (I think). I tried loading the webpage in GeoServer itself as suggested on StackExchange, but this produced the same result.

The WFS layer was only loaded on one occasion: when loading through an AJAX request. But in this case the features were loaded onto 'Null island' although specifying the correct CRS. I will give you the code for this as well:

var firstWFSsource = new ol.source.Vector({
format: new ol.format.GeoJSON()
})
$.ajax({
type: "GET",
url: 'http://localhost:8080/geoserver/postgis_test/ows?service=WFS&' +
'    version=1.1.0&request=GetFeature&typeName=postgis_test:clusterresult&' +
     'maxFeatures=500&SRS=EPSG:31370&outputFormat=text%2Fjavascript&'+
      'format_options=callback:getJson',
dataType : 'jsonp',
jsonpCallback : 'getJson',
success:function(data){
    // If response is valid
    var geojsonFormat = new ol.format.GeoJSON();
    console.log(data);
    // reads and converts GeoJSon to Feature Object
    var features = geojsonFormat.readFeatures(data);
    firstWFSsource.addFeatures(features);
}
});

var firstWFS = new ol.layer.Vector({
    source: firstWFSsource
});


texts: Rotating a vector layer in QGIS with qgsAffine (or other method)?

texts:
I would like to rotate a set of vector points in QGIS an arbitrary number of degrees around a central point (or arbitrary point).

This is similar to a recent question about creating a regular grid; it was suggested there to use the "Affine Transformation" tool (which I assume meant the plugin) to rotate or shift a grid of points an arbitrary angle or distance. I suspect I'm not understanding how it works, and haven't been able to get it to work.

I create a regular grid of points in QGIS, and ensure that the UTM zone is set correctly for both layer and project, enable editing for the layer, then open the plugin dialog (qgsAffine):

I select 'whole layer' and then, wanting to rotate the entire field of points by 15°, put 15 in both 'rotation' boxes (which may be where things are going wrong). The operation results in rotating the points somewhere off-planet!

Is this the right tool for the job? I would like to rotate a set of points about their common centre, ideally.

Update: qgsAffine is just a thought; if we can do this in any QGIS tool I'll be happy!

Update 2: qgsAffine is usable IF you know the right numbers to plug in (see answer below, thanks Mike!). Spreadsheet / calculator works fine, or here's the R function to get the numbers directly:

## Compute correct affine numbers for qgsAffine plugin
affine <- function(originX, originY, rotAngle) {
  A <- rotAngle * pi / 180
  scaleX <- scaleY <- cos(A)
  rotX <- sin(A)
  rotY <- -sin(A)
  transX <- originX - cos(A) * originX + sin(A) * originY
  transY <- originY - sin(A) * originX - cos(A) * originY
  aff <- data.frame(scaleX, scaleY, rotX, rotY, transX, transY)
  return(aff)
}

So, to rotate a grid of points in northern Uganda (UTM 36N), affine(578988, 419210, 30) gives:

     scaleX    scaleY rotX rotY   transX    transY
1 0.8660254 0.8660254  0.5 -0.5 287174.7 -233330.5

... which, entered into the qgsAffine dialog, correctly rotates the points.


texts: Is snapping working when changing from NAD83(decimal degree) to NAD83 Zone UTM 11N(meters)?

texts:
I've used the following ESRI snippet  to Convert the display extents in Pixels (at the current map scale) and then return out the map units.

public System.Double ConvertPixelsToMapUnits(ESRI.ArcGIS.Carto.IActiveView activeView, System.Int32 pixelUnits)
{

  if(activeView == null)
  {
    return -1;
  }
  //Get the ScreenDisplay
  ESRI.ArcGIS.Display.IScreenDisplay screenDisplay = activeView.ScreenDisplay;

  //Get the DisplayTransformation
  ESRI.ArcGIS.Display.IDisplayTransformation displayTransformation = screenDisplay.DisplayTransformation;

  //Get the device frame which will give us the number of pixels in the X direction
  ESRI.ArcGIS.Display.tagRECT deviceRECT = displayTransformation.get_DeviceFrame();
  System.Int32 pixelExtent = (deviceRECT.right - deviceRECT.left);

  //Get the map extent of the currently visible area
  ESRI.ArcGIS.Geometry.IEnvelope envelope = displayTransformation.VisibleBounds;
  System.Double realWorldDisplayExtent = envelope.Width;

  //Calculate the size of one pixel
  if(pixelExtent == 0)
  {
    return -1;
  }
  System.Double sizeOfOnePixel = (realWorldDisplayExtent / pixelExtent);

  //Multiply this by the input argument to get the result
  return (pixelUnits * sizeOfOnePixel);
}

and then this to check if the snapping is within tolerance ITolerance which is in Pixels(Example 10 pixels) :

Double dblToleranceX;
dblToleranceX = snaptolerance(ITolerance);// call ESRI Snippet
dblProximity = dblToleranceX;
ouble dblToleranceY = dblToleranceX;
pInputpoint = gpMxDocument.ActiveView.ScreenDisplay.DisplayTransformation.ToMapPoint(intX, intY);
pRectangle.PutCoords(pInputpoint.X - dblToleranceX, pInputpoint.Y - dblToleranceY, pInputpoint.X + dblToleranceX, pInputpoint.Y + dblToleranceY);
pSpatialfilter.SpatialRel = esriSpatialRelEnum.esriSpatialRelIntersects;
IFeatureClass pFeatureClass = gpFeatureWorkSpace.OpenFeatureClass(strFeaturecls);
pSpatialfilter.Geometry = pRectangle;
pSpatialfilter.GeometryField = pFeatureClass.ShapeFieldName;
pFeatureCursor = pFeatureClass.Search(pSpatialfilter, false);
pFeature = pFeatureCursor.NextFeature();

if (pFeatureClass.FeatureCount(pSpatialfilter) < 1)
 {
    MessageBox.Show("Snapping is not within tolerance");
                    return;
                }

The code work perfectly on NAD83(decimal degree)  but does not work  when Changing the data frame's coordinate system  to NAD 83 UTM  Zone 11N (Meters).

Appreciate your help if you can figure out a way around this issue in order to make the code work for these kind of situation .


texts: LeafletJS Measurement of distance in deepzoom tilelayer

texts:
I have been trying to understand how LeafletJS determines distances between two points in a layer. I am trying to implement the custom scale plugin, (and have done so successfully) on a static tiled image using the leaflet deepzoom tilelayer plugin.

Added Information

Using CRS.Simple, the image zoom is adjusted by powers of 2 by the looks of it:

CRS.Simple seems to operate in powers of 2:

o.CRS.Simple=o.extend({},o.CRS,{projection:o.Projection.LonLat,transformation:new o.Transformation(1,0,-1,0),scale:function(t){return Math.pow(2,t)}})

Constants that I know:

DPI
tile size (256 x 256 pixels)
map size in pixels (bounds of the image)
leaflet zoom level
size of image in the real world

What do I change in Leaflet to allow me to make accurate scale bars?


texts: How to share business logic among multiple applications using different apis with arcgis?

texts:
I've been working with ArcGIS and Esri WPF API for a year, but now me and my job colleagues felt the need to build a new application using the Esri Javascript API. Our application using WPF API has a lot of business logic that the new one certainly will have too. Is there a way to share this logic among the applications, considering that they use different APIs and must coexist?

Is building GP Services a good option for sharing this kind of logic (like validating feature insertions)?


texts: Specify an alternative URL for a "tileLayer" to use in Leaflet?

texts:
Let's say I have a:

var myTileLayer = L.tileLayer('https://{s}.tile.SomeMapServer.com/{z}/{x}/{y}.png');

and if some tiles fail to load, I want them to be loaded from

  'https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png'

I have tried to achieve this using 'errorTileUrl', but that probably doesn't take "{}" placeholders.

var myTileLayer = L.tileLayer('https://{s}.tile.SomeMapServer.com/{z}/{x}/{y}.png',
{errorTileUrl: 'https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png'});