texts:
I have a large world map in a shapefile from which I need to crop a smaller region. I am planning to use in in java based simulation. How to do that?

I am completely new to GIS world so can anyone recommend a decent tutorial on basic manipulations to help me understand what is actually going on?


texts: Calculating polygon length and width using open source GIS?

texts:
How can I calculate the length and with of polygons using open source tools?
By "polygon length" I mean the length of the longest line within the polygon (update: to be correct, I just need the longest line within the polygon's convex hull) and the width as the longest measure perpendicular to the length measurement.

I'm trying to create a non-axis aligned minimum bounding rectangle for a polygon. The length of this minimum bounding rectangle would be the polygon's diameter.


texts: Invalid expression error 000358 when using in_memory workspace?

texts:
Make Feature Layer with an SQL Selection Expression against a geometry property like length fails when the source feature class is in an in_memory workspace. An otherwise identical model works just fine when the feature class source is in a regular workspace (file-gdb).

ERROR 000358: Invalid expression Shape_Length = 0
Failed to execute (Layer of Zero Length Lines).

I'm using ArcGIS Desktop 10.3.1 but don't believe the issue to be version specific.

Fail:

Success:


texts: Substracting rasters using ArcGIS Raster Calculator?

texts:
I need to adjust a DEM (raster) and I was wondering if the Raster Calculator is the right tool to use, and how?

I’m modelling floods in Xochimilco (Mexico) and via OpenStreetMap I have acquired the shapefile of water bodies. This shapefile I have converted to a raster with the same dimensions as the Digital Elevation Model: both rasters have the same projection and the same cell size (10x10).

My goal is to ‘burn in’ all bodies of water  1 meter into the DEM, so that the cells of the DEM with water become 1 meter lower than their initial value.

I figure that the I should assign the water raster a ‘value’  of 1, and then the formula in Raster Calculator would be [Original DEM – Water Raster = New DEM].

What steps should I take?


texts: Converting Northings/Eastings to Longitude/Latitude using R issue with Northern Ireland?

texts:
I have found a lot of questions about conversions of Northings/Eastings to Longitude/Latitude in R. I have been working with this tutorial from Alex Singleton which is really useful in explaining the process and giving sample code. I have adapted the code for my own CSV file, and worked through the steps.
My issue is, when i plot the data onto the map, it plots accurately in relation to each other, but does not plot over Northern Ireland, instead plots relational to eachother, but floating over wales.

Is there a simple fix for this to shift the plot?

setwd("~/R/maps")
ukgrid = "+init=epsg:27700"
latlong = "+init=epsg:4326"
holes <- read.csv("Pothole_Enquiries_2016.csv", header = TRUE)
map <- readRDS("GBR_adm2.rds")
plot(map, col = 'grey', border = 'darkgrey')
holes <- subset(holes, select = c("EASTING","NORTHING"))
holes$holes_ID <- 1:nrow(holes)
coords <- cbind(EASTING = as.numeric(as.character(holes$EASTING)),NORTHING = as.numeric(as.character(holes$NORTHING)))
holes_SP <- SpatialPointsDataFrame(coords, data = data.frame(holes$holes_ID), proj4string = CRS("+init=epsg:27700")
holes_SP_LL <- spTransform(holes_SP, CRS(latlong))
points(holes_SP_LL,col=2,pch=18)

I am presuming that this is something to do with the UKGrid codes, but i have only started to use spatial plots so this is a little new to me.


texts:
I'm using Arcpad 10.2 R5 installed on my new Trimble Juno. In my apm project I'm using axf files instead of shp. During editing axf I can not select any point/vertex, although I can add new point, fill forms, but the point is unselectable. I click on a point, but it doesn't mark. While of this I can edit and select shp points.
I've tried reinstalling Arcpad, but it doesn't help. The same axf file copying to my second Trimble Juno works perfect and it is fully editable/selectable.
Any idea for this problem?


texts: i calculated some Pointlayers with Qgis and i want to find them in the field. Can i load these points including the GPS-data on mysSmartphone and use it as a GPS-device? Are there any Apps?


texts: Missing styles when adding data from WFD layer from the EA in QGIS 3.16.4

texts:
So I have been trying to open data I have previously saved to our PostGIS DB about a couple of years ago. Specificaly these layers relate to the data from the EA WFD River Basin Districts Cycle 1 and 2.
Now when I add these to QGIS the style disappears but the data table is there and visible just not the polygons.
These will however display via GeoServer as a WMS but not WFS which has the same issue.
See image below;


texts: Is Postgres and Disk IO limiting pgRoutings performance?

texts:
I am doing this Q&A style because I have done my google research and I could not find any thrustworthy and 'benchmarked' data. So I did it myself.

CPU vs DISK as 'the limitting' factor on performance?
How to improve performance through config tuning?
How to improve performance through query tuning?


texts: Import in QGIS a XML with metadata and a reference to the attribute information ISO 19115

texts:
I have created an XML file containing metadata for a specific layer. The XML is
in ISO 19115/ISO 19139. I have also added attribute information (based on ISO19110) as a reference.

<gmd:contentInfo>
 <gmd:includedWithDataset>
   <gco:Boolean>0</gco:Boolean>
 </gmd:includedWithDataset>
 <gmd:MD_FeatureCatalogueDescription>
   <gmd:featureCatalogueCitation xlink:href="http://########:####/layers/path/feature_catalogue"/>
 </gmd:MD_FeatureCatalogueDescription>

I am trying to upload this metadata file in my layer on QGIS using the
MetaTools plugin. Most of the metadata (title, abstract etc.) are uploaded
correctly but all the attribute information described in the "feature_catalogue" are not loaded.

As it can be seen in the screenshot the Content Information section is empty.

Is there a way to upload metadata in a layer in QGIS and include also the
feature's attributes?


texts:
I am building an app in the ArcGIS Web AppBuilder.
I am trying to pass the geometry resulting from the drawbox dijit to a geoprocessor. Its an extent.
But when I attempt it and log the parameter to console from the geoprocessing script its always empty.
I get the geometry from the drawbox dijit like so :-
_onDrawComplete: function(graphic) {
  var geometry = graphic.geometry;
  this._bathyList(this.bathyList, geometry)
},

I then want to pass this geometry into a geoprocessor. I tried passing it in directly and it complained it needed to be a gpfeatureset.
So I have attempted to create a gpfeatureset which seems to have got round the above error.
The problem is that when I log the param out on the geoprocessing end it is always null.
var features = [];
  features.push(geometry);

  var featureSet = new FeatureSet();
  featureSet.features = features;

  console.log(geometry)
  console.log(featureSet)

  // setup gp stuff and params
  var params = { "boolList":boollist, "poly": featureSet };
  var gp = new
  Geoprocessor('https://arcgispreprodent.admiralty.co.uk/server/rest/services/Jim/BathyClip/GPServer/BathyClip')

  gp.submitJob(params, this._processresults, this._jobstatus, this._joberr)
},

I am not sure how to debug this now, I assume I must be building an invalid gpfeatureset but it doesnt error.
I can change the datatypes, I can chacge the gp script, I just dont know how to get it to pass the polygon.
Heres my gp script
boolList = arcpy.GetParameterAsText(0)
poly = arcpy.GetParameterAsText(1)

arcpy.AddMessage('Param 0 : ' + boolList)
arcpy.AddMessage('Param 1 : ' + poly)

Output_Zip_File = os.path.join(arcpy.env.scratchFolder, 'output.zip')

tiflist = [u"Whole_UK_Mosaic", u"Whole_UK_COMORG"]

def cliptif(clippee, clipper, output):
    desc = arcpy.Describe(clipper)
    extent = desc.extent
    arcpy.Clip_management(clippee, str(extent), output, in_template_dataset=poly, clipping_geometry="ClippingGeometry")

for i in range(len(tiflist)):
    if boolList[i]:
        outTif = os.path.join(arcpy.env.scratchFolder, 'output'+str(i) + '.tif')
        cliptif(tiflist[i], poly, outTif)

zipper = zipfile.ZipFile(Output_Zip_File, 'w')
for i in range(len(tiflist)):
    zipper.write(os.path.join(arcpy.env.scratchFolder, 'output'+str(i) + '.tif'), tiflist[i]+'.tif')
zipper.close()

for i in range(len(tiflist)):
    if boolList[i]:
        outTif = os.path.join(arcpy.env.scratchFolder, 'output'+str(i) + '.tif')
        arcpy.Delete_management(outTif)

arcpy.SetParameter(2, Output_Zip_File)

I am not sure why my second parameter is always null when logging from the geoprocessor, I might be falling foul of asynchronous stuff somewhere or building feature set wrong.


texts: Layers fail to draw when generating tile cache in AGS

texts:
I recently created an MXD-based map service on AGS 9.3.1 SP2. The map seems to render fine as a dynamic service, but when I convert to a cached service and attempt to generate tiles, the operation fails:

ERROR 000569: Failed to generate cache. All server contexts failed to cache map.
Failed to execute (ManageMapServerCacheTiles).

The following lines show up in the server error log:

<Msg time='2011-06-30T17:27:25' type='ERROR' code='103800' target='service.MapServer' methodName='Map.Draw' machine='hostname' process='3248' thread='2756'>Streets 1-2: The operation was attempted on an empty geometry.</Msg>
<Msg time='2011-06-30T17:27:25' type='ERROR' code='10003' target='service.MapServer' machine='hostname' process='3248' thread='2756'>Failed to draw layer Streets 1-2</Msg>

Frankly I'm not sure where to go from here. As I said earlier, the service renders fine in a dynamic configuration. Any ideas?


texts: How can I render heatmap contours from QGIS into tiles?

texts:
I have a heatmap and contours in QGIS, and I want to export these heatmaps as tiles. This works well with tiff, but I want to include contours into images, to make them look better. In QGIS, I can make a contour layer, but it can't be exported into raster format, only in vector.

How can I make a raster countour layer or vector tiles with QGIS or linux-based gis tools?


texts: Selecting certain fields from feature class to create new feature class in ArcPy?

texts:
I am stuck on how to create multiple features classes from one feature class using a select few fields.
I have the following attribute table for a shapefile.

I would like to select the 'USAF', 'COUNT', 'ADMIN_1' fields and append one year field from 2000, 2001, 2002...2014), such that the attribute fields of the new feature classes will be:

'USAF', 'COUNT', 'ADMIN_1', '2000'

'USAF', 'COUNT', 'ADMIN_1', '2001'

'USAF', 'COUNT', 'ADMIN_1', '2002'
. . .

Till 2014.
I have been trying to work it out using modelbuilder and a python script (separately) and haven't been close to successful. Here is a script I wrote to test on the first year, but the resultant shapefile is identical to the input
 #import the arcpy site-package
 import arcpy
 from arcpy import env
 env.workspace = "Z:\\GEOG517\\siu850437305\\Project\\Results" #set current workspace
 env.overwriteOutput = True #allow outputs to be overwritten

 #List variables
 in_file = "Z:\\GEOG517\\siu850437305\\Project\\shapefiles\\sahel_ppt.shp"

 arcpy.MakeFeatureLayer_management(in_file, "stationslyr.shp")# Make a layer   from the feature class
 stations = "stationslyr"

 #Create field list iterate through using for loop
 fields = arcpy.ListFields(stations)
 for field in fields:
     for field.name == '2000':
         if field.name == "2000":
             arcpy.SelectLayerByAttribute_management(stations, "new_selection" )
             arcpy.CopyFeatures_management(stations, "yr2000")# Write the selected features to a new featureclass
             print ("Created yr2000")
             break

