texts: SetSpatialFilter worked unexpected when using OGR by Python

texts:
I am using OGR Python binding to read the shapefile, this is the code:

.....
bbox="POLYGON ((20.18164062 30.63867188,20.63867188 30.63867188,20.63867188 30.18164062,20.18164062 30.18164062,20.18164062 30.63867188))";
print "fs:before:",datalayer.GetFeatureCount()
datalayer.SetSpatialFilter(bbox)
print "fs:after :",datalayer.GetFeatureCount()
for feature in datalayer:
    geometry = feature.geometry().Clone()
    if not geometry.Intersect(bbox):
        continue
    # do my job

However the above codes will give me output like:

fs:before: 156396
ERROR 1: Attempt to read shape with feature id (157823) out of available range.
..................................
ERROR 1: Attempt to read shape with feature id (160730) out of available range.
More than 1000 errors or warnings have been reported. No more will be reported from now.
fs:after : 0

It seems that OGR try to read a feature no exist in the shapefile. Since there are no features with the id in the output. Neither for id "157823" or "160730".

And if I comment line:

datalayer.SetSpatialFilter(bbox)

No errors occur.

And this is a screenshot of the data and the bbox used as filter:

http://pbrd.co/1kAn4R3

The green features are my data, and the red square is the bbox used to filter the features.

What's the problem?

Update the complete codes:

def renderTile(coord,items):
    bbox=create_bbox(coord,0.15)

    for item in items:
        datalayer,datasource=_open_layers(item)
        _get_features(coord,bbox,datalayer])

def _open_layers(path):
    datasource=ogr.Open(path)
    datalayer = datasource.GetLayer(0)
    return datalayer,datasource

def _get_features(coord,bbox,datalayer):
    #print "fs:before :",datalayer.GetFeatureCount()
    datalayer.SetAttributeFilter('{0} <= {1}'.format('FSCALE',coord.zoom+1))
    datalayer.SetSpatialFilter(bbox)
    #print "fs:after :",datalayer.GetFeatureCount()

    for feature in datalayer:
        #print feature.GetFID()
        geometry = feature.geometry().Clone()
        if not geometry.Intersect(bbox):
            continue

        geometry = geometry.Intersection(bbox)
        handle_result(feature,geometry)

coord =(10,12000,2000)
items= ["/home/kk/data/shps/points.shp","/home/kk/data/shps/roads.shp"......];
renderTile(coord,items)

The above codes are served by werkzeug. The coord is create by HTTP parameters, the items are read from configuration file.


texts: QGIS - PostGIS offline editing with constraints

texts: I note that QGIS offline editing (QGIS 2.28.9) of a PostGIS database does not seem to bring the PostGIS database constraints with it to SQLite (lookup tables etc).  Is this being considered?


texts:
I've tried several different approaches to getting the simple geoprocessing package that I created to execute from my C# project . . .

        ESRI.ArcGIS.Client.Tasks.Geoprocessor _localGPService = null;
        string gpkPath = ConfigurationManager.AppSettings["EsriRuntimeGeoPackageLocation"];
        string geoObjectPath = ConfigurationManager.AppSettings["EsriRuntimeGeoProcessorObject"];
    LocalGeoprocessingService.GetServiceAsync(gpkPath, GPServiceType.Execute, (gpService) =>
    {
        if (gpService.Error != null)
        {
         MessageBox.Show(gpService.Error.Message);
        }
        _localGPService = new ESRI.ArcGIS.Client.Tasks.Geoprocessor(gpService.UrlGeoprocessingService + geoObjectPath);
        DataContext = this;
    });
    List<ESRI.ArcGIS.Client.Tasks.GPParameter> parameters = new List<ESRI.ArcGIS.Client.Tasks.GPParameter>();
    _localGPService.ExecuteAsync(parameters);

I've also tried the following . . .

 private LocalGeoprocessingService _localGPService = null;
    string gpkPath = ConfigurationManager.AppSettings["EsriRuntimeGeoPackageLocation"];
    _localGPService = new LocalGeoprocessingService(gpkPath, GPServiceType.Execute);

    _localGPService.StartAsync((callback) =>
    {
        if (callback.Error == null)
        {
            CoverButton.IsEnabled = true;
        }
        else
        {
            MessageBox.Show("Error starting gp service");
        }
    });
    string geoObjectPath = ConfigurationManager.AppSettings["EsriRuntimeGeoProcessorObject"];
    ESRI.ArcGIS.Client.Tasks.Geoprocessor gp = new ESRI.ArcGIS.Client.Tasks.Geoprocessor(_localGPService.UrlGeoprocessingService + geoObjectPath);
   List<ESRI.ArcGIS.Client.Tasks.GPParameter> parameters = new List<ESRI.ArcGIS.Client.Tasks.GPParameter>();
   gp.ExecuteAsync(parameters);

When it gets to the last line (execute), I don't get an error, but it also doesn't do anything.  I've been tracing through the code, and the geoprocessor object does get created and invoked with a URL, so I don't know why it does nothing. Tracing through, it appears as if everything is OK.  The tool that the gpk was created from runs OK when I execute it in ArcMap.
Any ideas as to what might cause this? I just don't know what to look for anymore.
Thanks!


texts: Can I interactively type multiple string records as input for GP Widget?

texts:
I am running a GP script tool where a required input is set up as a RecordSet for the user to type multiple string values that will be used by the tool for further analysis. However, when I publish my tool as GP service and consume it within Web AppBuilder 2.3, the input RecordSet becomes a JSON set of fields:

{"displayFieldName":"","fields":[{"name":"OBJECTID","type":"esriFieldTypeOID","alias":"OBJECTID"},{"name":"Name","type":"esriFieldTypeString","alias":"Name","length":50}],"features":[],"exceededTransferLimit":false}

What I would like, instead, is an input parameter where the user can interactively type multiple string values that will be passed to the task. How can I accomplish this? Should I not use RecordSet?


texts: How to calculate percentage of selected area from DEM?

texts:
I am working with elevation data and performed a local operation to convert it from meters to feet.

I am interested in finding the area (expressed in percentage) that is above 4,000 feet, but I am having difficulty recalling how to do this.


texts: What is the procedure to publish a WMS using Geoserver?

texts:
I wanted to be able to load a network connection in Google Earth (or Google Maps) to display daily updated maps using Geoserver. I want something very simple but I have never published web maps before.

Can anyone explain the steps required to do this?

Update:

Meanwhile i managed to learn most part of it. Now i only need to know how to perform spatial analysis and publish the results all automatically.


texts: Transformation of user-defined projected coordinatesystems

texts:
I want to transform from one projected coordinate system (pcs) to another using ArcGIS 10. Unfortunatelly the starting system is not predefind by ESRI and thus I´m not sure the software is doing what I want. The parameter of the starting system are the the following:

Europe_Lambert_Conformal_Conic
Authority: Custom

Projection: Lambert_Conformal_Conic
False_Easting: 0,0
False_Northing: 0,0
Central_Meridian: 10,5
Standard_Parallel_1: 48,667
Standard_Parallel_2: 53,667
Latitude_Of_Origin: 51,0
Linear Unit: Meter (1,0)

Geographic Coordinate System: GCS_WGS_1984
Angular Unit: Degree (0,0174532925199433)
Prime Meridian: Greenwich (0,0)
Datum: D_WGS_1984
  Spheroid: WGS_1984
    Semimajor Axis: 6378137,0
    Semiminor Axis: 6356752,314245179
    Inverse Flattening: 298,257223563

... and this should be transformed into: ETRS_1989_LAEA (WKID: 3035). What I have done is:
1. unprojecting to GCS WGS 1984
2. transformation to GCS ETRS 1989 (using the transformation option: ETRS_1989_To_WGS_1984)
3. projecting to ETRS_1989_LAEA

Now my question: does ArcGIS recognize the different userdefined standardparallels of the Lambert Conformal Conic projection as there is a predifined Lambert Conformal Conic projection already included in te software while unprojecting. The parameters of the predifined psc are the following:

Europe_Lambert_Conformal_Conic
WKID: 102014 Authority: ESRI

Projection: Lambert_Conformal_Conic
False_Easting: 0,0
False_Northing: 0,0
Central_Meridian: 10,0
Standard_Parallel_1: 43,0
Standard_Parallel_2: 62,0
Latitude_Of_Origin: 30,0
Linear Unit: Meter (1,0)

Geographic Coordinate System: GCS_European_1950
Angular Unit: Degree (0,0174532925199433)
Prime Meridian: Greenwich (0,0)
Datum: D_European_1950
  Spheroid: International_1924
    Semimajor Axis: 6378388,0
    Semiminor Axis: 6356911,946127947
    Inverse Flattening: 297,0

