texts: Sentinel 5P Extraction through Google Earth Engine : scale

texts:
I've been trying to extract Sentinel 5P imagery through GEE and I am confused on which scale to extract using Export.image.toDrive function. Usually, we give 30m resolution for LANDSAT imagery as the resolution is 30m. But what should I specify the resolution for Sentinel 5P?
My project's aim is to propose green measures in a city that has highly concentrated pollution pockets.
Though I have downloaded GeoTIFF with 500m/pixel resolution and 5000m/pixel, I don't think they are right to perform a fruitful scientific analysis.
PS: I am new to Satellite imagery processing.


texts: GeoServing thousands of rasters, is there a limit on the number of images?

texts:
I have a web application on which the users do some processing that outputs rasters. Those rasters are transferred to geoserver for later visualization through WMS.

The GeoServer is running on TomCat6, this geoserver is feeded daily with hundreds of ~2mb rasters.

A few days ago it reached ~3000 datastores (and layers) and it stopped responding every time a new raster is added, then it needs to be rebooted.

I followed the optimization guideline on the geoserver documentation and changed to Oracle Java and set the container options (more java memory and throughput garbage collector).

That solved the problem and it's running very fast.

But the number of rasters continues to increase and my question is if there is a limit on the number of raster stores that geoserver can handle so that I should expect to have the same problem in the future.


texts: Tools to compare unversioned copies of geospatial data

texts:
My office has about 100 GB of geospatial data of varying forms (imagery, feature classes in geodatabases, LIDAR, tables, etc) that have gone somewhat unmanaged for the decade during which they were collected. I have now been put in charge of inventorying and restructuring the data. The data is stored on the filesystem and is unversioned, but many datasets and geodatabases have been copied to new places without additional metadata describing what or why.

