texts: PyQGIS script for loading several .xyz files into QGIS gives NameError

texts:
I found a script for loading several csv files into QGIS and modified it for .xyz Data. Yesterday on one Computer it worked and today on a different Computer it's not working anymore.

Is there something wrong with the code?

I tried to use the following code in the python console

import glob, os

path = "C:\Users\J\Desktop\xyz_data"

os.chdir(path)

for file in glob.glob("*.xyz"):
    uri = "file:///" + path + file + "?delimiter=%s&crs=epsg:4647&xField=%s&yField=%s" % (" ", "field_1", "field_2")
    name = file.replace('.xyz', '')
    lyr = QgsVectorLayer(uri, name, 'delimitedtext')
    QgsMapLayerRegistry.instance().addMapLayer(lyr)

After the "path" line I get the following error message:

  Traceback (most recent call last):
    File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\code.py", line 63, in runsource
      code = self.compile(source, filename, symbol)
    File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\codeop.py", line 168, in call
      return _maybe_compile(self.compiler, source, filename, symbol)
    File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\codeop.py", line 99, in _maybe_compile
      raise err1
    File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\codeop.py", line 87, in _maybe_compile
      code1 = compiler(source + "\n", filename, symbol)
    File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\codeop.py", line 133, in call
      codeob = compile(source, filename, symbol, self.flags, 1)
    File "", line 1
  SyntaxError: (unicode error) 'unicodeescape' codec can't decode bytes in position 2-3: truncated \UXXXXXXXX escape
  os.chdir(path)
  Traceback (most recent call last):
    File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\code.py", line 90, in runcode
      exec(code, self.locals)
    File "", line 1, in
  NameError: name 'path' is not defined

When I enter the last line of the Code into the python console nothing happens. No layers are uploaded and I don't get error Messages.


texts:
I am trying to find a good source of roads data, the idea is mapping an entire state of a country in Brazil, so what I need is really large and good data.

Can I get all the roads of Google Earth in a specified area in a easy way so I can used it in ArcGIS?

If that is not possible, what source instead of OpenStreetMap can I use to do that task?