Then tried to append the years on using:
field_names = []
fields = arcpy.ListFields(stations)
for field in fields:
     field_names.append(field.name)

but that failed too!
I am not even sure there is a tool to do this. I am new to python scritping for geoprocessing tasks so I am not looking for a really complex method just a simple and straightforward way to do this.


texts:
I am trying to create a pattern in QGIS that could look like a wall with several bricks but I'm not quite sure I understood the custom dash pattern tool.

When I try to do something like this, the pattern never gets perfect, see the image below. Is there anyone who could advise me something?


texts: How to use St_intersects with different geometry type

texts:
I try to run this really simple spatial query with PostgreSQL+Postgis extension.

CASE 1:

SELECT a.id, b.id
FROM PolygonLayer1 a, PolygonLayer2 b
WHERE St_intersects(a.geom,st_centroid(b.geom))

CASE 2:

SELECT a.id, b.id
FROM PolygonLayer1 a, PolygonLayer2 b
WHERE St_intersects(st_centroid(a.geom),b.geom) -- This time I get the centroid of PolygonLayer1

Both PolygonLayer1 (~8000 features) and PolygonLayer2 (~10000 features) are polygons (or Multipolygons).

The case 1 run in less than 10 seconds
The case 2 never ends

I really don't understand the difference between the 2 cases.

So the question: What is the difference between the 2 cases ? How can I get a working query for both cases ?

Specification:

If I run: EXPLAIN CASE1, I get a cost of ~20'000 "operations"

"Nested Loop  (cost=0.00..21878.84 rows=196 width=8)"
"  Join Filter: _st_intersects(st_centroid(a.geom), b.geom)"
"  ->  Seq Scan on superposition_bruit a  (cost=0.00..6089.04 rows=8804 width=728)"
"  ->  Index Scan using affectation_com_old_geom_gist on affectation_com_old b  (cost=0.00..1.53 rows=1 width=591)"
"        Index Cond: (st_centroid(a.geom) && geom)"

========================================================================

If I run: EXPLAIN CASE2, I get a cost of ~40'000'000 "operations"

"Nested Loop  (cost=0.00..40530509.25 rows=196 width=8)"
"  Join Filter: ((a.geom && st_centroid(b.geom)) AND _st_intersects(a.geom, st_centroid(b.geom)))"
"  ->  Seq Scan on superposition_bruit a  (cost=0.00..6089.04 rows=8804 width=728)"
"  ->  Materialize  (cost=0.00..3813.98 rows=13332 width=591)"
"        ->  Seq Scan on affectation_com_old b  (cost=0.00..2744.32 rows=13332 width=591)"


texts:
I am trying to use the orthophotos as my background in project. I have used gdal2tiles to prepare the tiles for xyz. Then in the next step, I have uploaded the folder in our server. Then I make a new connection to xyz tiles with the following example link--->   https://xxxxx.xxxx.solutions/planning_storage/xxxxx/{z}/{x}/{y}.png
The xyz is made but I can not read it in my QGIS.
I am using QGIS 3.16 in Windows 10.


texts: Seeking native spatially-enabled SDK to deploy GIS data on iOS and Android mobile devices?

texts:
I am only interested in iOS and Android options that allow me to build native solutions to render/collect data. Response time is very important, since I want to be able to draw many geometries - not just points. Extra points for automatic clustering options.
Licensing is very important - if it makes my resulting product fully Open Source (for examples, forces me to GPL my resulting code) then it is not an option.

Route-Me (iPhone/iPad): Slippy Map to show tiles.
Spatialite (Android/iPhone/iPad) GIS spatially-enabled backend.
Couchbase Mobile (Android/iOS) Library to do synchronization (disconnected editing) with CouchDB NoSQL Databases. GeoCouch provides spatially-enabled backend.
ArcGIS for Smartphones SDK  (Android/iOS) All-encompasing library.
Cloudmade SDK (Android/iOS) displays their own layers with some minor options to overlay your own. I wasn't sure I should list it, but this is pre-emptive.
Every platform has their own limited options (allow to put some basic markers) out of the box. This includes MapKit and the Android Maps.

I did not list gvSIG Mobile because of the licensing (GPL) and that the iOS portion is under development. Also, OpenLayers and Sencha Touch are web-based, so not exactly what I am looking for since IMHO, the UX is still lacking for current-generation Smarphones. And any of the web-to-native compilers like Titanium Appcelerator and similar are usually just very simple wrappers to the native options.

What good options that are not listed do I have? Is there a comprehensive list anywhere?


texts:
I tried to  work on TIMESAT but always facing a problem.
Here is the method I have used. Seems correct, but I get a warning at the end.

I have created a NDVI file from ENVI 4.8 >> Transform >> NDVI >> saved Output data type as Floating Point. Next I use Basic Tools >> Stretch Data >> Stretch Type Linear >> Stretch range by value -0.036 to 0.833 which is my NDVI range to Output DATA Range to 0 to 255. and saved the DATA Type as Integer as Timesat requires 8 bit integer(0-255).

I did the same procedure for my 5 NDVI files, except I used a different Stretch range by value X to Y, based on their own range.

Now I try to open them in TIMESAT >> TSM GUI, I created a txt file with the absolute path like these 5:

C:\timesat31\NDVI\trial\ndvi_1
C:\timesat31\NDVI\trial\ndvi_2
C:\timesat31\NDVI\trial\ndvi_3
C:\timesat31\NDVI\trial\ndvi_4
C:\timesat31\NDVI\trial\ndvi_5

They belong to 1 year data.
Number of Row 1461 Number of Col 2255

Row to process 500 to 500

Column to process 500 to 500

but I get a following warning:

  Error while evaluating unicontrol Callback

Could anybody help me out?


texts:
I'm working a set of 30meter DEMs (~40 DEMs) for New England. I've seen a lot of debate on whether to work with (A) the individual DEMs first, and then mosaic the final raster calculations together, or (B) mosaic the original DEMs and then run all other processing on the new mosaic. I'm trying to determine what the best approach is for my specific project.

My eventual goals are as follows:

Calculate slope for all of New England
Calculate hillshade for all of New England
Calculate aspect for all of New England
Compare slope, hillshade, and aspect with 3 other regional layers (not derived from DEMs)

Would it be best to calculate slope, hillshade, and aspect for each DEM and have separate rasters that then can be mosaicked?

Or mosaic first and calculate on one giant DEM mosaic?

Because the spatial reference of the original DEMs is in UTM, I'm also running into some projection issues (mostly resolved in previous question: mosaic dem with multiple coordinate systems). My data spreads across 2 zones; UTM 19N and UTM 18N.

I have a VERY basic understanding of GIS (using ArcGIS 10.2).

I've mosaicked first and run into some issues with changed values between the original DEM and mosaicked DEM. I just did a small mosaic (2 30m DEM tiles) to see what's happening with the Python code. The code is pasted below:

arcpy.MosaicToNewRaster_management(input_rasters="'F:/MyBook Backup/All Data Files/GIS Data/Geodatabases/Unprojected2.mdb/Pned30m41073';'F:/MyBook Backup/All Data Files/GIS Data/Geodatabases/Unprojected2.mdb/Pned30m41074'",output_location="F:/MyBook Backup/All Data Files/GIS Data/Geodatabases/Unprojected2.mdb",raster_dataset_name_with_extension="Pned41073n074",coordinate_system_for_the_raster="PROJCS['USA_Contiguous_Albers_Equal_Area_Conic_USGS_version',GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Albers'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['central_meridian',-96.0],PARAMETER['Standard_Parallel_1',29.5],PARAMETER['Standard_Parallel_2',45.5],PARAMETER['latitude_of_origin',23.0],UNIT['Meter',1.0]]",pixel_type="32_BIT_FLOAT",cellsize="#",number_of_bands="1",mosaic_method="BLEND",mosaic_colormap_mode="MATCH")


texts: Convert ArcGIS network dataset 10.1 to 10

texts:
I am using ArcMap 10 and have a network dataset made with ArcGIS 10.1.

I can not open the network dataset for the version of my Arcmap.

Is there a tool or way to convert the network dataset from 10.1 to 10 without having to redo the network dataset?


texts:
I am trying to use an ArcGIS (10.2) personal database MDB as a container for data to pass to a QGIS user. I formerly used QGIS 1.8 for Mac OSX, and after some fiddling around, I was able to open a personal geodatabase and load one or more features from it.
I just upgraded to QGIS 2.0 and upon trying to check my MDB test file before sending to my colleague, I am unable to load the geodatabase. In fact, personal geodatabase as a file type no longer appears in the available types QGIS is prepared to open in the File Type drop down box in the load vector data operation.
Is there some setting I need to tweak or an additional plugin I need to install in order to regain this functionality?


texts: I know programming in Java and I am familiar with GIS concepts as well. I do GIS programming using javascript and PHP. I want to now do it using Java. How should I go about it? PLease let me know where I should begin with? Links to tutorials or simple Java based GIS project works will really be appreciated!


texts: For every point, distance to nearest polygon of type X?

texts:
In ArcGIS Desktop, can you calculate as an attribute value the distance to the nearest polygon of type X?

Assume we have a point layer of wildlife observations and we want to know how far each of them is from the nearest wetland.

Your answer should mention whether it is specific to ArcGIS Desktop v9 or 10, and the necessary license level.


texts: Generate vertices on line with respect to point on line

texts:
I have a road layer in line and point layer snap over line i.e. road layer.
I had use QGIS 3.4 "snap geometry to layer" but no vertices are generated on the line layer.
Now I want to generate vertices on line with respect to point feature.
How to do it in QGIS?


texts:
I am using arc 10.2. Every time I add an array both python 2.7 and as a tool in Arc I am getting a "TypeError: Abstract Class". I have also tried changing my script multiple times to make the array differently but every time I call on the arcpy.Array() I get this error.

I have run the script on other computers with Arc and it works perfectly so I can't see an error with the script.

I tried uninstalling and reinstalling ArcGIS. I tried reverting back to 10.1 and 10.0 and I still get the same error.


texts: Subdividing square buffer in sub-units using QGIS?

texts:
I've been trying for a while now to create a grid inside a square buffer, so that it can be divided into a pre-defined number of sub-units. Creating a grid that covers the whole area is not ideal because the grid's squares will not be properly aligned with the square buffers and I need that in order to be able to count how many grid's squares are occupied by houses inside the buffer. I have 150 buffers in a single layer, and i'd like to have the grids only inside them.

I am still quite new to GIS.


texts: Transfering ArcGIS custom symbology to another machine

texts:
I created custom symbols on my PC at home using:

Layer Properties dialog box > Symbology tab > Advanced > Rotation

Now I want to transfer the symbology settings to another machine on campus.

How can I do this?

I can't locate a file which contains the setting.


texts:
We have a map rendering engine (renders from MapInfo Interchange Format). This engine uses Mercator Projection. Now we have to insert a raster layer (orthophotos of a city), we have GeoTIFF files with the .tfw files.

The idea: we create ~200x200px png files for each zoom level and the client will ask a WCF service for the images with parameters: zoom level, top left coordinate, bottom right coordinate.

My problem is: how to make the vector and raster layers work together? I don't really want to change the aspect ratio of the raster images.

Thanks in advance, Mart


texts: Updating polygon attribute based on mode value from cluster of points using PostGIS?

texts:
I'm updating a PostGIS table of polygons denoting property boundaries with some values that form part of a code for sites scattered across the properties. Neither dataset is 100% accurate so I can't rely on a direct attribute transfer based on overlap; instead I need to find a way to take the most common value from the points and apply them to the polygon. Attaching an image to clarify:

That image also shows an example where a polygon needs to be defined for EY-183 but that's another matter.

Any ideas on how to do this, either in PostGIS?


texts: can we update imagemosaic layer with geoserver rest api?

texts:
I am new to GIS and geoserver. I am trying to add and modify imagemosaic layers to geoserver with geoserver rest api. Actually I could add the layers. But I have modify them with rest api while adding or after adding. The parameters I have to set is:

Layer name, srs handling, declared rss, background values, use jai image read

Can I set these values with rest api? If there is documentation about it i will read but I couldn't find any. And if we can set these values I have to solve out how?

Thanks,,


texts: How can I access components from a spatial list and create buffer (width =0) to fix invalid geometries in R?

texts:
I'm creating some concentric buffers up to 70 km away (every 10 km) around polygons of a SpatialPolygonDataFrame.
illustrative picture

I used this code to create the buffers:
#I created a numeric vector with the buffers widths that I need#

     buffering<-c(seq(10000, 70000, by=10000))

     #I creted a vector of characters

     buff_labs <- c(seq(1,10, by=1)) %>% paste0("ring", .,"km",US)

      #I associate the character vector with the numeric vector

      names(buffering) <- buff_labs

      #cal my data
      setwd("~/prodes/min")

        minas <- readOGR("minas.shp")

        # this function generates a list with the seven buffers that I need (70km, 20km, 30km, 40km,
          50km, 60km e 70km).

        buffers <-mapply(FUN = rgeos::gBuffer,
                          width = buffering,
                          MoreArgs = list(spgeom =minas, byid=TRUE),
                          SIMPLIFY = FALSE,
                          USE.NAMES = TRUE,
                          )

Until here the code works very well for all my shapesfiles
But I only need the 10 km wide ring for each distance and I used this code to remove the buffers centroids:
        centers <- c(minas, buffers[1:(length(buffers)-1)])

        rings <-mapply(FUN = raster::erase,
                       x = buffers,
                       y = centers,
                       SIMPLIFY = FALSE,
                       USE.NAMES = TRUE)

This code works very well in seven of eight shapefiles that I have, but in one generates this error:
RError in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, unaryUnion_if_byid_false, :
TopologyException: Input geom 0 is invalid: Self-intersection at or near point -669458.50954959029
2325717.8159672269 at -669458.50954959029 2325717.81596...

But, I found a possible solution.
After creating the buffers, I have to create a new buffer (width = 0), in each buffer, because this error is usually related to invalid geometries generated when I made the buffers.
But I don't know how to put that in my code because, when I crieating the buffers I generate a vectors list with seven components. And I need to modfy one by one with (width =0) to corect the geometries.
And now, I really don't know how to do this.


texts: What does the ST_Union of the geometry column of two tables produce?

texts:
I'm working in Postgres 9.5. I have two tables, each with a column geometry(MultiPolygon,4326).

I want to creat a single MultiPolygon that is the union of every polygon in the two tables.

I have tried running ST_Union, but it isn't doing quite what I expected. If I do this:

SELECT
ST_GeometryType(ST_Union(table1.geom, table2.geom))
FROM table1, table2;

This produces lots of MultiPolygon results, instead of a single one, which suggests to me that the result of a simple ST_Union is not a single polygon.

What is ST_Union(table1.geom, table2.geom) actually giving me, conceptually? Why are there multiple polygons?

What do I need to do to get a single multipolygon - something like ST_Union(ST_Union(table1.geom), ST_Union(table2.geom))?


texts: Measure my Cats movements , Accuracy of Mobile phone GPS at metre scale . how to improve it

texts:
I wanted to get some information related to an idea I have
I have a lot of cats and a big garden.

The cats generally roam through the garden and this location is fixed and Known

If i have known area , in this case my garden , does this make getting location information GPS , more accurate

Can i track my Cats accurately if i hook them up to a mobile phone.(Picture 1)

I have found the thread here about GPS accuracy What is the maximum Theoretical accuracy of GPS?

How accurate will this be , to 2 metres is this relatisic with mobile GPS

If i had more than one cat with a mobile , Could these devices talk to one another to give an more accurate picture ( picture 2 )

If i was able to put some fixed mobile reciever points in Could be combined with cats mobiles to give more accurate read out  ( picture 3)

The outputs of what i like to be able to Capture

xy information per second for my Cat

How accurate would this be .
Could this be done , by only using mobile phones. Are there better options

Picture 1

Picture 2

Picture 3


texts: Using PHP and JavaScript to dynamically create shapefile or KML

texts:
While I have used GIS software (ArcGIS and just started playing with QGIS at home), I have never looked into the code or file structure of shapefiles.
I am working on a project for my work, a  web-based application running on a LAMP server and written in PHP and JavaScript. The idea is a search box that will be used to search the "geographical hierarchy" of a specific location. This would allow my users to type in a city, and it would return the state, country, and continent in which the city is located. Google can provide this information and more, but this will be used on a closed network so that's not an option. I am using a MySQL database containing all continents down to the city level.
I still have a little way to go, but I know what needs to be done to finish this part.
I recently thought it would be nice to add the ability for my users to enter a geo-coordinate and get the same result as if they had entered the city the coordinate is located in. Right now I do not have any coordinates in my database, but I know I have to reference them somehow. I believe there is an SDE database on this network, but i have not tried to access it or see what it would take to pull data from it. I would like to just reference it to get the geo data rather than duplicating the data in my database, if possibly (to keep my database as small as possible).

Another feature I would like to provide (now, if it is not too difficult, or as a future update) is the ability for my users to dynamically build a KML (as a lot of users would be using GoogleEarth) or a shapefile based on the users' specific request (single country/multiple countries, states or cities within countries, etc).
To build the KML, I know that I would need a list of decimal lat/long for the shape of the polygon. I am not sure what would be needed for the shapefile.
My reasoning for this feature is to allow GoogleEarth users the ability to quickly/easily generate a KML and, for GIS users, my tool could just be another front end to the SDE to generate specific shapefiles.

So, specifically:

What do I need to do to be able to look up a location based on a coordinate located within that location?
What is the best way to dynamically create a KML?
How can I dynamically create a shapefile keeping the data from the original database?

Since this is a closed network, I am limited to the libraries and scripts that other people have already added. I am hoping that i can code most of this by hand or transcribe it, and I would be most comfortable using PHP and/or JavaScript. However, if there is a specific script or library that I would need, there is a process to get stuff moved to this network, but it is not a simple, and then it has to be approved, which is its own process.

As this is a work in progress, I am open to suggestions to improve what I have done so far or what my thinking for the rest of my program.


texts: ogr2ogr kml to shp conversion created 184 shapefiles?! What am I missing?

texts:
I'm trying to convert this KML dataset of ice arenas to a shapefile using ogr2ogr (GDAL 1.8). After some initial troubleshooting to get the script running, it's creating 184 different shapefiles (736 unique files including required siblings)!

My ogr2ogr script is essentially this, with some extra flags. The -append and -skipfailures flags are necessary (not sure if precision or geometry are doing anything, but I'm pretty sure they don't hurt). I'd prefer not to use -skipfailures, as clearly it means losing data, but without it the script won't finish:

ogr2ogr -f "ESRI Shapefile" "E:\4_GIS\HockeyArenas\shp\ice_rinks.shp" "E:\4_GIS\HockeyArenas\doc.kml" -lco PRECISION=false -nlt "geometry" -append -skipfailures

Ultimately I want to move this data into PostgreSQL/PostGIS, but if it's making 184 shapefiles, I don't want to litter my PostGIS db with 184 tables..

Anyone know how to get one shapefile, preferably without using -skipfailures?

Thanks, community.


texts:
I would like to use QAD plugin to emulate AutoCAD, especially to draw some dimensioning. There are two possibilities to do so: click on 2 points or select an object.

I want to use the object selection on a line with more than 2 vertices, I would like something like that, but with the actual length of the line (currently the length is the one between the start and end points of the line, and I did it manually with the first method) :

And when I try to use the object selection, I get a result limited by the only two vertices of the section I clicked, as such (I clicked on the almost horizontal section):

Would someone know how I could proceed to solve that?


texts:
I have a world map which I assumed used a Mercator projection but when I tried to plot points on it using a formula for the Spherical Pseudo-Mercator projection (from here: http://wiki.openstreetmap.org/wiki/Mercator#Python) but the Y coordinates seem off. How can Identify which projection is it using? Can someone provide pseudocode/Python for transforming lat,long pairs to x,y pixels for it?
If it may help, this image comes from a KML file.

  <GroundOverlay>
    <name>1986-2010</name>
    <altitude>0</altitude>
    <altitudeMode>clampToGround</altitudeMode>
    <Icon>
      <href>KG_1986-2010.png</href>
    </Icon>
    <LatLonBox>
      <north>90</north>
      <south>-90</south>
      <east>180</east>
      <west>-180</west>
    </LatLonBox>
  </GroundOverlay>


texts: Is inverting slope map same as inverting surface in ArcMap?

texts:
Is inverting a slope map the same as inverting a surface?
I found code on this website https://support.esri.com/en/technical-article/000006694
The code is following:
(([elevgrd] - Z_Max) * -1) + Z_Min

Can I use slope map instead of elevgrd ,maximum slope value instead Z_Max and minimum slope value instead Z_Min?
I want to do something like least cost path analysis in ArcGIS but my main goal is to find the hardest path between two points.


texts: Legend styling in QGIS: How to create overlaying points based on size

texts:
I have a categorized point layer, were the size of the points is based on a feature. Normally in the legend, there would appear multiple points with different sizes. However they consume a lot of space. I would like to have the points one on top of another (see picture below).
Is that possible without having to do it all manually?


texts: What is the syntax for gdal_calc.py to return a certain value?

texts:
I'm trying to return the values from a slope raster where the values = 0.

I am getting 'invalid syntax' returned. This is what I'm trying

gdal_calc.py -A slope.tif --outfile=result.tif --calc="A=0"

I'm trying gdal_calc.py because the internal QGIS calculator is only returning NaN values no matter what I do.


texts: Tweaking location of zoomControl (cartodb.js)

texts:
I have succeeded in creating a header outside of the CartoDB UI, but the header is blocking the zoom controls, as can be seen below.

Is there any way to modify the zoom control location to adapt to the existence of this header?


texts: Use Proj (Proj4) to convert UTM coordinate to Cartesian coordinate relative to an Earth tangent plane with its origin define as lon/lat?

texts:
I may have asked in the wrong place here, https://stackoverflow.com/questions/61147814. Then I discovered http://gis.stackexchange.com, so I decided to ask again here:

I have a UTM coordinate (zone, hemisphere, northing, easting, height) where height is in meters.

I want to convert this coordinate into a custom format (or vice versa) where I choose an arbitrary plane that is tangent to the Earth. The tangent point of the plane has a longitude/latitude with height of zero (it "touches" Earth at that point).

The X axis of the tangent plane points North, Y points West, and Z points up. I refer to this coordinate system as "North/West/Up" or "NWU".

Once I have the tangent plane, I define all points relative to the plane's origin (tangent point) in NWU coordinates of meter units.

An NWU coordinate basically looks like (lon, lat, X, Y, Z) which can also be referred to as (lon, lat, N, W, U) where (lon, lat) describes the touch point of the tangent plane and the origin from which all NWU points are relative to.

This is what the tangent plane looks like:

Given a UTM coordinate (zone, hemisphere, easting, northing, height), how can I use PROJ to convert it to/from the NWU coordinates of (lon, lat, N, W, U) where (lon, lat) is predefined? I.e. What would the PROJ strings for this be?

If the (lon, lat) are hard-coded into the strings, that's fine. I can template them depending on which tangent plane I choose.

I don't know if the latitude is geocentric or geodetic. For this case assume the same Lon/Lat as what Google Maps and most other mapping systems use.


texts: Unsplit/Dissolve multiple touching lines in Stream Network using ArcGIS Desktop?

texts:
I am trying to dissolve multiple lines for stream networks into individual features. The dataset contains many stream networks so a complete dissolve of the whole dataset wont do as I would like to only dissolve streams that are in the same network (illustrated in the pictures below).

I've tried using the Unsplit Line tool and the Dissolve tool in ArcGIS 10. The tool will only join two lines together, i.e. lines that share a common end point. This is illustrated in the screenshots below.

Here is the orignal dataset with the lines split, a single line feature is highlighted

This is the same feature highlighted after using the unsplit/dissolve tool. As you can see it has only joined two lines together rather than all of the lines that are touching.

The features in the same network do not have a common attribute either.


texts: Adding __init__.py in eclipse python interpreter' libraries

texts: I am trying to Make eclipse understand the API in of QGIS. the __init__.py module is available in this address: C:\Program Files\QGIS Pisa\apps\qgis\python\qgis\__init__.py but when I try to add it in python interpreter.libraries tab, it can not be found and I just can add the following address: C:\Program Files\QGIS Pisa\apps\qgis\python\qgis I want to know is it right or I should find __init__.py to add it in Libraries?


texts: QGIS standalone application in python can't access 3d:tessellate processing tool

texts:
I'm developing a standalone QGIS application and I want to run the tessellate (3d:tessellate) processing tool using the following code:

import processing
from qgis.analysis import QgsNativeAlgorithms
from qgis.core import QgsApplication
from processing.core.Processing import Processing

# Supply path to qgis install location
QgsApplication.setPrefixPath("C:\\OSGEO4~1\\apps\\qgis", True)

profile_folder = '.'
# Create a reference to the QgsApplication.  Setting the second argument to False disables the GUI.
app = QgsApplication([], False, profile_folder)

# Load providers
app.initQgis()

Processing.initialize()
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())
processing.run("gdal:aspect", {'INPUT': 'D:/OneDrive/Personnel/Daniel/QGIS/Kingston/Kingston_sub9.gpkg|layername=Road', 'OUTPUT': 'TEMPORARY_OUTPUT'})