Without going into what we should have done or how to structure it in the future (I'm working on it), I'm trying to figure out if there is a decent way to compare the feature class data and geodatabases when I encounter a copy to determine what changed between the two copies. I can use Beyond Compare to compare many of the types of data.  I came across a few GIS.SE questions on the topic of comparing specific items, but I don't know that they would work - at least not without further scripting, which isn't desirable but is an option.

Getting differences between old and new versions of ArcSDE database?
How can one compare two geometries in arcpy successfully? (this one could work with finessing, but I'm hoping for something better)

I realize that the comparison is a pretty difficult operation given all of the different types and properties of GIS data, but I'm hoping for something that can at least handle differences in spatial location, associated data tables, and ideally metadata. Commercial software, on the one end, and relatively barebones code on the other, are both ok for me.

Does this even exist?

Edit: Just to be clear - we're using ArcGIS 10.0 and 10.1 - I'm happy to use solutions that involve other open source or commercial providers too, but any ESRI solutions would need to be compatible in 10.0. Also, our data are not currently in a true database (SDE, postGIS, etc), but I'm willing to move the stuff I want to compare into another format if it means I can compare it.


texts: How to do a spatial search without select() using PyQGIS?

texts:
Using QGIS 1.9.0-master, I want to do a spatial search for features in a vector layer without selecting them, i.e., without using QgsVectorLayer.select(), QgsVectorLayer.selectedFeatures() et al.

More specifically:

I want to get the feature IDs inside a rectangle. That could be accomplished with something similar to QgsSpatialIndex.intersects() if I could gain access to the spatial index of the layer's dataprovider.

Unfortunately, I didn't find a way to do it. I was only able to find how to create the index using QgsVectorDataProvider::createSpatialIndex()

Does anyone know if this is possible?


texts: Changing aspect ratio of legend in ArcMap?

texts:
I'm making some good ol' maps at work (Using 10.x) and I'm having difficulty getting my legend to fit with the amount of descriptive text I have and the amount of space available to me (as well as to make it look... you know... decent. Which can be difficult a lot of the times with how clunky layout view in ArcMap is).

Is there even a way to theoretically get ArcMap to un-gray-out the checkbox that has "Preserve Aspect Ratio" checked? Nothing I have tried has worked and my google - fu ended with someone who had a similar question and a response that was basically, "Sorry dude, ArcMap just be like that some times."

Ideally I would rather not convert the legend to graphics since you can't make changes easily to the legend after the fact, so is there a better way to do this without converting to graphics?

Here's a picture if it helps refresh your brain on where this is in the legend properties:


texts:
I have two rasters, one depicts current water depth and the second depicts the original pond depth.  I wish to calculate the volume of the region between the two.
-be gentle, I'm a noob.


texts: Does CartoDB support the SQL MERGE command?

texts:
I am developing a mechanism that uploads records to CartoDB - first by importing the records to a new table, and then insert the records from the new table to an existing table.

The problem is that the existing table may already contain some of the records (each record has a unique record_id) - in this case I would like to override the existing values with the values of the new record.

This could be done in two steps, by deleting the existing records (whose record_id is in the new table) and then insert all the records from the new table.

The SQL MERGE command should allow to do that in a single command. Is this command supported by CartoDB? Can you give an example with the exact the syntax? (I am getting a 'Syntax error at or near "MERGE"').


texts:
I have to amend the following code so that on the map it's possible to draw a second circle graphic as well as a line, polygon and point. I've tried to add an else if statement much like the polygon type but that hasn't worked.

Can anyone advise?

function addSearchGraphic(geometry) {

    drawToolbar.deactivate();
    map.enableMapNavigation();
    var type = geometry.type; //can only ever be a point or polygon

    if (type === "point") {
        $('#point').button("option", "label", "Draw Another Search Point");
        $('#polygon').button("disable");
        $('#line').button("disable");
        $('#circle').button("disable");
        symbol = drawToolbar.markerSymbol;
        searchgraphicsLayer.add(new esri.Graphic(geometry, symbol));
    } else if (type === "polygon") {
        $('#point').button("disable");
        $('#line').button("disable");
        symbol = drawToolbar.fillSymbol;
        if (esri.geometry.polygonSelfIntersecting(geometry)) {
            alert('Problem! The area you have drawn has a boundary that crosses over itself.  Please redraw your area');
        } else {
            $('#polygon').button("option", "label", "Draw Another Search Area");
            searchgraphicsLayer.add(new esri.Graphic(geometry, symbol));
        }
    } else if (type === "polyline") {
        $('#line').button("option", "label", "Draw Another Search Line");
        $('#point').button("disable");
        $('#polygon').button("disable");
        $('#circle').button("disable");
        symbol = drawToolbar.lineSymbol;
        searchgraphicsLayer.add(new esri.Graphic(geometry, symbol));
    }
    bufferSearchGeometries();
}


texts: PostGIS/pgrouting pgr_createTopology return an empty table

texts:
I am new of Pgrouting, I have created a new database with PostGIS and PgRouting extension, than I have populated my DB using the osm2po sotware and now I have my graphs. However I need to create the vertices table in order to make the network routable, I run the following command:

SELECT pgr_createTopology('grafo_padova', 0.00001, 'the_geom', 'id');

Postgres returns

NOTIFICA:  PROCESSING:
NOTIFICA:  pgr_createTopology('grafo_padova', 1e-005, 'the_geom', 'id', 'source', 'target', rows_where := 'true', clean := f)
NOTIFICA:  Performing checks, please wait .....
NOTIFICA:  Creating Topology, Please wait...
NOTIFICA:  -------------> TOPOLOGY CREATED FOR  0 edges
NOTIFICA:  Rows with NULL geometry or NULL id: 0
NOTIFICA:  Vertices table for table public.grafo_padova is: public.grafo_padova_vertices_pgr
NOTIFICA:  ----------------------------------------------

Successfully run. Total query runtime: 77 msec.
1 rows affected.

And the result is an empty table. The table "grafo_padova" has a geometry column named "the_geom" and it is filled.

Can someone help?


texts: Implementing topology rule Line Feature Class should not overlap with Polygon Feature Class in ArcGIS Desktop?

texts: Is there a way to setup such a rule in ArcMap 10 Topology that Line Feature Class features shouldn't overlap with Polygon Feature Class features?


texts: Arcpy ERROR 000732 when ApplySymbologyFromLayer

texts:
I want to split a shapefile by certain attribute and then apply a symbology layer, finally output it as a KMZ file.

R1 = outpath + 'R1.lyr'

arcpy.MakeFeatureLayer_management(fc, R1, "WHERE CAUSE")

arcpy.ApplySymbologyFromLayer_management(R1, symbologyLayer)

arcpy.LayerToKML_conversion(R1,outpath+'R1.kmz')

It raised error: arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Input Layer: Dataset R1 does not exist or is not supported
Failed to execute (ApplySymbologyFromLayer).


texts: Why can't I upload shapefiles into my SpatiaLite database

texts:
When I try to import a layer into my SpatiaLite Database, I get this error message.

Why is this happening now ?


texts: Can you programmatically change the button image for a ESRI.ArcGIS.Desktop.AddIns.Button?

texts: I'm new to the AddIn customization model in ArcGIS 10 and I'm trying to figure out how to change a button image for a button. I'm looking for the equivalent to the BaseCommand UpdateBitmap method and haven't had any luck.  Is this a limitation of the AddIn model?  Or did I miss something obvious?


texts:
I am looking for a rather large DEM --- on the trillion cell scale --- to run some tests on.

So far, I have found:

200 billion cells: The SRTM Global Dataset 30m DEM
120 billion cells: The National Elevation Dataset 10m DEM (higher resolutions have incomplete spatial coverage)
70 billion cells: The PAMAP Lidar 1m DEM

Are there any publically-available DEMs, which are not stored in an Arc Geodatabase format, which meet or exceed the 1 trillion cell threshold?


texts: Merging Personal Geodatabase (MDB) data into ArcSDE using FME?

texts:
Im a novice in FME and I have a MDB file that i want to write/insert to an SDE database using FME. I created a simple translation that generally reads the data from the MDB file and writes it to the SDE database, although the translation is successful and it prints out the number of written rows, when i checked the SDE database, it wasn't able to right any data on it.

i then checked the log and i found this specific warning:

  MS SQL Server (Spatial) Writer: Attribute data ' too large to fit in
  specified data type for columnSHAPE'

The shape column holds the geometry data in the MDB file. FME sets the SHAPE datatype to 'oleobject' in the Reader and I noticed that the SHAPE column in the writer was automatically converted to 'text' datatype.

Is there anything else i need to configure to make this work?

UPDATE: I am using FME Database Edition 2014 (20140120 - Build 14235 - WIN64)

Here are some screenshots of my FME translation:


texts: Tick boxes for a qgis plugin written in python

texts:
in a python plugin for qgis I want the user to be able to select one or more of several options via tick boxes.

I imagine it to look like e.g. Qgis 2.2 tool "Vector" -> "Data Management Tools" -> "Join attributes by location". In the "Attribute Summary" section of this tool the user is able to tick several functions on and off:

At the moment I use a string input (##options=string "option1,option2,option3"). The user has to enter all options separated by a comma. This works but it does not look very pretty.

Is this possible with Pyhton? Any suggestions? Thanks in advance.


texts: Divided road segments from XTools Pro need unique, informative IDs

texts:
I can divide segments with my road using XTools Pro, but I want to give a useful identifying number to each segment.

For example: If D100 road has been divided into 20 segments, the segments must be numbered 1,2,3,4,5,....20, based on the kilometer location of segments.

How can I do this?

I use ArcGIS 10.2.2. Let me try to explain with an example:

I divided road segment A that is 60 km long with XTools, and each of the resulting segments is 1 km. In total I have 60 subdivided segments. But XTools adds same ID for all subdivided road segments (all 60 subdivided parts of segments A have same ID).

It is hard to detect the kilometer of subdivided parts of segments. Considering that I have a lot of segments, it is impossible to know the location of any (randomly selected) subdivided segments.


texts:
I have 5 python scripts in the following folders :

.\Signal_Slot_Example\main.py
.\Signal_Slot_Example\ProgressDialog\init.py
.\Signal_Slot_Example\ProgressDialog\Progress.py
.\Signal_Slot_Example\SliderDialog\init.py
.\Signal_Slot_Example\SliderDialog\Slider.py

Main.py when run calls the other 2 scripts and creates a simple slid that interacts with a progress bar and when I run it in Idle it all works fine.

When I run it in the QGIS python console I get :

Traceback (most recent call last):
  File "C:\OSGEO4~1\apps\Python37\lib\code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "<string>", line 1, in <module>
  File "C:/OSGEO4~1/apps/qgis-ltr/./python\qgis\utils.py", line 672, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
ModuleNotFoundError: No module named 'SliderDialog'

Why won't it work in QGIS but runs fine in Idle ?

All the example code I have run in QGIS python console and in idle worked the same the only difference here is that I am calling classes from other .py scripts.


texts:
Is it possible to average a stationary position with a GNSS (GPS, GLONASS...) with a centimeter accuracy without external corrections (like differential or RTK)? If so, how long does such averaging have to take?

My problem is I have an RTK setup that is working quite nicely - in consists of two u-blox NEO-M8T evaluation boards and I was able to get an RTK fix (static mode) with RTKLIB very quickly (1-2 minutes). I have one problem though - I would like to gather precise location data in unknown area where I have no precise known location around that I could use for the base station so I am thinking of averaging the position for base station from a single receiver myself... so far I used the position from u-center for RTKLIB - I guess RTK can do good relative precision from rover to base but rover's absolute precision goes back to the base station precision I guess...


texts: How do I assign point value ranking system to multi-ring buffers?

texts:
I created a multi-ring buffer of different distances, and I want to assign a point ranking to each distance as part of my analysis. So for example: farms within the 0-25 mile buffer would get 4 points, farms within the 25-50 mile would get 3 points, etc. How do I do this?

I created a spatial join to get counts of the outputs, and now I want to integrate the values with the counts. Like, there were 4 farms within the 0-25 mile buffer, so the point ranking value would be 4 times 4.


texts: Putting result of selection of plots in frame of QGIS?

texts:
I have made a plugin for QGIS in pyQt. Then I can make a query from a PostGIS database and I get the plot or bunch or plots results of that query in a layer.
This layer is overlapped on the top of a layer who show all the plots.
But sometimes I get only one of two plots, and they "lose" between more than 5000 plots, then I would like to put the plot or bunch of plots in a frame.

Any idea how could I do that?


texts:
From a WFS (with EPSG: 6247 displayed in the window) I added several layers and changed the CRS to EPSG: 3116 (the same I'm working on my current project). The layers  are not overlapping  even if "on the fly transformation" is enabled. Assuming that is a projection problem I made a reprojection (save as ....) and change the coordinate system but the situation remains the same.

I tried different CRS on all layers but the situation does not change. Even I downloaded  the .prj again and replaced it in some layers without success.

The only difference I find is that the  "extent" changes in the metadata. What I understand is that this change with the reprojection, is not independent of change with respect to CRS.

I'm I missing something so obvious ?

I followed previous posts without success:
- Same crs, different spatial reference system units, shape files / QGIS
- Adjusting extents so that same-CRS shapefiles overlap in QGIS?
-Why do my Shapefiles not align correctly?

.


texts: Getting week number from date using Arcade with ArcGIS Online

texts: How can I get the week number (week position in year) from a date with ArcGIS Online’s Arcade?


texts: Why does the undefined "ghost classes" arise during supervised classification in ERDAS IMAGINE?

texts:
After a supervised classification in ERDAS IMAGINE I have some undefined "ghost classes" in the attribute table. Also have a look at this youtube video at about 7:12.

I'm wondering why they arise and how they affect the classification result and the accuracy assessment?

My first consideration was that they arise within the merging in the signature Editor.
Or is it to pixels, that could not be categorized during the classification process because the possibility of a membership to one ore another class is not high enough?For instance I used the maximum likelihood classifier which categorises the pixels depending on the possibility of being a member of a specific class. In the latter case should they be marked in some way? They would also have to be taken into account in the accuracy assessment because the pixels are still there but were not assigned to a class.


texts:
Can anybody tell how to completely clean up your desktop computer from QGIS files ( Windows Registry, user folders and program folders)?

It's for QGIS version 2.4 and 2.10.

'Ok' and 'cancel' buttons are missing from dialog boxes in QGIS project, and despite uninstalling QGIS the same error appear when installing new version, even 2.8.5 (Long term release).

Therefore a complete clean up of QGIS is needed.


texts:
I am trying to calculate my End_Streets field within two joined tables. I want to remove the redundant SECONDARYN value from the STREETS field.

" & ".join(!IntersectAndEnds2.Streets!.split(" & ").remove(!StreetSplit2.SECONDARYN!))

Here is a portion of my table:

Here is the error I am receiving when I run the expression in the field calculator.

Error typed out:

  ERROR 000539: Error running expression: " & ".join(u"ACORN CT & LARAN RD".split(" & ").remove(u"ACORN CT"))
  Traceback (most recent call last):
  File "", line1, in  TypeError Failed to execute (CalculateField)


texts: What prior art can be cited to refute Apple's claim to have invented Interactive Map?

texts:
In U.S. Patent Application 20130339891 Apple Inc. claims "An
interactive capability enables a user to dynamically adjust the
content of an electronic map" as well as "Geospatial applications can
be integrated with the map to provide information pertaining to the
area displayed on the map, and to refine search results to those that
are relevant to the area." Apple seeks to patent these ideas based on a number of claims.

The claims the patent application makes seem overly broad and applicable to many GIS capabilities; starting with "A method for displaying a map on a computing device, Comprising: storing information to be displayed on the map in a memory of the computing device, the stored information comprising a plurality of different layers of information, wherein each layer contains a respective type of information."

Surely there is prior art in many GIS systems that can be cited to refute this claim. Please do so.


texts: Benchmarking Compact Cache v1 using ArcPy

texts:
I'm trying to benchmark the read performance of different cache types. One of the cache types I am benchmarking is the ESRI Compact Cache v1.
If there is not a function in the ArcPy library which can read in a single tile from a compact cache, given x, y, and z coordinates, does Esri provide a lower level library for such functionality?


texts: R - not the same number of raster cells after reprojection

texts:
I have a raster projected in EPSG:21781 (available here) that I'd like to display with mapview. Raster values correspond to the grid id, and need to be preserved.
If I understood correctly, mapview is calling leaflet::addRasterImage() to do so, and reproject the raster into EPSG:3857 (pseudo-mercator) through raster::projectRaster().
The first problem is that the reprojection is done through a 'bilinear' method (by default) whereas I need the nearest-neighbor method. To overcome this problem, I have reprojected the raster manually using the following command:

grid <- raster::raster("~/grid100.tif")

grid2 <- projectRaster(grid,crs="+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs",method="ngb") # corresponding to epsg 3857

summary(values(grid))

Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
1 1035406 2070572 2074309 3111142 4157490 3554337

summary(values(grid2))

Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
1 1056410 2097896 2093243 3133844 4157490 3917853

table(is.na(match(values(grid),values(grid2))))

FALSE TRUE
6206341 1488539

Conclusion: 1'488'539 cells are dropped after reprojection...
Any idea of what is going on here?
Thanks


texts: Converting lat, long and values to raster doesn't project correctly


texts: Doing mosaic colour correction in ArcMap?

texts:
I have created a mosaic dataset of two landsat 8 raster tiles.  Since there is a non uniformity and a colour variation i am not able  to do supervised classification properly.I want to do a mosaic colour correction.  After opening the Mosaic colour correction window it is greyed out and I am not able to use the function.

Should I activate any extension to use it?
Is there any other way to solve this problem other than using mosaic color correction tool.


texts:
I'm writing a script to automate post and reconcile a target and edit version. I consistently get the warning message "The reconcile process was not performed. The specified target version has no edit versions to reconcile with." I double checked that in the Administer Geodatabase tree view that the edit version is a child of the target version.  I'm also the owner since I create those test versions.

I've tried googling the message and found someone asking the same question on Geonet in June 2015 with no answers.


texts: QGIS, I uploaded an SVG file as a label. I set the size at 7 however I have two labels appearing massive

texts:
Here is a photo below. I am using QGIS 2.18.14

How can I solve this issue and possibly have the labels spread out throughout the rest of the lines and not limited to the ones above?


texts: I have a lot of photos in the geodatabase, I want to export the photos to a file. How do I export the photos from geodatabse? What tools are used?


texts: How to implement subtype functionality in Oracle Spatial?

texts:
I am planning to shift my entire database from arcgis to oracle.
I will not be using ArcSDE services also.
I need subtitutes for ESRI's functionalities like domain, subtype, relationships in oracle.
I have some how worked out two them i.e. domains and relationships, but no alternatives for subtypes have striked me.


texts: Subtracting columns (that include null values) in ArcMap Field Calculator

texts:
I’m desperately trying to run a very simple command where I subtract one column with others. However the problem is that some rows have NULL values, therefore the result of subtracting is also NULL.
I've tried to use the script found on this forum, but I guess I didn’t write it correctly, as it is not working as it should .

def stack(item1,item2):
  itemList = [item1,item2]
  myList = [item for item in itemList if (item != None)]
  return (item1-item2)

Is there a standard function for subtracting two or multiple fields in ArcPython?


texts:
QGIS has a plugin called OpenLayers. From that plugin, we can get a very pretty map. I have downloaded some shapefiles for the region (a city in my case) I am targeting.

When I put the shapefiles on top of the OSM map, everything looks perfect (playing a little bit with transparency...)

Now I would like to cut the OSM map area not contained in my polygon shapefile while still keeping all the OSM attributes (names, roads, colors etc...) and keeping the zoom-able ability of the OSM map.
THEN I'd like to export the result as a SLD. I'll use that into GeoServer.

How would I go and do that or is there another way for doing that?


texts: How to extract polygon location data (County) to a point shapefile (only has DMS)

texts: I have imported an excel spreadsheet to my map and the points have fallen into certain locations. How do I add the location of the polygon that the point has landed in, to the attributes table of the point without hand jamming it all in one point at a time? For instance a table full of school locations has no county data in the table. How would I add the county to the schools attribute table?


texts: Using field calculator unable to keep time and date format

texts:
I was looking to keep the date/time as is when using field calculator to keep the format as I'm not able to maintain the time in the new field "da". This is what I get when trying to do so. I tried searching for solutions, but to no avail so far.

Update based on comments:

I'm using field calculator to solve this.
I'm calculating into a date field from string field.
Data source is a shapefile


texts: What is layer file (*.lyr) in ArcGIS Desktop?

texts:
My organization uses *.lyr files consistently as a means for our users to view the data and build their mxds.

It works great for our purposes, but we have run into an issue where we don't think the users actually understand what a .lyr file does or why its used, and we spend a lot of time dancing around an issue because of miscommunication in definition/terms.

How do you explain what a .lyr file is? Assume the user has light/moderate GIS knowledge?


texts: Data loss and misplaced values after Append

texts:
I attempted to append 10 shapefiles after carefully verifying the schema. I keep losing large amount of data and getting misplaced values (in wrong columns/fields).

I am using ArcMap 10.2.2 and total appended dataset size is 2mill plus observations (parcel polygons and attribute data).

Any ideas of what could be happening? I double/tripled check the schema in all 10 files.


texts: Object based image classification using ArcPy or Python

texts:
I am looking to classify several Landsat 5 TM scenes using object based image analysis within ArcMap either using a built in tool, writing script in ArcPy, or writing a separate Python code (this is less ideal since I am not sure how to run this code in ArcMap). I have been running pixel level classifications on these scenes but my images are not as accurate as they should be. I am fairly new to writing Python script but feel more comfortable using this language than eCognition Developer (which has poor documentation regarding utilizing the program as a beginner and running this type of analysis).

Do tools exist within ArcMap to separate the scene into objects?

Can ArcPy code be written to segment the scenes?

If so, are there any effective resources out there regarding this?


texts: Make a line from a set of GPS points close to a specific GPS point with PostGIS

texts:
I have to develop a map-matching algorithm. At the moment I only match GPS points with the closest road segment. I'm looking for a way to compare the directions between a GPS point and a road to improve the matching. I have already looked at many subjects and I found a way to compare the directions but I've an issue for getting the direction of each GPS point.

I use PostGIS to store and manipulate my data, and I want to compute the azimuth of each GPS point.  My problem is to have the set of points which can be used to build a line and calculate the azimuth of this line with:

ST_Azimuth(ST_StartPoint(lineBuilt), ST_EndPoint(lineBuilt))

My idea was to take all GPS points inside an area around the current GPS point to build a line and finally give the direction of this line to the current GPS point. My problem is to get the actual GPS point, currentPoint in the query below:

SELECT degrees(ST_Azimuth(ST_StartPoint(new_line),ST_EndPoint(new_line)))
FROM
(
SELECT DISTINCT ON(gps.gid) gps.gid, ST_MakeLine(gps.geom ORDER BY gps.gid) AS new_line FROM gps WHERE ST_DWITHIN( currentPoint?, gps.geom, 50)
) AS query;


texts: How to style feature in mapserver mapfile based on values from database?

texts:
We have a table "Shapes" in PostGIS that consists of (among other things) a geometry column (the_geom) and an ID column (type_id) which is a foreign key to another table called "ShapeTypes".

A given row in the Shapes table may have a Type ID value or it may be NULL. If it does have a Type ID, then we want to style the feature based on that. The issue is that we don't know how many different values there may be for Type ID, so we can't just add a fixed number of CLASS entries into the mapfile. Our goal for these shapes (polygons) is to set the fill color... and we can estimate that there will be between 3 and 10 – a color range would be fine. We can also do a SQL JOIN to get the name for each Type, so our query results are essentially the geometry, the type ID, and the corresponding type name. The name is important to the end-user who don't see the IDs elsewhere in the application.

Is it possible to color the shapes based on that Type ID and have the legend display the type name - even though we don't know in advance how many types there will be?


texts: Matplotlib error when installing Molusce or lifemapper in QGIS Brighton

texts:
I am installing plugins via the QGIS desktop internal interface, and have had twice now messages associated with matplotlib.  The first error was for something called lifemapper and now for MOLUSCE.  This second plugin is quite popular, so I am not sure if something failed to install properly when I installed the Brighton QGIS software today, or if the plugin is not configured to work with Brighton.

The specific error I get with MOLUSCE is:

  No module named matplotlib.backends.backend_qt4agg.


texts:
I'm trying to save a shapefile that came in EPGS:3005 coordinate system, so all the shapefile polygon points have 6 or 7 digits before the decimal point (e.g. 1223282.8231), and I want to save/convert the shapefile to EPGS:4326 or any other CRS that has 2 or 3 digits before the decimal point (i.e. longitude and latitude to be processed and uploaded to Google My Maps.

I've tried all the methods shown here, and it appears like the layer properties have been saved as EPGS:4326, but when I run my python script (that creates the KML file for Google Maps) all the coordinates are still in the 6 or 7 digit format. How do I fix this?


texts: SQL Query Syntax gives ERROR 000358: Invalid expression?

texts:
I'm writing a script that will cycle through a shapefile of building footprints for an entire state, create a folder for the individual counties in that state, and then intersect the building footprints file with a county polygon stored.

Because I'm working with multiple states and hundreds of counties, I've created two variables that will take the state and county names of the current iteration- I keep running into an issue when it comes to writing a query that uses those variables.  My goal is to create the following query:

"ST_NAME = 'Colorado' AND "POLYGON_NM" = 'Weld'

I'm not familiar with writing sql queries in arcpy and the syntax is throwing me.  Here's what I currently have as my query:

dictionary = set(zip(statenames, countynames))
for item in dictionary:
    try:
        stname = str(item[0])
        coname = str(item[1].replace(r' ', '_'))
    except:
        stname = str(item[0])
        coname = str(item[1])
    if stname == 'Colorado':
        coco = arcpy.CreateFolder_management((basefile + '/Colorado/Counties'), coname)
        query = '"ST_NAME" =' + "'" + stname + "'" + " AND" + '"POLYGON_NM" =' + "'" + coname + "'"
        fl = arcpy.MakeFeatureLayer_management(counties, 'in_memory/counties_lyr')
        coselect = arcpy.SelectLayerByAttribute_management(fl, '', query)
        fc = arcpy.FeatureClassToFeatureClass_conversion(str(coselect), 'in_memory', coname)
        arcpy.Intersect_analysis([Colorado, fc], (str(coco) + '/' + coname))

I keep getting the following error for the Make Feature Layer tool:

ERROR 000358: Invalid expression

In the past I used a different syntax, but it only produced empty shapefiles, if no errors:

query = "ST_NAME" = ('"' + stname + '"') AND "POLYGON_NM" = ('"' + coname + '"')


texts: Excel XLSX data with coordinates displayed as points

texts:
I am preparing a course for QGIS newbies, which includes handling Excel spreadsheet data in QGIS.

It's so easy: Just drag an XLS, XSLX or CSV from Windows Explorer to QGIS!

For data with no coordinates, I teach them 1:1 joins to shape files to obtain positions, e.g. a place name. For data with coordinate attributes in the CSV, I teach them how to use the "Add Delimited Text Layer" function, and the same goes for Excel XLS files using XY Tools plugin.

The problem is XLSX files, which can easily be dragged and dropped into QGIS and displayed as attributes. But I haven't found a way to plot points from (correct) coordinate information - without previous conversion to XLS or CSV. I would like the operation to be smooth and inviting for the QGIS newcomers.

Question: How can you drag-drop an Excel XLSX file and plot the coordinates using the fewest steps?


texts:
Is it possible to arrange multiple attribute tables in tabs in QGIS?

It annoys me that when opening a new attribute table it appears beside the old ones and even shifts other parts of the GUI around  if there is a lack of space.


texts: Free software to create a map server that can be installed directly in a subdomain without a server install?

texts:
I want to create a web map page using my own shapefiles but all the map server software that I’ve seen need to be installed directly in the server.

Is there any software that allows me to create a map server just by putting the files in some folder and "execute" them (like php or js, for example)?


texts: Point to raster conversion: point not in the center of raster cell

texts:
I created a point dataset by adding a x/y coordinates from an Excel table. Here's my point shapefile:

I am converting it to a raster using cell size of 1 using the point to raster in Conversion Tools. I do not understand why the points are not located in the center of each raster cell. In fact for some cells (on the lower left hand side), the points are in the center while for some it is not. I think this is also leading to the white strip in the middle.

Why are the points not in the center of the raster cells, and how do I fix it?


texts: Delete zero values for Lat-Long columns on a batch of files

texts:
I have a script to delete 0 values from a shapefile attribute table which I have converted into a csv.

setwd("C:/Users/Florian/Documents/R/Evan")
data = read.csv("AF_486_2005_culled_geo_Events_shp.csv", sep=",") #read data as matrix (arrays)

nrow(data) # test how many rows
Not0 <- which(data$LATITUDE == 0) #output which rows = 0
data <- data[-Not0,] # new data = old data with rows != 0
nrow(data) # test how many rows
write.csv(data, file = "AF_486_2005_culled_geo_Events_shp_-00.csv") #output new file

Only problem is I have hundreds of these csv files. How do I apply this script on all of them? I realize this is a novice question, and from my research I know that functions like 'lapply' exist, I am just inept at using them.


texts: How to convert to NoData in Google Earth Engine?

texts:
I have the next code:

I am trying to make a mask (only one shape that I would like to use as an analysis region) with: all elevations less than 4.100 meters (from DEM), all area with water (NDWI) and all slopes less than 15 degrees (from DEM).

I could mask it but I don't know how to join the three layers and I dont know either how to convert the 255 value to NoData.

How can I do it to insert all this variables in one layer and then use that layer as a region?

var region= MNDWI_DPS

//NDWI

var NDWI = ee.Image('MODIS/MCD43A4_NDWI/MCD43A4_005_2016_05_16')

var visualize = NDWI.visualize({bands: ['NDWI'], min: -1, max: 0}).clip(region);

//DEM+SLOPE

var MDE = ee.Image('USGS/SRTMGL1_003')

var terrain = ee.Algorithms.Terrain(ee.Image('USGS/SRTMGL1_003'));

var visualize1 = MDE.visualize({bands: ['elevation'], min: 0, max: 4100}).clip(region);

var slope = ee.Terrain.slope(MDE);

var visualize2 = slope.visualize({bands: ['slope'], min: 0, max: 15}).clip(region)

Map.addLayer(visualize.clip(region))
Map.addLayer(visualize1.clip(region))
Map.addLayer(visualize2.clip(region))


texts:
I have some computation on a numpy array and want to display the result in QGIS. So far ok, but loading the result into QGIS fails. Here my code fragment:

zData2 = np.subtract(myarray,zData2) #do some computes
output_raster = gdal.GetDriverByName('MEM').Create('',ncols, nrows, 1 ,gdal.GDT_Float32)
output_raster.GetRasterBand(1).WriteArray(zData2)   # Writes my array to the raster

newLayer = QgsRasterLayer(output_raster,'my new raster') #this failes...
QgsMapLayerRegistry.instance().addMapLayer(newLayer)


texts: Is there a pgRouting-like library for MySQL?

texts: As described in the OsmDB.pm page I have created a schema in MySQL from an OSM file. Now I have to find the length of the shortest path between two given points (preferably using the A* algorithm). How do I do this? I know pgRouting is the library for PostGIS/PostgreSQL. Is there a similar library for MySQL or can I achieve this using an SQL query?


texts:
I have this function for labeling layers using arcpy
This function works totally fine but I have the field hardcoded into the function

def label(wah):
    for lyr in arcpy.mapping.ListLayers(mxd):
        if lyr.name == wah:
            if lyr.supports("LABELCLASSES"):
                for lblClass in lyr.labelClasses:
                    lblClass.expression = "[field]"
                    lblClass.showClassLabels = True
                    lyr.showLabels = True

i want to be able to change the labeling expression for every layer. I have tried this way below, but the labeling does not show up at all.

def label(wah, field):
    for lyr in arcpy.mapping.ListLayers(mxd):
        if lyr.name == wah:
            if lyr.supports("LABELCLASSES"):
                for lblClass in lyr.labelClasses:
                    lblClass.expression = "'["+field+"]'"
                    lblClass.showClassLabels = True
                    lyr.showLabels = True

Can someone help me fix my code. i am aware many programmers use .format instead of just using string with the + sign. what is the advantage of that?

when i call the function should it be

label("layer", field)

or

label("layer", "field")


texts: Is it possible to do supervised image classification in QGIS? I require it for the estimation of cover management factor (C factor) in USLE


texts: Getting WMS from GeoServer with OpenLayers

texts:
I have got Layers on geoserver. my goal is to show that layer using OpenLayers.

My script is below and returns blank map. How would I fetch and show layer on web?

 <script type="text/javascript">
       var layers = [

       new ol.layer.Tile({
           extent: [-13884991, 2870341, -7455066, 6338219],
           source: new ol.source.TileWMS({
               url: 'http://myIP/geoserver/myWorkspaceName/wms',
               params: { 'REQUEST':'GetMap','BBOX':'437762.766104487,4376882.69760637,522118.932880649,4461344.83731443', 'LAYERS': 'myWorkspaceName:myLayerName', 'TILED': true  ,'SRS':'EPSG:5255','FORMAT':'application/openlayers','WİDTH':'767','HEİGHT':'768'},
               serverType: 'geoserver'

           })
       })
        ];
        var map = new ol.Map({
            layers: layers,
            target: 'map',
            view: new ol.View({
                center: [-10997148, 4569099],
                zoom: 4
            })
        });
    </script>


texts: Classifying GPS locations according to Land area classification using ArcPy?

texts:
I have a vector layer for the GPS locations and a raster for the landuse/landcover. I've realized a supervised classification to classify the categories in the Lulc and a maximum likelihood classification (even though I do not know if this last step was useful).

The ideal would be to have a table with the number of location points found for each categories. Like this I can have a percentage, in which categorie ( ex forest) we found 56 locations so 21% for example.

Can anyone help?


texts: Is it possible to force uppercase in a geodatabase string field?

texts:
I have certain string fields in my geodatabase tables with which I would like to enforce an all-uppercase rule. These fields are ones that do not lend themselves to domains (such as address fields).

Is it possible to create a geodatabase string field that supports only uppercase characters inherently? I would like to be able to use the built-in ESRI tools within ArcMap to edit my tables and have them converted to uppercase on the fly.

I am using ArcSDE and MS SQL-Server. Answers that utilize front end validation on ESRI's stack are preferred, but a generalized SQL solution that parses inputs on the RDBMS side with something like an INSTEAD OF trigger is acceptable as well.

Edit: As the accepted answer consists of SQL editing on an RDBMS backend of ArcSDE, Devdatta Tengshe has noted that it can be very dangerous to do this if you do not know what you are doing. It could possibly corrupt your data.

If your situation is conducive, you should use BradHard's answer instead (creating a database view as your data source).

ESRI does not actually discourage the usage of SQL triggers on the backend to maintain data. However, ESRI instructs in their documentation that you should not use SQL on your RDBMS to edit the following ArcSDE data:

  Do not use SQL to edit the following:

  ArcSDE geodatabase system tables
  ObjectIDs (or row IDs) that are maintained by ArcSDE or the geodatabase
  Tables or feature classes that participate in the following geodatabase functionality:

  Annotation feature classes and feature-linked annotation
  Archived data
  Cadastral or parcel fabrics
  Domains
  Geodatabase replication
  Geometric networks
  Network datasets
  Relationship classes
  Representations
  Schematic datasets
  Subtypes
  Survey datasets
  Terrains
  Topology

As long as the fields that you want to force uppercase do not fail any of the requirements listed above, then you can safely use Devdatta Tengshe's SQL trigger method below to force uppercase on whatever individual string fields you desire.


texts: Calculating SCS curve numbers using DSM instead of DTM

texts:
Are the curve numbers (calculated by SCS Curve number method) through DSM (Digital Surface Model) valid? Most of the SCS curve number studies focus on calculating curve numbers through a DTM (Digital Terrain Model).

As the DSM has heights of man made (buildings, roads, siphon structures etc.) and natural features (rocks, trees etc.) is it better to model bridges on it as well?


texts: Converting non equidistant point data to raster?

texts:
I have a csv file in which first column is longitude, second is latitude and third is temperature. However, the points are not at all equidistant i.e. for each pair of lat-lon I have a value. I want to convert it to 2D raster.

I tried with matlab (triangulation and griddata interpolation), but it was giving me 3D raster (surface).

Also, in ArcGIS I tried using XY data as x-lon, y-lat and z-temp but it was not showing georeferenced points. Is there any way ?


texts: PostGIS ST_Transform seems to treat the vertical datum but not the equivalent pyproj Transformer, even so they both rely on Proj

texts:
I've narrowed down a problem I met when processing some points having altitude values to the following:
It seems that PostGIS and pyproj do not treat the vertical datum in the same way but both are actually based on Proj.
Detailed explanation
Here's a simple reproducible PostgresSQL query and its response:
SELECT
  ST_AsText(
    ST_Transform(
      ST_SetSRID(
        ST_MakePoint(8.37909, 47.01987, 1000),
      4326),
    2056)
  ) AS geom;

Here, the altitude of 1000m seems to be evaluated as an ellipsoidal height (relative to the WGS84 ellipsoid I guess?)
Hence =>

geom (text)

"POINT Z (2671499.8913080636 1208075.1135782269 951.4265527641401)"

And here is the same transformation in Python using pyproj:
import numpy as np
from shapely import ops as sp_ops
from shapely.geometry import Point
import pyproj
from pyproj import transform, Proj, Transformer

project = Transformer.from_crs('EPSG:4326', 'EPSG:2056', always_xy=True)
point = Point(8.37909, 47.01987, 1000)
geom = sp_ops.transform(project.transform, point)

np.set_printoptions(suppress=True)
print(np.array(geom.coords))

which results in:
[[ 2671499.8913080636   1208075.1135782297    1000.  ]]

As you can see here, the altitude stays the same after the transformation.
I also tried using the EPSG code for the true 3D WGS84, namely EPSG:4979 according to: Vertical Datum transformation using Pyproj and https://epsg.org/crs_4979/WGS-84.html?sessionkey=00zgez8kt5 but the result is exactly the same.
At this stage, you can notice that PostGIS' result is actually "accurate" according the local geoid-ellipsoid separation which is roughly 48.21m at this position, so 1000 (ellipsoidal_height) - 48.21 (separation) ≅ 951.79 (orthometric_height).
So, why does the altitude stay the same when using pyproj and not when using PostGIS when 1) both actually rely on Proj itself and 2) the EPSG codes I used are the same in both cases?
May it be that PostGIS is a little to lenient?
It is said on the ST_Transform documentation that:

Requires PostGIS be compiled with Proj support. Use
PostGIS_Full_Version to confirm you have proj support compiled in.

and when I actually check for its full version, it's aligned with that statement:
SELECT PostGIS_Full_Version();

"POSTGIS=""3.1.0 5e2af69"" [EXTENSION]
  PGSQL=""130""
  GEOS=""3.7.1-CAPI-1.11.1 27a5e771""
  PROJ=""Rel. 5.2.0, September 15th, 2018""
  LIBXML=""2.9.4""
  LIBJSON=""0.12.1""
  LIBPROTOBUF=""1.3.1""
  WAGYU=""0.5.0 (Internal)""
TOPOLOGY"

On the other hand, my pyproj version is:
>>> pyproj.__version__
'3.0.0.post1'

But I cannot figure out any reference to altitude on vertical datum in the pyproj documentation https://pyproj4.github.io/pyproj/stable/api/transformer.html
I would be very grateful if someone had a clear explanation as to why things are currently the way they are. And how could I properly handle the case in Python, i.e. the same way PostGIS does?


texts: Displaying "World From Space" projection pole-to-pole?

texts:
The "Earth From Space" projection automatically displays the planet with a slight tilt to it. While I am sure this is more accurate, I need to display it un-tilted. I've tried adjusting projection specifications, but I can't seem to make it sit straight.

Does anyone know if this can be done?

There must be a way to show a planetary view from another angle.


texts: ArcGIS treats two similar projections definition differently

texts:
I am about to recreate a global grid in ArcMAP, based on the Behrmann projection I though.
The map I want to recreate has this definition (in the .prj file):
PROJCS["WGS_1984_Cylindrical_Equal_Area",
GEOGCS["GCS_WGS_1984",
DATUM["D_WGS_1984",
SPHEROID["WGS_1984",6378137.0,298.257223563]],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]],
PROJECTION["Cylindrical_Equal_Area"],
PARAMETER["false_easting",0.0],
PARAMETER["false_northing",0.0],
PARAMETER["central_meridian",0.0],
PARAMETER["standard_parallel_1",30.0],
UNIT["Meter",1.0]]

Which I though was the exact same as ArcMAPs Behrmann projection:
PROJCS["World_Behrmann",
GEOGCS["GCS_WGS_1984",
DATUM["D_WGS_1984",
SPHEROID["WGS_1984",6378137.0,298.257223563]],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]],
PROJECTION["Behrmann"],
PARAMETER["False_Easting",0.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",0.0],
UNIT["Meter",1.0]]

But ArcMap treats them slightly differently and can project one to the other. And they don't exactly line up. What is the difference?


texts: How to limit input of Python script parameter to polyline only when using Layer file?

texts: When the parameter is feature class I can filter the input to polyline only. How to do this when the input is layer file at CURRENT MXD? I'm looking for something that will operate before the script starts running. I'm using Layer file so the user can choose the input from the data frame at the current MXD.


texts: Calculate Latitude and Longitude from Range, Azimuth, and Elevation

texts: I was wondering if there is a way to calculate the Latitude and Longitude of a given point provided you know the range, azimuth, and elevation of said point?


texts: ArcPy - SDE Connection Causes "pythonw.exe" Stops Working

texts:
I came across this problem when tring to access some layers in the enterprise database (Oracle). The codes are written as:

cwd = sys.path[0]

# Process: Create Database Connection
DBConnFile_Name = "Enterprise.sde"
DBConnFile = os.path.join(cwd, DBConnFile_Name)

arcpy.CreateDatabaseConnection_management(
    cwd,
    DBConnFile_Name,
    "ORACLE",
    "<instance name here>",
    "DATABASE_AUTH",
    "<username here>",
    "<password here>",
    "SAVE_USERNAME",
    "",
    "SDE",
    "TRANSACTIONAL",
    "SDE.DEFAULT",
    ""
)

# Get the enterprise database layers
in_Council_Pg = os.path.join(DBConnFile, "DISTRICT_COUNCIL")
in_CrimePrev_Pg = os.path.join(DBConnFile, "DISTRICT_CRIME_PREVENT")
in_ZIP_Pg = os.path.join(DBConnFile, "ZIP")
in_Neighborhood_Pg = os.path.join(DBConnFile, "NEIGHBORHOOD")

I want to use those layers in enterprise database for further processing using other ArcPy functions. When I use the layers saved at local directory as arguments, the entire program runs perfectly fine. But when I try to use the layers in SDE as arguments, Windows tells me that "pythonw.exe has stopped working". One of the functions is written as:

# Process: Spatial Join - Council Districts
arcpy.SpatialJoin_analysis(
    CrimeIncid_Pt_WGS,
    in_Council_Pg,
    CrimeIncid_Pt_Council,
    "JOIN_ONE_TO_ONE",
    "KEEP_ALL",
    "",
    "WITHIN",
    "",
    ""
)

print "Stops here."

The programs fails at the function before printing "Stops here".

Do anybody have similar issue?


texts: Are there any known problems using QGIS on Mac OSX Mountain Lion?

texts:
A few days ago Apple released a new version of the OSX (10.8) and I'm thinking about upgrading.

So, I was wondering if anyone had any problems with GIS applications (QGIS) and other similar software after the upgrade?

Thanks,


texts: How to obtain to Max zoom levels for Google layer using OpenLayers

texts:
I was trying to obtain maximum zoom level ( level 22) for Google street map using Openlayers.layer.google script.

The following code does not seem to work. It works to the maximum zoom  level when the zoom variable is set to 15, but below or above 15, it sets centre to world extent.

Any idea what went wrong? we tried using Google maps API v2 and v3, same result for Maxmimum zoom, except for value = 15.

var position = new OpenLayers.LonLat(lon, lat).transform(fromProjection, toProjection);
var zoom = 22;
map.setCenter(position, zoom);


texts: Map Service Update - clearing browser temp files

texts:
I'm noticing whenever I make a change to an ArcGIS map service (add layer, remove layer...etc) I always have to clear server cache and my internet browser temp files in order for my flex viewer app to not give an error.  This can be problematic for any user viewing the application before and after an update has been made to the service.  What is the strategy to prevent these errors from occurring on users accessing an app that may need periodic updates?

Changing the application url (I'm guessing there must be something better than this)?


texts: How to remove last drawing feature in OpenLayers 3?

texts:
I have already found Undoing last point when drawing linestring in OpenLayers 3? (with working Plunker) but it can't undo misclicks when drawing polygon. So what I want is remove last drawing.

Does anyone know how to do that?


texts: Finding minimum-area-rectangle for given points?

texts:
As you see in the figure,  the question is:

How to find the minimum-area-rectangle (MAR) fitted on the given points?

and a supporting question is:

Is there any analytical solution for the problem?

(A development of the question will be to fit a box (3D) to a cluster of points in a 3D point cloud.)

As a first stage I propose to find the convex-hull for the points which reforms the problem (by removing those points are not involved in the solution) to:
fitting a MAR to a polygon.
The required method will provide X (center of rectangle), D (two dimensions) and A (angle).

My proposal for solution:

Find the centroid of the polygon (see Finding center of geometry of object?)
[S] Fit a simple fitted rectangle i.e., parallel to the axes X and Y

you may use minmax function for X and Y of the given points (e.g., polygon's vertices)

Store the area of the fitted rectangle
Rotate the polygon about the centroid by e.g., 1 degree
Repeat from [S] until a full rotation done
Report the angle of the minimum area as the result

It looks to me promising, however the following problems exist:

choose of a good resolution for the angle change could be challenging,
the computation cost is high,
the solution is not analytical but experimental.


texts: How to quickly turn regular spaced points into a grid?

texts:
I'm looking for a quick way in ArcGIS Desktop or Pro to turn a set of regular spaced points into a grid like shown here

Currently I am creating a fishnet first and then conduct a spatial join but I have to manually manipulate the fishnet coordinates to make the points shown as centroids. Is there a tool that creates a full grid like this? I am not looking for the Point to Raster tool as this does not create a full grid as shown here.


texts: Setting feature label in HTML5 Viewer of Geocortex Essentials?

texts:
While I can edit layer properties in Geocortex Essentials 4.4 and change its label name and other properties, feature label never appears on my feature layer. Feature layers are loaded in Geocortex Essentials Site from ArcGIS Online account.

Can any one please guide me with some basic understanding. I m newbie to Geocortex Essentials and the Geocortex Viewer for HTML5.


texts:
This question has been marked as a duplicate multiple times now.  Please read the entire question before flagging it as a duplicate.  The previous "duplicates" for questions about copy parallel do not produce the same result as copy parallel tool.

I am looking for an arcpy process that EXACTLY replicates the copy parallel editor tool in ArcGIS 10.3.  This is not the same thing as simply offsetting lines with arcpy.  That is simple, but it does not replicate the copy parallel tool.

I am using ArcMap with a STANDARD licence.  Advanced licence functions won't work as a solution (buffer with flat endcaps)

As an example, I have run the script from this answer https://gis.stackexchange.com/a/181545/56638.  I used an offset of x = 20 and y = -20
which produces the red lines in the image below (yellow line is original line).

Then I used the copy parallel tool and set the offset to 20.  That produced the green lines in the image.  These lines are exactly 20m from the original line.

The problem is that when offsetting with the arcpy solution, the angle of the line is not considered.  Therefore, the red lines (produced using arcpy method) are more than 20m from the original line.  If you want the offset line to be a certain distance from the original line, you would need to know the angle of the line and calculate the distances to offset x and y based on trig.  I am struggling to come up with a solution for this.

Another example with the arcpy solution, with offset set to x = 20 and y = 0.  This demonstrates the problems that bends cause and highlights the need to take into account the angle or bearing of the original line.

This is the code that I have tried so far:

tmpwrkspace = r'C:\Users\jasbal\Documents\ArcGIS\Default.gdb'
tempFeatureClass = tmpwrkspace + "\temp_rline"
tempFeatureClass2 = tmpwrkspace + "\temp_lline"
fc_line = r'C:\Users\jasbal\Documents\ArcGIS\Default.gdb\RWP_Rev0_4'
if arcpy.Exists(tempFeatureClass):
    arcpy.Delete_management(tempFeatureClass)
arcpy.CopyFeatures_management(fc_line, tempFeatureClass)
if arcpy.Exists(tempFeatureClass2):
    arcpy.Delete_management(tempFeatureClass2)
arcpy.CopyFeatures_management(fc_line, tempFeatureClass2)
xOffset = 20
yOffset = -20
xOffset2 = -20
yOffset2 = 20
with arcpy.da.UpdateCursor(tempFeatureClass, ["SHAPE@XY"]) as cursor:
    for row in cursor:
        cursor.updateRow([[row[0][0] + xOffset, row[0][1] + yOffset]])
del cursor
with arcpy.da.UpdateCursor(tempFeatureClass2, ["SHAPE@XY"]) as cursor:
    for row in cursor:
        cursor.updateRow([[row[0][0] + xOffset2, row[0][1] + yOffset2]])
del cursor

My idea is to calculate the angle of the line and to use that in a trig function to calculate x and y offset based on the distance I want to be from the original line.

Something like this:

distance = arcpy.GetParameterAsText(0)
angle = "Angle Field" obtained with search cursor
xOffset = distance/cos(angle)
yOffset = distance * tan(angle)

I am not sure if I have my trig correct.  I need some advice on this part.

I am not looking for a complete code snippet for copy parallel.  The only part I need to figure out is how to calculate the offsets correctly for a set distance from the original line based on the line angle.


texts: How to obtain a live route on a carto map?

texts:
My main goal is to build an app that send lat lon from my smartphone to my carto (postgis) database and show the route "live" (in real time) as a line and not as a lot of points.

i am able to save every 40 seconds lat and long in the database as point using SQL API

example:

Cartourl/user/myusername/api/v2/sql?q=INSERT INTO my_table (description, the_geom) VALUES ('myroute', ST_SetSRID(ST_Point(lat, lng),4326))&api_key=myapykey

now i think i have to use ST_MakeLine to obtain a line from points.
but can i do this every time that i submit a point? has to be done before?
which query i have to use?

or i can do something like:

Cartourl/user/myusername/api/v2/sql?q=INSERT INTO my_table (description, the_geom) VALUES ('myroute', ST_SetSRID(ST_MakeLine(lat, lng),4326))&api_key=myapykey

or

SELECT
  ST_MakeLine(the_geom_webmercator ORDER BY ST_X(the_geom), ST_Y(the_geom) ASC) as the_geom_webmercator,
  min(cartodb_id) as cartodb_id
FROM
  my_table


texts:
I am new to python, so bear with me.
Trying running reclassify in pyscripter returns the error:

Message File Name   Line    Position
Traceback
           21
    Reclassify  C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\sa\Functions.py   5236
    swapper C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\sa\Utils.py   47
    wrapper C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\sa\Functions.py   5230
        C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\geoprocessing_base.py    498
ExecuteError: ERROR 999999: Error executing function.
Failed to open raster dataset
Failed to execute (Reclassify).

My code is simple, and it works in the python window in arcmap, so I am guessing it has something to do with how i set pyscripter up. I successfully tried SplineWithBarriers, so it is working:

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

# Set environment settings
env.workspace = "V:/USER/b/Processing"

# Set local variables
inRaster = "barr7"
reclassField = "VALUE"
remap =  RemapValue([[1,9],[2,8],[3,1],[4,6],[5,3],[6,3],[7,1]])

arcpy.CheckOutExtension("Spatial")

# Execute Reclassify
outReclassify = Reclassify(inRaster, reclassField, remap, "NODATA")


texts: How to validate OpenLayers area calculation using PostGIS?

texts:
I am working on a GIS API that work on top of OpenLayers.

I have a feature that calculates the area of a drawn polygon. I have tested it on the browser and the areas I got, appear to be correct. For example, I have measured some country areas and they matched the real areas approximately.

Surely, I have to run some unit tests to prove that the area calculations are correct. I use PostGIS to generate the area value, but I get completely different results than the ones from OpenLayers.

In PostGIS, I do this:

SELECT
st_astext(st_transform(st_geomfromtext('POINT(50 0)', 4326), 900913)),
st_astext(st_transform(st_geomfromtext('POINT(-60 -5)', 4326), 900913)),
st_astext(st_transform(st_geomfromtext('POINT(-30 -30)', 4326), 900913)),
ST_Area(
    ST_GeomFromText(
        'POLYGON(
            (
                50 0,
                -60 -5,
                -30 -30,
                50 0
            )
        )'
    ), false
) as area

I get these results:

"POINT(5565974.53966368 -7.08115455161362e-10)";
"POINT(-6679169.44759641 -557305.257274577)";
"POINT(-3339584.72379821 -3503549.84350437)";
19582793771831 // area

I introduce these three points in my test, like this:

var points = [
    new OpenLayers.Geometry.Point(5565974.53966368, -7.08115455161362e-10),
    new OpenLayers.Geometry.Point(-6679169.44759641, -557305.257274577),
    new OpenLayers.Geometry.Point(-3339584.72379821, -3503549.84350437)
];
var ring = new OpenLayers.Geometry.LinearRing(points);
var polygon = new OpenLayers.Geometry.Polygon([ring]);

Then, I used the getArea() method from OpenLayers.Measure, but I got an area of: 17050040174394.607 sq. meters

According to PostGIS, the area of this polygon should be: 19582793771831 sq. meters.

It is a big difference. what could I be doing wrong here?

Thanks

PS: I can't publish the whole code here, because it is work related.

EDIT 1 - CLOSE THE LINEAR RING

var points = [
    // Point 0 - Start point
    new OpenLayers.Geometry.Point(5565974.53966368, -7.08115455161362e-10),
    // Point 1
    new OpenLayers.Geometry.Point(-6679169.44759641, -557305.257274577),
    // Point 2
    new OpenLayers.Geometry.Point(-3339584.72379821, -3503549.84350437),
    // Point 3 - End point (the same as start point)
    new OpenLayers.Geometry.Point(5565974.53966368, -7.08115455161362e-10),
];

I tried to close the Linear Ring, but I still get: 17050040174394.607 sq. meters


texts: Export from SDE to SHP in ArcGIS 9.3.1 - problem: not exporting all fields

texts:
I'm exporting a layer from SDE to SHP by python script, and I'm doing it in 2 phases: the first includes select by location (for picking the relevant features) and the second phase includes the Select (Analysis) tool which exports the features to the output SHP.

PROBLEM:

The Select(Analysis) tool omits fields - in my case, OBJECTID and GLOBALID - which need to be in the output. (I know it's the tool's fault because when I do a manual data>export data, it exports the fields all fine and dandy.)

Does anybody know how I can fix this?


texts: Display feature count in atlas legend only for features included on atlas page?

texts:
I am working with QGIS 2.16, and I am wondering if there is a way to include a feature count only for the features displayed in each atlas page. For example, San Francisco only has 120 of a certain feature, while the total count for all metro areas in the dataset is 1500. Is it possible for the atlas page legend to only display the 120 and not the 1500?

Apologies for any lack of clarity. The question is simple but a little hard to explain.


texts:
I made a Join between a shape and a dbf, and now I want to calculate some data from the table that came as a result from that Join. Everytime I try to run the Field Calculator it says there's an error. I read I should re-save that table as a different shape so I can edit the attributes, but:

I don't know if that correct and
I don't know how to do that.

I've using ArcGIS for Desktop for a very short time.


texts: GDAL 2.1 DelaunayTriangulation Give bad results (sometimes)

texts:
I am using the latest stable release of GDAL 2.1.  When I add vertexes to a polygon in a loop sometimes it does not return the correct results.  If I put those same values in hard coded, the right result is returned.

OGRPolygon     poly;
OGRLinearRing  ring;
OGRLineString  line;
for(pVert = Vertexes.begin(); pVert != Vertexes.end(); pVert++)
    ring.addPoint(pVert->X, pVert->Y);
poly.addRing(&ring);
poly.closeRings();

OGRMultiPolygon* poOutPolys = (OGRMultiPolygon*)poly.DelaunayTriangulation(0.0001, 0);
int NumOut = poOutPolys->getNumGeometries();

This results in the following values in the debugger watch window, but only 2 triangles returned (there should be 3 with a 5 vertex polygon):

-       ((*(OGRSimpleCurve*)(&(*(OGRLineString*)(&ring))))).paoPoints,5 0x07b1a330 {x=2771.5700000000002 y=4720.0959999999995 } OGRRawPoint *
+       [0] {x=2771.5700000000002 y=4720.0959999999995 }    OGRRawPoint
+       [1] {x=2774.0639999999999 y=4722.5940000000001 }    OGRRawPoint
+       [2] {x=2589.1100000000001 y=4722.5940000000001 }    OGRRawPoint
+       [3] {x=2590.6970000000001 y=4721.0050000000001 }    OGRRawPoint
+       [4] {x=2640.1439999999998 y=4720.0490000000000 }    OGRRawPoint

If is use the following code, 3 triangles are returned.

OGRPolygon     poly2;
OGRLinearRing  ring2;
ring2.addPoint(2771.570, 4720.096);
ring2.addPoint(2774.064, 4772.594);
ring2.addPoint(2589.110, 4722.594);
ring2.addPoint(2590.697, 4721.005);
ring2.addPoint(2640.144, 4720.049);
poly2.addRing(&ring2);
poly2.closeRings();
OGRMultiPolygon* poOutPolys2 = (OGRMultiPolygon*)poly2.DelaunayTriangulation(0.0001, 0);
int NumOut2 = poOutPolys2->getNumGeometries();

Debug watch window:

-       ((*(OGRSimpleCurve*)(&(*(OGRLineString*)(&ring2))))).paoPoints,5    0x07b1a0e8 {x=2771.5700000000002 y=4720.0959999999995 } OGRRawPoint *
+       [0] {x=2771.5700000000002 y=4720.0959999999995 }    OGRRawPoint
+       [1] {x=2774.0639999999999 y=4772.5940000000001 }    OGRRawPoint
+       [2] {x=2589.1100000000001 y=4722.5940000000001 }    OGRRawPoint
+       [3] {x=2590.6970000000001 y=4721.0050000000001 }    OGRRawPoint
+       [4] {x=2640.1439999999998 y=4720.0490000000000 }    OGRRawPoint

What is going on?  Am I doing something wrong?


texts: Merging multipolygon data into a polygon

texts:
I have a MultiPolygon representation of some shapes in POSTGIS, I am trying to convert them to Polygon data.

I tried, st_geometryn(r.geom, 1) but that didnt work because the way polygons are embedded I think. Just returns a partial.

Tried ST_ConcaveHull(ST_Collect(r.geom), 0.99), no luck returns nothing

Tried:

SELECT f.gid, ST_AsGeoJSON(ST_Union(f.geom)) as singlegeom FROM (SELECT gid,
       (ST_Dump(geom)).geom As geom FROM public.nyc_boroughs ) As f GROUP BY gid"

but getting back to Multi Polygon,

I feel like a little bit of direction could help me figure out.

Data is here http://www.sendspace.com/file/2savyw


texts: How to nest GeoJSON features based on their geometry in javascript?


texts: QGIS Several overlap / intersect geoprocessing

texts:
First time asking here, I hope you guys can help me:

So I have this polygon layer with multiple overlapping buffer polygons (these buffers are based on a point layer).

I also have another layer (let's call it "forest") from which I need to extract polygons that are covered by the buffers (using either Clip or Intersection).

My problem is if two buffers overlaps, nor the clip tool or the intersect tool can duplicate the polygons so each buffer has "their own" polygons extracted from "forest".

That might not be clear, so here's a situaton: Let's say that Buffer B1 and Buffer B2 completely overlaps. I need polygons from "Forest" covered by B1 so I can run stats on them, but I would need copies of the same polygons (since B1 and B2 completely overlaps) associated with B2 so I can run the same stats for B2.

I could do it the long way by selecting each buffer individually and clip "Forest", so I would have one layer by buffer, and then merge them all togheter, but I'd like to save time ...

Any clue?


texts: Unable to retrieve treeID after segmentation

texts:
I encountered a problem that treeID feature can not be retrieved after implementing Dalponte2016() chm-based segmentation method. The segmentation results can correctly show in multiple colors by using
plot(las_res, bg = "white", size = 2, color = "treeID")

but these treeIDs are actually not available (NA), and I cannot export the treeIDs in both las/laz and .csv formats.


texts: Nearest point does not intersect line in 'sf' R package

texts:
This is a similar question to my recent questions on how to find the nearest point along a polyline and how to split a polyline by an intersecting point using only the sf package.

My question has mainly to do with the task in the first question I asked but addresses an unresolved problem: given 378 or more points, find the nearest points that intersect a line. I know what functions to use in sf, the problem is that many of the nearest points are still some small distance (e.g., 10^-11 units away) from the line.

Using sf_0.6-4 the st_nearest_feature and st_nearest_points functions will locate the nearest point along a polyline. However, due to the high precision of the coordinates the nearest point is still a very small distance from the reach and does not intersect it. According to the sf author setting the precision may help, but I've had no success in my own attempts to do intersects using a precision of 1000 (three digits following the decimal).

I've also tried using st_snap with a tolerance of 1e-09 but it does not consistently snap the point; in fact this tool should be unnecessary if the precision is used correctly (which evidently I am not doing). Regardless of messing around with precision settings or snapping the points, the points do not consistently intersect the line.

Does anyone know how this can be done? I've included all sites and nearby rivers and a script below demonstrating my basic approach with sf, along with examples of how I've used st_set_precision and st_snap. I've uploaded the input data here.

library(sf) # requires version 0.6.4
library(tidyverse)

pts <- st_read("./inputs", "sites", stringsAsFactors = F)
rivers <- st_read("./inputs", "rivers_buffers", stringsAsFactors = F)

pts.near <- lapply(pts$SiteId, function(i){
  # Select site, site buffer, and reaches within the site buffer
  pt <- filter(pts, SiteId==i) # get individual site
  pt.buffer <- filter(pt.buffers, SiteId==i)

  rivers.buffer <- st_intersection(rivers, pt.buffer)

  # # Set precision
  # p <- 1000
  # pt <- st_set_precision(pt, precision=p)
  # rivers.buffer <- st_set_precision(rivers.buffer, precision=p)

  # Select nearest reach
  reach.nearest <- st_nearest_feature(pt, rivers.buffer)
  reach.nearest <- rivers.buffer[reach.nearest,]

  # # Locate nearest point along nearest reach (Returns line)
  pt.nearest <- st_nearest_points(pt, reach.nearest)
  pt.nearest <- st_cast(pt.nearest, "POINT")[2]

  # Snap point
  pt.nearest <- st_snap(pt.nearest, reach.nearest, tolerance=1e-09)

  pt.dt <- data.frame(SiteId = i,
                      ReachId = reach.nearest$FID,
                      Intersects = st_intersects(pt.nearest, reach.nearest, sparse=FALSE)[1,1],
                      Distance = st_distance(pt.nearest, reach.nearest),
                      X = st_coordinates(pt.nearest)[1],
                      Y = st_coordinates(pt.nearest)[2]
                      )
})

result <- do.call(rbind, pts.near)

PS; if I can't solve this in R I will commit heresy and use ArcGIS.


texts: I've been working in the alpha of 1.9 for quite a while (because of doing bug reports) and I like it a lot. It's more than a year since 1.8 was released, so I wonder when 1.9 is expected out.


texts: What QGIS features are available for shapefile layers but not PostGIS layers?

texts:
Our office has a shared Network drive with limited storage and shared access to a PostGIS DB. We're trying to migrate away from using shapefiles in QGIS to pulling layers from PostGIS.

Are there any editing/styling features that wouldn't work with a PostGIS layer that do with a shapefile layer?

One thing I can see is that PostGIS layers would have to be stored de-normalized to enable editing.

e.g. this

streets:
gid|street_name|geom
---|-----------|----
1  |'Main St.' |line
2  |'Main St.' |line

not this

streets:
gid|street_name_id|geom
---|--------------|----
1  |1             |line
2  |1             |line

street_names:
street_name_id|street_name
--------------|-----------
1             | 'Main St.'


texts: Measuring distance between two points using digital elevation model (DEM) and ArcMap

texts:
I have a number of points mapped on a digital elevation model (DEM). The points were mapped on the DEM based on their GPS location (latitude, longitude).
What is the best way to measure the distance between two points on the map?
The measuring tool?
Can I use the DEM to calculate the distance?
If so, how?


texts:
I have successfully manage to make the workshop from of pgrouting (wrote for v1) with pgrouting v2.
Now I want to use pgrouting 2 on my postgis database created with osm2pgsql.

The problem is first to create the segments and this blog post was very helpfull for this.
However,I get:

gis=# select pgr_createTopology('network', 0.01, 'the_geom', 'gid');
ERROR:  We can not determine the srid of "network.the_geom"!

Is there a way to set the srid properly so that pgr_createTopology finds it.
Or is there an alternative way?
I can't believe osm2pgrouting is the only importer for pgrouting.
Thanks


texts: How to read geographic coordinates when the Shapefile has a projected spatial reference

texts:
I have a polygon Shapefile and want to read the polygon's vertices using OGR and Python. The reference system of the file is projected. However, for further processing I need the geographic coordinates of the vertices.

This is my code:

import ogr
path_to_infile = '/path/to/shapefile.shp'
driver = ogr.GetDriverByName('ESRI Shapefile')
infile = driver.Open(path_to_infile)
layer = infile.GetLayer()
for area in layer:
    area_shape = area.GetGeometryRef()
    area_polygon = area_shape.GetGeometryRef(0)
    no_of_polygon_vertices = area_polygon.GetPointCount()
    for vertex in xrange(no_of_polygon_vertices):
        lon, lat, z = area_polygon.GetPoint(vertex)

Coordinates (lat & lon variables) will be read according to the spatial reference. In my case, projected coordinates are displayed. Is there a way to directly extract geographic coordinates in OGR?

I know that .prj files of projected reference systems also store information on the corresponding geographic coordinate system. In arcpy for instance, geographic coordinates of a Shapefile with projected spatial reference can be accessed using the GCS property of arcpy.Describe(dataset).spatialReference:

import arcpy
# accessing projected coordinates
spatial_ref = arcpy.Describe(dataset).spatialReference

# accessing geographic coordinates
spatial_ref.GCS

When spatial_ref.GCS is defined, I get geographic point coordinates for a polygon's vertices.

Is there a similar property for OGR available? Reprojecting is not my first choice, as I'm working in a custom reference system.


texts: Multiple memory transfers Python Script using Gdal (gdalwarp) w/ GPU on JPEG2000 Datasources

texts:
I'm looking at writing a python script that will open a jpeg2000 datasource and decompress it on the GPU, then warp that data using gdalwarp with opencl (also on the GPU)
My question is, will the dataset be copied to ram from the video memory after it's opened, then back to the GPU for the warping, then back to the RAM? or will the data stay in the GPU memory between opening being that both operations happen on the gpu?

If gdal does this in the more efficient way, then I'll be wanting to buy a license or Kakadu jpeg2000 driver but I don't want to fork up the money for it unless there is a good indicator that this is going to work as expected.


texts: Selecting features in attribute table using maximum aggregate function

texts:
I have a table containing distances from one ID (an address) to another ID (bus stops).

There table contains distances from all addresses to all different bus stops.

Now I want to select the feature where the distance between the certain address and the bus stop is the shortest.

In QGIS I've tried to do something like this:

where:

"meter" is the distance column;
"fromID" is the address.

There are 16 different addresses and a little more than 6000 features in the table. What I need is a table with only 16 features.

You can, of course, use minimum function as well.

Can anybody help out here?


texts: Is it possible to get all layers and/or metadata from a Mapbox GL source?

texts:
In a mapbox-gl map, given the source:

map.addSource('myTileSource', {
    type: "vector",
    tiles: "url to source tiles.pbf",
});

is it possible to get all the layers and/or metadata for the tiles inside that source? This is possible in node with node-mbtiles, but wondering if it is possible directly in mapbox-gl via javascript? Something like map.getLayersFromSource would be awesome.


texts: Error when calculating final polygon after subtracting all intersecting polygons using shapely

texts:
This question is a follow-on question from Calculating final polygon after subtracting all intersecting polygons using shapely which was answered.
I have a list of centroids that I buffer to 5000 meters radius. I need to take each resulting circle and run it against the remaining centroids(also buffered to 5000 m) and find the final shape after subtracting all intersections.
In other words, I need the circle shape A to iterate over the dataset, find all circle shape(s) that intersect with it (B and C) and then create a column with a geocoordinates of the polygon, that is resulting after all intersecting circle areas are subtracted (see image below) and if it does not intersect with any other circle(s), just add a geocoordinate for the entire circle:

Here is the code I have:
df.head()

# Output:
#                cell_centroid  pop     id
#0  POINT(21.25602 55.903539)   9004.0  0
#1  POINT(21.119156 55.538406)  517.0   1
#2  POINT(21.249535 55.890204)  10369.0 2
#3  POINT(21.359934 55.904164)  800.0   3
#4  POINT(24.433042 56.104589)  800.0   4

df['cell_centroid'] = df['cell_centroid'].apply(wkt.loads)

gdf = gpd.GeoDataFrame(df, geometry="cell_centroid", crs = 4326).to_crs(crs=3857)

gdf["buffer"] = gdf.buffer(5000)
gdf["difference"] = None

gdf = gdf.set_geometry("buffer")

for i, row in gdf.iterrows():
    others = gdf[gdf["id"] != row["id"]]
    diff = row["buffer"].difference(others.cascaded_union)
    gdf.loc[i, "difference"] = diff

gdf = gdf.set_geometry("difference")

The error I am getting now is:
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-122-ad757ced6ebf> in <module>
      3     diff = row["buffer"].difference(others.cascaded_union)
      4 #    gdf.loc[i, "difference"] = gpd.GeoSeries(diff)
----> 5     gdf.loc[i, "difference"] = diff
...

ValueError: Must have equal len keys and value when setting with an iterable

What's the right way to assign the resulting polygon (diff) to the GeoDataFrame?
I have found GeoPandas setting geometry for a row with .loc and Multipolygons for a similar issue, but its solution does not work in my case.


texts: arcpy Make Feature Layer using a where clause

texts:
From a table of lines and polylines I am I'm building a seperate table of points (nodes) with their X and Y coordinates stored as attributes. I want the nodes to be unique, so I'm checking to see if a node with the same coordinates has already been inserted by making a feature layer from the nodes table with a SQL clause to selecting any rows with the same coordinates. If the feature layer has zero rows then I insert the new node otherwise I do nothing and move to the next line/polyline.

The problem I have is that the SQL clause in the Make Feature Layer seems to be selecting rows even where the coordinates do not match. It works fine for the first one but then after that it picks 1 row each time when most or all of the time it should pic zero rows.

This is making my head hurt. Can any one help please?

Here is the script..

import arcpy
arcpy.env.overwriteOutput = True

fc_Nodes = 'C:\Users\Rob&LP\Documents\GIS General\DRN Configuration\DRN Sample\LinkNodes_point.shp'

arcpy.DeleteRows_management(fc_Nodes)

fc_Links = 'C:\Users\Rob&LP\Documents\GIS General\DRN Configuration\DRN Sample\DRN_SevUp_14Recs_polyline.shp'

shapeName = arcpy.Describe(fc_Links).shapeFieldName

rows = arcpy.SearchCursor(fc_Links)
InCurs = arcpy.InsertCursor(fc_Nodes)

Ins = 0

for row in rows:

 feat = row.getValue(shapeName)
 ANode = feat.firstPoint

#check to see if node has already been inserted
 arcpy.MakeFeatureLayer_management(fc_Nodes, "fl_P1", """"X" = ANode.X""")
 if int(arcpy.GetCount_management("fl_P1").getOutput(0)) == 0:

  Ins = Ins + 1

  print "NumRows: " + str(arcpy.GetCount_management("fl_P1").getOutput(0))
  print "Inserted: "+ str(Ins)
  print "ANodeX: " + str(ANode.X)
  print ""

  newNode = InCurs.newRow()
  newNode.SHAPE = ANode
  newNode.setValue("CHECK", 1)
  newNode.setValue("X", ANode.X)
  newNode.setValue("Y", ANode.Y)
  newNode.setValue("ID", Ins)

  InCurs.insertRow(newNode)
  del ANode
  arcpy.Delete_management("fl_P1")
 else:

  print "NumRows: " + str(arcpy.GetCount_management("fl_P1").getOutput(0))
  print "NOT Inserted : " + str(ANode.X)
  print ""
  del ANode
  arcpy.Delete_management("fl_P1")

arcpy.Delete_management("fl_P1")
del InCurs
del rows

print "--------------End---------"

`


texts:
Loading a DEM dataset with .asc extension in Arcmap 10.7 (here the dataset http://www.sardegnageoportale.it/webgis2/sardegnamappe/?map=download_raster - select "DTM passo 10"; here the file http://webgis.regione.sardegna.it/scaricocartografiaETL/modelliDigitali/DTM10m/478.zip)
one of the elements of the mosaic is not displayed correctly. In fact, it appears in a single shade of gray and shows a range of elevation values that do not correspond to the real values.

The same item displayed in qgis 3.10 does not show the same problems and displays correctly.
I tried to open the file with a text editor but without being able to identify the problem by comparing it with the other files of the dataset that are displayed correctly.
How can I solve this problem?


texts: Uncaught TypeError: Cannot read property 'toString' of undefined from ArcGIS JavaScript API?

texts:
I updated my application to 3.12 and am getting the following error:

Uncaught TypeError: Cannot read property 'toString' of undefined

It is being triggered by the following line:

coordMenu.addChild(new dijit.layout.ContentPane({
    content:"

"
}));

Any ideas?  It was working fine in 3.11


texts: Filter on distance point field and choice field

texts:
I have model shop

shop

class Shop(models.Model):
CATEGORY = (
    ('BookShop', 'Bookshop'),
    ('Bakery', 'Bakery'),
    ('Dairy', 'Dairy'),
Shop_category = models.CharField(max_length=200, null=True, choices=CATEGORY)
Distance = models.PositiveIntegerField(null=True,)
Shop_location = models.PointField(null=True)

I want to set filter on choice field and distance.

 The point field get location of user and calculate distance between shop and user

 how i can set so that if the distance calculated between shop and user less then are equal to distance given by shop then shop will print and one more condition choice filed it is compulsory to match?


texts: Calculating area of pixels in continuous raster using ArcGIS Desktop?

texts:
I have a raster with continuous values between 0 and 1.

I want to calculate the area (can be a number of pixels too) in this raster which has values in the specific range - say between 0.1 and 0.2 and 0.3-0.4.

In other words, the goal is to see which interval occupies the larger area.

What tool could be used for that?


texts: Postgres/PostGIS Spatial Query Performance Problem

texts:
I have a query in Postgres/PostGIS that I am a bit confused about. The spatial intersection happens very fast, ~ 2 seconds for 10,000, but the update takes > 15 seconds for the same number of points.

a_table contains ~ 66,000,000 points

b_table contains ~ 3,000 polygons

Both are spatially indexed (gist and clustered) and the b.wanted_value column is indexed as well.

Any suggestions for improving the speed of this query would be sincerely appreciated.

UPDATE a_table a
SET a.needed_value = (
  SELECT b.wanted_value
  FROM b_table b
  WHERE ST_Intersects(a.the_geom, b.the_geom)
);


texts: How can I copy columns from layer 1 to layer 2?

texts:
I have two layers, the first layer is pipe: that has some fields (column) title like : length, id, material. Also I have another layer that is empty (without any fields). How can I copy fields (column) title from layer 1 to layer 2?

I am new to QGIS.


texts: Identifying polygons that lie between two polygons using ArcMap

texts:
I am working on a single layer which includes an island constructed of multiple polygons (used to represent territories).
I would like to look at dispersal between territories, specifically, how many territories are traversed during dispersal. E.g. on the attached schematic island below, dispersal between territory 10 and 6 would result in 1 territory traversed.

I am using ArcMap 10.6.


texts:
I have a raster that represents a classification, this means each pixel can contain a value in a small set (for instance [0,1,2,3,4,5,6,7,8,9,10]). My file is around 2.4 MB which is pretty small considering that it is 5000x5000 pixels. When I create a naive copy using python gdal bindings the output file is 100 MB, around 40 times bigger! (by naive, I mean, I open the raster, extract the data array and write it to a new raster, using gdal_translate yields same result). Looking closely to both files using gdalinfo I see something interesting, in the original file:

Band 1 Block=256x256 Type=Byte, ColorInterp=Gray

in the new file:

Band 1 Block=5000x1 Type=Byte, ColorInterp=Gray

I guess that somehow the tif format can optimise the size when it finds redundant information, for example 256x256 pixel blocks with the same value. The question is, how can I write in 256x256 blocks instead of 5000x1 blocks?


texts: Keeping ids and attributes of source features when creating QgsGraph?

texts:
I followed the example in QGIS documentation for Network Analysis library to create a graph from a vector layer.  When building the graph, the features of the line layer are added as arcs (edges) to the graph. I am looking for a way to assign the arcs with either the original feature ids or their attributes. I have tried two methods but none has worked so far.

Method 1

When building a graph a QgsArcProperter is used to assign costs to the arcs. The QgsArcProperter has a requiredAttributes method, which based on the documentation is called by QgsGraphDirector for fetching attributes from source layer.

When using the QgsDistanceArcProperter, that inherits its methods from the QgsArcProperter, the requiredAttributes method returns an empty list and the properties of the arcs return a list with the cost of the arc.

# build graph with length cost
# vl (a vector layer)
director = QgsLineVectorLayerDirector(vl, -1, '', '', '', 3)
properter = QgsDistanceArcProperter()
director.addProperter(properter)
builder = QgsGraphBuilder(crs)
pStart = QgsPoint(475869, 4203007)
tiedPoints = director.makeGraph(builder, [pStart])
graph = builder.graph()

properter.requiredAttributes()
# >>> returns []

# get the properties of an arc with 0 index
graph.arc(0).properties()
# >>> returns the cost of the arc e.g. [10.89]

When using a custom cost properter as shown at Properter to get travel time as cost for network analysis, then the requiredAttributes return the indices of the given columns. However, when accessing the properties of the arcs, only the cost is returned. From my understanding, the properties should return a list values of the attributes as specified by the requiredAtttibnutes method.

# custom cost properter
class CustomCost(QgsArcProperter):
    def __init__(self, costColumIndex, defaultValue):
        QgsArcProperter.__init__(self)
        self.cost_column_index = costColumIndex
        self.default_value = defaultValue
    def property(self, distance, feature):
        cost = float(feature.attributes()[self.cost_column_index])
        if not cost or cost <= 0.0:
            return self.default_value
        else:
            return cost
    def requiredAttributes(self):
        return [0, self.cost_column_index]

# build graph but replace properter with custom cost
properter = CustomCost(13, 0.01)

properter.requiredAttributes()
# >>> returns [0, 13]

# get the properties of an arc with 0 index
graph.arc(0).properties()
# >>> returns the cost of the arc e.g. [10.89]

Method 2

I have also tried to loop through the polylines and the multi-polylines of the layer to reproduce a sequential index based on their segments. I store the information in a dictionary with arc indices as keys and feature ids as values. Every feature id is added twice as I am using a bidirectional graph.

indices = {}
idx = 0
for f in vl.getFeatures():
    if f.geometry().wkbType() == 2:
        for i in f.geometry().asPolyline()[1:]:
            indices[idx] = f.id()
            idx += 1
            indices[idx] = f.id()
            idx += 1
    elif f.geometry().wkbType() == 5:
        for pl in f.geometry().asMultiPolyline():
            for p in pl[1:]:
                indices[idx] = f.id()
                idx += 1
                indices[idx] = f.id()
                idx += 1

However, if tied points are used when building the graph then this solution does not work. If a point is tied to the graph, then the segment that is closest to it is split and therefore the sequence of the index breaks.


texts: Deleting Service Area Layers using ArcPy?

texts:
Using arcpy I can make a service area layer in Network Analyst like this:

outSAResultObject = arcpy.na.MakeServiceAreaLayer(inNetworkDataset, outNALayerName,
                              impedanceAttribute, "TRAVEL_FROM", "10",
                              "DETAILED_POLYS", "NO_MERGE",
                              hierarchy = "NO_HIERARCHY")

However, where is this created? in "inNetworkDataset" and the GDB associated, but it's not visible in the GDB or ArcMap. It is there somewhere: when I try to run the script again it says

  Failed to execute. Parameters are not valid. ERROR 030036: A layer
  with this name exists. Failed to execute (MakeServiceAreaLayer).

I change the outNALayerName and it works. I would mainly like to know how to remove or delete old service areas, but am also vaguely curious if it is possible to access them using ArcMap or Catalog


texts: Choosing projection for routes that fall within two UTM zones?

texts:
I'm creating charts for routes that fall within two UTM zones; 30S and 31S. The offshore routes are in the Alboran / Balearic Sea - connecting Spain to Algeria.

What projection should I use?


texts: How to build a shortest spanning trees that respects existing connections using QGIS?

texts:
I have

a) a point layer of connection points to an existing network, i.e. an existing network can be accessed here.

b) a point layer of sites that should be connected newly to the existing network

I want to generate the shortest spanning trees to connect the point of layer b) to any of the points of layer a).

With the ReconstructLine plugin I can build the shortest spanning tree between all point of both layers, but this tree is more than I need as it builds a new network between every points. I would need a spanning tree plugin that respects the already existing connections between the points of layer b.


texts:
I have a line table and a points table in PostGIS. I need to split the line table with the points table and return all geometries.

Tried...

CREATE TABLE lines_split AS (
SELECT
a.seg_name, ((ST_DUMP(ST_SPLIT(a.geom,b.geom))).geom) as geom
FROM lines_table a
INNER JOIN points_table b
ON ST_INTERSECTS(a.geom, b.geom)
);

This runs successfully with unusual results. It returns a subset (pink lines) of the lines which are not split by the points...

I have this now...

I want this...

How do I perform this task in PostGIS?


texts: Creating a GRASS mapset from a bash script

texts:
This is a follow up to a previous question on the creation of GRASS locations from a bash script.

Having successfully created a location, I would now like to create a new mapset within this location. The script looks like this:

#!/bin/sh
# GRASS session variables
GRASSDB="/tmp/grassdata"
GRASSLOC="4326" # This is the EPSG code of the CRS used
GRASSMAP="comparison"

# path to GRASS binaries and libraries:
export GISBASE=/usr/lib/grass74
export PATH=$PATH:$GISBASE/bin:$GISBASE/scripts
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$GISBASE/lib

# Create GRASS database folder
mkdir -p "$GRASSDB"

# Create a new GRASS location
echo "==== Creating new GRASS location ===="
grass74 -c -e "epsg:$GRASSLOC" "$GRASSDB/$GRASSLOC"

echo "==== Creating new mapset ===="
grass74 -c -e "$GRASSDB/$GRASSLOC/$comparison"

It produces the following output:

$ sh comparison.sh
==== Creating new GRASS location ====
Cleaning up temporary files...
Creating new GRASS GIS location/mapset...
Cleaning up temporary files...
==== Creating new mapset ====
Cleaning up temporary files...
Creating new GRASS GIS location/mapset...
ERROR: Failed to create new location. The location <4326> already exists.
Exiting...

Unlike what is described in the wiki, the -c argument always attempts to create a new location, and not just a mapset. How can I create just a mapset within an existing location?


texts: Importing and exporting KML into and from DashDB for Geospatial Analysis

texts:
I am a recent convert to IBM's DashDB, and I considering proposing to use it at my work. My case would greatly be bolstered if I could show good, easy integration for geospatial analytics data, namely loading and performing SQL filtering on geodata currently in .shp or .kml formats. If it would be possible to also export the filtered data into a KML as a result that would be AMAZING.

So, to give a practical example, say I have an .SHP file with all the zipcodes in the US, I want to export that into the database, then select the shape for the 02138 zip code and send it to the query-sender in KML format.

Does anyone have experience with that?


texts: How to not display a layer, if its below a other semi transparent layer?

texts:
I have two polygon layers (PostGIS). One of them is semi transparent. If the other layer is below the semi-transparent one, it should not be displayed. My situation at the moment is:

The grey buildings outside the red (semi transparent) area should be displayed, the ones inside should not.

I can not physically edit the layers and have to find a "display-solution" with rules or something.


texts: Is there a way to force QGIS to use scaleable markers for map prints?

texts:
The implementation of scaleable vector graphs in QGIS is very cool, however, when svgs are used for point pattern fill of polygons, these markers are obviously rasterized, blowing up the filesize of pdf-prints tremendously..

This even happens with simple fill for polygons, i.e. polygon filling with diagonal lines, which will be rasterized for pdf-prints. The pdf-filesizes change a 10-fold when using these fill types.

Can anybody tell me if there is an option or workaround to force QGIS to use scaleable markers for map prints?


texts: Using projections to obtain perspective views

texts: I understand that a perspective view is essentially a type of projection. Is it possible to create projection definitions for a typical GIS system that will mimic a perspective view from a particular location? If so, how feasible would it be to use common 2D GIS tools to create perspective renderings of data, similar to what can be achieved in 3D tools such as Arc Scene and Google Earth?


texts: Disabling Torque autoplay + resetting UI

texts:
I’m building a map in CartoDB.js that includes multiple layers including a Torque layer. Because there is other information on the map and a sidebar to read, I don’t want the Torque layer to automatically play.

There seem to be a few easy ways to do this – torqueLayer.stop() / torqueLayer.pause() / torqueLayer.toggle() all work to prevent it from playing. And I saw this previous answer: CartoDB torque chart stop from playing on page load

But they don’t reset the controls UI. With this solution, the icon to the left of the default slider shows two bars for pause, rather than triangle for play. It seems the UI only responds to clicks.

How can I change the UI to match the state of the animation?


texts: I have received multiple tif. files like roads, forest_dis_scaled, grass_dis_scaled. I am trying to create a useful map for the study area but I am not sure if it is even possible to convert these files to shapes and use it on a map. I have tried raster to polygon but it didn't the results I was looking for.


texts:
I'm looking to have AIS data ( Carbon + from True Heading) in real time on QGIS.

From AIS I got the data as following:

!AIVDM,1,1,,A,35?7p455h53k<LH>;dR`wbgoP000,0*7F

After using a specific software to read it that's how it come back:

466004000;under way ;128°'; 0.0kt;25.020450N;51.645125E;306.9°;511°;59s; 150324 000000;serial#1(A)[1]

It seems that i can't load the txt file as coordinates are with "N" & "E".

Did someone tried already to import data from AIS in real time on QGIS?

Likely, someone will know how to transform the data directly from AIS, without using a specific software to read it.