As you may recognize I`m not not really sure what ArcGIS is doing in the transformation process of of one psc to another, in contrast the transformation from gsc to another gsc is understandable.


texts: Selecting polygon as Input in QGIS Modeler?

texts:
I need to create a model where I select two polygon and the first one is clipped with the second one. Difference function is doing it. But asks for two input layers.

How can I add as Input a Polygon to be selected from canvas?

All the inputs I try ask me for a whole layer and then I can filter by selected features in it.

I have both polygons in the same layer so I click "Only selected features" in both input layers (same input layer in both) and it's not working.

Clipper plugin is not working for me. Maybe because layer is Multipolygon.


texts: Tiling Vector data without losing geometrical and topological information

texts:
I am trying to figure out a way to sub-divide a vector based map into smaller blocks. I can tile the map based on different strategies such as q-tree. However, it is unavoidable to cut some polylines and polygons that cross the boundaries of the blocks. By doing this, we will lose their geometric and topological attributes.

The tiled data can only be used for display but not for spatial analyses because they have partly lost their geometric and topological attributes. Do you have any suggestion or experience with sub-dividing the vector data and still maintaining their us usability for spatial analyses?


texts: GDAL osgeo for python treating Z dimension differently

texts:
I have a space defined by a polyhedral surface, and I would like to check if a point is inside this space, and I would like to do this in python with gdal.

After trying the default installation of gdal with conda, there were some functions that did not work that I wanted regarding sfcgal, so I compiled gdal from source with sfcgal (as this seemed to contain extra 3D functionality) and geos,for python3. I have tried the following in python:

>>> from osgeo import ogr
>>> ps2 = 'POLYHEDRALSURFACE Z (((0 0 0,0 0 1,0 1 1,0 1 0,0 0 0)),((0 0 0,0 1 0,1 1 0,1 0 0,0 0 0)),((0 0 0,1 0 0,1 0 1,0 0 1,0 0 0)),((1 1 0,1 1 1,1 0 1,1 0 0,1 1 0)),((0 1 0,0 1 1,1 1 1,1 1 0,0 1 0)),((0 0 1,1 0 1,1 1 1,0 1 1,0 0 1)))'
>>> geom2 = ogr.CreateGeometryFromWkt(ps2)
>>> geom2.CloseRings()
>>> point3 = 'POINT Z (0.5 1.5 0.5)'
>>> point3_g = ogr.CreateGeometryFromWkt(point3)
>>> point3_g.Within(geom1)
False
>>> # This is correct behavior
>>> point2 = 'POINT Z (0.5 0.5 1.5)'
>>> point2_g = ogr.CreateGeometryFromWkt(point2)
>>> point2_g.Within(geom1)
True
>>> # This is unexpected behavior!

Is there anything I am missing for gdal to treat the Z dimension like one of the other dimensions?


texts: Using SHAPE field in a covering index for SQL Server

texts: Can I use the SHAPE field in a feature class as one field in a covering index? In addition to SHAPE, we almost always see reads on COLUMN2. I cannot create an index like this (SHAPE, COLUMN2) in ArcCatalog, but I could in SSDS. Is this advisable?


texts: Difference between Landsat-8 data in Earth Engine data catalog and Cloud collections?

texts:
Recently, I am trying to export some Landsat-8 surface reflectance data from Google Earth Engine, and the values of these data are around thousands, which can also verified by example code provided by Google Earth Engine (min:0, max:3000 for visualization):

However, it is not very convenient for Google Earth Engine to export a huge size of data, so I found that Google Cloud Storage has public bucket (cloud.google.com/storage/docs/public-datasets/landsat) for Landsat-8 data according to the suggestion from @spatialthoughts.

When I downloaded some Landsat-8 data from Cloud collections and open them in ENVI, I found DN values of these data are about 10 times bigger than those data from Google Earth Engine. Plus, I also found FAQ of Google Earth Engine mentioned that:

Therefore, I am wondering is there any preprocess procedure for Landsat-8 data in Google Earth Engine? I really need many Landsat-8 full images exactly same as Google Earth Engine.


texts: Open link in new window when clicking on marker

texts:
I have my map showing my different markers - from info pulling from MySQL. But I need to let the marker open html when clicked on.
At this moment the popup is just the marker label and does not take you to the link provided.
<!DOCTYPE html>

    <title>AJC Tracking</title>
    <meta charset="utf-8" />
    <link rel="stylesheet" href="https://d19vzq90twjlae.cloudfront.net/leaflet-0.7/leaflet.css" />


texts: UnicodeDecodeError: While installing gdal in a new Python Virtual Environment

texts:
I was trying to create a Python Virtual Environment for some initial GIS development. I installed pip and virtualenv and then tried the following commands from the Linux-Terminal to accomplish the task

$ cd Test/virtualenvs_for_test
$ virtualenv --no-site-packages postgis-cb-env
$ source postgis-cb-env/bin/activate
$ pip install gdal

But found the following error message:

  error: command 'gcc' failed with exit status 1

  Cleaning up...
  Command /home/gislinux/Test/virtualenvs/postgis-cb-env/bin/python -c "import setuptools, tokenize;file='/home/gislinux/Test/virtualenvs/postgis-cb-env/build/gdal/setup.py';exec(compile(getattr(tokenize, 'open', open)(file).read().replace('\r\n', '\n'), file, 'exec'))" install --record /tmp/pip-h8eZSu-record/install-record.txt --single-version-externally-managed --compile --install-headers /home/gislinux/Test/virtualenvs/postgis-cb-env/include/site/python2.7 failed with error code 1 in /home/gislinux/Test/virtualenvs/postgis-cb-env/build/gdal
  Traceback (most recent call last):
    File "/home/gislinux/Test/virtualenvs/postgis-cb-env/bin/pip", line 11, in
      sys.exit(main())
    File "/home/gislinux/Test/virtualenvs/postgis-cb-env/local/lib/python2.7/site-packages/pip/init.py", line 185, in main
      return command.main(cmd_args)
    File "/home/gislinux/Test/virtualenvs/postgis-cb-env/local/lib/python2.7/site-packages/pip/basecommand.py", line 161, in main
      text = '\n'.join(complete_log)
  UnicodeDecodeError: 'ascii' codec can't decode byte 0xe2 in position 42: ordinal not in range(128)

I have no sure idea what this error is exactly. Although initial web-hunting suggested a need to change encoding from ascii to utf8, but how?


texts: gdalwarp does not use lanczos resampling on image with alpha channel

texts:
I want to warp Geotiff image of shaded relief (RGBA) in projection EPSG:8353 to the projection EPSG:3857 to use it in Mapnik web tiled map. I want to use lanczos resampling, but the result image looks "jagged", like nearest resampling has been used. I am using this command:
gdalwarp \
  -of GTiff -overwrite -r lanczos  \
  -tr 1.194328566968441 1.194328566968441 -tap \
  -t_srs EPSG:3857 shading-sjtsk.tif shading-web.tif

The target resolution in georeferenced units is set to be the same as zoom level 17 in web mercator projection of tiled map.
Please could you advise me, why the lanczos resampling is not used? I tried also other methods (like cubic), but image always seems to be indistinguishable from nearest.
Input image:

Output image:

Original tif images are shared here.


texts: Continuous buffer from individual polygons in QGIS

texts: I have a stream segment for which I am trying to find a buffer around with criteria as 5m of elevation. So for each straight line within a stream segment, I have found buffer (polygon) for it. Now I need to merge this buffers in such a way that they form a continuous single polygon i.e. without any gaps and overlaps between them. Is there any way in QGIS or any python code to achieve the same?


texts: Draw a line of desired length and bearing in OpenLayers?

texts:
How can I create a straight line when I know the coordinates of the first point, and I know the distance and angle to the end point but I don't know the coordinates of the end point?

Example:

Start point 1: -44.14984703 -19.82517271
End point 2: (unknown lat and long) distance => 200km (note: at 75 degrees from point 1)

Note: if there is also a way in which I can discover the lat and long based on the distance from point 1 and the degree (75º), I think it would work too, I just don't know how to apply.

(Translated from Portuguese using Google Translate.)


texts: Reprojecting LiDAR data with libLAS, error: "latitude or longitude exceeded limits"?

texts:
I have a lidar file (file.las) with the following spatial reference:

Spatial Reference:
PROJCS["NAD83 / UTM zone 11N",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.2572221010002,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-117],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","26911"]]

Then I reproject it to WGS84 with las2las:

las2las --a_srs EPSG:26911 --t_srs EPSG:4326 -i file1.las -o output.las

And I get:

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

But when I try to revert the transformation: from WGS84 to NAD83/UTM 11N, I get an error:

las2las --a_srs EPSG:4326 --t_srs EPSG:26911 -i output.las -o wgs2utm.las
ERROR 1: latitude or longitude exceeded limits
error: Could not project point for ReprojectionTransform::latitude or longitude exceeded limits0.

Why I got this error, if I am just reversing the transformation (originally from NAD83 to WGS84, and then from WGS84 to NAD83)?.


texts: Splitting table column and using to update another table in ArcPy

texts:
I have a table (T1) which I need to populate (using python for ArcGIS) the fields Val1, Val2 and Val3

ID Val1 Val2 Val3

0

1

2

with data from another table (T2):

ID Data

0 V1

1 V1

2 V1

0 V2

1 V2

2 V2

0 V3

1 V3

2 V3

so the original table (T1) becomes:

ID Val1 Val2 Val3

0 V1 V2 V3

1 V1 V2 V3

2 V1 V2 V3

I currently do this using a SearchCursor nested in the UpdateCursor.

I thought a quicker solution would be to split T2 into three columns or a temporary table, and then update T1 using the update cursor.

Can anyone help me with this?

My code currently:

rows = arcpy.UpdateCursor(Flow_Paths)
for row in rows:
    X_A = row.X_Cord_A # First X Cord
    Y_A = row.Y_Cord_A # First Y Cord
    rows_i = arcpy.SearchCursor(point_shp)
    for row_i in rows_i:
        if X_A == row_i.xCentroid and Y_A == row_i.yCentroid: #Find upstream path lengths and extract maximum lengths, elevations, and min area
            ElevationMax = row_i.DEM
            AreaMin = row_i.CatchArea


texts: Are my reflectance values of bands in range and what does DRA in GIS 10.1 do to bands?

texts:
I have done the TOA correction to red and near infrared bands using raster calculator and calculated NDVI from lansat 8 image.

The NDVI range came as : upper range=184 and lower range=(-235).
When I click on the DRA check box it comes in a range of -1 to +1.

So my question is what does DRA do to the bands? Are my corrected reflectance range value correct?


texts: How to get a feature that lies closest to a point in openlayers?

texts: Im trying to perform routing with openlayers, geoserver and mysql. In order to get the start and end points, I need to be able to get a feature of the node layer (which has point features) that lies closest to an arbitrary point I click on the map. Does openlayers provide any functionality for this. If not, kindly suggest any alternative methods for achieving this.


texts: How can I get the correct mean from focal() with a circular window?

texts:
I am trying to perform focal statistics on a raster in R using a circular window (defined with focalWeight) and the focal function. In my specific case, I want each cell to be the mean of all cell values within a 100 m radius. Because I have NAs in my dataset, I followed the advice of Robert Hijmans in this post, which is to give weights a value of 1 and use fun = mean rather than the default sum function.

This seemed a perfect solution; however, the values produced are wrong. (This post describes a similar problem but unfortunately never received an answer.)

After much reading and testing I have come to realise that the problem is linked to the denominator used for the calculation of the mean. My question is, can I achieve the correct result without resorting to a clunky, slower custom function?

Here is a small example of the problem (demonstrated without NAs to keep it simple). Let's create a raster where all cells have a value of 1.

# Create a raster of 9 x 9 cells (res 1 m when projected in GB national grid)
r<-raster(extent(0,9,0,9), ncol = 9, nrow = 9, vals = 1, crs = "+init=epsg:27700")   # all the cells have value of 1

# Set circular window of 2m radius
fw <- focalWeight(r, 2, "circle")

Because there are no NA values, we can use the sum function and the weights. As expected, all cells (except edges) receive the correct mean value of 1.

focal_r <- focal(r, fw)

summary(values(focal_r))  # all values of 1

I would expect that changing the weights to values of 0 (outside circle) and 1 (inside circle) and using the mean function would yield the same result (as suggested in R. Hijmans' answer). However:

fw[fw > 0] <- 1   # replacing weights by 1

focal_mean_r <- focal(r, fw, fun = mean)

summary(values(focal_mean_r))  # now highest value is 0.52

The highest value we can achieve is 0.52, which I found out is the ratio of "1" cells to all cells in the focal window:

ncell(fw[fw > 0])/ncell(fw)

In other words, the function sums the raster values for cells falling in the "1" circular window, but divides the result by the total dimensions of the weights matrix. But surely this is not a desirable behaviour? The 0 cells (and NA cells for that matter, when na.rm is selected) should not count in the denominator. So is this an issue with the method? Or am I doing something wrong?

I even tried to coerce the focal window 0s to NA as suggested in method #4 of this blog post, but it doesn't work and I get a blank raster.

So my workaround is a custom function where I multiply the weights of 1 by the cell values, sum them all, and divide by the number of non-zero cells of the focal window:

custom_r <- focal(r, fw, fun = function(x){sum(x * fw, na.rm = TRUE)/ncell(fw[fw > 0])})

summary(values(custom_r))  # correct results with all values of 1

For the visual folks, you can see artifacts arising from the wrong calculation (left) when I view my raster in QGIS, while the output raster looks normal and smooth using my custom function (right). (NB: My large raster currently has only values of 1 or NA, although values could in practice range from 0 to 2.)

If anyone can confirm that this is the only way to go, I'll live with it. But the custom function runs about 3x slower than the built-in mean method, so if I'm missing something I would really like to know!


texts: Saving Bing QuickMapServices satellite layer without losing image quality

texts:
I'm using QuickMapServices to add a Bing satellite layer. but when trying to save what's in the Map view to a raster layer, I notice degradation in quality as if it's doing a compression process on the image. How can you save it without losing quality?

I'm using the 'Project > Save as Image'.

As I can't use any of the Raster options like Translate or Clipper etc as they require selection of a raster first.


texts: Get progression bar into PyQGIS 3 standalone script

texts:
When I run a treatment inside QGIS with Processing toolbox (native algorithms) I get a progression bar :

Is there a way to use this percentage inside a standalone PyQGIS script? I searched in the description of processing.run() but it seems nothing is linking to this progression bar.


texts:
I am trying to setup the Data Frame in QGIS to a minimum scale of 1:0 like ArcMap does. Zooming in goes down to 25.000:1

Any suggestions?


texts:
While adding items in the ArcGIS Online platform, an error appears with the message:

  Invalid URL

The file I am trying to upload is a shapefile I receive from a client. It had no defined projection so I exported the shapefile using ArcMap in WGS84, to make it compatible with the ArcGIS online platform. I have tried to upload it in Firefox, Chrome and Explorer, and the same message appears: Invalid URL.

The file is a .zip with the shapefile extensions (.cpg/.dbf/.prj/.sbn/.sbx/.shp/.shp.xml/.shx).

I have tried to "Add item" inside the Content page, and the Invalid URL message appears. When I try to upload it on My Map - "Add Layer From File", the message is "Unable to import Shapefile"

Any solutions?


texts: Creating buffers that will not overlap in QGIS

texts:
I am dealing with a .shp with a number of points in it. The distances between points range from ~300m to 5km.

I wish to create a 1 km buffer around each point with as little overlap as possible. Each buffer should have an ID. Where 1km buffers overlap with each other, I would like to give that buffer a different ID.

So lets say I have 50 points, and I create a 1 km buffer around every point. Owing to the short distances between my points in the Shapefile, overlapping between buffers are inevitable. Where there is overlap, I would like to give one of those buffers a particular ID. I would like to replicate this throughout my entire dataset until there are no buffers of the same ID overlapping.

I am sure this is possible in QGIS with a combination of the vector manipulation tools, I am just struggling on the best approach.

I am using QGIS 2.0.1 Dufour.


texts: QGIS Server 2.18 GetPrint doesn't wait for XYZ tiles

texts:
I am trying to generate a PDF output with a Print Composer using QGIS Server 2.18. The map includes a couple of layers from local shapefiles, and also an XYZ Tile Server layer.

QGIS Server does not seem to wait for the XYZ Tiles to load when generating an output. The Composer template loads correctly and so do all the layers based off shapefiles on the server's file system. But the XYZ tiles from the external site simply do not render at all. This happens in both the GetMap and GetPrint requests.

If I spam refresh on the GetPrint request several times, the XYZ do load. I believe that under the heavier load, the XYZ tiles do load in time to be included in the output. And once they have been loaded once, they seem to cache and will load on all subsequent requests.

Is there anything I can do to ensure that QGIS Server will wait for the XYZ tiles to load every time before returning the result of the GetPrint request?


texts: Shortening code that provides a database description using ArcPy?

texts:
I need advice on how to reduce the length of the sample code below. I need to perform the same work on several more datasets.

import arcpy, os

arcpy.env.workspace = r"\\gis\Env.gdb\Water"

datasets = arcpy.ListDatasets(feature_type='All')
datasets = [''] + datasets if datasets is not None else []

for ds in datasets:
    for fc in arcpy.ListFeatureClasses(feature_dataset=ds):
         path = os.path.join(arcpy.env.workspace, ds, fc)
        print(path)
############
arcpy.env.workspace = r"\\gis\Env.gdb\Soil"

datasets = arcpy.ListDatasets(feature_type='All')
datasets = [''] + datasets if datasets is not None else []

for ds in datasets:
    for fc in arcpy.ListFeatureClasses(feature_dataset=ds):
         path = os.path.join(arcpy.env.workspace, ds, fc)
        print(path)
###########
# need to repeat for several more datasets

I also have plans to take this further. I just wanted to clean up this segment as I go along. My goal is to create a text file that identifies, within a file or SDE database, all feature classes, their fields, and the data type and length of each field.


texts:
This is my first time using QGIS(2.14.18), I construct this map, but I don't know how to measure the distances (D1, D2, D3, Dx ...) between the circles and the stars. So I have 2 layers (based in 2 CSV files), with thousands and thousands of lines, the image below is just a test.

The distances I need are only between stars and circles, never between star to star or circle to circle. One star can have infinite circles, but a circle can only have one star link. In the two CSV files exist a key ID that stars and circles share with each other, but I don't know how to use it in QGIS to solve this problem.

I have tried using the distance matrix in QGIS, but I don't know if I'm doing this right. I think the distance matrix is giving redundant results, like 2 times the same distance or calculating the distances between stars and circles with no link between them.

I need these distances in meters or kilometers, not degrees.


texts: QGIS Layer creation failed for shapefile

texts:
QGIS v3.4.4

New project and trying to add a shapefile (line) layer w/ manual data entry.

  Layer>Create layer>New Shapefile layer

Fill out the form w/ a file name, geometry type and I added a Date field to the list w/ the default id field
Click OK
"Layer creation failed. Please check message log for further information"
Message log has time stamp, "WARNING:  (blank)"

What am I missing?


texts: Performing Statistics operation on field in Attribute table?

texts:
I am working with the Shapefile consisting set of points. One of its field (let's call it X) contains numeric data, ranging from 0 to 1. In order to find out the average value of field X, I opened Attribute table, navigated to the field X and clicked with the right mouse on it. Showing me the list of options, the option Statistics is not active, meaning I cannot click on it.

It is a text field holding double. When I tried to perform some filtering using this attribute X, the operation looked like X >= '0.99' for example. But it worked, it indeed performed the operation correctly.

Has anyone ever faced this situation?


texts:
I have a CQL filter that only displays a certain number of points, however when I try to apply the CQL to my existing leaflet map through the options for a layer nothing displays. No errors are thrown in the console.

Any ideas?

http://SERVER.com/geoserver/Workspace/wms?VERSION=1.1.1&REQUEST=GetMap&SERVICE=WMS&LAYERS=Workspace:Layer&STYLES=Point&BBOX=3026,2004,657587,1244444&SRS=EPSG:27700&WIDTH=570&HEIGHT=1035&FORMAT=image/png&TRANSPARENT=TRUE&CQL_FILTER=NAME='Leeds'

Below is what I have tried

 var FilterWMSlayer= L.tileLayer.wms("http://Server.com/geoserver/Workspace/wms",
    {
        layers: 'Workspace:Layer',
        format: 'image/png',
        transparent: true,
        CQL_FILTER:NAME='Leeds'

    });

   var FilterWMSlayer= L.tileLayer.wms("http://Server.com/geoserver/Workspace/wms",
    {
        layers: 'Workspace:Layer',
        format: 'image/png'&CQL_FILTER=NAME='Leeds',
        transparent: true
    });


texts: Google Earth Engine: Join FeatureCollection and List

texts:
I imported a csv called PADDYPOINTS containing the headers below as a FeatureCollection:
["town","paddy_field","section","longitude","latitude"]
Then, I extracted NDVI values on each point using the .getRegion() function
var PADDYPOINTS = ee.FeatureCollection('users/earthengine/paddy_coordinates');

var dataset = ee.ImageCollection('LANDSAT/LC08/C01/T1_8DAY_NDVI')
  .filterBounds(PADDYPOINTS)
  .filterDate('2017-01-01', '2017-12-31');

var poiArray = dataset.select('NDVI').getRegion(PADDYPOINTS, 30);

However getRegion() only returns coordinates: ["id","longitude","latitude","time","NDVI"] It doesnt return the other columns
How do I join the coordinates from the FeatureCollection above with the getRegion() result to also export the ["town","paddy_field","section"] columns together as one result?
P.S.: I cant seem to join the original csv with the result using python-pandas as the coordinates are not the same anymore after processing


texts:
I have a vector layer with 2498 features (houses) that I want to rasterize/produce a raster image from but can't figure out how:

Column MARK_Z = is Meters above see-level for the base of each house
Column Z = Total meters above see-level for each house
Column BYGG_H = is each buildings total height from base

So ideally I would want code two of them into the new raster attribute, but Raster > Conversion > Rasterize (Vector to Raster) only let's me chose one, so I'm looking at different ways to solve this.
I have a PostGIS up and running with Qgis, How can I rasterize this with PostGIS through Qgis?
More information:
Download this exact vector layer and see if you can sort it out?
https://drive.google.com/open?id=0BwFHJFrNX33hZ2E0UlEtZDlUOXc
Running QGIS 2.18.0 'Las Palmas' but could change if another version works better?


texts: Transformation problem WGS84 (EPSG:4326) and ETRS89 (EPSG:4258) in QGIS

texts:
I want to transform vector points from EPSG:4326 to EPSG:4258 and vice versa for survey application. Currently the difference between the two is about 70 cm in xy direction. Using QGIS to transform does not yield correct results (both CRS are basically identical in proj4).

I read about this problem but have not found a workable solution for whole layers. In my layers I have usually up to ten ground control points that need to be transformed.

Is there a solution out there that you could point me to?


texts: Updating table with nearest node from another table for pgRouting?

texts:
I have a roads table that has been transformed into nodes for use with pgrouting.

I also have a table of geocoded addresses which contains lat, lon and geom. I am trying to write a query that updates the addresses table with the GID of the nearest node.

Something like:

update addresses set nearest_node_id=(SELECT gid FROM roads
    ORDER BY geom <-> ST_GeometryFromText('POINT(addresses.lat addresses.lon)',4326)
    LIMIT 1  );

The Select gid ... portion of this works fine with string values subbed for addresses.lat & lon. But I can't figure out how to get the whole thing working.


texts: Python: color Canadian provinces in a shapefile based on given values

texts:
The goal here is to create a choropleth map of Canada in Python. Suppose I have a dictionary with values referring to each Canadian province/territory:

values={'Alberta': 1.0,
 'British Columbia': 2.0,
 'Manitoba': 3.0,
 'New Brunswick': 4.0,
 'Newfoundland and Labrador': 5.0,
 'Northwest Territories': 6.0,
 'Nova Scotia': 7.0,
 'Nunavut': 8.0,
 'Ontario': 9.0,
 'Prince Edward Island': 10.0,
 'Quebec': 11.0,
 'Saskatchewan': 12.0,
 'Yukon': 13.0}

Now I want to color each province based on the corresponding value in values, using a continuous colormap (e.g., shades of red). How to do that?

So far I have only been able to plot the Canadian provinces/territory within matplotlib, but their shapes appear in a unique color, and I don't know how to change that according to the numbers in values.

This is where you can find the shapefile: http://www.filedropper.com/canadm1

And this is my code to date:

import shapefile
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
#   -- input --
sf = shapefile.Reader("myfolder\CAN_adm1.shp")
recs    = sf.records()
shapes  = sf.shapes()
Nshp    = len(shapes)
cns     = []
for nshp in xrange(Nshp):
    cns.append(recs[nshp][1])
cns = array(cns)
cm    = get_cmap('Dark2')
cccol = cm(1.*arange(Nshp)/Nshp)
#   -- plot --
fig     = plt.figure()
ax      = fig.add_subplot(111)
for nshp in xrange(Nshp):
    ptchs   = []
    pts     = array(shapes[nshp].points)
    prt     = shapes[nshp].parts
    par     = list(prt) + [pts.shape[0]]
    for pij in xrange(len(prt)):
     ptchs.append(Polygon(pts[par[pij]:par[pij+1]]))
    ax.add_collection(PatchCollection(ptchs,facecolor=None,edgecolor='k', linewidths=.5))
ax.set_xlim(-160,-40)
ax.set_ylim(40,90)

This is the image I am getting so far:


texts: OSR Python testSR.ImportFromEPSG(4326) returns empty spatial reference

texts:
This is empty using Python .ImportFromEPSG(4326):

import osr
#this fails:
testSR = osr.SpatialReference()
testSR.ImportFromEPSG(4326)
print testSR.ExportToPrettyWkt()

#why did the import from EPSG fail above?
testSR.SetWellKnownGeogCS("WGS84")
print testSR.ExportToPrettyWkt()
GEOGCS["WGS 84",
DATUM["WGS_1984", and so on...


texts:
I have a very large shapefile with over a million parcels within it.

The problem I am having is that the polygons representing each parcel are very jagged.

What I am looking for is a simple and fast way of essentially "simplifying" these polygons and thus getting cleaner lines out of things.

Naturally, in doing so I need to ensure all parcels are still linked together so there are no spaces between the parcel boundaries.

What's the best way of achieving these objectives?


texts: Can a Carto web map auto-pan to follow point animation?

texts: Setting up an animated map in Carto Builder that follows a point across the Pacific. In order to see local route details the map must be scaled in, but then not all of the route is visible at once. Can Carto auto-pan to follow animations?


texts:
I'm trying insert a feature into a FeatureClass, but when I try to assign the geometry to featurebuffer the app return the next Error:  The remote procedure call failed. (Exception from HRESULT: 0x800706BE)

IWorkspaceEdit workspaceEdit = (IWorkspaceEdit)this._getWorkSpaceFromArcSDE();

workspaceEdit.StartEditing(true);
workspaceEdit.StartEditOperation();

IFeatureClass fc = this.getFeatureClassFromArcSDE(layerName);
IFeatureCursor pFeatureCursor;
IFeatureBuffer pFeatureBuffer;
IFeature pFeature;

QueryFilter qf = QueryFilter)_pServerContext.CreateObject("esriGeodatabase.QueryFilter");
qf.WhereClause = "";
pFeatureCursor = fc.Search(qf, false);
pFeature = pFeatureCursor.NextFeature();

if (pFeature != null) pFeature.Delete();

pFeatureBuffer = fc.CreateFeatureBuffer();
pFeatureCursor = fc.Insert(true);

**pFeatureBuffer.Shape = geometry;** //Here return the error

pFeatureCursor.InsertFeature(pFeatureBuffer);
pFeatureCursor.Flush();

workspaceEdit.StopEditOperation();
workspaceEdit.StopEditing(true);

Marshal.ReleaseComObject(pFeatureCursor);
Marshal.ReleaseComObject(pFeatureBuffer);
Marshal.ReleaseComObject(fc);


texts:
What is the difference between a delta table and a base table?
I am confused by these two words.

Delta tables

  The adds and deletes tables for a dataset are collectively referred to
  as the delta tables because they store changes made to the dataset.

Base tables

  The base table is the core table of a feature class. It contains all
  the nonspatial attributes and, if using a SQL geometry type, the
  spatial attribute as well.

  The term base table is used to differentiate this core table from
  other side tables, such as the delta tables, ArcSDE XML tables, or the
  f and s tables used by the sdebinary geometry storage type.

  When you look at a feature class through your database management
  system's user interface, you see the base table. For example, if your
  geodatabase contains a versioned feature class named prj_sites,
  you will find a table named prj_sites in your database. That table
  is the base table.

  Base tables are also called business tables.


texts: QGIS: why doesn't the "graduated" classification work for my column?

texts: I'm new to QGIS so bear with me. I need to show something on a map using the "graduated" classification. The thing is the program doesn't give me this option. I believe it's been reading the column as text instead of numbers, but how can I change that?


texts: R in QGIS Processing without installing GRASS

texts:
I have installed QGIS on Mac, using Homebrew:

  brew install qgis2 --with-r

from osgeo/osgeo4mac. I did not want Grass ( although I have laso tried --with-grass . In the processing toolbox, I went to activate the rscripts. So far so goo, but when hitting ok, QGIS complains that the path to Grass 6.4 is not found ( not sure why it does that). I cannot activate my R. ANy hints how to go about this? There is no option to delete the grass integration, and despite not activating grass, this complains.


texts:
I am pretty new to ArcGIS for Server and ArcSDE but cannot find a solution in the Esri Documentation or Forums.

I don not know how to use sdemon or even where to start this command line tool.


texts:
I want to enable read support from Erdas ECW files in gdal on windows.

I have GDAL binaries from http://www.gisinternals.com/sdk/ and i also have a self compiled version.

when compiling myself I have found the options in nmake.opt and could add support myself, but if its possible I would rather add it to the precombiled distribution from gisinternals.  Is this possible? or will I need to compile it my self?


texts: Averaging MODIS 8-day LST for each pixel in a region to generate annual mean LST

texts:
I have downloaded 8-day MODIS LST products (MOD11A2) for the whole year of 2011  and is there any method in R to extract those pixels in clear sky conditions with  emissivity error < 0.01 and LST error < 1 K and average LST values for each pixel to generate only one raster map for 2011? My current code is not complete as below.

    lst <- lapply(list.files("D:/LST2011/Surf_Temp_8Days_1Km_v6/LST_Day_1km",
                     pattern = "LST_Day_1km", full.names = T), raster)
    csd <- lapply(list.files("D:/LST2011/Surf_Temp_8Days_1Km_v6/Days_Clear",
                     pattern = "Days_Clear", full.names = T), raster)
    mean_lst <- temporalComposite(lst, csd, fun = function(x) mean(x, na.rm = TRUE))


texts: QGIS installation has stopped running when startup?

texts:
My laptop is Windows 10 and I have installed QGIS 2.14 64bit twice; because due to some problems I uninstalled the software and reinstalled it. Now: "qgis has stopped running when startup"

I have searched for solutions on GIS SE and here is what I have done:

I delete the .qgis2 but it fails;
I tried to run qgis-ltr-bin.bat as "gis.stackexchange.com" guides but it fails;
I tried to uninstall and delete everything about QGIS in C: drive as "gis.stackexchange.com" guides but it fails;
I tried to install other software like 2.18 64bit but it also fails.


texts:
i have posted once a question about not getting the "get more" icon within the plugin manager tab.
I solved my problem by specifying the proxy server that supports the network to which my PC (windows7 OS) is linked.
However, a very odd thing happened while i was giving a training course on QGIS. Several of the attendees used their own laptop (all windows) and none of them got the get more icon within the plugin manager.
All of them downloaded QGIS Dufour from the QGIS website and followed the straightforward installation process. Why this bizarre thing is happening? How can i solve the problem?


texts: Basing an 'ELSIF' expression on several fields

texts:
I am still fairly new to coding so please excuse what may be simple mistakes.

I am working in ArcMap 10.5 and have an Overhead Conductor feature layer that I need to label based on what field  (Conductor A, Conductor B, Conductor C, Conductor N, Construction Statuse) says what.

For example:

If Conductor A contains a value and Conductor B and C contain 'N/A' then the label for Overhead Conductors should display what is in the Conductor A field.

If Conductor B contains a value and Conductor A and C contain 'N/A' then the label for Overhead Conductors should display what is in the Conductor B field, and so on.

HOWEVER, the twist comes in when 'Construction Status' field contains a 'Proposed' status instead of an 'Installed' status. When this field has a 'Proposed' status I would like the label to be displayed in a larger, bold font.

I have included the original labeling expression and the expression of when I added what I hoped would make the font appear larger and bold when a 'Proposed' status was displayed in the 'Construction Status' field. I have also included a screen shot of how the fields appear in the attribute table.

The label expression format I am using is done in VBScript.

If anyone knows of an easier way to accomplish this in a different script I am open to suggestions but I am unfortunately stuck due to my lack of coding experience.

ORIGINAL EXPRESSION

Function FindLabel ( [conductora], [conductorb], [conductorc], [conductorn] )
dim a
if [conductora] <> "" and [conductora] <> "NA" then
a = a &  [conductora]
elseif  [conductorb]  <> "" and [conductorb] <> "NA" then
a = a &  [conductorb]
elseif  [conductorc]  <> "" and [conductorc] <> "NA"then
a = a &  [conductorc]
elseif  [conductorn]  <> ""and [conductorn] <> "NA" then
a = a &  [conductorn]
end if
FindLabel = a
End Function

MODIFIED EXPRESSION

Function FindLabel ( [conductora], [conductorb], [conductorc], [conductorn] )
dim a

if [conductora] <> "" and [conductora] <> "NA" then
a = a &  [conductora]
elseif  [conductorb]  <> "" and [conductorb] <> "NA" then
a = a &  [conductorb]
elseif  [conductorc]  <> "" and [conductorc] <> "NA"then
a = a &  [conductorc]
elseif  [conductorn]  <> ""and [conductorn] <> "NA" then
a = a &  [conductorn]

elseif [CONSTRUCTIONSTATUS] = "Proposed" AND [conductora] <> "" and [conductora] <> "NA" then
a = a &  [conductora]
FindLabel = "<BOL><FNT size = '12'>" + a &  [conductora]  + "</BOL></FNT>"
elseif  [CONSTRUCTIONSTATUS] ="Proposed" AND [conductorb]  <> "" and [conductorb] <> "NA" then
a = a &  [conductorb]
FindLabel = "<BOL><FNT = '12'>" + a &  [conductorb] + "</BOL></FNT>"
elseif [CONSTRUCTIONSTATUS] ="Proposed" AND [conductorc]  <> "" and [conductorc] <> "NA"then
a = a &  [conductorc]
FindLabel =  "<BOL><FNT = '12'>" + a &  [conductorc] + "</BOL></FNT>"
elseif  [CONSTRUCTIONSTATUS] ="Proposed" AND [conductorn]  <> ""and [conductorn] <> "NA" then
a = a &  [conductorn]

end if
FindLabel = a

End Function


texts: Problem loading GeoJSON file on geojson.io

texts:
I am newbie in GeoJSON file. I am having difficulty to load my GeoJSON file on geojson.io. Here is my GeoJSON file:

https://gist.github.com/anonymous/0abb8f926977d577b52d.

The source file is from:

http://data.gov.sg/common/search.aspx?cs=1&q=Region+Census2010&a=AVA,HPB,IDA,IPOS,IE,MINDEF,SCDF,SPF,NHB,PUB,SLA,SSC,STB,WDA,URA&da=CPF,CAAS,EDB,EMA,HSA,IRAS,LTA,MPA,MCCY,MOE,SING,MFA,MOH,ICA,SPS,MSF,DOS,MOT,MAS,NYP,NAC,NLB,NPAR,RP,SP,SPIB&c=12,03,02,07,09,04,08,05&dc=01,10,06,11&f=SHP&df=CSV,KML,PGDB,TXT,OTHERS,XLS,XML&page=1

named RegionalCensus2010.

It is a shapefile. I convert the shapefile to GeoJSON using QGIS. I don't have any issue when I am using GeoJSON file from other resources like US government file, but somehow I just can't load my GeoJSON file. I have checked the structure of the file which seems to be correct but I am not sure why I can't load map layer on geojson.io.


texts: Unknown Coordinate system for San Francisco LiDAR LAS Files! How it is possible to assign a coordinate system to LAS file?

texts:
I have downloaded all the publicly available LiDAR data from https://viewer.nationalmap.gov/basic/ for San Francisco area. Files cannot be converted since the coordinate system is unknown in all the data covering at least the Berkeley area.

How I can assign the coordinate system (which theoretically is not possible because I have to first convert them and then assign the coordinate system) or solve the problem?


texts:
I am a beginner with computer programming  and I trying to use Python pack, sentinel-sat 0.12, to download images of Sentinel 2, but I get back a response I don't understand:

  sentinelsat.sentinel.SentinelAPIError: HTTP status 500 Internal Server
  Error: IOException : Server returned HTTP response code: 400 for URL:
  http://localhost:30333//solr/dhus/select?q=beginposition:[2017-08-01T00:00:00Z%20TO%202017-08-30T00:00:00Z]%20cloudcoverpercentage:[0%20,%20100]%20platformname:Sentinel-2%20footprint:"Intersects(POLYGON((16.6085%2038.8957,16.6221%2039.4901,17.2206%2039.4802,17.2020%2038.8860,16.6085%2038.8957)))"&wt=xslt&tr=opensearch_atom.xsl&dhusLongName=Sentinels+Scientific+Data+Hub&dhusServer=https%3A%2F%2Fscihub.copernicus.eu%2Fdhus%2F&originalQuery=beginPosition%3A%5B2017-08-01T00%3A00%3A00Z+TO+2017-08-30T00%3A00%3A00Z%5D+cloudcoverpercentage%3A%5B0+%2C+100%5D+platformname%3ASentinel-2+footprint%3A%22Intersects%28POLYGON%28%2816.6085+38.8957%2C16.6221+39.4901%2C17.2206+39.4802%2C17.2020+38.8860%2C16.6085+38.8957%29%29%29%22&rows=100&start=0&sort=ingestiondate+desc

The script I have used is:

from sentinelsat.sentinel import SentinelAPI, read_geojson, geojson_to_wkt
api = SentinelAPI('**********', '*********','https://scihub.copernicus.eu/dhus')

# search by polygon, time, and SciHub query keywords
footprint = geojson_to_wkt(read_geojson('KR_REFEFA_POLY_extent_4326.geojson'))
products = api.query(footprint,
                     date = ('20170801', '20170830'),
                     platformname = 'Sentinel-2',
                     cloudcoverpercentage = '[0 , 100]')

The GeoJSON file SRC is in WGS84 (EPSG:4326).

I hope can you help me with this problem.


texts:
I tried combining the code in the example below, with my own code but it does not work: the timeslider continues to count up from 0 and does not display the text alongside the number.

How to Change Torque Map Slider Text

My code is shown below.

What I really want is for my time slider to show the values of one of the columns in the data uploaded to Carto as it plays, preferably with the label "Years B.P.".

But for now all I want to know is why the code from the example does not work with my code?

<!DOCTYPE html>

    <title>javascript_slider_testing_V2</title>
    <meta name="viewport" content="initial-scale=1.0, user-scalable=no" />
    <meta http-equiv="content-type" content="text/html; charset=UTF-8"/>
    <link rel="shortcut icon" href="http://cartodb.com/assets/favicon.ico" />

    <link rel="stylesheet" href="http://libs.cartocdn.com/cartodb.js/v3/3.15/themes/css/cartodb.css" />


texts: Evaluate minimum points required to define polygons based on an error threshold

texts:
I need to build a model or write an algorithm to do this.

here is my cadastral problem: I have plenty of polygons and I need to make a report of coordinations of their vertices. This is pretty easy: I can extract vertices using ArcGIS toolbox and add XY coordinate to each point again using ArcGIS toolbox. However, since each polygon has 20-30 vertices, I need to condense the information. I need to make reports of the most important vertices which define a polygon.
The problem here is that these are polygons with irregular shapes. Some can be defined by using 4 points and some require 10 for example. Moreover, my datasets are huge and I definitely need an automated approach to this.

Here is what I thought I have to do:

1) Polygone to vertices > add XY to vertices

2) for each set of polygon vertices use an iterative approach as follows:

2.1) start with 4 vertices

2.2) calculate area

2.3) evaluate area error percentage (area error/correct area = abs(real area-calculated)/correct area)

2.4) if area error < desired threshold (say 0.05)

then end step 2, go to next polygon

else  add one vertex, go to 2.2

3) Make a report of selected vertices

However, I have a problem with this. How should I select the correct vertices in step 2.1 and 2.4. For example, a polygon having 20 vertices, may not be defined using 4 of its vertices, while it can be very well defined using another 4 vertices.

I have access to Arcinfo+python, Xtools, and Matlab+Mapping toolbox. I don't know if I have to use Matlab for this problem or not, since I doubt the efficiency of ArcPy in calculation of new polygon areas.


texts: Is it possible to change legend symbols in the QGIS print composer legend?

texts:
I have some layers with very complex data-defined symbols.
If I use a legend in the print composer these symbols are not drawn correctly.

Is there a way to change the legend-symbols with pyqgis, so I can use custom PNG or SVG images instead of the default legend symbol for these layers?

I know how to add a Push Button to the Print composer and how to connect it with a function.
I want to add a button to the Legend-Settings to allow me to replace the automatically generated legend icon with a custom image.
So what I still need is the information of how I can access the legend-symbols with pyqgis/pyqt, and how to replace them with a QImage on a QLabel or something like that?

Very basic mockup of the button used to change the symbol:

Automatically generated Legend:

Legend with custom legend symbol:

I already found out how to access the items of the legend in the print composer but not yet how to access the symbol itself:

import qgis
from PyQt4.QtCore import *
from PyQt4.QtGui import *

activeComposer = iface.activeComposers()

for item in activeComposer:
    if item.composerWindow().windowTitle()=='test':
        for i in item.items():
            if isinstance(i,QgsComposerLegend):
                #print i
                #print i.model()
                legend = i
                for i in xrange(legend.modelV2().rowCount()):
                    posteleg=legend.modelV2().index(i, 0)
                    print posteleg.data()

Edit 1:

I have also found out how to Access QIcon-Objects of the legend-tree but cannot swap them yet:

def run(self):

        activeComposer = self.iface.activeComposers()
        #print(self.resolve('icon.png'))
        for item in activeComposer:
            if item.composerWindow().windowTitle()=='test':
                for i in item.items():
                    if isinstance(i,QgsComposerLegend):
                        legend = i

                        layerIcon = QIcon(os.path.join(os.path.dirname(__file__), "icon.png"))

                        for i in xrange(legend.modelV2().rowCount()):
                            posteleg=legend.modelV2().index(i, 0)
                            posteleg.model().iconGroup().swap(layerIcon)
                            print posteleg.data()

Here is one real life example where you can see a symbology combined of many symbol layers:

This will end up in the legend like this:

As I need the proper symbol in the legend I would like to make a screenshot of my symbol, crop it and use it as an image in my legend.

I know I could just overlay a seperate image on top of my legend which covers the automatically generated symbol but I would like to have a "cleaner" solution that allows me to replace the symbols in the legend with custom images.

Edit 2:

In the meantime I have found another way to gain access to the legend entries:

from qgis.core import QgsLegendRenderer, QgsComposerLegendStyle

compDict = {}
for comp in iface.activeComposers():
    # workaround to get name: read it from window title
    compDict[comp.composerWindow().windowTitle()] = comp.composition()
if "mycomposername" in compDict:
    itemLlegend = compDict["mycomposername"].getComposerItemById("mylegend_id")
    if itemLlegend:
        print itemLlegend

tree_layer_layer =  itemLlegend.modelV2().rootGroup().children()
for item in tree_layer_layer:
        if item.layerName()=="MyLayername":
            print "match"
            QgsLegendRenderer.setNodeLegendStyle(item, QgsComposerLegendStyle.Group)

This allows me to access the QgsLayerTreeLayer objects and I can switch the legend Style (Group,Subgroup,Hidden).
But I still have no idea how to access the legend-symbol.
Any ideas?


texts: Finding number of trips in each segment of road

texts:
I have street line layer for a city. I also have a lines travelled by the vehicle, which was created from the GPS points collected in the study (OBVIOUSLY, the path of vehicles are not accurate, so they dont fall exactly in the road line). I want to know the number of times the road is traversed, or, i want to know the number of trips for any segment of road.

I am using Arc GIS. I do not have knowledge of python, but could learn and edit if there is code available for similar tasks.

For example, if one road segment between two blocks is used twice, then i want one column on attribute table of the road segment to have a value of 2.


texts: Versioning with Oracle Workspace Manager in QGIS?

texts:
I'm working with a workspace enabled oracle database. Through the database connection properties in QGIS meanwhile it is possible to set the workspace or version. So I can edit my spatial data in a given version.

Is it also possible to write a new workspace/version with QGIS?

I want to edit my spatial data and save them into a new workspce/version of my oracle database.


texts:
i use this example for doing my own wfs-T example. . And my code are same under link. My only problem is i don't change polygon default style. it gives error when i use openlayers style map. How can i change polygon style , i want to change color ,opacity,etc of polygon , please help?

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


texts: Building CSW query to discover available products in GeoNetwork

texts:
My task is to implement data discovery and download starting from a CSW catalogue based on GeoNetwork.  I got into the GetCapabilities of it but I do not understand how to move from the metadata catalogue of CSW to the download operation.

Which is the CSW type of query to build to discover the available products?


texts: Improving display speed of a virtual layer in QGIS

texts:
Following on this question.
I have a set of reference lines and a non geometrical table of events. This table contains the mileage of the event along the reference lines, amongst other attributes.
I upscaled my application to use real world data. It works, but the display is slow.
I have read about spatial indexes, but I am not sure about how to make it work with a non spatial table.
My virtual layer query is the following :
SELECT
    e.fid /*:int*/,
    e.'LE_idRiver' AS RouteID /*:int*/,
    r.'LB_idRiver' AS RivName /*:text*/,
    e.'LE_Value' AS val /*:int*/,
    ST_Line_Substring(r.geometry, e.'LE_FROM' / st_length(r.geometry), e.'LE_TO' / st_length(r.geometry)) as geom /*:linestring:2056*/
FROM  "Linear_events"as e
Join    "base_lines" as r ON (e.'LE_idRiver' = r.'fid')

The baselines are quite long (up to 50km) and some events are short (~10m), so there are quite a few of them. I was thinking about splitting the baselines every km and use the processes described here but I would rather have a single geometry if possible.
Moreover, my hypothesis is that the whole layer is redrawn even when I zoom in, meaning that this approach would lead to an improvement at a lower scale, but the gains would vanish at a larger scale, leaving the problem unsolved.
Is there a way to improve the display speed of my layer at all zoom levels in this use case?


texts: How to Remove Multiple Labels for Stacked Points? ArcGIS 10.2.2?

texts:
I have a point feature class, in certain areas there are stacked points. These stacked points are accurate. However, when I label the points I am getting a label for each point, which is okay for each unique attribute that is being displayed. However, there are some points that have the same attribute information being displayed, where I would only like one label.

Example:

How can I display one label for all stacked points if their attribute information is the same? So that it would look like:

Any thougts? I am using a point feature class from ArcSDE 10.2.2.

I cannot alter that data in any way and I am required to used this dataset. Also, I will be publishing it as a dynamic mapping service and documentation has suggested that Maplex Label Engine is not recommend since it slows down performance.

I am looking for a work around taking account into these requirements.


texts: ArcMap modelbuilder CalculateValue precondition not working

texts:
I am using ModelBuilder in ArcMap 10, this is my model:

Cell Statistics shall be calculated for only those rasters dating in summer months (April - October), input raster name format is yYYYYmMM, e.g. y2007m10. Therefore my precondition as boolean value:

The model is working but calculates cell statistics for all input rasters, not only for the summer months. Probably my code is incorrect, or are the elements wrongly ordered?


texts: Show elevation profile from linestring's z values

texts:
I have linestring geometries with z values and I want to display them in a nice elevation profile (with x/z scaling functions, several profiles to overlay, and so on).

Is there any QGIS function or plugin for it?

I've seen tools for digitizing elevation profiles from DEM's, e.g. Profile tool. These tools have such a height profile display, but they seem to require the height digitization I don't need.

I could go to Excel, yes. But I'd like to stay inside QGIS in order to avoid another interface.


texts:
I have a task I would like to solve with GRASS GIS 7, if somehow possible.
I want to create a cross section from a landslide - and would like to calculate the area of that cross section/profile. Any way to do this easily without using ArcGIS?


texts: Summarizing Fields based on content in QGIS?

texts:
already read many answers for similar questions, but i didnt find a solution for my problem. Dont know how to start.

Problem: shp with about 50 fields. Each field contain NULL, 1 or 2. Now i need two new fields. One containing the number of fields with "1", the other with "2".

I am using QGis 2.16.2 and 2.14; OS Win10

I'm trying @Joseph's solution, add the hint @kyle reed and make some modifications to start as an expression:

@qgsfunction(args='auto', group='Custom')
def R1func(feature, parent):
    lists = []
    layer = qgis.utils.iface.activeLayer()
    for feat in layer.getFeatures():
        f = feat.fields()
        num = f.count()
        for i in range(num):
             if feat[i] == 2:
                 lists.append(i)

    return len(set(lists))

and also try a short version

@qgsfunction(args='auto', group='Custom')
def R2func(feature, parent):
    lists = []
    layer = qgis.utils.iface.activeLayer()
    for feat in layer.getFeatures():
        lists.append(feat.fields())

    return lists.count(2)

Version 1 returns for every row 60 (for number 1) or 37 (for number 2)
Version 2 returns 0

Maybe i make a mistake.


texts: How to check if a layer exists in a geodatabase through ArcToolBox or Python?

texts:
I have got a python to update a file geodatabase (*.gdb).

Every time the code executes, it needs to check if a certain layer has been already there in the previous geodatabase; if so, delete it and create a new one; otherwise, overwrite it with a new one.

I am just wondering if I can achieve this by ArcToolBox or a python code?

Thanks

Alex


texts:
I am sorry if this is a relatively dumb question, but I was wondering if there is a way to store a SketchUp file (.skp) in a PostGIS database?  I was thinking about as a blob, but I wasn't sure if anyone has had success with this.  The record would be a point location.  Right now I am doing a file path link to the file on the directory, but I was hoping there was a way to store the file in the database.  I would like to use this for data storage, but also as a means to retreive and view the data as well.

Thanks...


texts: Keys for merging US state/county ANSI's with GADM Unique identifiers?

texts:
I've got county-level data from USDA NASS, and I want to merge it with gridded climate data.  A preliminary step is the associate the counties with the coordinates of their centroids.  If I merge by names, about 5% of the counties are left out, because of different representations of names.  For example, the city and county of Chesapeake, Virginia is listed as "Chesapeake" in GADM and "Chesapeake City" in the NASS data.

The NASS data has ANSI codes, while the GADM data has its own scheme.  Is there some sort of concordance table or dictionary that allows me to link these two?

Here is an example of the NASS data:
http://quickstats.nass.usda.gov/results/E03CA5B5-AE87-3C9D-B8E7-662F85AD3CC0

And here is an example of the GADM data:

2> library(raster)
2> adm <- getData('GADM', country='USA', level=2)
trying URL 'http://biogeo.ucdavis.edu/data/gadm2.7/rds/USA_adm2.rds'
Content type 'unknown' length 13936198 bytes (13.3 MB)
==================================================
downloaded 13.3 MB

2> head(adm)
  OBJECTID ID_0 ISO        NAME_0 ID_1  NAME_1 ID_2  NAME_2 TYPE_2 ENGTYPE_2
1        1  244 USA United States    1 Alabama    1 Autauga County    County
2        2  244 USA United States    1 Alabama    2 Baldwin County    County
3        3  244 USA United States    1 Alabama    3 Barbour County    County
4        4  244 USA United States    1 Alabama    4    Bibb County    County
5        5  244 USA United States    1 Alabama    5  Blount County    County
6        6  244 USA United States    1 Alabama    6 Bullock County    County
  NL_NAME_2 VARNAME_2
1        NA      <NA>
2        NA      <NA>
3        NA      <NA>
4        NA      <NA>
5        NA      <NA>
6        NA      <NA>

The county ANSI's don't appear to correspond to the GADM's ID_2 field.


texts: How to create a 3D shapefile from a raster?

texts: I´m trying to make a 3D shapefile from a raster feature but not working. Tried already to run the command in OSGeo4W but nothing happen. I'm using the Qgis 1.8.0 Lisboa. The purpose is to get a final 3D dxf file! Can anyone help me?


texts: Leaflet error: h.fn is undefined when using layer.on('loading') event

texts:
In my Leaflet application (version 1.0.1 1.0.3) I connect to ArcGIS Server via the Esri Leaflet plugin. I load the layer as dynamic layer:

var layerUrl = "http://localhost:6080/arcgis/rest/services//service1/MapServer";
var layer1 = L.esri.dynamicMapLayer({
   url: layerUrl,
});

I only load this layer when a checkbox is clicked:

if (checkBox.checked){
    map.addLayer(layer);
}

Now my goal is, to show the cursor as 'wait' during the layer is loading. So I add this line:

layer.on('loading' , cursorWait());

The cursorWait function is this:

function cursorWait() {
    document.body.className = 'wait';
}

CSS:

.wait, .wait * { cursor: wait; }

Which works great, but now I get error messages on my console:

TypeError: h.fn is undefined  leaflet.js:5:4965
    fire                     file:///C:/folder1/Leaflet/LeafletLib/Leaflet/leaflet.js:5:4965
    _renderImage             https://cdn.jsdelivr.net/leaflet.esri/2.0.0-beta.7/esri-leaflet.js:5:10474
    _requestExport/<         https://cdn.jsdelivr.net/leaflet.esri/2.0.0-beta.7/esri-leaflet.js:5:17673
    _createServiceCallback/< https://cdn.jsdelivr.net/leaflet.esri/2.0.0-beta.7/esri-leaflet.js:4:24163
    bound                    self-hosted:916:17
    createRequest/httpRequest.onreadystatechange https://cdn.jsdelivr.net/leaflet.esri/2.0.0-beta.7/esri-leaflet.js:4:9620

My question is: does someone know why these error messages occur?
What really bugs me, that they appear only if I add the layer.on('loading',..) line, but then it appears also when I zoom in or out, although I am calling this only when the checkbox is checked. Also it seems to print an error message for each tile (?).
Am I using the layer.on('loading',..) event wrong? Or is it a bug? Should I just ignore the error messages?


texts:
I'm trying to read in a shapefile and keep getting a File Not Found error and this exception:

  ShapefileException: Unable to open /Users/name/Documents/shapefile.shp

Here's what I've got:

import shapefile

# read the shapefile
reader = shapefile.Reader('/Users/name/Documents/shapefile.shp')

I've tried several variations of that and still get the error. The shapefile exists and it works. It's also in the same folder as the dbf and shx files. Any ideas?


texts: How to request group of layers from Geoserver with turn off/on capabilities? Geoserver layer organizer

texts:
I would like to organize my layers by topics, and for each of them (if chosen on the client-side), I would like to request a single-tile image, which contains the predefined layers in the correct rendering order (e.g. if one topic has 5 layers, it renders them from the bottom to top like this: water, roads, highways, railways, poi). But, I also want to be capable to request just four of it, or three, etc. in any combination but still in the correct order.

My first idea was Geoserver's "Layer group", but as I've learned that is really more suited for group of layers which are not meant to be turned on/off seperately. (Please let me know if I am wrong about this, though.) And I also kinda lose Getfeatureinfo capability with this.

After some research I figured out, that in the GetMap request, if I add all my layers' name in the correct order (from bottom to top-->LAYERS=water,roads,highways,railways,poi), I get the same result, and I can also narrow down the list to the required ones.

The GetMap solution seems to by right, but in that case I need to fill up the list for each topic on the client side.

So my main question: is there a solution to get a list of layers from Geoserver by topic in the right rendering order? Or is there an extension for organizing layers this way? (I haven't seen any extension for that purpose)

My idea was that if I organize all layers of a topic to a workspace, I can request a getCapabilities for the workspace WMS endpoint, parse and read out the layers, but unfortunately it gives back the layers in a pseudo-random order, so not really good. Maybe a custom getCapabilities XML can solve this?

Please share your ideas, what's the best practice to do this?

I'm on Geoserver 2.4.2.


texts: Getting List of Selected Rows for Specific Layer using Layer Object in ArcPy?

texts:
What is the fastest/cleanest way to get a collection of rows that are selected in a certain layer? For some background, I am trying to hide selected features in a specified layer.

The solution involving the use of SearchCursor uses the layer name, but my setup may involve duplicate named layers
Another solution I found copies the selected rows into a new layer, which may be time costly for my setup.

Is there a way to get the selected rows by using a layer object which is provided by a script parameter?


texts: Including weighted values in zonal statistics using stand-alone PyQGIS script

texts:
I am having trouble to include the (spatial) weight of a pixel when using zonal statistics in a stand-alone python script (jupyter notebook 6.0.3 on Ubuntu 18.04.5).
My data contains of a .tiff (values between 0 and 1) and a shapefile overlapping several pixel. When computing the "mean", zonal statistics only accounts for pixel whose centroids are within the shapefile. However, I would like the pixel values contribute to the mean according to their weighted area instead of this binary approach ("in or out").
I did some research, this is exactly what I am looking for. Unfortunately, it is meant for PyQGIS in QGIS and I have not managed to implement it in my stand-alone script (I am stuck with "ModuleNotFoundError: No module named 'qgis'").
Has someone successfully applied this approach or found any way outside of QGIS?
If necessary, I can provide my code and screenshots.
EDIT 1:
Assuming that the two right pixels in my example .tiff have the values 1 and 0.75, the two left ones 0.25 and 0. Zonal statistics calculates the "mean" like this:
(1 + 0.75 + 0.25 + 0)/4 = 0.5
as all centroids are within the red shapefile (frame). Let's assume that the two right pixels are only 75% within the shapefile. I would like then the mean to be computed as
(0.75x1 + 0.75x0.75 + 0.25 + 0)/4 = 0.390625.
How can I accomplish this within Python? I have numerous shapefiles with different outlines, so I would like to automatize the mean computation process.


texts:
I'm looking for something exactly like GeogTools, but the download link is not working.

GeogTools is an Excel add-in that provides specialized functions for:

Converting latitude and longitude to UTM coordinate
Converting UTM coordinates to latitude and longitude
Converting U.S. State Plane coordinates to latitude and longitude*
Converting to and from Minnesota and Wisconsin County Coordinate Systems

* Conversion from U.S. State Plane coordinates to latitude and longitude (using NAD83) uses an adaptation of the National Geodetic Survey's SPCS83 conversion program.


texts:
In PCI Geomatic, I have two raster layer.raster1 identified like %1 and raster2 like %2.

if (%1 <= 10) then
   %2 = 60
elseif (%1 > 10 and %1 <=50) then
  %2=80
elseif (%1 > 50 and %1 <=90) then
  %2=120
else
  %2=0
endif

Please, How can I do this in QGIS?


texts: calculating the difference between 2 angles using st_azimuth and dot product

texts:
the function below tries to calculate the difference between 2
angles. I use the mathematical reasoning is as follows :

1) Obtain the 2 angles :

angle1 = ST_Azimuth(ST_StartPoint(l1), ST_EndPoint(l1))
angle2 = ST_Azimuth(ST_StartPoint(l2), ST_EndPoint(l2))

2) Obtain vectors (of lenght == 1) from 2 angles :

V1 = (cos(angle1), sin(angle1))
V2 = (cos(angle2), sin(angle2))

3) Compute the dot product :

V1 * V2 = V1[x] * V2[x] + V1[y] * V2[y]
=
cos(angle1) * cos(angle2) + sin(angle1) * sin(angle1)

Since the lenght of the vectors are both 1, it follows that :

the angle between the two vectors = ArcCos(dotProduct)/(length(v1)*legnth(v2)

= ArcCos(dotProduct)

The problem in the function below is that the dot product is yielding
values greater than 1, which should not be possible
if the math reasoning is correct, and that is my question :

is the reasoning correct ?

CREATE OR REPLACE FUNCTION angleDiff(l1 geometry, l2 geometry)
RETURNS FLOAT AS $$
DECLARE angle1 FLOAT;
DECLARE angle2 FLOAT;
BEGIN
  SELECT ST_Azimuth(ST_StartPoint(l1), ST_EndPoint(l1)) into angle1;
  SELECT ST_Azimuth(ST_StartPoint(l2), ST_EndPoint(l2)) into angle2;

  RETURN acos(cos(angle1) * cos(angle2) + sin(angle1) * sin(angle1));
END;
$$  LANGUAGE plpgsql


texts: How to obtain X & Y values of a projected point shapefile?

texts:
I know how to add X and Y coordinate values of a shapefile to the attribute table.

Now I would like to know how to achieve the coordinate values (in m) from a projected shapefile?

I am working with QGIS but I also have ArcGIS at my disposal.


texts:
I need to run my code through a lot of rasters that i've downloaded with "Download Raster" function, but some rasters have no information, so my script fails. I wonder if there's any way to use an if/else statement to solve my problem, in other words, if there's statsitics in the image, run the code, else, ignore the raster and select the next one.

I tried this, but it's not working...

for tiff in rasters:

    NodataResult = arcpy.GetRasterProperties_management (tiff, "ANYNODATA", "")
    Nodata = NodataResult.getOutput

    if Nodata = 1

        CalculateStatistics_management (tiff)

    else

        arcpy.MakeRasterLayer_management(tiff, "Z:\\04 - Temp\\Charles\\NDVI_tempo\\NAL\\temp_rast1.tif", "", "", "3")
        arcpy.MakeRasterLayer_management(tiff, "Z:\\04 - Temp\\Charles\\NDVI_tempo\\NAL\\temp_rast2.tif", "", "", "4")


texts: Georeferencing Cadastral map in DWG using AutoCAD Map 3D?

texts:
I want to georeference a cadastral map in dwg format. I have the lat-long of 3+ points on the map but do not have the scale.

I tried "Set Location" in Autocad Map 3D but it references using the lat/long & direction for only one point and then one has to scale the map. I do not have a scale.

How do I georeference the map without a scale?


texts: Querying geodesics (shortest paths) using sDNA

texts:
I am querying geodesics with sDNA 4.0.2.

Is there a way to output a set of geodesics based on specified origin-destination pairs (specified with polyline/link IDs)?

Right now when I use the "Geodesics" tool and enter the Origin IDs and Destination IDs, it will output geodesics of all possible combinations of the input origins and destinations. This is not very efficient computation- and storage-wise, because I only need to know a subset of those geodesics.

For each geodesic, I would like to get its length and the IDs of the polylines/links that comprise the geodesic. I think the program has already computed all the shortest path between all possible origins and destinations, but I am not sure if that information is available for the end users to play with in the ArcGIS tool interface.


texts: Polygonizer plugin for QGIS halts in Windows?

texts:
I've had some success and some problems using the Polygonizer plugin in the past.  I'm using Polygonizer 2.1 in QGIS 1.8 on Windows 7.  I'm aware there are sometimes problems using Polygonizer in Windows, but I was wondering if any solutions or work-arounds had come up since I researched the issue the last time.

The linear layer I want to convert at present is fairly complex and consists of 18693 lines.  When I use the "new method", it reaches 44% completion and just stops.  I've had it running for about 3 hours now at 44%.  Only rarely do I get an error message.  When I use the old method, it is extremely slow, completing about 1% per hour.  Last year I was able to run Polygonizer successfully on an older computer with Windows 2000, but that machine has since died!  Does anyone have any suggestions?

The rarely appearing error message reads:

  An error has occured while executing Python code:

  Traceback (most recent call last):
    File "C:/Users/George/.qgis//python/plugins\Polygonizer\PolygonizerDialog.py", line 103, in threadFinished

  msg = QMessageBox.question(self, 'Polygonizer', 'Polygonization finished in %03.2f seconds. \n %d polygons were crested. \n Load created layer?' % ((self.t2 - self.t1), polyCount), QMessageBox.Yes | QMessageBox.No, QMessageBox.Yes)

  NameError: global name 'polyCount' is not defined

  Python version:
  2.7.2 (default, Jun 12 2011, 15:08:59) [MSC v.1500 32 bit (Intel)]

  QGIS version:
  1.8.0-Lisboa Lisboa, 6416f38


texts: How to use ArcPy to copy selected features to an existing feature class?

texts:
Is there a way to copy a selected feature into an existing feature class?

Sort of like append would do, but only the selected features.

I want to iterate through a table and select features with a matching attribute in a feature class and copy that selected feature into an existing feature class.

Basically, it'd be like selected a feature in ArcMap and copying and pasting it into a different feature class.

Is this possible?


texts:
I'm new to using GIS.

What I'm trying to do is to create a table that shows how far are the schools I have currently highlighted (using the select by location feature) away from the Major arterial Roads and National/state highways in 3 criteria (close 0-300, medium 301-600m and >600m). I've been trying to create a new layer based on select by location- the same way I usually do for select by feature but it doesn't work.
Any advice?
I've included a screenshot for reference.


texts:
I've got a numpy array and need to sum the values in the rows. Then use the Update Cursor to update a field with this value. However the update cursor is updating the field with the same values.

tbl = arcpy.da.FeatureClassToNumPyArray(dg, ave_fields)

with arcpy.da.UpdateCursor(dg, "Weekly_Ave") as ucursor:
    for urow in ucursor:
        for i in range(len(tbl)):
            urow[0] = sum(tbl[i])
            ucursor.updateRow(urow)


texts: Difficulty with nested if statments in Survey123

texts:
I'm new to this and having great difficulty with the below nested if statments. What I want it to do is:

If the user selected Periodic and Dropoff time or Collection time
then produce the integer that the user chose from the
periodic_speed_limit_low list.

If the user selected Periodic and Off peak time then produce the
integer that the user chose from the periodic_speed_limit_high list.

If the user selected Fixed  then produce the integer that the user
chose from the fixed_speed_limit list.

if(selected(${speed_limit_type},'Fixed'),int(${fixed_speed_limit}),if selected(${speed_limit_type},'Periodic') and selected(${time_day}, 'Dropoff time') or selected(${time_day}, 'Collection time'),int(${periodic_speed_limit_low),if selected(${speed_limit_type},'Periodic') and selected(${time_day}, 'Offpeak time'),int(${periodic_speed_limit_high))

Can someone tell me why Connect isn't accepting the above?


texts: How to calculate the distance between the centroid and the nearest point of the polygon in QGIS?

texts: I have a shapefile as shown in the picture. I have estimated the centroid points of these polygons. How to estimate the distance between between the centroid and the minimum edge of the polygon in QGIS?


texts: Uploading unprojected shapefiles using GeoServer REST API

texts:
I have a number of shapefiles that lack the .prj file. I'm able to upload them to GeoServer using the REST API. However, I'm unable to serve them after upload. Which brings me to the question as to whether I need to reproject before uploading to GeoServer. If I need to reproject, then which options are available?

I'm using GeoServer 2.1.3.

UPDATE: Layers created from the shapefile lack projection information. When I update the Declared SRS field, the layers are served without issues.


texts: Getting list of QGIS field names, not including hidden fields using PyQGIS

texts:
Currently, for a python plugin for QGIS 2.8.3, I get the list of fields for a layer using pendingFields as mentioned in this thread:

fields = self.layer.pendingFields()

I now have a few hidden fields (set via the fields tab in layer properties). How can I get a list of fields that are not hidden?

I see the possibility to get Edit Type, but am not sure where to go from here :

EditType {
  LineEdit, UniqueValues, UniqueValuesEditable, ValueMap,
  Classification, EditRange, SliderRange, CheckBox,
  FileName, Enumeration, Immutable, Hidden,
  TextEdit, Calendar, DialRange, ValueRelation,
  UuidGenerator, Photo, WebView, Color,
  EditorWidgetV2
}


texts: Rectangle for Clip_management using arcpy

texts: I want to draw an interactive rectangle on screen which should be used as the extent for clipping geometry for raster datasets. Any suggestion?


texts: My identify window no longer appears when I identify a feature in ArcMap.


texts: Mobile offline mapping and editing of shapefiles?

texts:
I've got a client that wants to create a mobile app that can view and edit shapefiles stored locally on the device to enable it to work offline. There is a possibility of converting the shapefiles to KML or other format if needed.

It would need to work on Android, iOS and Blackberry so a HTML5/Phonegap solution is preferred.

Does anyone have any recommendations?


texts: Translating ArcGIS field calculator python expressions in to QGIS field calculator python expression?

texts:
I'm transitioning to QGIS from ArcGIS.  I'm having trouble getting my python field calculator expressions to work in the QGIS advanced python field calculator.

Here is a screenshot of the working script in ArcGIS:

Here is what the Advanced Python field calculator window looks like in Q.  I've made some changes based on posted Q&A's that I've seen here and in other places:

I've gotten some syntax errors that I think I've worked out, but then I'm also getting some "NoneType" errors, which I have no idea what they mean.  I'm assuming the "Global Expression [Optional]" block in QGIS is the same as the optional code block in ArcGIS?

What does a "NoneType" error mean?

Below I'll post a snippet of the code.

    def getvalue( Exemption_Code, Ag_Value, Appraised_Value, TaxRate_Municipal, Exemption_Local_Homestead, Exemption_Local_O65, Exemption_Local_Disability, Exemption_Vet_3, Exemption_Vet_4 ) :
      if (Exemption_Code == '000') :
        value (Appraised_Value / 100) * TaxRate_Municipal
      elif (Exemption_Code == '002') :
        value (Ag_Value / 100) * TaxRate_Municipal
      elif (Exemption_Code == '004') :
        value 0

    value = getvalue( <Exemption_Code>, <Ag_Value>, <Appraised_Value>, <TaxRate_Municipal>, <Exemption_Local_Homestead>, <Exemption_Local_O65>, <Exemption_Local_Disability>, <Exemption_Vet_3>, <Exemption_Vet_4> )


texts:
While opening the software the error message is poping up.

The procedure entry point ucrtbase.terminate could not be located in the dynamic link library.api-ms-win-crt-runtime-I1-1-0.dll

Installed version is QGIS 3.10.8


texts: Calculating distance between two points with QGIS 2.14.0

texts:
I have two shapefiles with the following information:

List A:
ID 1 | xA | yA | ID 2
List B:
ID 1 | xB | yB | ID 2

{x;y} = {lat;lon}

I want to calculate the distance between two points, e.g. ID 1: distance from {xA,yA} to {xB,yB}.

The KBS is WGS 84 and I need the result in meters. Can anyone tell me how to do that?


texts: mapserver vrt et mapfile : file not found or unsupported format

texts:
I try to create a mapfile with data in a csv file. I have failed to draw the layer with a Virtual Format and a OGR connectiontype directly in the mapfile. So I have decided to store the connection in a different file like here http://mapserver.org/input/vector/VirtualSpatialData.html#steps-for-display
mapfile :

LAYER
CONNECTION "C:\blabla\basol_auv.ovf"
CONNECTIONTYPE OGR
DATA "test"
METADATA
    "wms-srs" "2154"
    "wms-title" "test"
END
NAME "pollution"
TYPE POINT
CLASS
    NAME "pollution"
    STYLE
        SYMBOL "pollution"
        SIZE 10
    END #----fin style
  END #-----fin class
END #-----fin layer sols pollués

basol_auv.ovf :

<OGRVRTDataSource>
   <OGRVRTLayer  name='test'>
    <SrcDataSource relativeToVRT="0">C:\OSGeo4W\Apache\htdocs\data\test.csv</SrcDataSource>
    <SrcLayer>test</SrcLayer>
    <GeometryType>wkbPoint</GeometryType>
    <GeometryField encoding='PointFromColumns' x='A' y='B'/>
    </OGRVRTLayer>

I have a new error message :

  msDrawMap(): Image handling error. Failed to draw layer named 'pollution'.
              msOGRFileOpen(): OGR error. Open failed for OGR connection in layer `pollution'.  File not found or unsupported format.

  I dont'know what to do to verify if it is a matter of path or format.