# Stop QGIS appllication
app.exitQgis()
app.exit()

I received the following error message saying it cannot find the algorithm:

C:\OSGeo4W64\apps\Python37\python.exe C:/Users/berge/PycharmProjects/qgis_geo_sim/runner_chordal_axis.py
C:\OSGeo4W64\apps\qgis\python\qgis\utils.py:744: DeprecationWarning: the imp module is deprecated in favour of importlib; see the module's documentation for alternative uses
  mod = _builtin_import(name, globals, locals, fromlist, level)
Traceback (most recent call last):
  File "C:/Users/berge/PycharmProjects/qgis_geo_sim/runner_chordal_axis.py", line 19, in <module>
    processing.run("3d:tessellate", {'INPUT': 'D:/OneDrive/Personnel/Daniel/QGIS/Kingston/Kingston_sub9.gpkg|layername=Road', 'OUTPUT': 'TEMPORARY_OUTPUT'})
  File "C:\OSGeo4W64\apps\qgis\python\plugins\processing\tools\general.py", line 106, in run
    return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)
  File "C:\OSGeo4W64\apps\qgis\python\plugins\processing\core\Processing.py", line 125, in runAlgorithm
    raise QgsProcessingException(msg)
_core.QgsProcessingException: Error: Algorithm 3d:tessellate not found

But I can run the same processing tool code in the QGIS python console and it's working fine:

processing.run("3d:tessellate", {'INPUT':'D:/OneDrive/Personnel/Daniel/QGIS/Kingston/Kingston_sub9.gpkg|layername=Road','OUTPUT':'TEMPORARY_OUTPUT'})

When I run the following code in the QGIS python console I can "see" what are the runnable processing tools and of course 3d:tessellate is there

for algo in QgsApplication.processingRegistry().algorithms():
    print(algo.id(), "------", algo.displayName())

3d:tessellate ------ Tessellate
gdal:aspect ------ Aspect
...passing 950 lines...
saga:zonalmultipleregressionanalysispointsandpredictorgrids ------ Zonal multiple regression analysis (points and predictor grids)
saga:zonalrasterstatistics ------ Zonal raster statistics

When I run the almost same code in a QGIS stansdalone script I have the following output and do not see the 3d:tessellate processing tool.

gdal:aspect ------ Aspect
gdal:assignprojection ------ Assign projection
... passing 949 lines...
saga:zonalmultipleregressionanalysispointsandpredictorgrids ------ Zonal multiple regression analysis (points and predictor grids)
saga:zonalrasterstatistics ------ Zonal raster statistics

What is wrong? Why I cannot "see" and use the tessellate "3d:tessellate" processing tool in a QGIS standalone application?

I used the following for this test:

OS: Windows 10
Version: QGIS 3.12
Installation: OSGeo4W64


texts: Javascript Api's ArcGISTiledMapServiceLayer won't show map

texts:
I am currently trying to show a map that has been tiled. I am using version 10 of the Arcgis server, and version 2.5 of the Javascript api.

The odd thing I am running into is that the map will not show when using the constructor arcGISTiledMapServiceLayer but it will show when I'm using the alternate constructor arcGisDynamicMapServiceLayer.

Here is a very simple code snippet i'm using to test out my functionality

dojo.require("esri.map");
dojo.require("esri.tasks.geometry");
var map = null;
var gsvc = null;
var pt = null;

function initialize() {
  map = new esri.Map("map");

  var layer = new esri.layers.ArcGISTiledMapServiceLayer("http://server.arcgisonline.com/ArcGIS/rest/services/ESRI_StreetMap_World_2D/MapServer");

  alert(map.toString());

  map.addLayer(layer);

}

dojo.addOnLoad(initialize);

As you all have no doubt realized, this example code is using arcgisonline.com instead of my actual service that I'm running.

There is one more very important piece of information, this code above, works it shows the map but the same code with my url doesn't leading me to believe something is wrong with my server setup, but I am fairly new to arcgis so I'm not leaving any stones unturned.

If there are any querstions please feel free to ask.


texts: Is an MPK file with Symbology from a TrueType Font supported on different machines?

texts:
I have exported an ArcMap 10.7.1. MXD map document to mpk for use on another server. The Font has been installed on the other System and can be viewed. However, the map symbols are not displayed correctly.
Do I need to restart the system on the other server for the layer to be rendered correctly or is another problem at play?


texts: How can I join Excel sheet with point shapefile?

texts:
How can I join an Excel sheet with shapefile (Point) and common field is station code?

In Excel one station code has multiple rows and I want to attach those values with shapefile.


texts: How do "join field" and "target field" work?

texts:
New user. This should be simple, but I can't figure out how to perform this basic join in QGIS. (It works in Cartographica.)

My shapefile (IL_BNDY_County_Py) has two fields/columns in its attribute table: COUNTY_NAM and CO_FIPS.

My source DBF file (IL_BNDY_County_Py-species-dbf4.dbf) has three fields/columns: COUNTY_NAM, CO_FIPS and SPECIES. (The file is a literal clone of the shapefile table, with the addition of "SPECIES".)

Both are available as layers in my project.