(OpenStreetMap has really good data, but there are some parts that are missing, some parts are really badly plotted and etc, and because of that it would take a long to adjust everything.


texts: SAGA grid tools not found on QGis 2.14.15

texts:
I am missing all grid tools from SAGA in QGis 2.14.15.

I know it's a question that has been asked before (there : SAGA grid tools missing in QGIS ) but unfortuately the answer given in this post did not fix my issue. I did have the GRASS folders not correctly configurated and fixed it, but it didn't have any influence on the SAGA tools (and I fail to see why it would have had any influence on it anyway since those are different plugins...).
I also tried to un-install and re-install QGis, but the tools still did not appear. I tried some other SAGA tools and they work perfectly. And QGis is in english so it's not an issue with the name of the algorithms either.

Any suggestion on how to fix this issue?

EDIT :
I am on Windows 10 and I have SAGA 2.3.2, which seems to be working fine on its own, but I'm not too familiar with the software on its own, I've always used it trhough QGis only... I can see the grid tools in SAGA's own interface by the way, they just don't appear in the QGis processing tool menu.


texts:
I was wondering if there is any good example of openlayers used in the HTML page with bootstrap3?

I searched around, and can only find bootstrap 2 examples, like this one: http://dev.openlayers.org/examples/bootstrap.html

I adapted my first map application from one of those examples. But when I tried to switch to bootstrap3, the layout was a mess. What used to be on my right half of the page now displays at the bottom beyond the normal viewport. Since my web page isn't that complex. I think it might be easier if I transfer the content to a bootstrap3 openlayers example.


texts: Click a polygon and have associated points be displayed

texts: I am mapping farms with the turnouts that deliver water to the fields. The farms (polygons) are parcels with a text field with 1-4 turnouts. Example- 2478,2479,4910,4911. The turnouts (points) have a unique ID in a number field. Ultimately, I want our users to be able to click on a farm and the turnouts that water that field are shown. Using ArcMap, how would I do that?


texts:
I'm trying to make a 5-year timelapse of monthly NDVI composites.
I apply a Hollstein mask, make composites and fill the masked area with mean NDVI pixels of the entire timeperiod.
Of course this method is questionable, and there are many ways, but I am mostly stuck on producing the timelapse video itself. I use the GIF production because I only use one band.
In GEE my composites look like this:

But when I produce a GIF timelapse, it looks like this:

I tried multiplying the NDVI by 100 to see whether that would induce more colour contrast, I tried to add more colours to the visualisation parameters, but the contrast remains the same.
Code: https://code.earthengine.google.com/dd50d25333a7ff8be05f76733bac6c3c


texts: How to consume an ArcGIS Server service from ArcObjects?

texts:
I have this end point to my ArcGIS Server:
http://myarcgisserver/arcgis/rest/services/prueba/GMRE/MapServer

I'm trying it with GisClient:

public void requestServer(string username, string password)
    {

        string url = "http://myarcgisserver/arcgis/rest/services/prueba/GMRE/MapServer";
        var gisServer = OpenConnection(url);
        MessageBox.Show(gisServer.ToString());

    }

private static IAGSServerConnection OpenConnection(string svcUrl)
    {
        //create a property set to hold connection properties
        var connectionProps = new PropertySet() as IPropertySet;
        //specify the URL for the server
        connectionProps.SetProperty("URL", svcUrl);

        //create a new ArcGIS Server connection factory
        var connectionFactory = (IAGSServerConnectionFactory2)new AGSServerConnectionFactory();
        var gisServer = connectionFactory.Open(connectionProps, 0);
        System.Runtime.InteropServices.Marshal.FinalReleaseComObject(connectionFactory);
        return gisServer;
    }

But I get a exception saying 'Invalid url'

My service in the browser:

How can I consume it from arcobjects?


texts:
My aim is to calculate the shortest distance between two point layers
One point layer contains Origin points, the other point layer contains Destination points.
Both layers contain more than 100 data points and that's why manual Point to point shortest path makes no sense to me.
Furthermore, the model contains a road layer on which all points are located and also another point layer.
This point layer contains all points where I have crossings of the streets.
Unfortunately, I have not found a way to calculate the shortest path for the Origin-Destination layers yet. The Model should match ID 1 of the origin layer to ID 1 of the destination layer and calculate the shortest path for it.
Is this anyhow possible in Qgis? I just found a bath process, but this just can run the whole layer and specific ids can not behead.
Continuing I want to save all point layers (crossings) that are on the shortest path. But this is not as important as the first goal.


texts: Loading layer from ArcSDE for SQL Server Database using PyQGIS?

texts:
I am use this code to load a layer from ArcSDE SQL Server database, but it only loads a table.

uri = QgsDataSourceURI()
uri.setConnection("localhost","", "arcsde", "sa", "****")
uri.setDataSource('dbo','Shelter',"geom")
vlayer = QgsVectorLayer(uri.uri(), "Shelter", "mssql")
QgsMapLayerRegistry.instance().addMapLayer(vlayer)

I try to use the follow code,it's can load layer from SQL server,but it can be edited,why?
uri = "MSSQL:server=localhost;database=arcsde;tables=dbo.SHELTER;trusted_connection=yes"
vlayer = QgsVectorLayer(uri, u"shelter", "ogr")
QgsMapLayerRegistry.instance().addMapLayer(vlayer)


texts: Running and scheduling processing jobs in QGIS 3?

texts:
Basically I want to write a script that would run on a daily basis on QGIS 3. This tutorial seems to have exactly what I want but the tutorial only works for QGIS 2.

QGIS Tutorials and Tips | Running and Scheduling QGIS Processing Jobs

How can I make the same thing in QGIS 3?


texts: Dynamically adjust the resolution of ol4 map to frame a feature within any size viewport?

texts:
I have a simple openlayers 4 map that calls an OSM base layer, plus a custom XYZ tile layer stored on my server:

var center = ol.proj.fromLonLat([-123.0, 44.95]);

var basemap = new ol.layer.Tile({
  source: new ol.source.OSM()
  });

var Salem = new ol.layer.Tile({
  source: new ol.source.XYZ({
    url: "http://ecotopia.today/cascadia/Tiles/PDX/{z}/{x}/{y}.png"}),
    minResolution: 0,
    maxResolution: 16000 });

var map = new ol.Map({
  layers: [ basemap, Salem ],
  target: 'map',
  view: new ol.View({
  center: center,
  resolution: 45 })
  });

By setting the center and resolution, I can make a particular feature from the custom layer fit perfectly within the viewport of a laptop with a 1366x768 screen (it's the city of Salem, Oregon, working example here). However, on a larger screen there will be lots more space around the feature, and on a smaller screen part of it will be cut, both of which I don't want. I realize I could set an extent or do a constrain center, but I want the user to be able to move freely through the whole map. It is only on the initial view that the feature should be correctly framed. Is there a way to set up a listener that gets the size of the user's viewport and then dynamically adjusts the resolution so as to correctly frame the feature when it first loads? Or is there any other strategy that achieves the same effect?

I made a jsfiddle here. Since my skills are limited, it would be great to see a complete piece of code or a working example.


texts: How to find coordinates from large CSV within multiple shapefile polygons

texts:
I have a CSV file of around 230 million rows, which contains a latitude column, a longitude column, a Date stamp, and an ID column. It look like this (I added the header in this example to help explain, the file does not have a header):

ID,                               Date,                     Lon,     lat
46d4089a713082d85452a2af64571644, 2016-11-30T12:57:11.000Z, 53.4529, -2.287
46d4089a713082d85452a2af64571644, 2016-11-30T12:57:26.000Z, 53.4521, -2.2859
46d4089a713082d85452a2af64571644, 2016-11-30T12:57:59.000Z, 53.4522, -2.2878
46d4089a713082d85452a2af64571644, 2016-11-30T12:59:01.000Z, 53.4547, -2.284
a6af7b30dc3ffa0ee7ecea02a2981b7d, 2016-11-30T13:03:01.000Z, 53.4457, -2.2774
693316c8b95e9fb07207400414714180, 2016-11-30T21:18:16.913Z, 53.2887, -2.6687

I also have two ESRI polygon shapefiles I created. What I want is a list of all the IDs that have at least one coordinate in both of the shapefiles.

Since the file is so large and could not be loaded into memory, my approach was to break the CSV file into 117 CSVs of 2 million rows. I then planned on using geopandas to read the smaller CSVs as a GeodataFrame, and use sjoin to find all the waypoints in shapefile 1.

I would then take the ID column of the waypoints within the shapefile as a list, take all waypoints in the large csv with an ID in the list, and do the same again finding the points in shapefile 2.

I tried to use spatial indexing to speed up the process, as it is a large file

I tried to test this using the first CSV of 2 million rows, using this code:

import pandas as pd
import geopandas
from geopandas.tools import sjoin
from shapely.geometry import *

waypoints = pd.read_csv(r'largefile_1.csv',sep=',', names=['TripId','lat','lon'],usecols=[0,2,3])

waypoints['geometry'] = waypoints.apply(lambda x: Point((float(x.lon), float(x.lat))), axis=1)

point = geopandas.GeoDataFrame(waypoints, geometry='geometry')
polygon  = geopandas.GeoDataFrame.from_file(r'Shapefile1.shp')
point.crs = polygon.crs

spatial_index = point.sindex
possible_matches_index = list(spatial_index.intersection(polygon.bounds))
possible_matches = point.iloc[possible_matches_index]
precise_matches = possible_matches[possible_matches.intersects(polygon)]

Is there a more efficient way of doing this, in another program or with different plugins?


texts: Quantifying similarity between two independent areas?

texts:
I'm assessing areas for suitability for species relocations, so essentially I'm looking to quantify how similar a possible relocation site is to a source site in several metrics. The species may travel 10km for example, so I have a point location for a source site, surrounded by a 10km radius. Within the radius I have values for every pixel (arbitrary size), for say altitude, or other numeric variables on a continuous scale that represent the habitat. At the relocation site, I also have a central point and 10km radius, but this may or may not overlap slightly with the source site radius. The pixelated areas within the two radii are not necessarily the same size however, as they can contain coastal areas, so different coastlines means different land area within the radii, and no data for sea area. EDIT: Spatial variability is not so important as habitat variables change more gradually over space, and not like a chess board/mosaic.

I'm looking for a way to quantify the overall similarity between the two areas. Essentially, taking the values of every pixel in each radius, I have two vectors of values of different respective lengths. I've found the Jaccard Similarity Index. EDIT: something like the Jaccard Index would be good but this is inappropriate here.

Would this be the best thing to do or is there a method more commonly used for this sort of spatial comparison? Preferably this would be done in R, otherwise QGIS.

https://www.statisticshowto.datasciencecentral.com/jaccard-index/


texts: Possible to convert only a single label in a shapefile to annotation?

texts: In ArcMap 10.3, is there a way to convert only a single label in a shapefile to annotation? I have a shapefile of 41 points (labeled 1 to 41). I only have one label that needs to be moved slightly to clarify which point it is associated with and converting to annotation would then make it where I could not easily change the placement properties etc for the entire shapefile. So, was wondering if anyone knew a way to move only the one label?


texts: Using Makefile when building QGIS plugin

texts:
I'm trying to learn how to build QGIS plugins (using Erik Westra's book Building Mapping Applications with QGIS), but am stuck on how to use make and Makefile.  I'm running this on Windows.

I've created the following Makefile:

PLUGINNAME = geometryInfo
PY_FILES = __init__.py geometryInfo.py
EXTRAS = metadata.txt icon.png
UI_FILES =
RESOURCE_FILES = resources.py

default: compile

compile: $(UI_FILES) $(RESOURCE_FILES)

%.py : %.qrc
    pyrcc4 -o $@ $<

%.py : %.ui
    pyuic4 -o $@ $<

deploy: compile
    mkdir -p $(HOME)/.qgis2/python/plugins/$(PLUGINNAME)
    cp -vf $(PY_FILES) $(HOME)/.qgis2/python/plugins/$(PLUGINNAME)
    cp -vf $(UI_FILES) $(HOME)/.qgis2/python/plugins/$(PLUGINNAME)
    cp -vf $(RESOURCE_FILES) $(HOME)/.qgis2/python/plugins/$(PLUGINNAME)
    cp -vf $(EXTRAS) $(HOME)/.qgis2/python/plugins/$(PLUGINNAME)

clean:
    rm $(UI_FILES) $(RESOURCE_FILES)

Westra's book says to run it by 'typing make deploy in the command line'.  I've worked out that I need to download a 'make' application for this, and have downloaded GnuWin, and added it to my path as C:\Program Files (x86)\GnuWin32\bin.

But this throws errors.

First, I get

pyrcc4 -o resources.py resources.qrc
process_begin: CreateProcess(NULL, pyrcc4 -o resources.py resources.qrc, ...) failed.
make (e=2): The system cannot find the file specified.
make: *** [resources.py] Error 2

I can overcome this if I create resources.py manually (but isn't that what the Makefile is supposed to do)?

Then, I get:

mkdir -p /.qgis2/python/plugins/geometryInfo
The syntax of the command is incorrect.
make: *** [deploy] Error 1

I can overcome that by switching the slashes around to \, and taking out -p (which Windows doesn't recognise).  Then it will create the directory (but only once - if I try a second time, it complains that the directory already exists).

The next line is cp -vf $(PY_FILES) $(HOME)/.qgis2/python/plugins/$(PLUGINNAME).  It doesn't like that ... perhaps I should change the code to Windows-like 'copy' instead of 'cp'?

But if I'm completely re-writing the Makefile and running half the processes manually, I can't help thinking I'm on the wrong track.

How should I be running the make deploy command when using Windows?


texts:
I am trying to integrate wind maps into an existing map and I would like to download the latest wind data for a local area like zyGrib does and then attach it into the map.

I am aware that there are sources where I could get this data for free but I have been unsuccessful in finding a source. Could someone point me to a free service that I could use to download the GRIB data preferably in bzipped format if possible.


texts: Why Stacking rasters return all values as NAs?

texts:
I have several rasters that I want to do some calculations ,basically calculating the moving average.

    saf
   class       : RasterStack
   dimensions  : 720, 1440, 1036800, 601  (nrow, ncol, ncell, nlayers)
   resolution  : 0.25, 0.25  (x, y)
   extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
  coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs

I am gratful to anyhelp


texts: How to get projection data (.prj) of raster image using GDAL/OGR?

texts:
I start with a boundary shape-file (think U.S. states) in lat/long that I convert to some projection, and then rasterize the result as a .tif.  What I need now are the raster coordinates of the original boundary points, as another shapefile.

With a .prf file for the raster projection, I can simply ogr2ogr the input shapefile with the .prj and get a shapefile where the points are now in raster-space; this is my initial inclination.  If there's a better way to do this, I'm open to suggestion.

tl;dr: I want to get the CRS/projection info for a raster coordinate system from a .tif file.

I would like to use command line tools.  ogr2ogr, gdal, etc.

Looking through the documentation of gdal/ogr2ogr, it seems I need WKT format, EPSG, etc for the projections.  Really, all I need to do is perform bilinear interpolation on the coordinates of the shapefile with the extents of the raster image.  Is there a quick tool to do something like this?

i.e., apply a function to every coordinate in a shapefile.


texts: QGIS Processing Script - Add field to output

texts:
My objective is to add a new calculated field to the output layer from a processing script.

This layer is created as follows:-

source = self.parameterAsSource(
            parameters,
            self.INPUT,
            context
        )

(sink, dest_id) = self.parameterAsSink(
            parameters,
            self.OUTPUT,
            context,
            source.fields(),
            source.wkbType(),
            source.sourceCrs()
        )

I am then adding a feature to the sink layer from source in an iteration loop

 for current, feature in enumerate(source.getFeatures()):
      value = 1
      #I want to add value to feature as an attribute, how??
      sink.addFeature(feature, QgsFeatureSink.FastInsert)

How do I add the value attribute to the sink layer?

Should it be done when declaring sink?


texts:
I would like to link a database with non-spatial data (it is just text and numbers) with a shapefile in QGIS. (Note: it is Microsoft SQL Server Management Studio) How do I do that?

If I am able to connect the database can I join both a specific table from SQL with the .shp file using a common ID?


texts: QGIS: Create polygon shapefile of ellipses from points using Geometry by expression tool

texts:
In QGIS 3.16.6 I need to create a polygon shapefile of ellipses from points.  I am doing the following:
Processing Toolbox -> Vector Geometry -> Geometry by expression

Input layer - points_in.shp

Output geometry type - polygon

Geometry expression - make_ellipse(transform(make_point("Lon-DD","Lat-DD"), 'EPSG:4269', 'EPSG:26913'), 635, 190, 90)

Modified geometry - ellipses_out.shp

Where "Lon-DD" and"Lat-DD" are fields in the point shapefile, the input CRS is GDD NAD83, the output CRS NAD83/UTM Zone 13N meters, and the semi-major and semi-minor axis are 635 meters and 190 meters and the azimuth is 90 degrees.
The Geometry by expression tool runs without complaints but the output shapefile is empty.
What am I doing wrong and how do I fix this?


texts: What are parameters and context in pyQgis processAlgorithm function?

texts:
I am working on a python scripts which are used in a Qgis module. In every script, there is the processAlgorithm function calling parameters and config. I can't find what are theses variables and where they are defined.
My purpose is to extract the code of theses scripts to avoid using qgis software.
Anyone knows how can I define my own parameters and config variables ?


texts:
I'm writing a plugin in QGIS and I'd like to generate a TIN surface using a vector layer with x,y and z (interpolation attribute).

The idea is to visualize the same result like in the Interpolation plugin but optimize it in my plugin, to facilitate users.

Does such a function exist?

Can I communicate with the Interpolation plugin to generate the result?

Thank you very much!

Ciao
Steffy


texts: Could not install Esri.ArcGISRuntime SDK package with Visual Studio 2015

texts: I am trying to install the ArcGIS Runtime SDK for .NET using Visual Studio 2015 (I see that it is not supported yet for VS 2017, so I installed 2015). I am follwing the ESRI documentation but having issues. I downloaded and installed the Runtime SDK (100.30) on my local and then opened my Visual Studio 2015. I am not seeing any of the Esri/ArcGIS templates under projects so I opened an existing solution; an AddIn, and tried using nuGet to download the package (I searched for ArcGIS Runitme SDK). I got the following error: "Could not install package Esri.ArcGISRuntime 100.30. You are trying to install this package into a project that targets '.NETFramework, Version=v3.5' but the package does not contain any assembly references or content files that are compatible with that framework." What am I doing wrong here?


texts: Updating back to GeoServer/Oracle from OpenLayers?

texts:
I have followed the example to setup a WFS Transaction successfully:

http://dev.openlayers.org/examples/wfs-protocol-transactions.html

My environment is GeoServer 2.6.1/OpenLayers 2.12.1. I do get the 4 toolbar icons on the upper right hand corner. (draw, modify, delete, save).  I used a layer that is connected to an oracle table in geoserver as the editable layer.

Now, when I click on the "save" button, it does not add a record to that table.  Of course, I was not expecting it to because it does not even know what key value I want to put in (as customerid) when it creates a new record with feature data; nor does it know which record to update the geometry for it its an existing key value (customerid).  The editable layer I created based on the above link is:

// Customers Editable VECTOR Layer
            var saveStrategy = new OpenLayers.Strategy.Save();

            customers = new OpenLayers.Layer.Vector("Customers", {
                strategies: [new OpenLayers.Strategy.BBOX(), saveStrategy],
                projection: new OpenLayers.Projection("EPSG:2236"),
                protocol: new OpenLayers.Protocol.WFS({
                    version: "1.1.0",
                    srsName: "EPSG:2236",
                    url: "http://server01:8085/geoserver/sf/wfs",
                    featurePrefix: "sf",
                    featureType: "customers",
                    featureNS: "http://a.org/layers01",
                    geometryName: "GEOMETRY"
                })
            });

The oracle table itself has the following fields:

-customerId
-geometry
-date_updated

How do I perform inserts through openlayers/geoserver/all the way to the oracle table? (To included key value and geometry, etc)  I would like the above example to actually save the geometry back to a specific record in the table.  I am not sure how to proceed on this.

I am able to do inserts!  It creates a new customerId with correct geometry.  However, it increments the primary key by one, but ignores the value I place in added_feature.attributes.customerId for a new customerId key.  And, I cannot get it to update at all.

I did set the customerId in the FeatureAdded method and set the state to OpenLayers.State.UPDATE.  I also checked that the customerId I placed in the added_feature.attributes.customerId variable does exist in the table.  This is the error I get when I try an update:

java.lang.RuntimeException: Parsing failed for FeatureId:
java.lang.NullPointerException Parsing failed for FeatureId:
java.lang.NullPointerException


texts: Date format conversion brings an unexpected empty space between the date parts

texts:
I always used the formula (below) in Python Calculator to convert the date format from Julian to Gregorian format. But,now a weird behaviour is happening and I am trying to figure out what is wrong in the processing. I am using QGIS 3.12 and after processing the formula no errors rises, but when I opened the new field in the table a little problem is detected. There is a empty space (4 spaces) between Year and Month parts of the date and I'm not able to identify this error and remove It. The new format date generated is "2016-***[space]***01-12 23:45" and the correct and expected result should be "2016-01-12 23:45".
Following the formula I'd used:
timestampparts=<Date_Time_Local>.split(' ')
dateparts=timestampparts[0].split('-')
timeparts=timestampparts[1].split(':')
newdate='{0[2]}-{0[0]}-{0[1]} {1[0]}:{1[1]}'.format(dateparts, timeparts)
value=newdate

And the output of conversion process:

There is something I can do to fix It?


texts:
I'm trying to intersect a polyline and polygon and write out the resulting polyline shapefile using ogr2ogr

Here is the command I've been trying:

    ogr2ogr -dialect SQLITE -sql "SELECT ST_Intersection(A.geometry, B.geometry) AS geometry, A.*, B.* FROM polygon A, polyline B WHERE ST_Intersects(A.geometry, B.geometry)" test.shp input.vrt

Here are the contents of the input.vrt file:

<OGRVRTDataSource>
    <OGRVRTLayer name="polygon">
        <SrcDataSource>polygon.shp</SrcDataSource>
    </OGRVRTLayer>
    <OGRVRTLayer name="polyline">
        <SrcDataSource>polyline.shp</SrcDataSource>
    </OGRVRTLayer>
</OGRVRTDataSource>

The resulting output is a .dbf file called "test" with all the desired fields, but no geometries. There is no accompanying shapefile with the same name. I've tried modifications of the command, particularly switching the order of the intersection. I know the features do in fact intersect, as I've done the intersection in ArcGIS and it works, producing 60+ features.

caveat - the polyline has ~1500 features whereas the polygon has been dissolved, and only contains one feature. I don't think this should matter.

I know there are other libraries to accomplish this, but I have an existing backend system built on gdal/ogr and I am looking to reduce dependencies to other libraries.

UPDATE: Debugging based on comments from @User30184.Here are is the output for one feature with the following ogrinfo call:

GEOS error: TopologyException: Input geom 0 is invalid: Self-intersection at or near point 97.562017383948856 16.504292693814882 at 97.562017383948856 16.504292693814882
OGRFeature(SELECT):1
geometry (String) = (null)
DN (Real) = 1
osm_id (String) = 27167387
code (Integer) = 5112
fclass (String) = trunk
name (String) = (null)
ref (String) = 408
oneway (String) = F
maxspeed (Integer) = 0
layer (Real) = 0
bridge (String) = F
tunnel (String) = F

Looks like the geometry field is not being computed. Any ideas?


texts: Automatically loading PostGIS relations in QGIS

texts:
I am working with a PostgreSQL/PostGIS (version 12.3/3.0.1) database and I am using QGIS (version 3.10) to visualize and work with my data. In the database I have two schemas:

schema for spatial data tables, e. g. gebaeude
schema for documents, e. g. dokument, dok_gebaeude (without spatial information)

I set a m:n relation between gebaeude and dokument using the table dok_gebaeude.
My tables are:

dokument: id, bezeichnung, dateityp, pfad
dok_gebaeude: id, fk_dokument, fk_gebaeude --> fk_dokument references dokument(id) and fk_gebaeude references gebaeude(id)
gebaeude: id, geom, nummer, nutzung

In QGIS I loaded the relationship with the >project properties>relation manager
and I can use it to connect documents (stored in table dokument) with gebaeude-objects. All work fine. BUT...
When I open a new QGIS project, I always have to discover and add the relationship again in my >project properties>relation manager. In addition, I have to readjust each time the attribute form of the table gebaeude (I´m using the drag and drop designer and implement there the m:n relationship).
Is there any possibility to automatically discover and add the relationship by starting a QGIS project?


texts:
I want to project my data into another coordinate system.

I tried using the Project tool but when I opened the attribute table after ArcGIS said it had finished it was empty.

I also tried using the merge function but it says the data has to be in the same coordinate system - which I can't do because ArcGIS loses the attributes needed to merge.


texts:
I can't see my home address on Google Maps as it is a new neighbourhood.

How can I get my home address added into Google Maps?


texts: How do I iterate through every cell in a continuous raster?

texts:
See this link for more details.

The Problem:

I want to loop through a continuous raster (one that has no attribute table), cell by cell, and get the value of the cell. I want to take those values and run conditionals on them, emulating the map algebra steps detailed below without actually using the raster calculator.

Per request of comments below, I have added details providing background to the problem and justifying the need to implement a method as such in the section below called "The analysis needed:".

The analysis proposed below, while being relevant to my problem by providing background,does not need to be implemented in an answer. The scope of the question only pertains to iterating through a continuous raster to get/set the cell values.

The analysis needed:

If ANY of the following conditions are fulfilled, give output cell a value of 1. Only give output cell a value of 0 if none of the conditions are fulfilled.

Condition 1: If cell value is greater than top and bottom cells, give value of 1:

Con("raster" > FocalStatistics("raster", NbrIrregular("C:\filepath\kernel_file.txt"), "MAXIMUM"), 1, 0)

Where kernel file looks like this:

3 3
0 1 0
0 0 0
0 1 0

Condition 2: If cell value is greater than left and right cells, give value of 1:

Con("raster" > FocalStatistics("raster", NbrIrregular("C:\filepath\kernel_file.txt"), "MAXIMUM"), 1, 0)

Where kernel file looks like this:

3 3
0 0 0
1 0 1
0 0 0

Condition 3: If cell value is greater than topleft and bottomright cells, give value of 1:

Con("raster" > FocalStatistics("raster", NbrIrregular("C:\filepath\kernel_file.txt"), "MAXIMUM"), 1, 0)

Where kernel file looks like this:

3 3
1 0 0
0 0 0
0 0 1

Condition 4: If cell value is greater than bottomleft and topright cells, give value of 1:

Con("raster" > FocalStatistics("raster", NbrIrregular("C:\filepath\kernel_file.txt"), "MAXIMUM"), 1, 0)

Where kernel file looks like this:

3 3
0 0 1
0 0 0
1 0 0

Condition 5: If any one of the adjacent cells has a value EQUAL to the center cell, give the output raster a value of 1 (using focal variety with two nearest neighborhood calculations)

Why not use map algebra?

It has been noted below that my problem could be solved using map algebra but as seen above this is a grand total of six raster calculations, plus one to combine all of the rasters created together. It seems to me that it is much more efficient to go cell-by-cell and do all of the comparisons at once in each cell instead of looping through each individually seven times and utilizing quite a bit of memory to create seven rasters.

How should the problem be attacked?

The link above advises to use IPixelBlock interface, however it is unclear from ESRI documentation whether you are actually accessing a single cell value itself through IPixelBlock, or if you are accessing multiple cell values from the size of the IPixelBlock you set. A good answer should suggest a method for accessing the cell values of a continuous raster and provide an explanation of the methodology behind the code, if not apparently obvious.

In summary:

What is the best method to loop through every cell in a CONTINUOUS raster (which has no attribute table) to access its cell values?

A good answer need not implement the analysis steps described above, it needs only to provide a methodology to access cell values of a raster.


texts: How can I filter a MultiPoint Z / xyz points against a 2d polygon

texts:
I would like to filter point data from LIDAR scans (3d) by parcel plot lines (2d).

I began with loading the lidar data into microsoft sql and performing a STIntersects between the MultiPoint and the Polygon. This provides a result, but loses it's 3rd dimension.

I then have tried loading the lidar .shp file into the fiona python library, calling .filter() to narrow down the features, and traversing points checking .within() against the plot line polygon. This seems to work but is super slow.

Next I loaded the lidar .shp into the rtree python library. However there doesn't seem to be the best 3d support. The methods that were obvious to me were limited by the inability to get a .bound on the shape that included the 3rd dimension. So when I added the shape to the index is complained the bound was 2d and the shape 3d.

Finally I attempted to load each point manually into the rtree index so that I could correct the 3rd dimension manually on the insert.

I have not validated that rtree will even work once I have a complete index.

Has anyone figured out a non brute-force method?


texts:
I'm trying to create 19 maps representing a highway construction. I need each map to show only the featured geometry and a  its classification
I´ve tried using expressions like this

"situ_Pista Nova" = 'Executado' AND @atlas_featureNUMBER = "Atlas_id"

The Atlas_id is a new column created to match the atlas_featurenumber. The atlas coverage layer is not the same as the layer I want to show

but every feature just falls and the rule else.


texts: I am looking for a way to create text on a MapInfo layer using Polylines.  So on a MapInfo layer, I need to create a line dataset, but I am after forming words - labelling items is not going to work for me.  Are there any tools that can help do this?  I have tried doing it freehand using a polyline, and a polyline with autonode mode selected, but frankly it looks rubbish


texts:
I need to get a GeoJSON geometry of this, then use a buffer for the TURF library.

Get get data from GeoJSON turf but would not admit it in any way.
any ideas ?

This is my code:

var feature_buff = select2.getFeatures();
var geojson  = new ol.format.GeoJSON();

for (var i=0;i<feature_buff.getLength();i++){
json.push(geojson.writeFeature(feature_buff.item(i)));
var obj =geojson.writeFeatures(feature_buff.item(i))

var geom = geojson.readGeometry(obj);

var buffer = turf.buffer(geom, 2000, 'kilometers')

this code fail. and this is the error:

Uncaught TypeError: (0 , ls[b.type]) is not a function


texts: CartoDB - Multipolygon not showing intersections

texts:
I'm having an issue with Multipolygons on CartoDB.

I have a very simple scenario: A multi polygon row on a table that needs to be displayed on the map. Fact is that when more than two polygons somehow colide, the MapView only display the portions that are not common between the polygons. It's very simple to demonstrate:

(Those are the original 3 polygons, manually rendered using Google Maps API)

(This is the CartoDB rendering. The overlapping area is missing!)

The craziest thing is that the renderer shows the portion covered by the three polygons. You can see the small area of Botafogo, this area is covered by the three polygons and is correctly shown on the CartoDB rendering.

Any idea why this is happening and how to fix it? I couldn't find anything useful on the documentation. The only workaround I could come up with is by doing some operations using Postgis (ST_Polygonize) but if I do this I'll lose the filtering options, since I'll be using a custom query. Also, to me, it makes no sense why CartoDB is not showing the overlapping area.


texts: How to calculate area of respective pixel?

texts: I wants to create NDVI Map, In which it shows an area of respected pixel value lies in map.I have NDVI value within -0.6 to 0.7. Suppose,I wants to calculate area of 0.5 NDVI value with scale of map.How it can be done? How to convert pixel value into area value?


texts: Creating JPEG images of geoprocessing output using ArcObjects?

texts:
I'm running batch processing on several datasets using a standalone EXE written in .NET. There are many output files for each dataset (>25 featureclasses/rasters). For every dataset I would like to grab a few specific outputs and create JPEG images of them for simple quality control. Ideally I would combine a few of the featureclasses together and create an image of the combination as if I were exporting an image of several layers as you can do using the Arc UI. I imagine that there will be a few steps involved in accomplishing this - can someone point me in the right direction?

For the sake of simplicity, you can assume that I have the full path to the files I want, which are stored in a FGDB on a drive. I can also derive the approximate extent of these featureclasses from the largest of them, so I can set a view extent.

For example:

I have just run the extension and created a series of output files which are stored in a FGDB. These featureclasses/rasters all cover the same discrete area and are projected.

Outputs (making these up for simplicity):

1) DEM of the area (represents the full extent).

2) Polyline road network.

3) Polygon building footprints.

4) Point featureclass representing manhole covers.