ogrinfo command gives the correct answer : ogrinfo basol_auv.ovf test


texts:
I have a field that contains various category names, lets say A, B, C, etc... I want to be able to know which categories are in the current data frame view.

How should I go about doing this?

I'm using ArcMap 10.5


texts:
I have got a point layer with about 1000 points. These point has a rule based symbology. So far so common.
In this layer I have also some actions implemented. If I unselect some rules, that only the objects of one rule are displayed. In my case ther are just 10 left. When I want to activate an action I click on an object, but it offers me several objects even those which are not displayed.
Is this a bug or is there a setting which I can set just to make a action to shown objects.


texts: Cannot use raster layers not defined as inputs in a model in QGIS 3.x

texts:
I have a problem in building raster models: both with QGIS 3.4.5 and QGIS 3.10.3, the raster calculator algorithm, used in a model, refuses to work if the requested layers are other than the defined raster layers as inputs.

The model is blocked with the following message:

An error occurred while performing the calculation
Error while running Raster calculator.

I'm looking for clear cuts in the forest with 2 satellite images n and n-1.

The model bug on the last algorithm (Raster calculator).


texts: gdal_rasterize is giving no data while converting a shapefile to raster

texts:
I have used next command:

gdal_rasterize -burn 1 -tr 30 30 fire_2010.shp fireRaster2010.tif