I want to add the SPECIES field to the shapefile.
Either COUNTY_NAM or CO_FIPS could serve as "key" field (to use terminology I'm familiar with from other software).

When I try to create a join, the Add vector join dialog shows the following options:

Join layer: IL_BNDY_County_Py-species-dbf4.dbf

Join field: (three menu options): COUNTY_NAM, CO_FIPS or SPECIES

Target field (two menu options): COUNTY_NAM or CO_FIPS

Is this what I should be seeing?

I can't figure out how to use this dialog to get the SPECIES field added to the shapefile. Part of the problem is that QGIS doesn't use the the term "key," so it's impossible for me to know which field QGIS wants to use to ensure the relationship. Given the options displayed, neither the "target" nor "join" field seems appropriate to serve as "key." How do I denote the "key" field?

If I choose SPECIES as the Join field, and CO_FIPS as the Target field, the results are duplicate columns for COUNTY_NAM and CO_FIPS with random mix of integer values and NULL values. SPECIES values are ignored and not copied over. I even tried creating a placeholder "species" column in the shapefile attribute table but still couldn't figure out how to get SPECIES copied over to species. What is the proper procedure for doing this? Is there another problem that I should be considering?


texts: Updating X and Y coordinates with UpdateCursor based on another shapefile using ArcPy?

texts:
I have two shapefiles: "broken_cities" and "us_cities."  "us_cities" contains point data representing cities in the US, and the x and y coordinates are correct.  "broken_cities" contains point data representing 10 cities in the US whose X and Y coordinates are incorrect.  There is no X and Y coordinate fields in the attribute table of either.

The goal is to update the broken cities X and Y coordinates using the correct coordinates from "us_cities."  I know I should use an update cursor, but I don't know how to combine the search cursor and the update cursor.  I am aware that my code does not work at all because I'm not sure of the logic of the workflow.

import arcpy

#set workspace
arcpy.env.workspace= r"Y:\Personal\jsmith\EGData.gdb"

fields=["OID@", "SHAPE@X", "SHAPE@Y"] #tokens to access

#create search cursors to look at values
broken_scur=arcpy.da.SearchCursor("cities_broken",fields)
cities_scur=arcpy.da.SearchCursor("us_cities", fields)

#OIDs of the 10 broken cities
broken_city_OIDs=[66, 68, 61, 74, 86, 80, 116, 94, 96, 97]

for row in cities_scur:
    for row1 in broken_scur:
        if row1[0]==broken_city_OIDs:
            row1[1]==row[1]
            row1[2]==row[2]


texts: Do not display specific map feature with CartoCSS

texts:
I have built an OpenStreetMap tile server following this guide and everything seems to work fine.

Now I want to customize the look of the generated tiles, and I started with something I tought it was simple: remove mountain peaks from the map.

I was only half succesful however, because I was able to remove the triangular peak marker, but not the peak name.

According to the CartoCSS documentation you can remove a specific symbolizer from the output using the syntax

symbolizer_name: none

So I started to hack the amentity-points.mss like this:

.points {
    //...

    [feature = 'natural_peak'][zoom >= 11] {
      marker: none;
      marker-file: url('symbols/peak.svg');
      marker-fill: @landform-color;
      marker-placement: interior;
      marker-clip: false;
    }

    //..
}

.text-low-zoom[zoom < 10],
.text[zoom >= 10] {
    //..

    [feature = 'natural_peak'][zoom >= 13],
    [feature = 'natural_volcano'][zoom >= 13],
    [feature = 'natural_saddle'][zoom >= 15],
    [feature = 'tourism_viewpoint'][zoom >= 16] {
      text: none;
      text-name: "[name]";
      text-size: @standard-font-size;
      text-wrap-width: @standard-wrap-width;
      text-line-spacing: @standard-line-spacing-size;
      text-fill: darken(@landform-color, 30%);
      [feature = 'natural_volcano'] { text-fill: #d40000; }
      text-dy: 7;
      [feature = 'tourism_viewpoint'] { text-dy: 11; }
      text-face-name: @standard-font;
      text-halo-radius: @standard-halo-radius;
      text-halo-fill: @standard-halo-fill;
      text-placement: interior;
    }

    //..
}

As I said, the peak marker does not get rendered, but its name does.

Am I missing something like another rule overriding my settings?


texts: Multiple tables from PostGIS to Leaflet with PHP and jQuery

texts:
I recently found this example:
https://northlandia.wordpress.com/2015/04/20/connecting-postgis-to-leaflet-using-php/comment-page-1/#comment-627

It is very helpful and it works fine for one PostGIS table.

My question is:
How can I fetch more than one table using this code?

This is my PHP script:

require('../php/connection.php');

$table = $_GET['table'];
$fields = $_GET['fields'];

$fieldstr ="";

foreach ($fields as $i => $field) {
    $fieldstr = $fieldstr. "l.$field, ";
}

$fieldstr = $fieldstr. "st_asgeojson(st_transform(l.geom,4326))";

$sql = "select $fieldstr from $table l";

//send the query
if (!$response = pg_query($conn, $sql)) {
    echo "A query error occured.\n";
    exit;
}
//echo the data back to the DOM
while ($row = pg_fetch_row($response)) {
    foreach ($row as $i => $attr){
        echo $attr.", ";
    }
    echo ";";
}

and this is my main app:

var map,
    fields = ["gid", "geom"],
    autocomplete = [];
$(document).ready(initialize);

function initialize() {
    map = L.map("map", {
        center: L.latLng(38.368858, 26.137269),
        zoom: 15
    });
    var tileLayer = L.tileLayer("http://{s}.tiles.mapbox.com/v3/mariosbatsaris.ilamhcgb/{z}/{x}/{y}.png").addTo(map);
    //next: add features to map
    getData();
};

function getData() {
    $.ajax("php/getData.php", {
        data: {
            table: "shelters_wgs84",
            fields: fields
        },
        success: function(data) {
            mapData(data);
        }
    })
};

function mapData(data) {
    //remove existing map layers
    map.eachLayer(function(layer) {
        //if not the tile layer
        if (typeof layer._url === "undefined") {
            map.removeLayer(layer);
        }
    });
    //create geojson container object
    var geojson = {
        "type": "FeatureCollection",
        "features": []
    };
    //split data into features
    var dataArray = data.split(", ;");
    dataArray.pop();
    //console.log(dataArray);
    //build geojson features
    dataArray.forEach(function(d) {
        d = d.split(", "); //split the data up into individual attribute values and the geometry
        console.log(d[2]);
        //feature object container
        var feature = {
            "type": "Feature",
            "properties": {}, //properties object container
            "geometry": JSON.parse(d[2]) //parse geometry
        };
        for (var i = 0; i < fields.length; i++) {
            feature.properties[fields[i]] = d[i];
        };
        //add feature names to autocomplete list
        if ($.inArray(feature.properties.gid, autocomplete) == -1) {
            autocomplete.push(feature.properties.gid);
        };
        geojson.features.push(feature);
    });
    //console.log(geojson);
    //activate autocomplete on featname input
    $("input[name=oid]").autocomplete({
        source: autocomplete
    });
    var mapDataLayer = L.geoJson(geojson, {
        pointToLayer: function(feature, latlng) {
            var markerStyle = {
                color: "#FFF",
                opacity: 0.8
            };
            return L.multiPolygon(latlng, markerStyle);
        },
        onEachFeature: function(feature, layer) {
            var html = "";
            for (prop in feature.properties) {
                html += prop + ": " + feature.properties[prop] + "<br>";
            };
            layer.bindPopup(html);
        }
    }).addTo(map);
};

function submitQuery() {
    //get the form data
    var formdata = $("form").serializeArray();
    //add to data request object
    var data = {
        table: "shelters_wgs84",
        fields: fields
    };
    formdata.forEach(function(dataobj) {
        data[dataobj.name] = dataobj.value;
    });
    //call the php script
    $.ajax("php/getData.php", {
        data: data,
        success: function(data) {
            mapData(data);
        }
    })
};


texts: Clipping land cover grid with shapefile in QGIS

texts:
One very common application of GIS to ecology is to layer a set of vector polygons over a land cover grid (or some other land feature raster) and clip it to determine the percentage of the cover classes under each polygon.  Thus, as a new user of QGIS I was kind of surprised that this operation was (apparently) not a readily possible in the QGIS program.

I did perform this operation years ago using ESRI software but I have been away from GIS for some time. For the QGIS program I have researched out several approaches that are claimed to work for this purpose and thought I would collect them here to see if more experienced users would comment on which is the better way to go.

The simplest approach I have found is this suggestion (Clipping raster with shapefile in QGIS?)  Assuming
that the raster data is in a format that QGIS can load, you can do
it inside QGIS with Raster -> Extraction -> Clipper. This is said to
work as long as the clipping polygon shapefile has the same CRS as
the raster?  This approach seems very simple but I cannot say if it
works as I find that I do not have the option
Raster\extraction\clipper on my version of QGIS (2.18.17)?  What is
more I am not sure this will return the number of pixels (or %land
cover class) under each polygon anyway?
One very promising approach seems to be to use the LecoS (land
cover analysis program) developed by Martin Jung. (https://conservationecology.wordpress.com/2013/04/01/another-lecos-update-and-more-to-follow-in-the-near-future/)  This seem to be
exactly the right operation but it might not something anyone can
download and use without some background information?  It appears
that the LecoS plugin requires some routines from the python
libraries including “scipy”, “numpy” and “pil (imaging)” in order to
run on Windows. To install these python routines for Windows one can
download the OSGEO4W Installer which I did (but I've also read that
you must choose the “advanced Install” and I think I used the
"routine" install when I did it) Is it possible to reload the
advanced install of OSGEO4W over top?

That question aside, it is not obvious to me is how proceed with downloading the libraries (python-numpy, python-scipy and python-imaging) then selecting them and thus allowing one to use LecoS to clip and return the information sought.  If anyone knows of a tutorial on how exactly to download and use this promising application I would love to read it.

I’ve read other suggestions but this is probably enough for now!  Land cover analysis is an operation I intend to perform regularly so I'd like to find the simplest approach.

I need to point out that this question has perhaps already been addressed here:Clipping raster with vector boundaries using QGIS?
I missed it earlier and therefore the first suggestion I made mention of (using the GDAL clipper tool) might be one solution to the problem.  I did not find this tool initially because I did not have the "GDAL tools" plugin checked on the list. My error.
However the clipper tool apparently only outputs the "clipped data" to another raster file.  Ideally I would like the output the data to something like a dbase file for further manipulation and analysis.


texts:
I created a plugin and I added my separate python file(added in plugin path itself) in the created qgis plugin. In that separate python file, I called a DLL file like below,

import sys
import clr
class pythonCall:
    def callMethod():
       clr.AddReference(r"D:\\K\\bin\\Debug\\ClassFile.dll")
       from ClassFile import DynamicCalc
       calc=DynamicCalc()
       print calc.getValue()

I called the separate python file in my qgis plugin python like below,

from callCSharp import pythonCall

I got an error as,

  ImportError: No module named clr

If I run the separate python file individually and its works fine. I got an error if I called it through QGIS.

I also tried adding the clr.pyd in qgis python path. it shows minidump error.
Kindly help me to fixed it.


texts: Perform sjoin in geopandas leads to:'AttributeError: 'GeoSeries' object has no attribute 'columns''

texts:
I try to find all the polygons in a dataframe thant intersect an other dataframe (made by buffering my polygons) usign sjoin.

First I did the buffers

gdf['geom_buf'] = gdf.geom.buffer(100)

Then I merged all overlapping polygons, then I transformed it into a geodataframe and exploded the geodataframe

union_gdf = cascaded_union(gdf['geom_buf']) +
union_gdf = gpd.GeoDataFrame(geometry=[union_gdf ])
exp_union_gdf = union_gdf.explode()

Then I tried to find all the buffer that the entities inside gdf intersect

results-intersect = gpd.sjoin(gdf["geom"], exp_union_gdf , how="inner", op='intersects')

I get this error

~\Anaconda3\lib\site-packages\geopandas\tools\sjoin.py in sjoin(left_df, right_df, how, op, lsuffix, rsuffix)
     46
     47     # due to GH 352
---> 48     if (any(left_df.columns.isin([index_left, index_right]))
     49             or any(right_df.columns.isin([index_left, index_right]))):
     50         raise ValueError("'{0}' and '{1}' cannot be names in the frames being"

~\Anaconda3\lib\site-packages\pandas\core\generic.py in __getattr__(self, name)
   4370         if (name in self._internal_names_set or name in self._metadata or
   4371                 name in self._accessors):
-> 4372             return object.__getattribute__(self, name)
   4373         else:
   4374             if self._info_axis._can_hold_identifiers_and_holds_name(name):

AttributeError: 'GeoSeries' object has no attribute 'columns'

I guess that exp_union_gdf is not a geodf but a geoseries.


texts: I have a custom CSV file with the columns UTM easting, UTM northing, red, green, blue. The latter values represent a RGB color of the pixel value. Now I want to convert this CSV to a raster GeoTIFF. I've found that there is a import function in QGIS (Layer -> Add delimited text layer) where I can import the CSV as vector layer. How can I transform it to a raster GeoTIFF?


texts: Trying to Save QGIS layer as a DXF; output is always really tiny

texts:
In multiple files on multiple computers, I've tried to export a layer in my QGIS file as a DXF so that I can open it in AutoCAD.  I've tried this for GPS located trails, contour lines, and a few other things.

I've tried right clicking the layer and using the Save As command and I've also tried the DXF Export from the dropdown menu.

It seems that no matter what I do, items that are scaled properly in QGIS get drastically scaled down when I save to DXF.  Lengths that should measure close to a mile read as .0617 (or something like that) in a drawing where 1 unit = 1 foot.

What really confuses me is that I've tried multiple installs of QGIS on multiple computers and it's the exact same every time, but yet I haven't yet found anyone who has heard of this issue.

I can upload the file to a host if someone would like to dive into it and see if they can figure it out.


texts: Publishing ArcGIS WMTS for non Esri clients?

texts:
We use Arcgis 10.4.1 products (Server, desktop). When we publish a service with
(Google - bing tiling scheme) caching capabilities we can use service in esri clients (JS SDK, Arcmap etc.) but can not use in clients like openlayers. Tiles are shown in completely wrong order.

What should we do to fix this?


texts:
I have way points of my field work, but some of points I saved are not related to my project.

How can i remove them in QGIS?


texts: How to test for previously joined layers before run a code

texts:
I've a function that join two added layers. The code works fine, but now I need to test for previously joined layers before run this function. The function code is right below:

def jointables(dict): # this dict provides [targetLayer:layerToJoin]
  for k, v in dict.items():
    target = QgsProject.instance().mapLayer(k)
    layerToJoin = QgsProject.instance().mapLayer(v)
    fieldToJoin = QgsProject.instance()
    symb = QgsVectorLayerJoinInfo()
    symb.setJoinFieldName('id_feature')
    symb.setTargetFieldName('id')
    symb.setJoinLayerId(layerToJoin.id())
    symb.setEditable(True)
    symb.setDynamicFormEnabled(True)
    symb.setUpsertOnEdit(True)
    symb.setPrefix('')
    symb.setJoinLayer(layerToJoin)
    target.addJoin(symb)

I've tried some things like this:

for k, v in dict.items():
    if QgsProject.instance().mapLayer(k).addJoin(QgsVectorLayerJoinInfo()) == True:
        break
    else:
        continue

or

for k, v in dict.items():
    if QgsVectorLayerJoinInfo().isEditable() == True:
       break
    else:
       continue

among others. I've missing something. Both conditionals has the same value True, doesnt matter if already exists an previously joined layer or not. How can I solve this?


texts: How do I display data from a django rest get request in a react-leaflet using onEachFeature?

texts:
I imported shapefiles into my PostgreSQL database with PostGIS extension and created a REST API using Django which has a GET response like this:
[
{
    "objectid": 1,
    "id_field": 5543,
    "county_nam": "B1",
    "const_code": 288,
    "constituen": "KAM",
    "county_cod": 47,
    "shape_leng": 1.55978672217,
    "shape_area": 0.057422949
},
{
    "objectid": 2,
    "id_field": 6634,
    "county_nam": "B2",
    "const_code": 3,
    "constituen": "KIS",
    "county_cod": 1,
    "shape_leng": 0.88558616055,
    "shape_area": 0.02332510618
},
{
    "objectid": 3,
    "id_field": 6600,
    "county_nam": "B3",
    "const_code": 10,
    "constituen": "KIN",
    "county_cod": 2,
    "shape_leng": 4.28418177601,
    "shape_area": 0.75826601303
},
{
    "objectid": 4,
    "id_field": 6731,
    "county_nam": "B4",
    "const_code": 13,
    "constituen": "KAL",
    "county_cod": 3,
    "shape_leng": 5.33308007101,
    "shape_area": 1.02533838151
}, ...

How do I display this response into a React Leaflet?
   async componentDidMount() {
    const rawcurrentuser = JSON.parse(localStorage.getItem("currentUser"));

    await axios.get('/api/gis/counties/')
      .then((response) => {
        this.setState({ county: response.data });
      })
      .catch(function (error) {
        console.log(error);
      });
  }

When I display as below, the data is not displayed.
<Map center={position} zoom={7} ref={this.mapRef}>
 {this.state.county && (
        <GeoJSON
          data={this.state.county}
        />
      )}
 </Map>

Displaying with Django view works fine using onEachFeature:
myLayer = new L.geoJson(data, {
  onEachFeature: function (feature, layer) {
          layer.bindPopup('<font size="4">' + feature.properties.county_nam)
        },
        style: countStyle
 }).addTo(count).addTo(map);

How do I use the same concept in my current React scenario?


texts: What are common ways of presenting features that share coordinates?

texts:
When I encounter features of the same class, but with identical coordinates I lay them out in a spiral arm formation.  I chose to pick the initial radius and angle at random (can't recall my actual motive).  My question is this:

Is there a standard way (or ways) of symbolizing features, from any WFS service, which would otherwise stack on top of each other in a 2D representation? For example, how does one 'best' display a transformer on a pole and the pole itself?


texts: Can I view a Google Maps map by bounding latitude and longitude coordinates

texts:
With Google Maps I can enter a single lat/lng coordinate and the map is drawn with that lat/lng in the center.

What I'd like to be able to do is enter two coordinate pairs, say the southwest coordinate and the northeast coordinate, and have a Google Map displayed with the appropriate zoom so that the map fills the bounds defined by those coordinates.

Is this possible? If so, how?

Thanks

And the Answer is... Go to Google Maps Code Playground!

Per Taylor's suggestion, I went to the Google Maps Code Playground and entered the following code and was able to get a centered map with the appropriate zoom level bounded by the coordinate pair.

function initialize() {
  var mapDiv = document.getElementById('map-canvas');
  var map = new google.maps.Map(mapDiv, {
    mapTypeId: google.maps.MapTypeId.ROADMAP
  });

  var southWest = new google.maps.LatLng(..., ...);
  var northEast = new google.maps.LatLng(..., ...);
  var bounds = new google.maps.LatLngBounds(southWest,northEast);
  map.fitBounds(bounds);
}


texts: Previous shapefiles do not show in re-opened projects

texts:
I have previously drawn a lovely shapefile in a project - all saved...good job. I open the very same project and, although the background raster maps project fine, my shapefile will not show. Zooming to layer gets me to the origin point of the system - co-ordinates 0,0! I check the CRS and they are different for the shapefile vs the raster map. I change them to read the same. No change. I make sure "on the fly transformations" is checked - no change. I highlight the raster map and click "Set project CRS from layer" and I get a distorted background map and still no shapefile.

Help! Why is it so difficult? Is it possible to turn off CRS once it is set? I'm never going to change it again once it is working properly, and all this hassle will disappear!

Doug


texts: How to copy a line from polygon to an other polygon QGIS 2.6.1

texts:
I have a polygon(state) now the state divide into two.
I have a new polygon with two state.

I want to copy only the new border(line) which divides the two state from new polygon to old.

I am trying to select a line from a polygon but it is not possible.


texts: Add image to bottom corner of MapBox.js map

texts:
I have a MapBox.js map up and running and I would like to add an image to the bottom left corner.  I tried adding in this code:
<img src="IMAGE URL" height=100 style=’position:fixed;bottom:0px;left:0px;z-index:999’/> with my image hosted on Dropbox but it did not appear.  I am using several different TileMill produced layers for this map.  How can I add an image that will always be at the bottom left corner of my map?


texts: Reprojected output data has different shape than input

texts:
I am using the python gdal library and the warp method (osgeo.gdal.Warp) and I have it set up like the below.

warped = gdal.Warp(
    destNameOrDestDS=destNameOrDestDS,
    srcDSOrSrcDSTab=srcDSOrSrcDSTab,
    srcSRS=srcSRS,
    dstSRS=dstSRS,
    dstNodata=float(srcNodata),
    format=fmt,
)

What I have noticed is that my output (x,y) shapes are different than my input.  Then when reverse the operation they change again, both times increasing.  Is this the expected behavior?  Am I missing some option that would keep these values stable?  I thought maybe I need to input a cell size or something, but I do not see that in the options.  I am a programmer, not GIS expert so maybe this is a lack of understanding on my part.

EDIT:

As a follow up to the below comment and answer, why can I transform a single point as in the python gdal/ogr cookbook?

from osgeo import ogr
from osgeo import osr

source = osr.SpatialReference()
source.ImportFromEPSG(2927)

target = osr.SpatialReference()
target.ImportFromEPSG(4326)

transform = osr.CoordinateTransformation(source, target)

point = ogr.CreateGeometryFromWkt("POINT (1120351.57 741921.42)")
point.Transform(transform)

print point.ExportToWkt()

What I had intended to do once I saw this behavior was transform the points individually instead of using the warp method.  What would have occurred in this instance?  I feel like warp is the better way to go I just want to be able to explain it.

Data Source:

https://www.nwrfc.noaa.gov/weather/netcdf/2020/20200330/QPE.2020033012.nc.gz
srcSRS = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'
dstSRS = "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 datum=NAD83 +towgs84=1,1,-1,0,0,0,0 +units=m"
dstNodata=-999.0

The above are the parameters that I use, if applicable.


texts: Leaflet: from postgre via PHP to geojson to website

texts: I'm still looking for a best way to lunch pgrouting extension on my website as a rookie .Can anybody tell me is it possible? I have OSM data into postrgesql and want to do some routing machine. I want to run it in Leaflet with PHP. And my question is...can I run it as a geojson file? I heard that I can turn php file into geojson which can be lunched in Leaflet as a layer. Maybe some example? Tutorials? Someone need to head up me on the right way, becouse I'm confused after milion options that I have found in last few days.


texts: Running models in python script from a list

texts:
For a while now we're running a script each night via task scheduler from which several custom models are being run.
In an effort to streamline the code and simplify adding more models in the future I'm attempting to create a list at the start of my script and calling them in a for loop.

import arcpy

Toolbox= r"S:\GEOINFO\05_Tekenkamer\Marc\Test1"
Models = ["TestDissolve","TestTerminals"]

#Import Toolbox
arcpy.ImportToolbox(Toolbox)

#Run models
for Model in Models:
    print ("Executing model " + str(Model))
    try:
       arcpy.TestDissolve_Test1()
       print ("Executing model " + str(Model) + "succeeded")
    except:
       print ("An error occured")

print("End of script")

The part where I'm stuck is the first line of the try block where I run the model. I'm trying to substitute the fixed part of 'TestDissolve' after arcpy. for the current Model in the list (Test1 is the alias of the toolbox).
The first thing I tried was phrasing it like 'arcpy.Model_Test1' but that gives an error that the model Model can't be found. I don't think you can compose it in a string and then run it, or at least I don't know how.

Another way to fix the issue of running multiple models could be just to run every model from the toolbox somehow but not all of them have to be run every day and some have to be run in a certain order so I didn't try that.

I'm not sure if I'm overlooking something obvious here or if it's even possible but I'm hoping anyone here can help me with this.

Some additional information: we're using ArcGis 10.3.1 (python 2.7) but we might be able to use python 3 for this on our server if that could solve the issue.

How can I call models with arcpy from a list of models?


texts: How to use a variable inside .expression() as conditional statements criteria in google earth engine?

texts:
I am trying to automate a segmentation process in google earth engine.

I have some parameters calculated outside of the .expression() function and needed to be used in defining different criteria inside the .expression() to segment images. An example of my code is as follow:

var meanV = normalizedTh.reduceRegion({
  reducer: ee.Reducer.mean(),
  geometry: boundry,
  scale: 30,
  maxPixels: 1e9
});

var stdV = normalizedTh.reduceRegion({
  reducer: ee.Reducer.stdDev(),
  geometry: boundry,
  scale: 30,
  maxPixels: 1e9
});

var mean =ee.Number(meanV.get('b7'));

var std = ee.Number(stdV.get('b7'));

var hT = mean.add(std);

var lT = mean.subtract(std);

print(mean,std,hT,lT)

//----Segmenting the image-----

var segmentedTh = normalizedTh.expression(
      "(b('b7') > hT ) ? 3" +
      ": (b('b7')  > mean ) ? 2" +
      ": (b('b7') < lT ) ? 1" +
      ":  0 "
          );

//----Visualizing the results----

var thpalette = ['FFA500','red', '1E90FF    ','blue'];
Map.addLayer(segmentedTh.clip(boundry), {min: 0, max: 3, palette: thpalette}, 'Segmented th');

I am getting this error:
segmented th: Layer error: Image.parseExpression: Expression contains undeclared variables: [hT, lT, mean].

Apparently, I need to convert the calculated parameter outside the expression() to some format which could be used in conditional statements inside it.
Would you please help me to resolve the problem?


texts: Constructing file path within Python for Feature Class within Feature Dataset?

texts:
Based on my understanding of the ArcGIS path structure, a feature class inside a feature dataset inside a geodatabase can be accessed by the following:

C:\Menu\Lunch.gdb\Sandwich\Cheese

Where Sandwich is the feature dataset and Cheese is the feature class.

In a python script, I want to construct this full path, so that I can check it using arcpy.Exists(x). I have the following variables (pseudocode):

WSpace = C:\Menu\Lunch.gdb\
FDSName = Sandwich
FCName = Cheese

Concatenating WSpace and FDSName gives the correct path to the feature dataset. What I need help with is tacking FCName onto the end. I would think the following code would work, but it renders a space instead of the backslash character:

fullPath = WSpace + FDSName + r"\" + FCName

Instead of r"\", I have also tried "\\", and "/", but they all render a space. This seems like a simple syntax problem.


texts: I have slope raster and that is categorized with slope in percentage(%) but my script not successful

texts:
I have slope raster and that is categorized with slope in percentage(%)
my condition is

if
             slope <= 1     then m = 0.2
         1 < slope <= 3     then m = 0.3
         3 < slope <= 5     then m = 0.4
         5 < slope          then m = 0.5

Here "slope" is raster data and m is other variable that change as per slope percentage.

Is used raster calculator in ArcGIS 10.1 with script shown as below:

Con((("Slope" <= 1),0.2)  |  (("Slope" <= 3  &  "Slope"  > 1),0.3) |
    (("Slope" <= 5  &  "Slope"  > 3),0.4) | (("Slope" > 5),0.5))

But it is not successful.


texts: Currently, I am working on SQL in PostgreSQL and I execute query of vector data by "ST_Union" function for the spatial data. Now, I want to access the result of query in QGIS 2.14, but I can't. Is there any one who can help me?


texts: Location-Allocation "Like" Problem using Polygon data


texts: Turning All layers on for all data frames in ArcMap using ArcPy?

texts:
I am using ArcMap 10.2. I used the following code to turn all the layers on for the Main Map data Frame.

dFrame = arcpy.mapping.ListDataFrames(current_mxd,"Main Map")[0]
    for lyr in arcpy.mapping.ListLayers(current_mxd, "", dFrame):
        lyr.visible = True

Does anyone know how to use ArcPy to turn every layer on for all the data frames?


texts: How to describe a MODIS Tile's distortion

texts:
When looking at the MODIS grid [1] and compare it to how the sinusoidal projection is distorted, I wonder how to both practically and formally describe the distortion factor for each MODIS grid tile.

By practically I mean, what's a center tile's (e.g. row 8, col 18) width and height in meters on ground compared to a tile of a region that is rather distorted on the sinusoidal projected map, like row 2, col 14?

My understanding is that, by example on the grid shown at [1] which should be a 1km resolution MODIS grid, each tile should be approximately 1200km * 1200km, but is that approximation nearly the same for mentioned tile regions?

And by formally I mean, how to mathematically describe the distortion factor for MODIS tiles?

[1] http://modis-land.gsfc.nasa.gov/MODLAND_grid.html
[2] https://en.wikipedia.org/wiki/Sinusoidal_projection#/media/File:Tissot_indicatrix_world_map_sinusoidal_proj.svg


texts: Update Cell based on other cell Fieldcalculator

texts:
I want to search for certain values in a column. If that value is found, I want to update another a cell in the same row, other column.

I tried

CASE WHEN "column1" = 9 THEN "column2" = 'SA' ELSE "column2" END

That did not work. How would I approach that?


texts: Knowing how many images are applied to a reduction

texts:
Is there a way to know across any type of reduction in Google Earth Engine, how many observations or images are in the stack that is being reduced, on a per pixel basis?

Here is an example of a reduction, but this would be accessible for any type of reduction: percentile, mean, linear fit, max/min, etc.

//This function masks clouds and adds quality bands to Landsat 8 images.
var addQualityBand = function(image) {
  return maskClouds(image)
    //NDVI
    .addBands(image.normalizedDifference(['B4', 'B3']))
    //time in days
    //.addBands(image.metadata('system:time_start'));
};

var ALPS_collection_2010= ee.ImageCollection('LANDSAT/LT05/C01/T1_TOA')
  .filterDate('2010-1-1','2010-12-31')
  .map(addQualityBand);
// Create a greenest pixel composite.
var ALPS_greenestPixelComposite_2010 = ALPS_collection_2010.qualityMosaic('nd');
// Display the results.
var vizParams = {bands: ['B3', 'B2', 'B1'], min: 0, max: 0.4};
var ALPS_2010_mask = ALPS_greenestPixelComposite_2010.clip(Alps);
Map.addLayer(ALPS_2010_mask, vizParams, 'Mask applied to 2010 composite');


texts: What do I need to create 3D model of real terrain?

texts:
I'd like to create 3D model of some real terrain (to be exact it's about some real roads and only as much as 80m to 100m of roadside - 40m to 50m of each side). As I don't have access to highly detailed local maps, my only source of topographic information of given terrain will be Internet.

Now my question - what would be the best approach to create such 3D model and make it as detailed as possible? I know I won't be able to get very detailed elevation data (AFAIK, DEM for Europe has only 90m resolution) but on the other hand it's still better than nothing. But also I don't know what would be the best source of detailed topographic data? Should I use Google Maps and work on its KML files or OpenStreetMap gives better results? Or maybe there are some other resources for getting the information I will need? BTW, would be some in-car GPS system of any help in this case?

I'd like to add that I want to model only the surface - the vegetation, buildings and other objects are not needed. Also I won't need any geographic co-ordinates of terrain as I need to have only right measurements in meters.

I'm totally new to this subject so any kind of help will be very much appreciated.

Thank you in advance for reply.

EDIT:
Thank you everyone for your comments. I'll try to add more information on this project.

Source of data - AFAIK, the most detailed DEMs for Europe (Western and South Europe and also North Africa are my areas of interest) can be obtained from ASTER and Landsat. Both have 30m resolution so it's as accurate as I need it to be.

3D model of terrain - I'd like to reconstruct some real roads that used to be parts of famous rallies or races in the 50's and 60's (Coupe des Alpes, Liege-Rome-Liege, Rallye del Sestriere or Tulpenrallye). Of course not whole roads as it would take ages for me, only some sections that can be about 7km long (that's game engine limitation).
So I need a 3D model of surface of some specific road so I could edit it in 3DS Max 5 program. Or some other 3D modeling program (FreeCAD or Pov-Ray) as I need 3DS Max 5 only for texturing, adding some objects like buildings, trees, people, sponsor banners and finally exporting max file to bmf file (it's a format that Rally Trophy uses). And the game looks like in the link below:

Rally Trophy - Game Intro

Now, my biggest concern is the right method of generating such 3D model of terrain - what kind of programs do I need? And also how to mark a road on it?

I'm sorry I cannot explain it any better. But please ask as many questions as you need to.


texts: Manage values from textbox widget in Google Earth Engine

texts:
I am a newbie in GEE and I am trying to build up an apparently simple widgets combination that allows to center the map around a point whose coordinates are entered by the users into text boxes.
I have tried the following code:

var lon = ui.Textbox({
  value: 1.0,
  placeholder: 'Enter longitude here...',
  onChange: function(value) {
    return(value);
  }
});
print(lon);

var lat = ui.Textbox({
   value: 1.0,
  placeholder: 'Enter latitude here...',
  onChange: function(value) {
    return(value);
  }
});
print(lat);

var Lo = ee.Number.parse(lon);
var La = ee.Number.parse(lat);
var button = ui.Button({
  label: 'Go to Location',
  onClick: function() {
    Map.setCenter(Lo,La,1);
  }
});
print(button);

which returns the following error:

Provided center object has invalid values for lat, lon, or zoom. Expected numeric values.

I am afraid that either I am not able to convert the text into number, or I am making a mistake in the structures defining "lon" and "lat"


texts: How to create points at line meeting points in ArcMap?

texts:
I have some countinous multilines(with 1 or more common points), and i would like to create points at the common points.
I tried with starting point and ending point creation, but it doesn't end or start at the wanted point.
Any idea?


texts: Defining number of cells in sample grid for IDW?

texts:
I'm trying to make IDW in R. In my tutorial, I found this approach which is easily understandable, expect of one point - why they used n = 6000 ?

How did they established the value to create a grid by spsample()??

library(maptools)
library(GISTools)
library(gstat)
require(deldir)
require(sp)

voronoipolygons = function(layer) {
  crds = layer@coords
  z = deldir(crds[,1], crds[,2])
  w = tile.list(z)
  polys = vector(mode='list', length=length(w))
  for (i in seq(along=polys)) {
    pcrds = cbind(w[[i]]$x, w[[i]]$y)
    pcrds = rbind(pcrds, pcrds[1,])
    polys[[i]] = Polygons(list(Polygon(pcrds)), ID=as.character(i))
  }
  SP = SpatialPolygons(polys)
  voronoi = SpatialPolygonsDataFrame(SP, data=data.frame(x=crds[,1],
                                                         y=crds[,2],
                                                         row.names=sapply(slot(SP, 'polygons'),
                                                                          function(x) slot(x, 'ID'))))
  return(voronoi)
}

library(gstat)
library(maptools)
data(fulmar)

fulmar.spdf <- SpatialPointsDataFrame(cbind(fulmar$x, fulmar$y), fulmar)

fulmar.spdf <- fulmar.spdf[fulmar.spdf$year == 1999, ]

fulmar.voro <-voronoipolygons(fulmar.spdf)

s.grid <- spsample(fulmar.voro, type = "regular", n = 6000)  ### WHY ?? !!!
idw.est <- gstat::idw(fulmar ~ 1, fulmar.spdf, newdata=s.grid, idp = 1.0)

# extract the unique xand y locations in the grid
ux<-unique(coordinates(idw.est)[,1])
uy<-unique(coordinates(idw.est)[,2])

# extract the predicted values and format var1.pred into a matrix of gridded values
predmat <- matrix(idw.est$var1.pred, length(ux), length(uy))

When I tried to rerun the code for my data, I found that the length(ux) and length(uy) should be multiplication of length(idw.est) to correctly fill my predmat matrix otherwise I have a warning msg :

Warning message:
In matrix(idw.est$var1.pred, length(ux), length(uy)) :
  data length [1003] is not a sub-multiple or multiple of the number of rows [49]

and my resulting interpolation is not correct.

Can you understand and explain me why did they used n = 6000?

Was it just by trial and error, or what trick can I use to set this number for any area?


texts: Importing raster in the same resolution/extent in SAGA

texts:
I'm using the Imcorr-Feature tracking in SAGA GIS for comparing displacements on  orthoimages and DSM. Every raster has to have exactly the same number of cells and map extent. This is where i'm struggling.

I used the tool "Grid Masking" and "Clip grid with polygon" in order to achieve this but the results were never exactly the same in terms of extent and number of cells. The results were very similar but not the same so I'm unable to run Imcorr...

I've also tried to create a new grid system, larger than my rasters in order to try to import it all it this new grid system. But can't import the raster in this new grid system.

Has anyone got any advice for having exactly the same map extent and number of cells in all my rasters?


texts: Split a large polygon with a layer of polylines or with a layer of polygons

texts:
I have a large polygon, of about a size of a municipality and i want to split it in about 12 - 13 pieces by using a related layer of either polylines or polygons (buffered polylines).

I tried using the Digitizing tools plugin but without success. I merge all the polylines or all the buffered polylines in order to have one feature that will be used to cut my municipality-polygon but it just does not work, without even have an error message. It just loads for a sec and then the municipality - polygon in still the same.

I also tried the erase function (difference) but i do now know if there is a way to split it after it is clipped.


texts:
The following expression enables snapping for PyQGIS 2.0 :

QgsProject.instance().setSnapSettingsForLayer(layer.id(), True, 2, 1, 10, True)

What is the equivalent expression in PyQGIS 3.0?


texts: ArcGIS Field Calculator returning error 000989?

texts:
I am trying to calculate a percentage of a field based off the value of another field.
if damage = x then return  the value of field !mrkt_imprv! * y

y = percentage value

my code is below its returning an error in the calculate field

def dam_cal():
    if !DAMAGE! == 'Affected':
        return !IMPRV_MRKT! * .1
    elif !DAMAGE! == 'Minor':
        return !IMPRV_MRKT! * .2
    elif !DAMAGE! == 'Major':
        return !IMPRV_MRKT! * .65
    elif !DAMAGE! == 'Destroyed':
        return !IMPRV_MARKT! * .85

the error is #000989 (line 2)
Parameters are not valid.


texts: Exporting a GeoTIFF from Google Earth Engine to Google Drive

texts:
I am trying to export my EVI analysis from GEE as a GeoTIFF and bring it into ArcGIS Pro. I have looked at many tutorials and have no idea why this won't work...I am not getting any error messages, it just will NOT export to my Google Drive from GEE... Here is my code. Am I doing something wrong?
var collectionID = 'LANDSAT/LC08/C01/T1_TOA'
var startDate = '2020-05-01'
var endDate = '2020-05-30'

var data = ee.ImageCollection(collectionID)
          .filterBounds(area)
          .filterDate(startDate, endDate)
          .filterMetadata('CLOUD_COVER', 'less_than', 10)
          .median();
print(data)

var FCC ={ bands: ['B5','B4','B3'],min:0, max: 0.3};
Map.addLayer(data,FCC,'Landsat8');

var evi = data.expression(
  '2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))', {
    'NIR' : data.select('B5'),
    'RED' : data.select('B4'),
    'BLUE' : data.select('B2')
});
print(evi)

Map.addLayer(evi, {min: -1, max: 1, palette: ['red', 'brown', 'yellow', 'green']}, 'EVI');

// Export the image, specifying scale and region.
Export.image.toDrive({
  image: collectionID,
  description: 'Y2020M05',
  scale: 30,
  region: area
});


texts: Fixing 'qodbc3:unable to execute statement' when querying SQL Server 2008?

texts:
While testing features of QGIS I tried to query each table using right click and query but I can't query sql server 2008 database :

Error:
An error occured when executing the query.
The data provider said:
[microsoft][odbc sql server driver][sql server] invalid object name 'table_test'
qodbc3:unable to execute statement

but it works fine with postgis.

Also, is there any why to join table and do the query?


texts: How to use Cross Origin Resource Sharing (CORS) from ArcGIS API Javascript Application?

texts:
I am getting the CORS error

  XMLHttpRequest cannot load http://services.arcgisonline.com/ArcGIS/rest/info?f=json.
  No 'Access-Control-Allow-Origin' header is present on the requested resource.
  Origin 'http://myserver.com' is therefore not allowed access.

when I use the ESRI  topo base map in my ArcGIS JSAPI application.
The application is hosted on IIS8 (Windows Server 2012 R2)

I am reading the ESRI ArcGIS API for JavaScript Samples CORS buffer
There is a suggestion to add the root URL for CORS enabled servers to

esri.config.defaults.io.corsEnabledServers

Can you provide example how to do that?

I am not sure if the problem is on my server/application or the ESRI map server


texts: Points and polygon with same coordinate system do not line up in ArcMap 10

texts:
1st, I added a shapefile with projected coordinate system SL_grid_99 to Table of contents.

2nd, I added my points (excel file) to the table of contents. I displayed XY data and exported point data as shapefile. This is the coordinate system for the points:

It seems that both, point and polygon, have the same coordinate system. However, they are not lining up. As can be seen in the map, the points are in the lower left-hand corner:

How can I make the points and polygon line-up?


texts:
I have a detailed map in qgis (vector, .shp) and I would like to split the map in four parts based on attributes to print each part of the map separately. Unfortunately, I didn't figure out how to do that, yet.

I can select the desired values (for one printout) in the attributes table and the region of the map gets highlighted. I could now save the selection as a new layer, reopen it, and print it. however, that doesn't seem to be the smart way, because i would end up with too many files.

I thought i could simply print only the selection. Is there a way to achieve this?

Thanks


texts: Finding Station order number along way using PostGIS

texts:
In QGIS and PostGIS I have two layers:

way (id, name, geom (multistring))
station (id, name, geom (polygon))

How do I find the station order along each way?

I tried using code of Finding all intersections of LineString and Polygon and order in which it intersects it using PostGIS
SELECT ROW_NUMBER() OVER() as id, station.nome as station, dp.path, dp.way_name, station.geom FROM s_602_ptrc.station, (
SELECT (ST_DumpPoints(geom)).path[2] as path,
       (ST_DumpPoints(geom)).geom as geom,
       way.name as way_name from s_602_ptrc.way
) as dp
WHERE ST_INTERSECTS(dp.geom, station.geom)
ORDER BY dp.way_name, dp.path;

with this result.

Only the stations that are nodes of the polyline are considered - But how to consider them all?
UPDATE
the first part of the code of @Geozelot goes all null values for the way_part:


texts: Multiple HDF5 files in single GDAL virtual raster

texts:
Does anyone know how to generate a single GDAL virtual raster linking multiple HDF5 files to each other? Each file represents different swaths of the globe with different lat/lon coordinates and data bands.

For now, I'm generating a single .vrt file for each HDF5 file, then using gdal_merge to collect them all together. However, I would like to do this in one process. The idea being to link all the files to generate a single global map.


texts: Loading of layer fails when using SPYDER

texts:
I configured SPYDER to work with QGIS by setting start file
and setting additionally some system paths

@echo off
SET OSGEO4W_ROOT=C:\PROGRA~1\QGIS2~1.18
SET QGISNAME=qgis
SET GRASSNAME=grass\grass-7.2.1
SET QGIS=%OSGEO4W_ROOT%\apps\%QGISNAME%
SET QGIS_PREFIX_PATH=%QGIS%
SET GRASS=%OSGEO4W_ROOT%\apps\%GRASSNAME%
SET SPYDER="C:\\PROGRA~1\QGIS 2.18\apps\Python27\Scripts\spyder3.exe"

CALL "%OSGEO4W_ROOT%"\bin\o4w_env.bat
call "%GRASS%"\etc\env.bat
path %PATH%;%QGIS%\bin
path %PATH%;%GRASS%\lib

set PYTHONPATH=%QGIS%\python;%PYTHONPATH%
set PYTHONPATH=%OSGEO4W_ROOT%\apps\Python27\Lib\site-packages;%PYTHONPATH%
set PYTHONPATH=%OSGEO4W_ROOT%\apps\qgis\python\plugins\processing;%PYTHONPATH%
set PYTHONPATH=%OSGEO4W_ROOT%\apps\qgis\python\plugins;%PYTHONPATH%

start "SPYDER aware of QGIS" /B %SPYDER% %*

import sys
print sys.path

sys.path.append('C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python')

sys.path.append('C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins')

sys.path.append('C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\\processing')

after those settings all imports are OK. However when i try to load a layer it fails, but the same code works in QGIS Desktop python environment. My goal is to use SPYDER (or another IDE e.g. Eclipse) to make life easier when debugging/designing QGIS scripts. Is it possible to integrate IDE with QGIS to acheve this ??

here is my code

from osgeo import ogr
from PyQt4.QtCore import QVariant
from qgis.core import QGis, QgsFeature, QgsField, QgsFields, QgsVectorLayer, QgsMapLayerRegistry
from processing.tools.vector import VectorWriter
import csv

from PyQt4.QtGui import *

wb = QgsVectorLayer('C:\\Users\\krzysztof.fajst\\Documents\\PYTHON\\GIS_LV networks\\os with roads feeders and substations.shp', 'os with roads feeders and substations', 'ogr')
layer = QgsMapLayerRegistry.instance().addMapLayer(wb)

also simple

layer = qgis.utils.iface.activeLayer()

fails saying qgis name not defined even if I imported iface with

from qgis.utils import iface

Krzysztof


texts: QGIS map scales and zoom levels of XYZ tiles

texts: I have a street name I want to see on a rendered tile at zoom level 15, but not level 14. How do I figure out what my minimum/maxixum scales for the rule-based labeling need to be?


texts: add additional data to geojson with leaflet

texts:
I'm using leaflet to create a world map with geojson. Now I need to add some additional information as popup for each country but I'm not sure how should I proceed it. Been finding solution for hours.

Should I modify geojson file to add 1 more field since it already listed coordinates of all the countries. Can I achieve this by php or js?
Should I just store the information in SQL ("US" -> "US's text", "Canada" -> "Canada's text,...) & use switch case to compare string in onEachFeature function?


texts: Using GPRecordSet as output parameter in Python Toolbox?

texts:
I would like to use a GPRecordSet in my python toolBox as an output parameter.

This geoprocessing will be share on my ArcgGIS server and in my ArcGIS portal.
I haven't see any documentation on how to use it.

Can someone give me a simple example to know how to use it?


texts: Downloading multiple datasets using the same reduceRegion in Earth Engine

texts:
I'm trying to get the average value within a week timeframe, for every band, for every dataset of Sentinel 5P on Google Earth Engine, for a 10km radius around a point.

The part I'm having trouble with is efficiently joining (currently with an inner join?) the various Sentinel 5P datasets together.

This is currently resulting in an empty dictionary as the result of .getInfo(), where as if I use reduceRegion on a single image (i.e. the mean of a single dataset filtered by date) I get a reasonable response. I've done this before by looping through the databases, but I'm assuming there is a much more efficient way.

How can I efficiently apply a reduceRegion across multiple earth engine datasets at the same time?

start_date = '01-01-2019'
next_date = '01-08-2019'
geometry = ee.Geometry.Point(lon=lon, lat=lat).buffer(1e5) #10 km around an arbitrary point

EARTHENGINE_DBS = [
    "COPERNICUS/S5P/OFFL/L3_NO2",
    "COPERNICUS/S5P/OFFL/L3_SO2",
    "COPERNICUS/S5P/OFFL/L3_CO",
    "COPERNICUS/S5P/OFFL/L3_HCHO",
    "COPERNICUS/S5P/OFFL/L3_O3",
]

# -- Make a date filter to get images in this date range. --
dateFilter = ee.Filter.date(start_date, next_date)

# --- Select imaged collection ---
ic = ee.ImageCollection(EARTHENGINE_DBS[0])
for db in EARTHENGINE_DBS[1:]:

    # --- Define the new data ---
    _new = ee.ImageCollection(db).filter(dateFilter)

    # --- Define an inner join. ---
    innerJoin = ee.Join.inner()

    # --- Specify an equals filter for image timestamps. ---
    filterTimeEq = ee.Filter.equals(** {
            'leftField': 'system:time_start',
            'rightField': 'system:time_start'
    })

    # --- Apply the join. ---
    ic = innerJoin.apply(ic, _new, filterTimeEq)

    # --- flatten ---
    ic = ic.map(lambda feature: ee.Image.cat(feature.get('primary'), feature.get('secondary')))

# --- aggregate ic ---
agg_ic = ic.mean()

mean_dict = agg_ic.reduceRegion(**{
    'reducer': ee.Reducer.mean(),
    'geometry':geometry,
    'scale': 1000,
    'bestEffort':True,
    'tileScale':16
})

mean_dict.getInfo()


texts: Kriging with uncertain observational data

texts:
I'm currently trying to use ArcMap 10.5.1 for this analysis:

Thus far, I have taken particulate matter trend estimates at monitoring sites across the U.S. and used Simple (Normal Score - due to slight non-normal distribution) Kriging to interpolate trends in between these sites. While I only have 156 sites, the Krige standard error prediction map shows small standard errors across the region. This suggests to me that the interpolation itself is fairly accurate to the observed data.

However, as was pointed out to me recently, this analysis does not take into account the uncertainty of the trend estimates. In other words, I take a site-specific estimate and then smooth those estimates across locations without accounting for uncertainty at the locations. This seems likely to underestimate uncertainty in the resulting Krige prediction. What I'm looking for is a way to incorporate my uncertainty (I have standard error values for each location) into the trend Krige.

Suggestions I've already received include:

Kriging the standard error at each site as a separate figure to show uncertainty across the region. Ideally, I'd like to find some way to include the standard error into my trend estimate so I only have one Krige prediction map.
Subset my data and cross-validate multiple times to confirm Krige prediction output. This still seems like I am not considering the inherent trend estimate uncertainty though.

Is there any way to include measurement standard error (which is different at each site) in the Krige prediction to reflect uncertainty?