What to accomplish:

Combine all of the files in the correct order (featureclasses overtop of the DEM) and save a JPEG file.

EDIT:

My extension runs as an EXE (not a toolbar) and doesn't actually open the Arc UI - it merely consumes the licenses and does the processing in the background. So, do I have access to ActiveView if I haven't "launched" the UI? Essentially, I don't have a TOC to populate with layers.


texts:
Does anyone know where I can find more information about the new Lightweight Runtime that will be included with ArcGIS 10.1?  This was discussed at the recent DevSummit and ESRI Federal User Conference.

I am looking for any general information about the product, but more specifically I am trying to find out if there will be any sort of licensing fees to use it.

More information can be found here:
http://esriaustralia.wordpress.com/2011/02/17/the-new-arcgis-runtime/

Thanks.


texts: Extracting the Number of Pixels within Buffers around Various Points in Google Earth Engine?

texts:
I have a Feature Collection of various points, each of which I've drawn a 2.5km buffer around using the following:

    var bufferBy = function(size) {
      return function(feature) {
        return feature.buffer(size);
      };
    };

After running a supervised SVM classification with only two possible output classes, either 0 or 1, I have determined the total number of pixels in both classes within the overlapping area of the buffers using the following:

    var sigP = classified_SVM_Park.reduceRegion({
      reducer: ee.Reducer.countEvery().group({
        groupName: 'signature',
      }),
      geometry: pBuffer.geometry(),
      scale: 30,
      maxPixels: 1e8
    });