to convert shapefile to raster. But I am getting no data in the raster. The same command works for another shapefile. The only difference I see is the polygons in this shapefile are small.

So, from the image it can be checked that fireRaster2010 has 0 0 [no values]. The fire2010 is a shapefile, and that is visible as blue colored patches.

here is the shapefile

here is the raster file- no data. The values 0 0 in the Layer Panel signifies that there is no data in the file. All the values are zero.


texts: Make ST_Polygonize in Postgis return dangles, cuts and invalid rings

texts:
Postgis allows using ST_Polygonize() from GEOS to make polygons out of lines. I have seen the same function available in Shapely library for Python, but in Shapely, the other method named polygonize_full() is also available.

In PostgreSQL 10.1 + PostGIS 2.3.0, I would do something like this:

SELECT (ST_Dump(foofoo.polycoll)).geom AS geomtextrep
FROM (SELECT ST_Polygonize(wkb_geometry) AS polycoll
FROM (SELECT wkb_geometry FROM the_table WHERE ST_GeometryType(wkb_geometry) = 'ST_LineString') AS foo) AS foofoo;

In Python3 + Shapely, I would do something like this:

result, dangles, cuts, invalids = polygonize_full(lines)

I have seen in the GEOS library that there are two functions:

GEOSPolygonize_full - Returns polygons, dangles, cuts, invalid rings
GEOSPolygonize - Returns only polygons