How would I go about extracting the number of pixels assigned 0 and 1 within each buffer--in other words, the number of pixels within the buffer around each individual point?

Link to the full code is here: https://code.earthengine.google.com/54d0806b66fd31476a7a314f2da11df6


texts: Automate creating an image pyramid store with Geoserver?

texts: Geoserver has a REST API to automate Geoserver commands, however I do not see any examples of how to automatically create an Image Pyramid store and associated layer. Is that possible? Can anyone provide me sample code?


texts:
I've been trying to list a ee.ImageCollection() in GEE but by using time indexing and I came up with this solution:
var geometry = ee.Geometry.Polygon([[-86.036445612425396,
11.175779794553685 ],
                                [ -86.036445612425396, 9.362885788210679 ],
                                [ -83.883633979893077, 9.362885788210679 ],
                                [ -83.883633979893077, 11.175779794553685 ],
                                [ -86.036445612425396, 11.175779794553685 ]]);

var start_date = '1979-01-01'
var end_date = '2019-12-31'
var Landsat = ee.ImageCollection('LANDSAT/LE07/C01/T1_TOA')
                    .filterBounds(geometry)
                    .filterDate(start_date,end_date);
var listOfImages = Landsat.toList(Landsat.size())

By doing this, I can index each image like this:
var firstImage = listOfImages.get(0)
Map.addLayer(ee.Image(firstImage))

However, this method indexes one scene of the collection each time. I mean, firstImage is the Landsat scene LANDSAT/LE07/C01/T1_TOA/LE07_015052_20000221.
I would like to know how to mosaic an image collection and index them by time. For example, getting the whole mosaic as a result of running this instruction: Map.addLayer(Landsat.median().select('B4')), in general, being able to get Map.addLayer(Landsat.get(nth-image).select('B4')).


texts: How to batch "gdal_merge" on Win7 in one command line command?

texts:
I want to create a mosaic of all tif files in a directory on windows 7.

How could I create a one-line command in the windows command line to loop through all tifs and merge them all together using gdal_merge.py?

I tried this;

for %i in (*.tif) do gdal_merge.py -of GTiff %i mosaic.tif

but it gives me an error about mosaic not being a supported dataset name??


texts: Leaflet marker icon does not change correspondingly

texts:
I am building an ASP.NET Core application using Leaflet.js library. On the map I would like to show containers(used for transferring building materials) as markers. These containers can have different statuses: take-able, old, normal

What I want to do is to show different marker-icon on the map for each status, correspondingly.
Let's see how far I could get:

function initMap(containers) {
    window.map = L.map('map', {
        center: [46.6726, 21.0856],
        zoom: 13,
        layers: [
            buildLayerGroupWithContainers(containers)]
    });

    L.tileLayer('http://{s}.tile.openstreetmap.de/tiles/osmde/{z}/{x}/{y}.png', {
        attribution: '&copy; OpenStreetMap contributors'
    }).addTo(map);
}

function buildLayerGroupWithContainers(containersInRawJSON) {
    var containerMarkers = [];
    var containersInJSON = JSON.parse(containersInRawJSON.replace(/&quot;/g, '"'));

    for (index = 0; index < containersInJSON.length; ++index) {
        var container = containersInJSON[index];
        if (container.IsTakeable) {
            containerMarkers[index] = buildMarker(container, takeableContainerIcon);
            alert("this is takeable: " + index);
        } else if (container.IsOld) {
            containerMarkers[index] = buildMarker(container, oldContainerIcon);
            alert("this is old: " +index);
        } else {
            containerMarkers[index] = buildMarker(container, containerIcon);
            alert("this is normal: " + index);
        }
    }

    return L.layerGroup(containerMarkers);
}

function buildMarker(container, cIcon) {
    return L.marker([container.Location.Latitude, container.Location.Longitude], { icon: cIcon })
        .bindPopup(
        "Konténer info:
" +
        "Id: " + container.CompoundId + "
" +
        "Tipus: " + container.Type + "
" +
        "Lerakás dátuma: " + container.LayDownDate + "
" +
        "Lerakta: " + container.LayedDownBy + "
" +
        "Tel.: " + container.PhoneNumber + "
" +
        "Megjegyzés: " + container.Others + "
" +
        "Öreg-e: " + container.IsOld + "
");
}

First of all you can see that there is a containersInJSON array. These containers come from the ASP.NET application, and the properties are set in C# code. So every property is set up properly.
In the for-loop firstly I check whether the container is takeable. If so, I build a marker for it with the corresponding container-icon. If it is not take-able, I check of it is old. If so I set the appropriate container-icon. If not, I build a marker with normal icon.
The setup of my markers are the followings:

var containerIcon = L.icon({
    iconUrl: '../images/markers/container.png',
    iconSize: [42, 42],
});

var takeableContainerIcon = L.icon({
    iconUrl: '../images/markers/container-with-green-flag.png',
    iconSize: [42, 42],
});

var oldContainerIcon = L.icon({
    iconUrl: '../images/markers/container-with-red-flag.png',
    iconSize: [42, 42],
});

Well, everything seems fine. But on the map, al the markers have the "takeableContainerIcon" marker-icon. At this point I was wondering of maybe my containers from c# are wrongly set. So that's why I put some alert methods after each marker-building, so for debugging purposes. The following line:

alert("this is old: " +index);

..is invoked once, which is the expected result. I have namely 12 containers, 11 are takeable and only one is old. But despite it, I can see only markers with takeableContainerIcon, and even the old marker is shown with takeableContainerIcon which is really interesting.

I could not get further.


texts: Reusing style and setting feature text label in OpenLayers?

texts:
In OpenLayers 5 I have a vector layer which can display 16K features in certain portions of the map. Each feature has a text label that displays a unique name.

If I create a new style for each feature so I can dynamically set the Text string I run into performance problems where the browser rendering crashes/hangs.

Alternatively, if I create a style cache and reuse the styles when the text label is set for one feature all the other feature text labels change.

Is there a way to reuse style object instances and have the text label be different for each feature?


texts: Create a center line between 2 polyline - QGIS

texts:
I want to create a center line for a river (polygon). I thought about simply creating a new line manually, but it would be a long and imperfect process. If there is a way to automate this process, that would be much better.

Can anyone help?


texts:
Running on Windows; 8gb memory, 500gb HD, i5 cpu; Raster file: 250mb

I am having a problem with zonal statistics advancing past 17%. I am trying to run statistics on a crop grown in certain counties. In the image below, you can see the counties of Washington state in the vector buffer. The raster data is dispersed throughout the state.

I have run zonal statistics on much larger raster files. The raster image comes from CDL NASS portal, which contains many more pixel values other than the one I have isolated. But it allows for a researcher to isolate the variable (crop) they are interested in and download that image without the rest of the data.

Any ideas?


texts: How to minimize size of Spatialite file by reducing excess projections, logging etc

texts:
I have a spatialite db with only 1 table. It contains a million simple polygons with 8 points each and is about 350 MB. It is intended for read-only use in a mobile app, so I want to reduce the size as much as possible. To that end I don't need 500 coord systems (described in 3 different tables!), activity logs etc.

I've tried dropping the excess coord systems with a delete command but it gives a foreign key constraint error.

delete from spatial_ref_sys where srid <> 4326

Simply saving as SQLite drops the size to 86 MB so that is the size of the attributes, primary key etc.
How do I trim ALL of the fat?


texts: How can I extract all distinct values from a PostGIS Raster?

texts:
I am working with an integer raster in PostGIS 2.1 and I would like to extract all distinct values in that raster. I was hoping to be able to do something along the lines of the following query, but I did not find a solution that works so far.

SELECT DISTINCT ST_Value(rast) FROM myrastertable


texts:
How do I find out if there are points in a multi-geometry feature as I suspect it could cause the error, especially the second one? When using ms_simplify on ms_simplify(drain[drain$H1 %in% unique(seq),], keep = 0.1), I get

Error: ms_simplify can only operate on (multi)polygons and (multi)linestrings

If I remove the %in% and subset with >, =, < signs it works. Also, the geometry column is labelled Shape for some reason and not geom as usually and is loaded from a .gpkg file drain st_read("water.gpkg"). Does anyone have ideas why these problems arise?

My code is quite long and I assume the error is specific to the data that I have. Therefore, no reproducible example this time :/


texts: Goodness of fit of variogram model using geoR

texts:
How can I find a measure of goodness of fit (let's say an R2) of a model fitted to a variogram by the function variofit using the geoR package?

The only value I have is the "minimised weighted sum of squares". Is it possible to obtain the R2 from this value? If so, how do I do that?

And, also, how do I find this same value in the case of a manual fit (using the eyefit function or manually specifying the parameters in the lines.variomodel function)?


texts: How to calculate position in a picture based on latitude/longitude

texts:
I have got a picture on which a user manually marks some points. Let's say for simplicity's sake these points always create a 4 point convex polygon. For each one of these points I know latitude and longitude.

Now the user enters latitude/ longitude for the 5th point and I have to show it on the same picture. What is the best way to calculate its screen position?

I have tried the following algorithm:

1) to calculate a geometrical center of the polygon on screen

2) to find an average longitude for the leftmost and  the rightmost
points and an average latitude for the topmost and the lowest points

3) to place a point at the same relation to the geometrical center
as its' geocoordinates are to this average point.

I have a test image, where I know exact geocoordinates for each point on the 10px grid, and I can say that this algorithm doesn't work.

What would be the better way?


texts: ArcGIS 10 Rest services regulary stop working

texts:
I have a problem with a new installation of ArcGIS Server.
The installation is ArcGIS 10 on Windows Server 2008r2 with IIS 7.
I have followed this guide to configure the IIS ArcGISServicesAppPool: http://support.esri.com/en/knowledgebase/techarticles/detail/36887
I did this because my cache wouldn't work without it.

The problem I have is that regularly the REST services stop working. The url to /arcgis/rest/services still lists all my services but when I want to open one, I get an error that the service is not available.
Also when I want to go to the ArcGIS/rest/admin/ I cannot log in.

I can solve this problem by recycling the ArcGISServicesAppPool in IIS. Then I can log in to the rest admin. But then I still need to clear the ArcGIS rest cache before my services work again.

One thing I've changed was to set the Idle time-out in IIS for the ArcGISServicesAppPool to 0, because I noticed that every time after that time out, my rest services didn't work anymore (same problem). So after this change, the services will run much longer, but still after a day or so, I need to recycle and clear the cache.

I've done some more checks and I notice that every time the ArcGISServiceAppPool automaticly recycles (via settings in IIS) the services stop being available.
The solution again is to manually recycle the application pool and then clear the rest cache again.

I've solved my problem by adding ArcGISWebservices user to the Administrators group. Don't find it a clean solution, but it's a solution. Can there be other problems with this? Security?


texts: Python: 'Select by Attribute' in ArcMap 10

texts:
I would like to get the search the field value

'SURFACE LOCATION%' in column 'RefName' in layer '*.dwg Annotation'

This is what I have so far...

# Find the Dataframe
df = arcpy.mapping.ListDataFrames(mxd, "Project Area")[0]

#Find the table named "*.dwg Annotation"
Table = arcpy.mapping.ListLayers(mxd, "*.dwg Annotation", df)[0]

whereClause = """RefName" LIKE 'SURFACE LOCATION%'""" % RefName

Anno = arcpy.SelectLayerByAttribute_management(AnnoTable, "NEW_SELECTION", whereClause)

Right now "Table" outputs the group layer name and the layer name.

However The link below shows that the the file path of the layer to work with SelectByAttribute.

http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//001700000071000000

Thoughts on how to get the file path of the layer?

and many thanks.


texts: Forest portals in the context of GIS and maps

texts:
Do you know any portals connected with forests and forestry in the context of GIS and maps?
The best for me will be portals where the data is served via WMS, WFS or WCS services.
The only portals i know are from Poland:

Regional Directorate of State Forests in Krakov

Regional Directorate of State Forests in Zielona Góra

I found a new one- it is the map of nature conservation in the forests of north-eastern Polish


texts: Finding a vertical kilometer between 2 points with specified distance in QGIS

texts:
I'm new to GIS software.
I want to find 2 points in selected area (lets say Czech Republic) with distance 5km or less and altitude difference 1000m or more.
I know that I can use GMTED datasheet with terrain data but then I'm little bit lost.

If it's a stupid question or it's been answered before, sorry for that.


texts: Google Earth - Force KML polygon on top of image overlay

texts:
In google earth:
Step 1.  Create polygon layer.
Step 2.  Overlay image file (mine is in jpg format).
Step 3.  Move polygon layer above and below image file in the "My Places" menu.
I cannot figure how to get the polygon layer to overlay on top of the image file.


texts: Is it possible to run two simultaneous search cursors?

texts:
I would like to search all feature classes in a geodatabase using a search cursor but I want to use two separate where clauses. I'm gathering information about a particular field and I need to know how it is different for each where clause.

The first where clause would be:

GIS_STATUS = 'Superseded'

and the other would be:

GIS_STATUS <> 'Superseded'

My suspicion is that I need to traverse the entire geodatabase twice using two loops. However, I wonder if there a way I can use two simultaneous search cursors? Parallel search cursors?

I'm using ArcGIS Desktop 10.4.1.

There seemed to be a bit of confusion about what I was trying to do, so I'll show what I had in mind.

Firstly, I define two sets.  One to hold the values in the METADATA field where GIS_STATUS = 'Superseded' and another where GIS_STATUS <> 'Superseded'.

superSet = set()
notSuperSet = set()

Then, I loop through all feature classes twice.  One for each where clause:

for dirpath, dirnames, filenames in arcpy.da.Walk(env.workspace, datatype = 'FeatureClass'):
    for filename in filenames:
        with arcpy.da.SearchCursor(filename, "METADATA", """GIS_STATUS = 'Superseded'""") as sCur:
            for row in sCur:
                superSet.add(row[0])

and

for dirpath, dirnames, filenames in arcpy.da.Walk(env.workspace, datatype = 'FeatureClass'):
    for filename in filenames:
        with arcpy.da.SearchCursor(filename, "METADATA", """GIS_STATUS <> 'Superseded'""") as sCur:
            for row in sCur:
                notSuperSet.add(row[0])

Then, I can use the set methods to determine how they are different.

superSet.difference(notSuperSet)

So, my question is: is this the best way to do this?  Do I need two loops?


texts: Creating custom offline maps for Collector on iPad without ArcSDE?

texts: Esri has informed me that you can only create custom offline maps for collector units if you have an ArcSDE set up with your GIS server.  They said ArcSDE is the only way to publish maps up to the ArcGIS Online cloud with "sync" and "offline" modes (which is required to pull the maps from the cloud to the collector unit).  I am trying find any workarounds so I do not have to purchase ArcSDE.


texts: Calculate Field based on field value in a second table

texts:
I have a table in a MSSQL Server Enterprise SDE database.  I need to join a table from a MS SQLExpress Workgroup SDE database and then calculate 11 fields in the Enterprise table based on values in the Workgroup Table then remove the join.

I have tried to do this in model builder add it takes forever to run.  I thought of selecting only those records that don't match before performing the calculate, but then thought maybe there were other ways to tackle this problems.

So now I am thinking that there are three alternatives that will do what I need in a quicker manor.  My question is can anyone give insight on which would be faster or easier to do.

Option 1
     Create a python code in ArcGIS 10.1 to perform the join select the mismatched records and perform the calculation then move on to the next field, when all done remove the join.

Option 2
     Use FME Desktop to perform the same steps outlined above.

Option 3
     Use SQL and perform the same steps above in Database

An input or suggestions would be appreciated.


texts: How to create a semi-transparent map collar in QGIS composer?

texts:
In QGIS 2.14, how can I create a semi-transparent "map collar" (white space around the outside of the page) in Composer?

I imagine the map collar would be some sort of graphic shape (not a feature) that is maybe 1 inch wide and covers the current data map.  Think of it like a donut:  white space that is semi-transparent white, while the interior is transparent.  The end goal is to create a "masking" effect that would focus the users eye on the interior of the map, but still allow them to see features around the edge of the map (and "hidden" under the collar).

A couple of ideas I had but cannot figure out a way to implement:

Buffer the map frame (to the outside only)
Create a graphic shape and "cut out" the interior


texts: Split linestring by selecting two points in OpenLayers?

texts:
I have a project where people can drag and drop a gpx file. This results in a multilinestring. People should be able to take a piece of their route and submit that. I was thinking of letting the user click two points and then taking the points in between.
First problem was making the multilinestring into a linestring because people can't select something when there is a gap betwee two linestrings.

I wrote a little function to take each coordinate from each linestring in the multilinestring.
I add this coordinate to a single linestring (lineSource) and as a point to another source (pointSource).

dragAndDropInteraction.on('addfeatures', function (event) {
if (event.features.length > 0) {
    let feat = event.features[0];
    if(feat.getGeometry().getType() === 'MultiLineString')
    {
        let lineStringArray = [];
        let lineStrings = feat.getGeometry().getLineStrings();
        for(let lineString = 0; lineString < lineStrings.length; lineString++){
            let linecoor = lineStrings[lineString].getCoordinates();
            for(let coor = 0; coor < linecoor.length; coor++){
                lineStringArray.push(linecoor[coor])
                pointSource.addFeature(new ol.Feature({
                    geometry: new ol.geom.Point(linecoor[coor])
                }))
            }
        }
        lineSource.addFeature(new ol.Feature({
            geometry: new ol.geom.LineString(lineStringArray)
        }));
        map.getView().fit(lineSource.getExtent());
    }
}
});

This is the result

So far, so good.
I also have a selectionLayer where the user can click on two points of the pointSource, with snapping (draw function). After they click two points, I was planning on looping over all the coordinates of the linestring and compare that coordinate with the coordinates of the two points.

As a test, I tried with the coordinates of one point

draw.on('drawend', function (e) {
  var coor = e.feature.getGeometry().getCoordinates();
  console.log(lineSource.getFeatures()[0].getGeometry().intersectsCoordinate(coor));
});

this always returns false.

I use snapping to add a point to my selectionLayer. How is that the feature from lineSource never contains the coordinate from the point in the selectionLayer?
No point in looping over the coordinates from the linestring if the linestring doesn't intersect the selectionLayer
What is the correct way to select points in a linestring? As far as I know, you can't add an id to a point/coordinate in a linestring


texts: How to set up a maximal zoom for polygons in Leaflet?

texts:
How to set up a maximal zoom for some layer, for example state polygons? I need that my US polygons which represent states have maxZoom 6, and when zoomed higher than 6, that US counties show up instead of states. Is that possible in Leaflet? All states are defined like Leaflet polygons.

Thanks in advance.


texts:
I have a table of points that looks something like this:

year,   latitude, longitude
2018,   42.02,    -87.68
2018,   41.91,    -87.78
2018,   41.64,    -87.68
2017,   41.63,    -87.56
2017,   41.89,    -87.67
2017,   41.67,    -87.80

Which I used Import XY to import from CSV as a set of points.

Now I want to perform some operation on it, but looping over the year attribute.

So far I got something like this:

import arcpy
imported_points = "imported_points.shp"

for yr in (2017, 2018):
  new_file = "poly_" + yr + ".shp"
  arcpy.SelectLayerByAttribute_management(imported_points ,
                                          "NEW_SELECTION",
                                          "[year] = " + yr)
  # Some operation, using this as an example
  arcpy.DirectionalDistribution_stats(imported_points , new_file, "2_STANDARD_DEVIATIONS", "#", "#")

However, when I try this I get the following error:

ExecuteError: Failed to execute. Parameters are not valid.
The value cannot be a feature class
ERROR 000840: The value is not a Raster Layer.
ERROR 000840: The value is not a Mosaic Layer.
Failed to execute (SelectLayerByAttribute).

What does this mean? What is the best way to do some operation like this?


texts: Image collections are not recognized as the same size for processing

texts:
I have two merged collections of images (LANDSAT and LANDTREND).

  var collection = Collection.select(['NDVI_to_ajust'],['LANDSAT']).merge(
    imageCollection.select(["NDVI"], ["LTGEE"]))

I need to store as property NDVI values calculated by each image collection, separately in array form, in 50 different samples (which are as a FeatureCollection called "antropVegPrima"). I used "reduceRegions" to do this.

var triplets = collection.map(
    function(image) {
    return image.reduceRegions({
      collection: ignoradoSamples50,
      reducer: ee.Reducer.first().setOutputs(image.bandNames()),
      scale: 30,
      tileScale: 4
    }).map(
      function(feature) {
      return feature.set({
        'imageID': image.id(),
        'date': ee.Date(image.get('system:time_start')).format('Y-MM-dd')
      })
    });
  }).flatten();

However, I have a problem that the size of the LANDSAT (124 x124) images is not the same as the LANDTREND (123 X12) image collection.

I was able to run this script for a sample collection determination, available at this link

However, for another collection of samples, the processing is not complete, as I demonstrated in this link


texts:
I create a script in order to make specific events according to the layers. I have already performed a sort according to the name of the layers. This is how I did it:

keep_layer = ['layer1', 'layer2']
for lay in QgsProject.instance().mapLayers().values():
       for k_layer in keep_layer :
           if lay.name() == k_layer :
               #Do something
           else:
               pass

This sorting allows me to make precise modifications on the layers with the same name as those entered in the keep_layer variable.
I also created a variable keep_group which will have the groups for which I want to make changes. Here is the type of subgroup I'm trying to select (screenshot) :

I thought of something like the previous sort but I can't select the name of the subgroups.
keep_group = ['group1', 'group2', 'group3']
for lay in QgsProject.instance().mapLayers().values():
        for k_group in keep_group :
            if lay.name() == k_group :
                for children in k_goup.children():
                    #Do something
            else:
                pass

How to select a subgroup according to its name in the LayerTree?


texts: Getting "missing FROM-clause" error when adding ST_Intersects constraint to PostGIS table

texts:
I have two tables in PostGIS: afrivs, which is a polyline of various rivers, and weirs2 which is an empty polyline layer where potential dam sites will be added.

Obviously dams have to cross a river to work so I would like to add a constraint to enforce this.  Following the examples here https://spatialdbadvisor.com/postgis_tips_tricks/127/how-to-apply-spatial-constraints-to-postgis-tables I tried adding the following constraint to check that any new dam polyline crosses a river polyline:

ALTER TABLE weirs2 ADD CONSTRAINT river_intersect
CHECK (ST_Intersects(afrivs.geom,weirs2.geom)=true);

However this gives me the following error message.

ERROR:  missing FROM-clause entry for table "afrivs"

I also tried framing it as a SELECT query to include a FROM clause then get
 syntax error at or near "SELECT".

Based on this https://dba.stackexchange.com/questions/149466/missing-from-clause-entry-with-constraint-on-value-from-type-in-postgresql, I tried putting the table names in () but then get:

ERROR:  column "afrivs" does not exist

What am I missing?

Using POSTGIS 2.4.4 r16526 on PostgreSQL 10.6.


texts: How to print a list of Tools from their respective Toolboxes with arcpy?

texts:
I would like to print a list of Toolboxes as listed in my Python code, however I would also like to print the tools in that Toolbox under neath the Toolbox name like so;

Analysis Tools(analysis)

Buffer_analysis

Clip_analysis

Data Management Tools(management)

DeleteRows_management

CopyRows_management

Here is the python code I have so far, but, when i run it it just prints a list of the Toolboxes without the tools that belong to them.
Here is my code;

 import arcpy
toolboxes = arcpy.ListToolboxes()
for toolbox in toolboxes:
print( toolbox)
toolbox = "Analysis Tools(analysis)" + \
"Data Management Tools(management)" + \
"3D Analyst Tools(3d)" + \
"Cartography Tools(cartography)" + \
"Conversion Tools(conversion)" + \
"Geocoding Tools(geocoding)" + \
"Geostatistical Analyst Tools(ga)" + \
"Linear Referencing Tools(Ir)" + \
"Spatial Analyst Tools(sa)" + \
"Spatial Statistics Tools(stats)" + \
"Coverage Tools(arc)"
startIndex = toolbox.find("(")
startIndex = toolbox.find("B")
tools = arcpy.ListTools("*_analysis" + \
"*_management" + \
"*_3d" + \
"*_cartography" + \
"*_conversion" + \
"*_geocoding" + \
"*_ga" + \
"*_Ir" + \
"*_sa" + \
"*_stats" + \
"*_arc")
print("/t + toolbox")


texts:
I have a CSV file with point cloud data of a terrain.
How do I use QGIS to calculate the volume of the terrain?

I tried using SAGA tool, natural neighbour, then use raster volume calculation.
However, I am pretty skeptical about the result.

I think the parameters in natural neighbour has not been set correctly. Can anyone advise me on this?


texts: Settings for point cloud in OpenDroneMap

texts:
I'm bit confused. When I process an image set in OpenDroneMap orthophoto creation works really good, also mesh creation (including textures). The only problem that I have is with point cloud. Always has missing data (see images below), I have been testing different parameter settings in settings.yaml file with no succeed.

I use LAS file to create DSM and DEM and this issue affects directly the final result. Does someone know how to solve this problem?

Orthophoto:

Point cloud (.las or .ply file)


texts: Finding the Geometry within which a Point falls

texts:
I have the following query, which takes on average about 4 seconds. How might one improve the query such that the same response comes back faster?

SELECT * FROM geometries g
    ORDER BY
      ST_Distance(g.wkb_geometry,
         ST_SetSRID(ST_Point(-77.8252923, 39.4328016), ST_SRID(g.wkb_geometry)))
    LIMIT 1

(I understand that this actually finds the nearest geometry to a point, which is the same as the geometry in which the point falls unless the point doesn't fall in any geometry. I'm open to seeing other approaches as well.)


texts: Changing the text of the Leaflet zoom buttons

texts:
I was wondering if there is a way to change the text of the Leaflet zoom buttons after they have been created.

I am aware that I could do the following:

document.getElementsByClassName("leaflet-control-zoom-in").innerHTML = 'text';

Or:

map.zoomControl.remove();
L.control.zoom({zoomInText: 'text'}).addTo(map);

But I was wondering if there is perhaps a more elegant way.

The docs do not mention a method but maybe I am missing something.


texts: Making elevation contours of raster smoother using QGIS

texts:
I am using QGIS 3.6.

I have raster of Europe map which I want to make a colored relief model.
After duplicating the layer, used hillshade on backward layer and colored the upper layer interpolation > discrete then 65 % transparency.
But the contour of color levels are not smoothed as the image shows.

Is there any way to do that?

Ex (photo):


texts:
I'm trying to build a simple multi-step selection (by attribute) with the QGIS Modeler. I created a simple sample layer to avoid data-driven errors. Now I want to query the layer and select features 1) by feature attribute "name" AND/OR 2) by feature attribute "age". The model executes the first query (name), picking up two features and then fails on the second, which should result in the selection of a single feature.

The log tells me "Error executing algorithm 1 local variable 'x' referenced before assignment See log for more details". The "History and log" (under Processing) ERROR section tells me to "See log for more details". So two logs tell me to look at a log. Is there a third log somewhere?

Another problem is that the selection comparison only allows one option (e.g. "=="). This is problematic because I want to be able to search for just "age" as well. Leaving the field for name blank selects nothing.

A few screenshots:

Using the modeler has also been a quite frustrating experience so far. Python errors out of the blue, not being able to delete algorithms before filling the input variable fields with garbage, nor being able to delete algorithms at all and QGIS just crashing. Is this "normal"?


texts: Merging shapefiles and saving them to file geodatabase using ogr2ogr?

texts:
I have multiple shapefiles. When it becomes one layer (one shapefile) the volume of shapefile gets over the limit of its capacity (I believe 4GB? is the limit of *.shp can have)

Therefore, I need to save them into FGDB using ogr2ogr

this is the code I came out with so far.

  for %f in (*disp_line.shp) do (if not exist test.gdb (ogr2ogr -progress --config shape_encoding "UTF-8" -f "filegdb" test.gdb %f -lco encoding="UTF-8") Else (ogr2ogr --config shape_encoding "UTF-8" -progress -f "filegdb" -update -append test.gdb %f -nln test_layer -lco encoding="UTF-8"))

when I see the result.. I get two layers.

t_sc1disp_line feture class has 20 rows and test_line has 286 rows
it is supposed to be 306 rows (20+286)..

What it means is that I have to merge the layers once again using QGIS or ArcGIS GUI tools.

Could you fix my weird OGR2OGR command?


texts: Read >1 geometry fields into QGIS via ODBC from non-spatial database

texts:
I have a non-spatial database in FileMaker Pro that has two geometry fields per record (line and point).  The database has a geometry_columns table and using the method outlined by @MerseyViking at Can QGIS read an ODBC connection?, I can create a live link between one of the geometry fields and QGIS.