Is there any way to replicate the behaviour of GEOSPolygonize_full with PostGIS? I would like to use the dangles, cuts and invalids to judge whether I have missed out on some lines in the input data, and flag these events so I can manually check the source files (AutoCAD drawings) and associated layers.

I would hope ST_Polygonize_Full() maybe existed in PostGIS, but as far as I can see, it does not.


texts:
Does ArcGIS 10 use parallel processing?

In which version of ArcGIS can we get this facility, if it exists?

Is ArcSWAT single-threaded ?


texts: Querying for underscore character in ArcMap?

texts:
For a standard LIKE query against an Oracle-based ArcSDE feature class, the underscore character represents a single character wildcard when used with a string.

I'm trying to impose a definition query to find a text string that starts with 4 digits exactly followed by an underscore character.

Does anyone know how I would specify an underscore character itself in a query or if/what the escape character might be?

MDHald's answer works for file geodatabases but my case is specific to Oracle. Falsely assumed that ArcSDE and file geodatabase query would function the same for this case.


texts: Closest facility analysis not offering shortest route in ArcGIS Network Analyst?

texts:
I want to get the shortest paths from some points to the nearest highway junctions. I extracted the highway junctions by doing the Intersect function from the network junctions and selected highways. Then I created the network database and did the Closest Facility Analysis (highway junctions for facility and points for incidents) to get the shortest path from the points to their nearest highway junctions. The problem is some of the obtained "shortest paths" are not the actual shortest paths since they skip the nearest highway junctions. It seems the highway junctions extracted from the network are not perfectly aligned on the network. I tried to use Integrate to integrate the highway junctions and network, and then recreated the network database, but it still didn't solve the issue. The following pictures show the problem. The red line represents the highway. The purple dots represent the highway junctions. The purple square represents the point needed to be analyzed to get shortest path.