Is there any way to get QGIS to read >1 geometry field per record via ODBC?  I want to create one layer in QGIS for each geometry field in the database, but as the geometry_columns table is static (i.e. the f_geometry_column value can't be set or chosen from QGIS), this does not seem to be possible.


texts: Copy and translate polygonz records according to attributes


texts: How are Google 3D objects made in Google Earth now?

texts: Since a few years, Google has stopped using SketchUp community Models to populate its 3D environment and switched to another technology. I'm wondering what GIS, LIDAR Image or 3D automatic method (it must be automatic) they are using to create this new realistic 3D environment all over the world (one can see for example underneath a bridge, under different trees and objects in 3D) ?


texts: How to plot elevation charts from a GPX file that has elevation data

texts:
Is there a service or package that converts a GPX file into an elevation chart?

I need either a CLI, a Web API, or a library/other software that could be wrapped up as a small service.

The GPX files already have embedded elevation -  I'm not looking for a service that returns elevation for a particular lat/lng.

Preferred output would be a PNG or SVG, though a javascript chart would also be good.


texts:
Can i create a map using plain image(jpg or png) with GeoServer and OpenLayers?

Scenario as i see it:

I upload image(somehow) to GeoServer and it converts it to a map layer
Then i create a simple html file with OpenLayers connected, and grab my layer from GeoServer

Any tutorials or documentation? I tried to find solution for this case by myself, but unsuccessfully.

Or may be i misunderstand something general in GeoServer?


texts:
I've been battling with adding layers to the TOC for some time now.  I have managed to single layers added but need to loop through a folder and add all layers.  Shapefiles I cannot get to add at all.  Only layer files.  Below is the code of where I'm at:

for layer in shp_List:
mxd = arcpy.mapping.MapDocument("CURRENT")
dataFrame = arcpy.mapping.ListDataFrames(mxd, "*")[0]

outlayer = lyr
layerfile = os.path.join(base_Folder, lyr + ".lyr")

arcpy.MakeFeatureLayer_management(layer, outlayer)
arcpy.SaveToLayerFile_management(outlayer, layerfile, "ABSOLUTE")
addlayer = arcpy.mapping.Layer(layerfile)
arcpy.mapping.AddLayer(dataFrame, addlayer, "BOTTOM")

#addLayer = arcpy.mapping.Layer(layer)
#arcpy.mapping.AddLayer(dataFrame, addLayer, "BOTTOM")
#arcpy.RefreshTOC()
#arcpy.RefreshActiveView()
del addlayer, mxd

The above code will add each layer to the TOC, then infuriatingly removes it and adds the next one.  When the whole script finishes there is nothing left in the TOC.  I've tried adding this script to a model and creating a derived output parameter which is then added as a model parameter with "add to display" checked.  I've also checked the geoprocessing options to make sure the box to add to display is checked.  I'm running Arc 10 (no service packs).  Can anyone help?


texts:
I draw a polygon on my map and I need to know if there are some features in it via WFS to get back a JSON. I guess I need to use CONTAINS property but when I read it it says that a contains b.

So I tried to write

  https://myurl.com/geoserver/wfs?service=wfs&format_options=id_policy:reference_no=false&version=1.1.0&request=getfeature&typename=biglayer:layer&outputFormat=application/json&CQL_FILTER=CONTAINS(POLYGON((465995.8321513263%205177213.444093402,525922.4623269045%205157645.564852396,434809.5246109744%205123401.776180638,465995.8321513263%205177213.444093402)),geom)

But I've got the error

  Binary geometry filter, but first expression is not a property name?

If I put geom in first position, it returns me 0 feature, and I don't get it neither...

EDIT : Well I finally understood that I have to use WITHIN but I've still got no result... Here's my url :

  https://myurl.com/geoserver/wfs?service=wfs&format_options=id_policy:reference_no=false&version=1.1.0&request=getfeature&typename=biglayer:layer&outputFormat=application/json&CQL_FILTER=WITHIN(geom,POLYGON((465995.8321513263%205177213.444093402,525922.4623269045%205157645.564852396,434809.5246109744%205123401.776180638,465995.8321513263%205177213.444093402)))

Any clue?


texts:
I just noticed that ogr2ogr.exe doesn't seem to work when executed with admin right ..?
I tested with the following parameters:

ogr2ogr.exe -f "GeoJSON" "C:\Users\root\Desktop\Plugin_new\Shapefiles_Karte\Landkreise\Bundeslaender.geojson" "C:\Users\root\Desktop\Plugin_new\Shapefiles_Karte\Landkreise\Bundeslaender.shp" -overwrite

I used a cmd started with admin privileges and one without. Both don't give any oputput when executed.
Executed without admin rights creates the geojson file. Executed with admin rights does nothing. No output created. No errors ..

Is this a normal thing?


texts: Reading coordinates from shapefile using shapely

texts:
I have coordinates and I am trying to find the roads they are on in a shapefile.
My code is as follows:
import shapefile
from shapely.geometry import Point # Point class
from shapely.geometry import shape # shape() is a function to convert geo objects through the interface

shp = shapefile.Reader(r'NBRoads/geonb_nbrn-rrnb_road-route.shp')

point_to_check = (-68.41856,47.28858) #example coordinate
all_shapes = shp.shapes()
all_records = shp.records()
print(shape(all_shapes[0]))
for i in range(len(all_shapes)):
    boundary = shape(all_shapes[i])
    if boundary.distance(Point(point_to_check)) < 1e-8:
       name = all_records[i][2]
       print("The point is in", name)

The problem is the LINESTRINGs coordinate system is not normal lat/long.
LINESTRING (2433182.259100001 7613044.969799999, 2433256.279800002 7613043.196899999, 2433348 7613041, 2433368 7613040, 2433388 7613040, 2433411 7613039, 2433455 7613038, 2433500 7613037, 2433542 7613036, 2433609 7613034, 2433618.4628 7613033.746199999, 2433618.502300002 7613033.745099999)
In my research I believe it has something to do with converting them using the .prj file from the shapefile but I cannot seem to make the link.
Here is the .prj info:
PROJCS["NAD_1983_CSRS_New_Brunswick_Stereographic",GEOGCS["GCS_North_American_1983_CSRS",DATUM["D_North_American_1983_CSRS",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Double_Stereographic"],PARAMETER["False_Easting",2500000.0],PARAMETER["False_Northing",7500000.0],PARAMETER["Central_Meridian",-66.5],PARAMETER["Scale_Factor",0.999912],PARAMETER["Latitude_Of_Origin",46.5],UNIT["Meter",1.0]]


texts: Transform an azimuth measured in a given projection system to its value in a separate projection system

texts:
I have a set of polygons in which one of the attributes is azimuth. This azimuth was measured in ESPG:3310 (California Albers). I want to transform this measure to ESPG:4326 (WGS 84). I know the values are different as when measuring the azimuth of a line in the two projections, I get a different value.

E.G. if I draw a line in QGIS and measure the azimuth between the first and last points using this code:

degrees(azimuth(start_point($geometry), end_point($geometry)))

I get a value of 179.721 when my layer is projected in ESPG:4326 and a value of 180.08 when projected in ESPG:3310.

I want to convert the value without the actual line being available to me. So I can't simply reproject the layer, but have to reproject the value.
Unsure if that would be necessary to have the location of the measured line, but if it is, we could use the centroid of the polygon.


texts: What's the difference between Lambert conformal conic projection and Lambert conformal conic alternative projection in PROJ4?


texts:
I am trying to plot points (latitude, longitude) inside a raster image in Python. First, I read the image using rasterio and plot it using the same package, after transform the image from uint16 to Byte (uint6):

#Transform to uint8
def scale8bit(image):
    scale = float(256) / (image.max() - image.min())
    return np.clip(np.round(np.multiply(image, scale)), 0, 255).astype(np.uint8)
#Read image and plot visible bands
src = rasterio.open(path)
array_8b = img_as_ubyte(img_as_float(src.read([3, 2, 1])))
a = scale8bit(src.read([3, 2, 1]))

Now, I have some WKB points which coincide --as far as I can tell from the data-- with the bounding box of the image. Hence, I am trying to run the following in order to visualize the images and the points (read using shapely):

from shapely import wkb
from rasterio import plot.show
for i, geom in zip(selected_images_paths[0:10], subset_images["geom"][0:10]):

  #Open file using rasterio and change to int8
  src = rasterio.open(i)
  src_8bit = scale8bit(src.read([3, 2, 1]))

  #Graph!
  fig, ax = plt.subplots(figsize=(30, 30))
  plot.show(src_8bit)

  #Open geom with shapely
  geom_wkb = wkb.loads(geom, hex = True)
  x,y = geom_wkb.xy
  ax.scatter(x, y, marker = 'o', s=50, c='r')

The units of the graphs are clearly on color intensity and not in long, lat, which are the units of the x,y tuple. I've tried different options, but I can figure out how to overlap both data sets.

Is there something wanting, or should I take another approach?

EDIT
I have been following the code described in this question (Plot shapefile on top of raster using plot and imshow from matplotlib). Nonetheless, the result is sub-optimal.

When I use ax = plt.gca() to set the axes in the subplot, the fig definition is completely ignored for some images and the resulting figure is small. I think the error must be on the ax or fig definition, but I've tried other suggestions in the documentation without success.


texts: Generating TIN from elevation data in Python

texts:
I would like to generate a 3d terrain model of an area, and I am stuck on creating the TIN (which I will convert to an STL). I know there are plenty of ways to use ArcGIS, QGIS, 3DEM, and various other programs to go through and generate a TIN for an area by hand and Stack Exchange answers explaining them. But that does not fit my use case, I need a solution in which I can script the whole process. Preferable in python, but if another solution exists I would not throw it out. I am also willing to sacrifice some accuracy for speed. The result will be printed and does not need to be perfect.
I have SRTM 1 arc-second data, which I downloaded from NASA. The first thing I attempted to do was generate a TRN out of it, which did work, however it was huge and unwieldy to use. There are way too many triangles. The plains have the same amount of triangle as the Himalayan Mountains. If the data were sparse, I would use scipy.spatial.Delaunay, but since it is a regular grid it doesn't simplify things any.
I also thought about just down sampling to make it smaller (or using coarser resolution data), but that doesn't actually address the real problem of needing the simplify the surface that was generated.
Is there a library I am missing to generate TINs in python?
Is there a better approach to what I am trying?
Possibly generate contour lines for an area, then triangulate them?
If that is possible, what library would I use to generate the contour lines?
If nothing else can someone point me to documentation, or a paper on how to generate a simplified TIN out of gridded elevation data?
If I have to write the library myself I will. The primary requirement I have is that it has to be something I can script. A cross platform solution would be ideal, I would like to run it on Linux, but I will take what I can get.


texts: How to fix the error between two coordinate systems?

texts:
How do you think the shift of line is due to morphological change of river? Or it is due to deviation between two coordinate systems?
Please point out !!


texts:
I am using osm2po to load the Europe dataset into a PostGIS database to use with pgRouting.

I find the default speeds in osm2po.config to work quite well, but when a maxSpeed tag is found this speed is used as-is to compute the cost.
Shouldn't it somehow be lowered by some ratio ? When used like this, osm2po underestimates the driving time IMHO.

i.e. in Belgium the maxSpeed on motorways is 120 km/h, but I'd rather use an average speed of 110 km/h.

Is it possible to either specify a ratio to apply to maxSpeed tags, or maybe some a conversion table of some sort ?


texts:
I was trying to execute this code:

# Load geojson files
country = gpd.read_file('svn_buffered.geojson')

train_shp = "large_training_aoi.geojson"
test_shp = "large_test_aoi.shp"
train_polyg = gpd.read_file(train_shp)
test_polyg = gpd.read_file(test_shp)

# reproject the polygons
target_crs = CRS.UTM_33N
country = country.to_crs(target_crs.pyproj_crs())
train_polyg = train_polyg.to_crs(target_crs.pyproj_crs())
test_polyg = test_polyg.to_crs(target_crs.pyproj_crs())

fontdict = {'family': 'monospace', 'weight': 'normal', 'size': 11}

fig, ax = plt.subplots(figsize=(20, 20))
test_polyg.plot(ax=ax, facecolor='w', edgecolor='r', alpha=0.5) #test area in red
train_polyg.plot(ax=ax, facecolor='w', edgecolor='g', alpha=0.5) # train area in green
country.plot(ax=ax, facecolor='w', edgecolor='b', alpha=0.5)
ax.set_title('Selected train AOI (in green) and test AOI (in red)');
plt.axis('off')

I tried to use the following syntax:

test_polyg = test_polyg.to_crs(epsg=32632)

and I got the error each time but only for the test data.

ValueError                                Traceback (most recent call last)
<ipython-input-20-323e17a794e3> in <module>()
    11 country = country.to_crs(target_crs.pyproj_crs())
    12 train_polyg = train_polyg.to_crs(target_crs.pyproj_crs())
---> 13 test_polyg = test_polyg.to_crs(target_crs.pyproj_crs())
    14
    15

1 frames
/usr/local/lib/python3.6/dist-packages/geopandas/geoseries.py in to_crs(self, crs, epsg)
   420         if self.crs is None:
   421             raise ValueError(
--> 422                 "Cannot transform naive geometries.  "
   423                 "Please set a crs on the object first."
   424             )

ValueError: Cannot transform naive geometries.  Please set a crs on the object first.

Any idea on how can I fix it?


texts: What is the maximum attachment file size in Esri file and ArcSDE geodatabases?

texts:
I am trying to figure out the maximum size of attachment in a file and ArcSDE geodatabases while working on a project specification details. I am not able to find these figures documented by Esri anywhere.

Running own tests, it seems as a file geodatabase has a limit of ca 600MB per file that is being attached. When trying to attach everything that is larger than 600MB, I am getting an error (where Points3D is the name of the feature class I am trying to attach a file to).

Just to see whether it is a file geodatabase only limit, I have run some tests with a SQL Server enterprise geodatabase and here I am getting stuck at somewhere between 700 and 800MB and either get the same error message (except that instead of the feature class name I have my SQL Server database filegroup name) or have ArcMap unexpectadly crashed and closed.

What is the official/documented maximum file size limit that one can upload to an Esri file geodatabase? I also welcome comments if someone has succeeded uploading anything larger than 1GB to a geodatabase as an attachment.


texts: How to save bookmarks with names in QGIS?

texts:
When I create a bookmark the name and project that I enter are not saved when I click Close.

The manual says there should be an OK button but there isn't. Using QGIS 1.8


texts: Create a grid over shapefile then save each section into GeoTIFF

texts:
I have a massive shapefile with 200,000+ building polygons. I've loaded these up in ArcGIS. I need to convert this into a usable raster (GeoTIFF) that can then be tiled using MapTiler.

My current process is as follows:

Click on the Full Extent icon in order to have all the polygons in
my current view. The scale is at 1:261,289.
Export this as raster using File > Export Map > TIFF > Write GeoTIFF tags (in Format tab). Set the dpi to 500.
When zoomed all the way into this GeoTIFF image, the result looks
too pixelated. You can barely tell they are polygons.

I need to the scale to be at 1:26,261 (or 1:25,000 for simplicity) in order to get decent looking polygons.

How can I create a grid over this vector file with 200,000+ polygons and then save each section into a GeoTIFF file at 1:26,261 scale?

I also tried this in QGIS. Similar issues, though somewhat better results in ArcGIS.


texts: Please can someone let me know how I can create a schema (using SimpleFeatureTypeBuilder()), read and write image data? Simple data types (Integer/Float etc) work fine but how do I handle image data?


texts: How to calculate the area of overlapping buffers within a single layer

texts:
I am very new to GIS, so don't know if this is EXTREMELY easy or impossible.

I have GPS coordinates for stores (about 15000) for which I have created a buffer zone around each store location obviously (or not) in a single layer. In metropolitan areas I now have overlapping buffer areas. I would like to calculate the area where buffers overlap.

I have established how to do this if each point were a separate layer, but with 15000 points, I just do not have the time or the energy to do this.

Please help.

I am working in QGIS.


texts: Label placement over other features and labels in ArcGIS Pro

texts:
I need place labels regardless of collisions with other layers. I use ArcPy and CIM. I need to place the label in the center of the polygon. At the same time, I need it not to change its position or disappear if it collides with another labels or feature. I try this code, but some labels are placed outside the polygons.
show_labels = m.listLayers("polygon23*")[0]
show_labels.showLabels = True

for lyr in m.listLayers("polygon23*"):
        if lyr.supports("SHOWLABELS"):
            if lyr.showLabels:
                for lblClass in lyr.listLabelClasses():
                    lblClass.expression = ''' "<FNT size='6'>" + $feature.label + "</FNT>"'''

label_show = m.listLayers('polygon23*')[0]
l_cim = label_show.getDefinition('V2')
lc = l_cim.labelClasses[0]

lc.maplexLabelPlacementProperties.polygonPlacementMethod = 'StraightInPolygon'
lc.maplexLabelPlacementProperties.placeOnlyInsidePolygon = True
lc.CIMStandardLabelPlacementProperties = 'allowOverlappingLabels'


texts:
I'm using a map for a QGIS project that has labels once you zoom in.

When I import the map onto Composer, I can't see the labels from my current zoom level, but they do show up once I export it as an image.

Does anyone know how I can prevent this from happening?

I would like QGIS to export the map exactly how I see it.


texts: How can I get QGIS Print Composer to print OSM map at desired zoom level?

texts:
Like the title suggests, I'm trying to print a map that uses OSM Stamen data from the OpenLayers plugin.

My map looks like this in QGIS:

And like this in the Print Composer:

But when I export the map to an image, it appears to be at a smaller zoom level.

Is there any way to actually make the map look like I'm desiring? I came across this thread, but it didn't lead me anywhere productive.


texts: Seeking most natural unit of analysis in between US county and state

texts:
Do you have any suggestions for the most natural unit of analysis in between US county and state?
I am looking at Core Based Statistical Areas but those don't seem to cover all counties.


texts: Create new attribute ranking polygons by size in QGIS

texts:
I'd like to create a new field in the attribute table which gives the rank of the polygon by size. For example, the polygon of the largest size would be rank 1 and the smallest rank N.

I was hoping sorting the table by shape size and then using $id would work, but that appears to be tied to the row ID #.

I am using QGIS.


texts: Project raster tool stays white at iteration in ModelBuilder?

texts:
I want to project a bunch of rasters (WGS84) to different UTM (WGS84) coordinate systems. The coordinate system depends on a field value of a point shape.

So on each raster lays a point with its spatial reference.

I scripted a little ModelBuilder toolchain for it but unfortunately the Project raster Tool stays white and doesn't execute.

Some detail information for the toolchain:

I extract the information wich point lays on the individual raster with the sample tool. Then I get for each raster a table with many no data points and one with the stored EPSG Information. This point is extracted by the table select tool and then I get the EPSG with the get field value tool and store it as spatial reference.

The spatial reference I use as Output coordinate system.
But the Project raster tool requires the "optional" Geographic Transformation although both coordinate systems use the WGS84 referece.

I tested the toolchain until the projection part and there is no error.
I've checked if the rasters lay in the UTM zone I want to project them to.
When I run the tool first, then disconnect the spatial reference and the Projection and then reconnect, the tool gets colorful. But when I run the toolchain again (without validation) only one Raster is projected several times to the same coordinate system until I cancel it.
at no point any error message occurs when running the tool
I've changed the data type of the field value to "Coordinate System" with no effect


texts: How can I multiply each band of multiband image by separate scalars using ArcGIS Desktop?

texts:
I have an image with 5 bands.
I need to multiply each band of this image by different scalars one time
(I'm lazy to do it one by one) and then get multiband image as a output.

For example (band1*1.1) (band*1.2) (band3*1.25) etc.
Then output return to an image with 5 bands.

How can I multiple each band of multiband image by separate scalars at the same time?


texts: Ignoring field types in Python list returned by arcpy.ListFields()?

texts:
I am generally just wondering if there is a nice and easy way in python to ignore all fields of a specifiy type (eg string) in a feature class.

At the moment I am doing it with

field_names = [f.name for f in arcpy.ListFields(featureclass)]
field_names = field_names[19:len(field_names)-2]

which ignores the first 19 and last 2 fields of my FC which are all string-fields but it would be better for me to have a code which just generally ignores all fields of a specific type


texts: Can I visualize a spatial movement (in a room) data using QGIS?

texts:
I am from a computer science background and I would say I have zero experience with GIS software and terminology. Can anyone confirm if I can use QGIS for the following problem:

I used a motion capture system to collect the movements of people in a room.

The system I used has its own coordinate system where the center of the room is (0,0)

The collected data is in the form:
timestamp, id, x, y
e.g. 1/27/2018 7:26:17 PM, user1, -1.62, 1.71

Where x and y represent the position with reference to the room center in meters. Look at the sketch:

So as now I have a list of points of the person movement, I want to create a heatmap of that. I used QGIS to import my data and I was able to create a heatmap but I have to zoom in to see a better view. I think this is related to CRS. Can I change the coordinate system to have (0,0) is the center of the view? and the total view box is 10 by 6 meters?


texts:
I am working with leaflet tiles loading tiles to a map div

http://jsfiddle.net/tBCd5/

I am now playing with some html5 canvas distortions like below

http://jsfiddle.net/UF4PK/

It would be cool to load base tiles to a canvas and then be able to use the above distortions on the base map.

Is it possible to load my map tiles into a canvas element rather than a div?

I have seen examples of creating canvas overlays but not an example where base tiles are loaded directly to the canvas element.


texts: How to not overlap cluster popups with FeaturePopups control?

texts:
FeaturePopups control resolve the problem of cluster popups, but there is a small design problem in my opinion: the event popup have the same anchor as the cluster popup and they partial or total overlap.

Is it possible to anchor the second (event) popup on the link on first (cluster) popup?

Thank you!


texts:
I have a DEM (saved as a .tif) and, using QGIS, I want to clip the DEM using a vector mask.

The function will not work.

Am I missing critical step?


texts: Faster, more efficient way to extract stream network from DEM

texts:
I am using QGIS and use r.stream.extract to extract the stream lines from a single DEM tile. However, this takes long (about 15 minutes) to complete, especially for a rugged terrain. I plan to extract streams for a much larger area with several merged DEM tiles (in geopackage format), but I am worried it will take even longer.

Is there a faster way to extract the streams (preferably using command line tools like ogr2ogr)?


texts:
Is there an average size for British postcodes? Or a method of determining the size for a particular British postcode? Thanks.
Considering creating a related program and this would be valuable information.

I'm referring to postcodes like:

BT28 1AW
LS9 8AH
L15 3DE

By size, I'm referring to the amount of area they cover (1 miles - 20miles)