Any suggestions?

I am using the ArcGIS Desktop 10.3.1 Advanced License.


texts: How to export tables to Excel from ArcGIS 10.1 using the xlsxwriter python module?

texts:
(I've edited this question to reflect some of the answers/insights from below...)

Original(ish) Question:
I'm wondering if someone can share a bit of insight (or actual code) on how exactly to export tables from ArcGIS 10.1 to an .xlsx spreadsheet (Excel 2010), specifically using the xlsxwriter python module.

(If you have a suggestion other than xlsxwriter, that's welcome, too, as long as it creates files with the .xlsx extension).

Some answers simplified from below:
If you have ArcGIS 10.2+ there are built-in Table To Excel tools.

If you aren't interested in the code, but just want an ArcToolbox script that works at the push of a button (if you aren't already using 10.2+), use this awesome Excel and CSV Conversion tool that includes a Toolbox in the .zip file.

Otherwise, if you're interested in the code itself, see below for some code specific to xlsxwriter that works. But note, "XlsxWriter is designed only as a file writer. It cannot read or modify an existing Excel file." (per the xlsxwriter FAQ's)

Check out some other Excel conversion modules for more options. Openpyxl seems like a good option for reading & writing.


texts: How to calculate difference between bounds in Google Maps?

texts:
Is it possible, after "bounds_changed" event, to get polygon coordinates covering the difference between new bounds and previous bounds?

Any idea how to do it? Any help would be greatly appreciated.


texts: QGIS3.10, Find all points that are located on a line

texts:
Inside my .shp there are two simple layer:

"Points_Example" (Type: Point, EPSG:25832)
"Lines_Example" (Type: Line, EPSG:25832)

Some of the points are located exactly on the line and I'm trying to filter out these points via Python, but for some reason QGIS is always saying there are no points on my line!
layerLines = QgsProject.instance().mapLayersByName('Lines_Example')
layerLine = layerLines[0]

layerPoints = QgsProject.instance().mapLayersByName('Points_Example')
layerPoint = layerPoints[0]

list = []
for line in layerLine.getFeatures():
    geom = item.geometry()

    for point in layerPoint.getFeatures():
        if point.geometry().touches(geom):
            list.append(point['id'])

I've also tried some other methods from "QgsGeometry", but nothing seems to work! How can I solve this?


texts:
I am working on a project that focuses on Scotland. When I use the print composer to add a scalebar the scale is wrong. It shows 50 km for a distance that is probably around 200 km. Following the advice in this forum, I changed my project projection to a UTM one, in my case, UTM zone 30 WGS 84 which should be suitable for Scotland. The on the fly (OTF) box is ticked as well. But I still have the same issue with the scalebar. I would be very thankful for any suggestions.

I am using QGIS 2.8.3 Wien

Sorry I can not find the post I read early, but it suggested to switch the project projection to UTM to ensure that the scale works properly.


texts: Populate Null Records with Last Non-Null in List

texts:
I'm trying to develop a script that will iterate through a list of values within a FBGDB table/field and insert the last non-null value that was found into the ones that are null.

Here's what the table looks like:

This is what I have so far:

import arcpy

updatecursor = arcpy.da.UpdateCursor("D:\Tools\Scripts\insert over null\Testing.gdb\Test_tbl", ["From_Value", "To_Value"])

for row in updatecursor:
    if row[0] is None:
        print "There are Nulls"
        row [1] = row[0]([-1])
        print row[1]
    else:
        row[1] = row[0]
        print "Has Value"
    updatecursor.updateRow(row)
del row, updatecursor

I know my issue is in row[1] = row[0]([-1]) but I can't figure out how to make it read backwards to the last non-null value which in theory would be the previous value in the list (assuming it's populating as it iterates through.


texts: Identifying feature showing all attributes are '0' on raster using QGIS?

texts:
I have been using the identify tool on QGIS fine all week, i went back to check a value and now everytime i identify the feature the attribute i get is '0'.

Has anyone encountered this issue or know of any possible fix.

See below image - for some background i am trying to see the flood depth from a raster which i can see isn't 0 based on the band colour but i can't get the exact figure.


texts:
I'm using geotools to get all layers from a Geoserver.

    WMSCapabilities caps = null;
    try {
        WebMapServer wms = new WebMapServer( sourceServer, 1000000 );
        caps = wms.getCapabilities();
    } catch ( Exception e ) {
        e.printStackTrace();
    }

All is fine, the code is working well when the capabilities XML is of an average size. The problem is when the XML is huge. I'm trying to play with the time parameter of new WebMapServer()but still getting the timeout error even with a big timeout value.

Failed to execute request http://10.5.115.22/geoserver/wms?REQUEST=GetCapabilities&VERSION=1.1.1&SERVICE=WMS
java.net.SocketTimeoutException: Read timed out
    at java.net.SocketInputStream.socketRead0(Native Method)
    at java.net.SocketInputStream.socketRead(SocketInputStream.java:116)
    at java.net.SocketInputStream.read(SocketInputStream.java:170)
    at java.net.SocketInputStream.read(SocketInputStream.java:141)
    at java.io.BufferedInputStream.fill(BufferedInputStream.java:246)
    at java.io.BufferedInputStream.read1(BufferedInputStream.java:286)
    at java.io.BufferedInputStream.read(BufferedInputStream.java:345)
    at sun.net.www.http.HttpClient.parseHTTPHeader(HttpClient.java:704)
    at sun.net.www.http.HttpClient.parseHTTP(HttpClient.java:647)

Is there any solution for this?


texts: Using multiple satellite constellations (e.g. GLONASS, GALILEO) by GPS unit?

texts:
I am trying to do research on GNSS/GPS units available and it seems like there is a large price jump for units with more satellite constellations (maybe I am erroneously making this correlation) that the unit is connecting to.

Specifically, I was trying to compare the arrow lite and arrow 100 (and other brands) and it seems like the only major difference but there is a price difference of 1k between the two.

What is the advantage of having multiple constellations?


texts:
First time using (QGIS on UBUNTU).
i am following some articles in order to publish qgis maps on the web.
My problem is (although i am the administrator) i am not able to save or copy my *.qgs in specific directories like /qgis-web-client/project.
I don' t even know where are some directories, since i am really confused with the ubuntu file structure philosophy.
Is there any way to change permissions? Is that a good idea?