texts:
Esri's telecom data model has 2 geometric networks, one for infrastructure (manholes and conduits) and one for fiber cables.

fiber cables are easily related to a section of conduit based on existing relationship class.

How can you label the conduit in a way that shows all the fiber cable_ID's that are related to that conduit?


texts: Proj.NET transform wkt to lat lon coordinate list

texts:
I am looking to transform from a WKT in 27700 into lat lon list of coordinates, but cannot quite work out how to do it using this library:

https://archive.codeplex.com/?p=projnet

I can easily do a single easting/northing point like this:

    public static double[] EastNorthToLonLat(double easting, double northing, CoordinateSystem fromCS)
    {
        var convertFrom = CS27700;
        switch (fromCS)
        {
            case CoordinateSystem.Epsg27700:
                convertFrom = CS27700;
                break;
            case CoordinateSystem.Epsg4326:
                convertFrom = CS4326;
                break;
            case CoordinateSystem.Epsg3857:
                convertFrom = CS3857;
                break;
        }

        var trans = new CoordinateTransformationFactory().CreateFromCoordinateSystems(convertFrom, CS4326);
        var result = trans.MathTransform.Transform(new[] {easting, northing});
        return result;
    }

I have seen this so I preume this has something to do with it, but I can't work out how to call the transform with the WKT, all it takes is the array of coordinates:

var fromCS = CoordinateSystemWktReader.Parse(wkt) as ICoordinateSystem;

My ulltimate is to get it into a format like this:

[[-0.165,51.185] , [-0.16,51.180], [-0.265,51.584]] etc


texts: How to unzip multiple folders with FME 2014

texts:
I have downloaded 21 zipped folders using the HttpFetcher transformer and I cannot figure out how to get FME to unzip these files.

I've seen some posts about using either the TclCaller or PythonCaller transformer, but I've tried that without luck.

Does anyone out there know how to unzip directories, using FME?


texts: Segmentation and Classification in ERDAS

texts:
I have segmented an image inn erdas and will want to use it in a classification exercise, how do i proceed?

My image is a low resolution 30m alos.


texts: Text in Google maps displayed as blocks in Chrome

texts:
Any idea why this is happening in Google maps while using Chrome? I can see the text display but then the text changes to blocks when it finishes rendering.

I'm using Chrome Version 54.0.2840.90 (64-bit) in Ubuntu 14.04


texts: Is there any way to perform spatial filter on multiple layers together?

texts:
I have a scenario where a featureclass[Child Version] is queried for some overlapping features on [SDE.DEFAULT] Version using ISpatialFilter, After getting few result, these features are looped and again queried the same way on [Child version]. (exactly opposite).
This takes lot of time. I still hope that the compexity can be reduced somehow. Is there any ways that we can run spatial filter in some kind of combination, so as to increase performance.
I have tried Spatial Cache, but because extent is not known it takes lot of time to fill the cache.

Any suggestions?


texts: How do I set mmqgis create grid options?

texts:
I am having a puzzle about using mmqgis create vector layer option. I want to use UTM CSR

this are my options. I want a grid of 20,000 meters in the east and north axis. Nevertheless i have an error screen stating: "invalid width/height 6.0/6.0". As I understand my settings I should have 6 lines separated by 20,000 meters
Am I understanding this in the wrong way? Thanks in advance


texts: QGIS 2.18.7: How can I change or add WMS parameters?

texts:
I'm trying to add open data from MétéoFrance (French meteorological organisation) as WMS with QGIS 2.18.7 on Windows 7.

The URL looks like this, I've added it to my WMS connections (with my valid token):

https://geoservices.meteofrance.fr/api/TOKEN/ClimNormalesWMORes40EssMetroWMS?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&BBOX=41.21598999999999791,-5.14082300000000014,51.48935699999999827,8.38489399999999918&CRS=EPSG:4326&WIDTH=515&HEIGHT=392&LAYERS=JANUARY_NORMALS_OF_MEAN_PRESSURE_REDUCED_TO_MSL&STYLES=DejaVuSans_circle_orange&FORMAT=image%2Fpng&DPI=96&MAP_RESOLUTION=96&FORMAT_OPTIONS=dpi:96&TRANSPARENT=TRUE

But if I add the layer to my empty project, nothing loads. The layer is in my list, but it doesn't show anything.
The support told me, that there is always styles=DejaVuSans_circle_orange or format image%2Fpng (%2F replaces /) missing in my request.

The layer source is indicated as

contextualWMSLegend=0&crs=EPSG:4326&dpiMode=7&featureCount=10&format=image/png&layers=JANUARY_NORMALS_OF_MEAN_PRESSURE_REDUCED_TO_MSL&styles=DejaVuSans_circle_orange&url=https://geoservices.meteofrance.fr/api/TOKEN/ClimNormalesWMORes40EssMetroWMS?VERSION%3D1.3.0%26BBOX%3D41.21598999999999791,-5.14082300000000014,51.48935699999999827,8.38489399999999918%26CRS%3DEPSG:4326%26WIDTH%3D515%26HEIGHT%3D392%26LAYERS%3DJANUARY_NORMALS_OF_MEAN_PRESSURE_REDUCED_TO_MSL%26STYLES%3DDejaVuSans_circle_orange%26DPI%3D96%26MAP_RESOLUTION%3D96%26FORMAT_OPTIONS%3Ddpi:96%26TRANSPARENT%3DTRUE

Another possible URL structure is:

https://geoservices.meteofrance.fr/services/ClimNormalesWMORes40EssMetroWMS?version=1.3.0&token=TOKEN&BBOX=41.21598999999999791,-5.14082300000000014,51.48935699999999827,8.38489399999999918&CRS=EPSG:4326&WIDTH=515&HEIGHT=392&LAYERS=JANUARY_NORMALS_OF_MEAN_TEMPERATURE&STYLES=DejaVuSans_circle_noir&DPI=96&MAP_RESOLUTION=96&FORMAT_OPTIONS=dpi:96&TRANSPARENT=TRUE

But the structure above is the one that their service gave me.

In my browser a png appears, it seems to be the right data.

Even if I choose png and the DejaVuSans_circle_orange layer from the connection list, it doesn't work.

How can I add this parameters to my request? Or is there another way to load the WMS in QGIS?


texts: Displaying a Digitial Elevation Model using GRASS GIS draped by an RGB Raster

texts:
I have the following:

Digital Elevation Model from SRTM
LANDSAT Imagery from the same area

I wish to display the DEM draped with the raster band combination (3,2,1) using GRASS GIS.

My attempts to use:

a. Create r.composite and then use it as an input to the color model for NVIS results in a B/W rendering of the raster surface.

b. Creating a VTK ASCII file using r.out.vtk is confusing as there is a confusion what is an input image and what is an rgbmap (created using the RED, GREEN and BLUE) bands of LANDSAT. There is of course no issue with the elevation data. However using paraview to display the output file data.vtk fails (does not work).


texts:
Ok, likely my bad for not making it to WhereCamp5280 and asking the FortiusOne geeks directly, but what's the process for creating GeoJSON 'tiles' for use in Polymaps?

The Population example over at polymaps.org notes that the data for the demo runs on Google AppEngine, and alludes to "rolling your own" vector tiles, but I can't seem to find any more information...

Anyone got some insight they can share?

Thanks!


texts: Should point data be equidistant projected when using ArcGIS IDW spatial interpolation?

texts:
I am wondering if anyone can clarify whether point data should be projected to an equidistant projection when using ArcGIS IDW spatial interpolation?

I am working on a dataset from Western North America spanning about 30 degrees of latitude. The data are currently in Lat/Long (NAD83).

Does ArcGIS "project on the fly" or somehow adjust for latitude when calculating distances to my sample points during the interpolation procedure or should I be supplying everything in a projection that preserves distances?


texts:
I am digitizing polygons from an ortophoto. In the snapping options, I enabled 'topological editing' and 'avoid intersection' of the polygon layer. I digitized some not-adjacent areas and then made a large polygon to fill the gaps between the layers. It worked perfectly and I do not see any overlaps, but when I use the Topology Checker it finds overlaps. I know how to clean them with GRASS, but I would rather prefere digitizing without overlapping....How does it works, or is there a bug?

I am using QGIS 2.14.1.


texts:
if (workspace=="x") {
        bbox = [3461431.90483,5555489.54044,3461896.81256,5555951.68634];
    } else
    if (workspace=="y") {
        bbox = [471730.87469,5524059.31339 ,472552.36397,5524962.18860];
    }
...
    map = new OpenLayers.Map("map",{
      projection : new OpenLayers.Projection("EPSG:31466"),
      units : "m",
      maxResolution : 2,
      maxExtent : new OpenLayers.Bounds(bbox)
  });

I work on a gis client that makes use of OpenLayers / Geoserver. Since I use the same code for several maps, it would be very useful if I did not have to hard-code the bounding box of the respective map in javascript, but to retrieve the coordinates from the geoserver. Is this possible?

-_________________________________

I know that I can receive the bounding box via GetCapabilities as a huge XML File that I would have to parse, but there must be a shortcut.


texts:
I encountered a problem in QGIS 2.10.1 when I tried to calculate the areas of polygons in a polygon layer.

I followed the instructions here (http://manual.linfiniti.com/en/vector_analysis/reproject_transform.html) to save the layer as one with selected CRS UTM zone 33S. It saved successfully; however, for some reason the saved layer with the new CRS won't show up on the map. Moreover, when I went to field calculator to try to calculate the area, it simply shows "nan" as results.

I have calculated areas in the same way before, and never encountered such problem. It seems that saving the original layer with UTM zone 33S causes the problem (since the saved layer won't even show up on the map). Any thought why this happened and how should I fix it?

Below is the link to the shapefile so that you can play with it to see the problem.
shapefile


texts: How to set checked "Auto Open Form" in identify results window, programatically?

texts:
I want to set check programmatically "Auto Open Form" check box in Identify Results window. Window is here and the check box is at the bottom-right:

The Identify Results window is not open at that time!


texts: Add base map to animation Google Earth Engine

texts:
I want to make an animation of fires using MCD64A1 with a base map something similar to this Sentinel 5-P animation:

What I could make until now using this code is:

// Define a collection
var col = ee.ImageCollection('MODIS/006/MCD64A1').filterDate('2019-01-01', '2020-12-31').select('BurnDate')

// Define a mask to clip the data by.
var mask = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017')
  .filter(ee.Filter.eq('wld_rgn', 'South America'));

print(mask);

// Add outline to the Map as a layer.
Map.centerObject(mask, 3);
Map.addLayer(mask);

// Define the regional bounds of animation frames.
var region =
    ee.Geometry.Polygon(
        [[[-86.88829713135783, 12.937849176033298],
          [-86.88829713135783, -30.715685432081326],
          [-30.81407838135783, -30.715685432081326],
          [-30.81407838135783, 12.937849176033298]]], null, false);

print('Region', region)

// Define RGB visualization parameters.
var visParams = {
  min: 0.0,
  max: 366,
  palette: ['FF0000']
};

// Create RGB visualization images for use as animation frames.
var rgbVis = col.map(function(img) {
  return img.visualize(visParams).clip(mask);
});

// Define GIF visualization parameters.
var gifParams = {
  'region': region,
  'dimensions': 600,
  'crs': 'EPSG:3857',
  'framesPerSecond': 5
};

// Print the GIF URL to the console.
print(rgbVis.getVideoThumbURL(gifParams));

which gives the following result:

I didn't find a way to do it.


texts: Convert a decimal degree to meters from the MySQL Spatial GLength

texts:
The MySQL Spatial function GLength returns the distance between two points, this distance is  in the same unit from the projection. Since I'm using the WGS84 (SRID:4326) the answer will be with decimal degrees. But I need to convert this to meters.

Searching I got the answer to my question that there is a proportion of 111200 meters (others say 111111 or 111044) to 1 degree.

Checking with QGis ruler with the geometrys generate I got a different number to that proportion: 80321 meters/1 dd

Which one is right?


texts: Extract all multiband pixel values inside the polygons from stars-proxy

texts:
I use stars and sf packages of R programming language for satelite image manipulation and data mining (with tidyverse packages).
I try to extract all pixel values from stars proxy object inside the many polygons. I use GeoPackage polygons as sf object and very large multiband GeoTiff image as stars proxy object. How can I decrease using RAM for extracting pixel values and joining it to polygons attributes? Finally I need to construct table with joined attributes from polygons and multiband image (as an picture).


texts: How do you create a CartoDB viz with a dynamic SQL statement?

texts:
I want to create a CartoDB viz where the where parameter in my SQL statement changes dynamically at runtime.

When I share my CartoDB viz from the dashboard, I can't help but notice the the SQL statement is passed as a URL var.

I thought that was an easy way to dynamically change the where clause.  After playing around  though, I found out that the parameter is ignored (and for good reason, to guard against SQL injection), though other parms like showing the title/desc and the ne/sw latlongs are processed.

So I'm now trying to create a viz using the CartoDB API.

But for some reason, I can't get the creating visualizations at runtime example to work.


texts: QGIS equivalent for ArcMap join for purpose of distances from points to polygon?

texts:
It is relatively easy to calculate distances from points to a polygon in ArcMap by doing a join and selecting the option for creating a field containing the minimum distance information. i.e. it generates a field containing the distance from each point to the boundary of the polygon.

Is there a relatively simple way to do the same thing in QGIS? I've tried NN Join, but it crashes each time I try to use it. I've also tried the Processing toolbox's 'distance to hub', but it gives the distance to the centroid of the polygon instead of the boundary.


texts:
I'm looking for some info about the possible limitations for the GLONASS navigation system, something similar to the GPS COCOM limits.
Is there any limitation on speed AND/OR altitude for this navigation system?

We are looking into embedding a GLONASS receiver into an high-altitude balloon project.


texts: ogr2ogr latitute longitude conversion formats

texts:
trying to convert shape files into mysql tables. Able to create table and load data successfully using OGR2OGR. Only challenge I see is latitude/longitude is not exactly the same format as one of my developer (no longer there) who did the conversion for me. Other than geometry shape, rest of the values are exactly matching.

For the same record my conversion shows
LINESTRING (3015615.1663366337 13849435.592079207,3015602.495049505 13849789.215841584)

My developer conversion shows
'LINESTRING(-95.69831901099997 29.789630642000077,-95.69832747399994 29.79060369800004)',3

what am I missing ? what should I do convert my shapefile into the format my devloper got it. Tried to convert into .csv as well but getting the same result.


texts: How to convert CRS projection of a map in QGIS?

texts:
I am totally new to GIS. I have a QGIS project. I am adding the screenshot.

As you can see the map is sort of rectangle in projection.  How can I project it like this (without the grids)

Thanks


texts: Percent coverage of a polygon in a raster grid

texts:
I have a wooded area polygon and a study area grid of cells with 300m resolution.

I would like to calculate a raster which shows the amount of each grid cell covered by the polygon. The density per km2 would also work, rather than just a binary 1 (has woods), 0 (does not).


texts: PostgreSQL COPY command shifting CSV values into wrong column

texts:
I am running the following query in PostgreSQL 10.1:

COPY addressbase FROM 'C:/Temp/ab_postgres/addressbase.csv' DELIMITER ',' CSV HEADER

It is returning the error:

ERROR: invalid input syntax for type numeric: "SY22 6XP"
CONTEXT: COPY addressbase, line 2, column y_coordinate: "SY22 6XP"
SQL state: 22P02

I fully understand that it's stating I can't put text in a numeric field. But my issue is that the value SY22 6XP is not actually in that field at all. So how come when I run this query, it moves the values in the CSV to different columns? See a snapshot of the data below - the value referenced here is in K2, not M2 which is where PostgreSQL is seeing it.


texts: Including Geometry in ModelBuilder Join?

texts:
I want to perform a join between a feature class and a table and then export the data as a feature class using ArcGIS ModelBuilder.

My current problem is that the geometry is not present in the output.

I have tried 2 methods so far:

1) Join Field: (which is supposed to work for Tables or Feature Classes)

2) Add Join: Creating a Table View and a Feature Layer from the associated data sources, and performing the Add Join

In each option I have tried the Feature Class (source) both as Input and Join.
The output never contains the Geometry - it is always a Table.

How can I include the Geometry in the output?


texts: How to give GPKG layer as output in processing algorithm?

texts:
I try to use qgis:creategrid processing tool. I want to give as output source gpkg layer which I have already create by defining it shape type, crs. However it does not work for me now. Here is my code:

parameters = {'TYPE': 2,
              'EXTENT': '%f,%f,%f,%f' % (xmin, xmax, ymin, ymax),
              'HSPACING': dx,
              'VSPACING': dy,
              'HOVERLAY': 0.0,
              'VOVERLAY': 0.0,
              'CRS': str(epsg),
              'OUTPUT': self.gpkgFile + "|layername=Grid_1"}

processing.run('qgis:creategrid', parameters)

self.gpkgFile is defined the gpkg source such as: D:\newProject\project1.gpkg

So 'OUTPUT' describes here: D:\newProject\project1.gpkg|layername=Grid_1

How can I use gpkg layer as 'OUTPUT'in such kind of processing tools?


texts: How does the radius work in Trilateration

texts:
I've been trying to understand the trilateration algorithm based on the wikipedia article. As of now, I have a system where I have the distance of 3 sensors wrt my position and i have the exact GPS coordinates of those 3 points.

How do I include the radius?


texts: Customize LabelLayer extent based on density/location using ArcGIS JavaScript API

texts:
I'm using a LabelLayer that can be seen only within a certain minScale - maxScale window. I'm trying to display labels of census tract IDs within Illinois, and the difference in label density between the Chicago area and the rest of the state is significant. Would there be a way to customize the minScale and maxScale settings based on the feature? I'm very new to JavaScript, so I'll just demonstrate what I'm looking for in pseudocode:

if ("County_nam" == "Cook") {
          minScale: 50000,
          maxScale: 0
    }
else if ("County_nam" == "Lake" or "County_nam" == "Dupage" or "County_nam" == "Will") {
          minScale: 100000,
          maxScale: 0
    }
else {
          minScale: 500000,
          maxScale: 0
    }


texts:
By clicking on a point on a map at a detected site, numerical values appear from a qc_id field automatically generated in the layer, in the attached SS showing the value 16.
I would like the value of another field to appear instead: could you show me how to operate in QGIS to bring up the values of the other airports_q field?


texts: QGIS count number of arguments in range of fields

texts:
I want to update an existing field with the number of values, in a range of preceeding fields.
I could do this easily in EXCEL with "=NUMBER ARGUMENTS", but I want to learn this in QGIS.

I want to update the field named "totaal" with the number of 'x' in the preceeding fields. Count the arguments in fields "5 jul" to "12 sep".


texts: Error while Importing data from .SHP files to spatialite using OGR2OGR

texts:
I was trying to import data from files provided here http://download.geofabrik.de/europe/france.html
 specifically http://download.geofabrik.de/europe/france/languedoc-roussillon-latest-free.shp.zip

I'm using OGR2OGR, I ran the command

ogr2ogr append -f SQLite db.sqlite gis_osm_landuse_a_free_1.shp -dsco SPATIALITE=YES

The error I have are

  ERROR 1: sqlite3_step() failed:   gis_osm_landuse_a_free_1.GEOMETRY
  violates Geometry constraint [geom-type or SRID not allowed] (19)
  ERROR 1: ROLLBACK transaction failed: cannot rollback - no transaction
  is active ERROR 1: Unable to write feature 2026 from layer
  gis_osm_landuse_a_free_1.

  ERROR 1: Terminating translation prematurely after failed translation
  of layer gis_osm_landuse_a_free_1 (use -skipfailures to skip errors)`

The corresponding spatialite tables are empty so it seems importing failed.

I have searched for solutions and it seems the Geometry column should be made lowercase.

I'm not sure about any solution, since I'm using SHP files. Is there a possible fix?

Note that many other SHP provided in the ZIP file linked above are converted without any issue. This issue happens for several SHP files.


texts: Exception raised when saving Geopandas DataFrame into a file

texts:
I am trying to apply the method to_file of a Geopandas DataFrame.

My DataFrame contains a column Polygon, a column Point, a string and float columns.

Here the output of the describe method:

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 1332 entries, 0 to 1331
Data columns (total 7 columns):
ID                    1332 non-null object
geometry              1332 non-null object
centroid              1332 non-null object
epsg4326_longitude    1332 non-null float64
epsg4326_latitude     1332 non-null float64
epsg3857_x            1332 non-null float64
epsg3857_y            1332 non-null float64
dtypes: float64(4), object(3)
memory usage: 72.9+ KB

The specific types are:

Polygon: shapely.geometry.polygon.Polygon
Point: shapely.geometry.point.Point
String: str
Floats: numpy.float64

When I try apply the method to_file i have the following exception:

geo_grid_df.to_file(filename='grid_50km.geojson', driver='GeoJSON')
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-78-d1eb792e5885> in <module>()
----> 1 geo_grid_df.to_file(filename='grid_france_50km.geojson', driver='GeoJSON')

C:\ProgramData\Anaconda3\lib\site-packages\geopandas\geodataframe.py in to_file(self, filename, driver, schema, **kwargs)
    411         """
    412         from geopandas.io.file import to_file
--> 413         to_file(self, filename, driver, schema, **kwargs)
    414
    415     def to_crs(self, crs=None, epsg=None, inplace=False):

C:\ProgramData\Anaconda3\lib\site-packages\geopandas\io\file.py in to_file(df, filename, driver, schema, **kwargs)
    109         with fiona.open(filename, 'w', driver=driver, crs=df.crs,
    110                         schema=schema, **kwargs) as colxn:
--> 111             colxn.writerecords(df.iterfeatures())
    112
    113

C:\ProgramData\Anaconda3\lib\site-packages\fiona\collection.py in writerecords(self, records)
    333         if self.mode not in ('a', 'w'):
    334             raise IOError("collection not open for writing")
--> 335         self.session.writerecs(records, self)
    336         self._len = self.session.get_length()
    337         self._bounds = self.session.get_extent()

fiona\ogrext.pyx in fiona.ogrext.WritingSession.writerecs()

fiona\ogrext.pyx in fiona.ogrext.OGRFeatureBuilder.build()

ValueError: Invalid field type <class 'shapely.geometry.point.Point'>

Any ideas?

pandas: 0.20.3 | geopandas: 0.4.0 | fiona: 1.7.10


texts: How to set workspace from input Feature Class?

texts: How to set workspace from input Feature Class ? is there an attribute of the source GDB for a feature class ? some times the FC is in a dataset and some time it isn't, How can i get the path of source GDB for setting workspace when the input is dynamic ?


texts: Loading Raster with SharpMap (GDAL) doesn't work

texts:
I have problems loading rasterfiles (.tiff) with SharpMap

I tried initializing a SharpMap GDAL Layer but nothing happens. It doesn't even throw and error when debugging after I set a path for the file which doesn't even exist.

SharpMap.Layers.GdalRasterLayer G = new
SharpMap.Layers.GdalRasterLayer("Raster", "C:\\*.tiff");

I have all the gdal.dll linked in my Visual Studio Project.

Since I also have a world file for my raster (.tfw) I also tried using a GdiImageLayer, which can optionally use a world file https://www.nuget.org/packages/SharpMap/

 SharpMap.Layers.GdiImageLayer H = new SharpMap.Layers.GdiImageLayer("Raster", "C:\\*.tiff");

but how can I set it?

Sorry for my bad English and being completely new to C# programming.


texts: What is the fastest version of gdal2tiles?

texts:
I have multi-CPU linux server, so I wonder what is the fastest parallel version of gdal2tiles.py ?

My Ubuntu version:

lsb_release -a
Description:    Ubuntu 12.04.5 LTS

How I install the package:

sudo apt-get install python-gdal

Package version:

apt-show-versions python-gdal
python-gdal/precise uptodate 1.10.0-1~precise2

And under /usr/bin/ I have gdal2tiles.py.

Does it make sense is to update python-gdal package > 1.10.0(maybe build newer version from source), where I can see changelog about gdal2tiles.py?

After some googling I found that there are some patched versions of original gdal2tiles.py and some parallel versions of gdal2tiles.py. So I wonder why this functionality is not included in official package?

If anyone point to stable and fast version of gdal2tiles.py I can run some benchmarks for community.

Here is my script that I use for benchmark using original single threaded gdal2tiles.py.

set -x # Dump each line of this script
set -e # Stop on first error

DATASET="my_dataset"
LEVELS_OF_DETAIL=19

time gdalbuildvrt -srcnodata "0,0,0" orthomaps/$DATASET/global.vrt orthomaps/$DATASET/orthophoto_*.tif

#Default gdal2tiles.py
GDAL_TO_TILES="/usr/bin/python /usr/bin/gdal2tiles.py"
time $GDAL_TO_TILES -z 0-$LEVELS_OF_DETAIL orthomaps/$DATASET/global.vrt default_gdal2tiles_$DATASET

What versions I have tried:

Original version of gdal2tiles.py (single threaded).
Patch from here to make original gdal2tiles.py parallel, produce missing tiles.
Version from here (also parallel), produce missing tiles.
Version from here (also parallel), I can't even run it.
Maptiler, works well, but cost money.
tiler-tools cause an error.

If creating a TMS tiles with parallel gdal2tilesis an issue, I can move to any other tile format which is supported in Cesium( by cesium tile provider) and have parallel tiler, any ideas?

Related questions:

Performance of google map tile creation processes

How to cut a 'tiled image' from a raster dataset using GDAL?


texts: How to show cell boundaries on a raster?

texts: I am trying to recreate a situation where each cell from a raster dataset has a black border as seen in this image. Any ideas?


texts:
I read this topic and I have same problem as others. When I use network analyst on one line, the module works successfully. But for multilines, this warning comes:

No route from location "Graphic Pick 1" to location "Graphic Pick 2".

What can I do? Could you explain more, please ?


texts: Re-ordering layers before exporting to web map?

texts: I have two layers and two groups that need to be listed above a google hybrid layer.  When I start qgis2web it reorders the layers placing the two groups (roads and places) under the google hybrid layer...thus obscuring the two groups.  How do I reorder the layers to place the google hybrid layer on the bottom before exporting?


texts: Avoiding No data driven pages were created error in ArcMap

texts:
I am trying to enable Data Driven Pages in my mapdocument but I keep getting an error that says:

  No data driven pages were created. Please verify that the index layer
  has a valid source and contains features for a page creation.

I know that that my layer "polytrial selection 5" exists.

What is happening?

I added some screenshots for more context.

.


texts: Smartest way to union different "layers" in PostGIS

texts:
I'm currently using this crude approach to combine several tables/layers:

CREATE TEMPORARY VIEW p_impd_800 AS (
SELECT st_union(geom) as geom FROM rt_201
UNION
SELECT st_union(geom) as geom FROM rt_202
UNION
...
SELECT st_union(geom) as geom FROM rt_207
    );

Is there any smarter and faster way you to achieve the same result?


texts:
I have a WFS Service I publish from ArcGIS Server.  I have used several clients to test and the service either does not display the geometry or it just throws an error.   The layer shows up in the TOC, I can access, and in some cases I can get to the attributes.  The ArcGIS MAP Service will load fine and I get geometry and attributes correctly.   Once I do that, I can find the geometry of the WFS, but it is off in a totally different area and rotated 90 degrees.   The MAP Service and WFS are the same, published with WFS capabilities set when I publish the .mxd.  MAP Service OK, WFS an issue...

When I load the service into ArcGIS Portal it seems to display fine.    Unfortunately, that is just a test and clients need to get at the service through WFS without ArcGIS Portal.

When I view in Bentley/GAIA/GRASS/QGIS, and others, I get a different error.   Again, the layers show up in TOC in all of the products I have tested but throw various errors or totally skewed geometry.

Here is an error I get in QGIS when I test it for example.

Layer xxxxxxx.GRID_SPACING_MAP: Download of features for layer
Bentley_SpacingMap:xxxxxx.GRID_SPACING_MAP failed or partially failed: Download of features failed:
Error downloading:
https://xxxxxxxxx/arcgis/services/xxxxxxx/SpacingMap/MapServer/WFSServer?SERVICE=WFS&REQUEST=GetFeature&VERSION=2.0.0&TYPENAMES=xxxxxxx_SpacingMap:xxxxxxxx.GRID_SPACING_MAP&SRSNAME=urn:ogc:def:crs:EPSG::4431&BBOX=1080020.71277404553256929,12232953.25763050653040409,1354911.23849608213640749,12477890.19503248669207096,urn:ogc:def:crs:EPSG::4431
 - server replied: Bad Request. You may attempt reloading the layer with F5'

The services pull data from ArcSDE, however I also tried from a File Geodatabase as a test and got the same results.

Ultimately I am trying to consume in Bentley.  I have been using other clients to try and debug, including QGIS, with more or less the same results.   ArcGIS Portal is the exception.   QGIS is version 2.18.

No complex features.   Just a simple polygon layer with around 300 features.   No transactions or editing.

Per ArcGIS Server documentation:
The WFS services you create are compliant with the WFS 1.1.0 and 2.0 specifications. They also support the WFS 1.0.0 specification on a read-only basis.


texts: If field is empty in attribute table, programatically populate subsequent field with that data

texts:
I have a CSV communicating directly with a series of points and lines. The CSV is a bit problematic, however. The hardware that provides the data only provides location data every few rows. That is, there is an attribute table with X and Y columns whose fields are only filled only at key times. The rest is other still useful data, but location is only updated intermittently.

I have my point feature structured so that only the bottom-most (most recent) attribute row is used to display one point. Naturally I need location data for that bottom-most row, otherwise the point disappears. Currently I'm not getting it but sporadically. My CSV -> point code looks like this currently (pyshp library used):

        cur_dir = os.path.dirname(os.path.realpath(__file__))

        file_location = os.path.join(cur_dir, "data", "mycsv.csv")
        out_file = os.path.join(cur_dir, "data", "mypoint.shp")
        idd, az, y, x = [], [], [], []

        with open(file_location, 'rb') as csvfile:
            r = list(csv.reader(csvfile, delimiter=','))
            reader = r[len(r) - 1:] #read only top row
            for i, row in reversed(list(enumerate(reader))): #reverse table
                idd.append(str(row[0]))
                az.append(float(row[2]))
                y.append(float(row[6]))
                x.append(float(row[7]))
                print "appended"

        w = shapefile.Writer(shapefile.POINT)
        w.field('ID', 'N')
        w.field('AZIMUTH', 'N', 12)
        w.field('Y', 'F', 10, 8)
        w.field('X', 'F', 10, 8)

        for j, k in enumerate(x): #write shapefile
            w.point(k, y[j])
            w.record(idd[j],az[j],y[j], k )
            print "recorded"

        w.save(out_file)

Here is a snippet of the CSV. You will see the chopped up columns on the right:

I have started by trying to add if x > 0: conditionals in various places below the multiple enumerators to possibly halt digitization until data is read, but I am given error min() arg is an empty sequence each time. I'm not even sure this is the right place to begin, as it would probably just put attributes out of sequence with each other. I really just need some sort of "if X or Y field is null, populate field with prior field's data" script.

I am using Python 2.7, developing out of QGIS 2.8 API, and running on Ubuntu 14.04,


texts:
I have a few hundreds of .tiff files that I converted as a raster stack. However, when I tried reading the file, it is classified as a RasterLayer. I was expecting it to be a RasterStack. The functions that I would use to my usual raster stack won't work i.e. subset(). Also when I convert it as.data.frame, it only give me the xy coordinates and the value of the raster. Did I make a wrong decision of stacking my rasters? Any hints on how to break down my rasters layers?

Btw, my raster looks like this:
class      : RasterLayer
band       : 1  (of  414  bands)
dimensions : 605, 973, 588665  (nrow, ncol, ncell)
resolution : 231.6564, 231.6564  (x, y)
extent     : 11267070, 11492472, 1355190, 1495342  (xmin, xmax, ymin, ymax)
crs        : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs
source     : /Users/manuel/UWThesis/RCodes/files/KTpostdam_evi.tif
names      : KTpostdam_evi
values     : -19030000, 79780000  (min, max)


texts: Cutline mosaicing areal imagery with open source software

texts:
I'm very new in gis, and i'm trying to find a fast solution to merge 200 historical areal images(black and white from 1944). I orthorectified them manually, using gcp-s. They have 1 to 10 meters error. So I cannot use simple mosaicing software. I want to use roads as boundaries. Is there any os software what can help me in this? Or shell i cut them by myself, and merge them together with qgis for example?

Bonus question: Is there any automated contrast and brightness balancer for this amount of images?


texts: How to determine the appropriate classification type in a graduated style choropleth map?

texts: I have to provide choropleth maps for for Energy savings, Emission savings, Cost Savings related to different administrative boundaries such as provinces and states, but I don't know which classification is more appropriate because my data for each province is separated from other parts and there is no link between the numbers. Any suggestion for how should I classify? I don't want to use equal intervals because the minimum is 20 and the maximum is 400 and there are few over 300 figures so equal intervals represent my data not in a good way, but about other methods such as jenks or deviation or quantile I have no idea!


texts: How to create 3D water depth with ArcScene?

texts:
I was trying to visualize flood depth in 3D using ArcScene. What I had for this purposes are:

high-resolution DTM of the valley
orthophoto and
Raster with flood depth information.

I have draped the orthophoto over the DTM and it looks like as follows:

Then I draped orthophoto over it:

Finally I overlayed flood map on top of orthophoto:

My question is , is there any way to visualize flood water in 3D with height. it looks like Flood depth raster only touches the base of the DTM. But I need to incorporate DTM + raster height in Flood depth map . How can I do it ?


texts:
I'm using ODL Studio to draw US zip codes, and though I can download the free database and get it to load, I can't figure out how to draw boundaries. The tutorial provided on their site explains how to do UK ones with a corresponding polygon shape file.

Is there a similar file for US zip codes?


texts: Configuring pop-up in ArcGIS Online web map of spatial intersection of located address and feature?

texts:
I have a simple web-map with a polygon layer that contains 5 features.  I would like to configure the web-map so that when the user geocodes an address that the the  information of the of the polygon feature appears in the pop-up of the located address.

Is there a way to accomplish this in ArcGIS Online?  If not, what is a viable and simple solution that requires little to no JavaScript modification?


texts: Displaying fields from geodata service using c#

texts:
I want to display my geodata stored in ArcGIS Server.
I did the connection to the server but the problem is the parameter versionName.  I will not use ArcSDE so I have ignored versionName.
Can someone help me to find the solution?
ESRI.ArcGIS.GISClient.IAGSServerObject so = GetSelectedServerProxy();

ESRI.ArcGIS.GeoDatabaseDistributed.IGeoDataServer gdso = (ESRI.ArcGIS.GeoDatabaseDistributed.IGeoDataServer)so;

//++ Queries for subfields if the Get Fields checkbox is enabled
ESRI.ArcGIS.Geodatabase.IQueryFilter qf = new ESRI.ArcGIS.Geodatabase.QueryFilterClass();

qf.WhereClause = "";

//++ based on the start index and count entered.
ESRI.ArcGIS.Geodatabase.IResultPortionInfo range = new ESRI.ArcGIS.Geodatabase.ResultPortionInfoClass();
        range.StartIndex = 0;
        range.Count = 500;

//++ Run a tablesearch based on the version, table name, fields selected, startindex and count
ESRI.ArcGIS.GeoDatabaseDistributed.IGDSQueryResultPortion res = gdso.TableSearch(versionName,"tabl", qf, range);


texts: Flow accumulation raster loading in QGIS 2.4.0 with null values (but flow direction loading ok!)

texts: I have been tasked with creating an open-source hydrology dataset for a community group. To do this I have exported the flow direction and flow accumulation rasters from ArcMap as TIFF files and then imported them into QGIS. However, a number of the flow accumulation rasters fail to import properly giving values of 0 rather than the range in the original ArcMap raster file. Can anyone help?


texts:
What is the smoothest route, the easiest user experience, to add gpx files to ArcMap?

Ideally, I'd like to interact with .gpx with the same facility as .shp, being able drag'n'drop onto the map canvas and so on. Ultimately the results will be feature classes in a file-gdb but shapefile is an acceptable intermediate.


texts: Free .exe for coordinate transformation lat-long - UTM

texts: I need to find a simple program for Windows which transforms lat-long to UTM. I know that there a websites to do this, but this is for use with the emergency services. The forestry commision use gps hardware which gives coordinates in Lat-Long. The software back in the Emergency central have to use UTM. The software doesn´t transform coordinates. They have asked me to find a solution to do a quick conversion. I wouldn´t like to put lives at risk by using a website which could go offline just when we need it, so I would prefer an installed program. Any good ones out there?


texts: How to read projection information from .dwg file

texts:
I have got an AudoCAD .dwg file which is based on some location of Sweden, but I don't know what is the internal projection of this file. How can I get projection information from this .dwg file?

Actually when I convert this .dwg file into a shapefile by using FME software, then it requires projection information (source file projection) for converting it into shapefile (the target format), and I don't know what is the projection of this .dwg file. Is there any way to get to know about the projection of source file .dwg ?


texts:
I am trying to get bike lane lengths from the following data set.

To do so I have tried using the add geometry attribute from the processing toolbox.

I have tried both calculate using Layer CRS and Ellipsoidal and gotten similar untrue results.

For instance, in running this it says one bike lane segment is only .08 km when it should be 1.18 km.

Any guesses for how to fix this (accuracy)?


texts: Remove Schema Lock after removing layer from TOC using python

texts:
I have a python addin tool that looks for layers in the TOC using ArcMap. I need to be able to replace these layers with new layers of the same name in the same workspace. I'm working with shapefiles. The lock file is not releasing when I remove the layer. I use arcpy.mapping.RemoveLayer(), then with the os module, go into workspace where the layer was located and use os.remove() to delete the files. Unfortunately this keeps failing due to the schema lock. Is there a way for ArcMap to release the lock. The lock is released if I right-click the shapefile and remove it from the TOC, but I need to wrap this functionality within the addin tool using python.

layerPath = layer.workspacePath
layerName = layer.name
arcpy.mapping.RemoveLayer(df, layer)
arcpy.RefreshActiveView()
for file in os.listdir(layerPath):
    if file.startswith(layerName):
        os.remove(os.path.join(layerPath, file))

How do I remove the schemalock?


texts: Get the Z values of existing lineStringZ with PyQGIS?

texts:
I can get the values of 3D QgsPoints with mypoint.get().z() When I iterate through a layer to get a specific polyline and then iterate through the polyline to get the start and end points- start_point = QgsPoint(geom[0]) I get a QgsPointXY and no Z values.   They are there until the iteration geom[0] and then poof gone.

Is there a way to include the z values with something other than QgsPoint(geom[-1])?

My related question is Changing the coordinate value of a Pointz object with python - which is actually how I got the points to create the polyLine

I am taking this route next
 Adding Z value from a Point field to a LineString using pyqgis

Below is the code I ran in the Python Console.   The output is also below.

layer = iface.activeLayer()
print(layer)
all_features = layer.getFeatures()
layer.startEditing()
for pipe in all_features:

        if pipe['name'] == '1212_1211':
            print('found it ********\n  ')
            pipe_id = pipe.id()

            geom = pipe.geometry()
            print(geom)
            new_feat = QgsFeature()
            geom = pipe.geometry().asPolyline()
            start_point = QgsPoint(geom[0])
            end_point = QgsPoint(geom[-1])
            print (start_point)
            print (start_point.x())

OUTPUT

<qgis._core.QgsVectorLayer object at 0x000001767C312EE8>
found it ********

<QgsGeometry: LineStringZ (26138560.91000000014901161 643574.7099999999627471 778.20000000000004547, 26138558.94000000134110451 643563.81000000005587935 0)>
<QgsPoint: Point (26138560.91000000014901161 643574.7099999999627471)>
26138560.91

UPDATE

I found a great chunk of code from this site that made this work for QGIS 3 - geom.constGet() is the key below. How to access shapefile m-values with pyqgis?

layer = iface.activeLayer()
print(layer)
all_features = layer.getFeatures()
layer.startEditing()
for pipe in all_features:

    if pipe['name'] == '1212_1211':
        print('found it ********\n  ')
        pipe_id = pipe.id()

        geom = pipe.geometry()
        line = geom.constGet() #line returns a QgsLineString object

        z = line.zAt(0) # z value of first vertex
        print(str(z))


texts: Error 000875 when using setting workspace with user input for slope spatial analyst tool

texts:
Why does the slope tool run fine until I use a user input to set the workspace?

Background: Below you will see the code that works fine at the top preceded by ##. When I change the workspace to be defined by a user input as shown beneath that, it no longer works. This is part of a larger script. In the larger script the workspace is set with the user input and it works fine with mosaic to new raster, project raster, and clip raster, but then it errors on the slope tool. I am now testing in the shorter script shown below so I don't have to wait for it to run each time.

I have tried changing drives, shortening the path, I switched everything from Esri Grid to tif because of the naming limitations, I tried a shorter name for the output raster.

##import arcpy
##from arcpy.sa import *
##arcpy.env.workspace = r'E:\My_Workspace\Processed_Products'
##
##
##
##if arcpy.CheckExtension('spatial') == 'Available':
##    arcpy.CheckOutExtension ('spatial')
##    outraster = Slope('DEM_Clip.tif',"PERCENT_RISE",1)
##    outraster.save('DEM_Slope1.tif')
##    print "The slope dataset has been generated"
##else:
##    print 'spatial analyst not available'

import arcpy
from arcpy import env
from arcpy.sa import *

arcpy.env.overwriteOutput = True

ws = raw_input('Please paste the workspace path here. ')
arcpy.env.workspace = "r" + "'" + ws + "'"

print " ws is " + ws

DEM_Clip = ws + "\\Processed_Products\\DEM_Clip.tif"

print "The file below has been created and it is in the Processed_Products folder"
print DEM_Clip

if arcpy.CheckExtension('spatial') == 'Available':
    arcpy.CheckOutExtension ('spatial')
    outraster = Slope(DEM_Clip,"PERCENT_RISE",1)
    outraster.save('DEM_Slope1.tif')
    print "The slope dataset has been generated"
else:
    print 'spatial analyst not available'

arcpy.CheckInExtension("spatial")

Output:

Please paste the workspace path here. E:\My_Workspace
 ws is E:\My_Workspace
The file below has been created and it is in the Processed_Products folder
E:\My_Workspace\Processed_Products\DEM_Clip.tif

Traceback (most recent call last):
  File "G:\GIS_Programming\488_Data\Lecture_10_Scripts\Lecture_10_4.py", line 33, in <module>
    outraster = Slope(DEM_Clip,"PERCENT_RISE",1)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\sa\Functions.py", line 6054, in Slope
    z_factor)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\sa\Utils.py", line 53, in swapper
    result = wrapper(*args, **kwargs)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\sa\Functions.py", line 6049, in Wrapper
    z_factor)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\geoprocessing\_base.py", line 506, in <lambda>
    return lambda *args: val(*gp_fixargs(args, True))
ExecuteError: ERROR 000875: Output raster: r'E:\My_Workspace'\Slope_DEM_Cl1's workspace is an invalid output workspace.
ERROR 000581: Invalid parameters.
Failed to execute (Slope).


texts: Framing of WFS GetFeature request using polygon intersect filter

texts:
I am trying to frame a Web Feature Service GET Request applying the INTERSECT filter but I am getting the following error. :

<ows:ExceptionReport xmlns:ows="http://www.opengis.net/ows" version="1.1.0"
language="en">
<ows:Exception exceptionCode="NoApplicableCode">
<ows:ExceptionText>Wrong coordinates.</ows:ExceptionText>
</ows:Exception>
</ows:ExceptionReport>

The request that I created is as follows :

"fluggs.wupperverband.de/secman_wss_v2/service/WFS_WV_Oberflaechengewaesser_EZG/guest?Request=GetFeature&typeName=secman_WFS_WV_Oberflaechengewaesser_EZG:Fliessgewaesser&Filter=<ogc:Filter><ogc:Intersects><ogc:PropertyName>Shape</ogc:PropertyName><gml:Polygon><gml:exterior><gml:LinearRing><gml:posList>51.293792 7.051833 51.309247 7.280830 51.218159 7.242377 51.205900 7.043937 51.293792 7.051833</gml:posList></gml:LinearRing></gml:exterior></gml:Polygon></ogc:Intersects></ogc:Filter>"

Where am I going wrong?


texts: How to update feature class from feature layer with ArcPy Search and Update Cursor?

texts:
I have a script where I select point features within a mile of major roads. If the feature falls within a mile of a major road then the field within_mile_road is updated with yes.  If not, it remains null.  I want to update the within_mile_road from the feature layers that I create back to the original feature class. How do I update the values that were written to the feature layer back to the feature class, here is my code below which does everything except update the original feature class.

import arcpy
arcpy.env.workspace = 'workspace'
    roads = 'roadsfc'
    fields = ['within_mile_road']
    for fc in arcpy.ListFeatureClasses():
        arcpy.MakeFeatureLayer_management(fc, fc+"_lyr")
        arcpy.SelectLayerByLocation_management(fc+"_lyr", "WITHIN_A_DISTANCE", roads, search_distance="1 Miles")
        with arcpy.da.UpdateCursor(fc+"_lyr", fields) as cursor:
             for rowFL in cursor:
                rowFL[0] == 'yes'
                cursor.updateRow(rowFL)


texts:
Is there a way to perform a 'reverse' of the clip function in ArcMap?

I am not exactly sure how to explain this so here is a diagram:


texts: Examples of self hosted cartodb implementations

texts:
I am very impressed with the ease of use of CartoDB and the fact that it is based on OGC standards. I think it would be very usefull for none-GIS specialists. However, I would like to know how well it is suited for internal use. Also because we already have many OpenGIS-services available and we would like to host our own data.

Therefore I'm searching for examples of self hosted implementations of CartoDB.

Please let me know if you know organisations who have done their own implementation of CartoDB.


texts: Character encoding for Oracle DB and Geoserver

texts:
I have a setup that uses Oracle database->GeoServer->OpenLayers->my JavaScript.

I use the UTF-8 character set in GeoServer and the JavaScript layers.
The Oracle instance is configured to use WIN1252 (NLS_CHARACTERSET=WE8MSWIN1252).

So, when I make a request from the JavaScript with the CQL_FILTER with characters such as ç, à, .. I get no results from the database because of the encoding.

What is the best way to deal with this problem of multiple character sets? My opinion is changing the database charset to UTF-8, but unfortunately is to risky, because I'm not the only one using it.

Any suggestions?


texts: PostGIS ST_Intersects or ST_Within for Point in Boundary

texts:
When I query with ST_Intersects or ST_Within I get the same result. I want to know which one to use for performance wise? Note: For points, there is no difference.
Using ST_Within:
EXPLAIN ANALYZE
SELECT id, name,short_name,type FROM locations WHERE ST_Within('SRID=4326;POINT(9.1011853 49.5994825)',bounding_polygon)

QUERY PLAN
"Index Scan using locations_bounding_polygon_idx on locations  (cost=0.29..19727.17 rows=1 width=76) (actual time=0.326..1.389 rows=8 loops=1)"
"  Index Cond: (bounding_polygon ~ '0101000020E610000001AC448FCE3322401D03B2D7BBCC4840'::geometry)"
"  Filter: st_within('0101000020E610000001AC448FCE3322401D03B2D7BBCC4840'::geometry, bounding_polygon)"
"  Rows Removed by Filter: 18"
"Planning Time: 0.126 ms"
"Execution Time: 1.405 ms"

Using ST_Intersects:
EXPLAIN ANALYZE
SELECT id, name,short_name,type FROM locations WHERE ST_Intersects(bounding_polygon, 'SRID=4326;POINT(9.1011853 49.5994825)')

QUERY PLAN
"Index Scan using locations_bounding_polygon_idx on locations  (cost=0.29..33.30 rows=1 width=76) (actual time=0.270..1.102 rows=8 loops=1)"
"  Index Cond: (bounding_polygon && '0101000020E610000001AC448FCE3322401D03B2D7BBCC4840'::geometry)"
"  Filter: st_intersects(bounding_polygon, '0101000020E610000001AC448FCE3322401D03B2D7BBCC4840'::geometry)"
"  Rows Removed by Filter: 18"
"Planning Time: 0.231 ms"
"Execution Time: 1.116 ms"

ST_Intersects: https://postgis.net/docs/ST_Intersects.html If a geometry or geography shares any portion of space then they intersect.
ST_Within: https://postgis.net/docs/ST_Within.html Returns TRUE if geometry A is completely inside geometry B


texts: ogr2ogr clip all shapefiles in a directory

texts:
I am new to ogr2ogr. Does anyone know if it is possible to loop through all shapefiles in a directory, clip and output them into an identical directory with a suffix "_clipped"?

I saw this similar question about converting all files in a directory. Problem here is that using this method I overwrite all my files.
Using ogr2ogr to convert all shapefiles in directory?

EDIT:

I want to make a .bat script to execute the code. This is how far I've come. It works if I output to the same folders as input, but I want it to output to some other folder. The user enters a new folder name, to where the clipped files will be saved.

Here is the returned errors:
Error 1: Failed to create file .shp file
Error 4: Failed to open shapefile 'correct path'

--
set parent=%~dp0
set ogr2ogrPath="C:\Program Files\QGIS 2.16.1\bin\ogr2ogr.exe"

set /p inFolder=Enter the name of the root folder containing the shapefiles to clip:
set /p outFolder=Enter a new name for the folder to save the clipped files:
set /p xMin=Enter xMinCoord:
set /p yMin=Enter yMinCoord:
set /p xMax=Enter xMaxCoord:
set /p xMax=Enter yMaxCoord:

for /R %parent%%inFolder% %%f in (*.shp) do %ogr2ogrPath% -skipfailures -clipsrc xMin yMin xMax xMax %parent%%outFolder%\%%~nf.shp %%f


texts: Exporting feature class to geometry objects and preserving attributes - ArcGIS, ArcPy

texts:
I'm using ArcPy and would like to copy features in a feature class to geometry objects and then perform a buffer on the geometry objects and copy them back to a feature class.

I am doing this in hopes that it will provide the performance and speed improvements that will allow a tool to run on larger datasets that are currently being hung up and/or failing in the normal buffer tool.

My question is, what is the best strategy for preserving an attribute field in this process? I need to preserve at least one ID field to join the original fields back to the buffered shapes after converting back to a feature class. Is this possible?

Here is my process so far for getting the geometries buffered:

import arcpy
i = 'PCPP_sites'
o = 'C:\buffer_test'
geometries = arcpy.CopyFeatures_management(i,arcpy.Geometry())
buffer_geom = []
for geom in geometries:
    b = geom.buffer(1)
    buffer_geom.append(b)
arcpy.CopyFeatures_management(buffer_geom, o)

How do I preserve attributes?


texts: I am using PostGIS and python3. I have a table with polygons objects, the polygons have year as an attribute. I need to identify the uniqueness of the polygons trough time. For instance I need to identify if polygon x from year y is most likely to be which polygon in another year. I realize that this is not an easy question, but any help or thought is more than welcome. I though on doing it by the percentage of interception for instance, but I am ot sure if it is the best analysis.


texts: Projection-related issue when extracting pixels from Landsat scenes?

texts:
I have some Earth Engine code that extracts surface reflectance pixels from Landsat scenes.
var imageCollection = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR");
var myfc = ee.FeatureCollection("users/adrianom/Meadow_17174");

// Function to cloud mask from the pixel_qa band of Landsat 8 SR data.
function maskL8sr(image) {
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = 1 << 3;
  var cloudsBitMask = 1 << 5;

  // Get the pixel QA band.
  var qa = image.select('pixel_qa');

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
      .and(qa.bitwiseAnd(cloudsBitMask).eq(0));

  // Return the masked image, scaled to reflectance, without the QA bands.
  return image.updateMask(mask).divide(10000)
      .unmask(-99).select("B[2-7]*")
      .copyProperties(image, ["system:time_start"]);
}

var ic = imageCollection.map(maskL8sr)
                        .select(['B2', 'B3', 'B4', 'B5', 'B6', 'B7'],
                                ['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2']);

var extractPixelsFromPolygons = function(feature){
  var ic_filtered = ic.filter(ee.Filter.calendarRange(2015, 2015, 'year'))
                      .filter(ee.Filter.calendarRange(12, 12, 'month'))
                      .filterBounds(feature.geometry());

  // Is this the right way to align???
  var proj = ic_filtered.first().projection();
  var latlon = ee.Image.pixelLonLat().reproject(proj);
  //

  var coords = latlon.select(['longitude', 'latitude'])
                 .reduceRegion({
      reducer: ee.Reducer.toList(),
      geometry: feature.geometry(),
      scale: 30
    });

  var lat = ee.List(coords.get("latitude"));
  var lon = ee.List(coords.get("longitude"));
  var point_list = lon.zip(lat);
  var pixelCoords = ee.Geometry.MultiPoint(point_list);

  var valuesAndCoords = ic_filtered.map(function(img){

    var pixelValues = img.reduceRegion({
      reducer: ee.Reducer.toList(),
      geometry: pixelCoords,
      scale: 30
    });

    // Fields derived from "myfc"
    var metaDictionary = ee.Dictionary({
      "Time/Date": ee.Date(img.get('system:time_start')),
      "UniqueID": feature.get('UniqueID'),
      "Area_m2": feature.get('Area_m2'),
      "Area_km2": feature.get('Area_km2'),
      "Landsat_ID": img.id()
    });

    var pixels_and_coords_dict = metaDictionary.combine(coords).combine(pixelValues);
    var new_feature = ee.Feature(null, pixels_and_coords_dict);
    return(new_feature);
  });

  return(valuesAndCoords);
};

var fc = myfc.map(extractPixelsFromPolygons).flatten().sort('UniqueID');
print(fc)

It filters the Landsat image collection by year, month, and by a geometry object. It then takes this image collection and grabs the first image's projection and uses this information to create a two-band image: one band representing latitude and the other longitude (I do this to essentially align all the Landsat pixels in different scenes together so that the coordinates are the exact same for each scene extraction). I use this image to get the pixel coordinates of my geometry object. I take these coordinates and use them to extract the pixel information from the Landsat scene. Finally, I export the coordinates, the pixel values, and some attributes from my polygon as a CSV. This works fine and dandy (feel free to try it yourself). However, I encountered what I believe to be a projection-related issue. When exporting the pixel values and their respective coordinates, take a look at the fourth feature in my feature collection that gets printed to the console. I noticed that I was extracting 254 pairs of coordinates, but only extracting 252 pixel values. They should be the same number: one pair of coordinates for one pixel value. Does anyone know how to fix this?


texts: Controlling features rendered with ArcGIS Server service?

texts:
Question: Using a mapping service set up with ArcGIS Server, how can I set up a parameter that limits which features are being rendered?  e.g. "/service?organization=megacorp"?

Background: I am working on a web mapping application which will be used by organizations.  They all have their own data, but it is stored in the same tables.

When they log in, they will be brought to a map which shows their data points.  Due to existing features, our preference is to store this data in a single table:

points_table
id | organization | name
------------------------------
1  | the state    | poor house
2  | megacorp     | tenements
3  | union carbide| bhopal

It looks like the proper way to secure this service will be to use windows integrated authentication, and although I'm not clear on the details, use that to control which features are being rendered.


texts:
Any chance somebody is working on plugin for QGIS or know about any project trying to get Timelapse Google Earth Engine (imagery) into QGIS?

There was a major update recently and it looks really useful for many implementations. If you don't know what I am talking about:

https://earthengine.google.com/timelapse/

"We sifted through about three quadrillion pixels—that’s 3 followed by 15 zeroes—from more than 5,000,000 satellite images … We took the best of all those pixels to create 33 images of the entire planet, one for each year. We then encoded these new 3.95 terapixel global images into just over 25,000,000 overlapping multi-resolution video tiles."


texts:
I am working on a raster, but I can't set a correct region for it.
In the Set region to match this raster map field, I select my raster map to change current region, but its resolution does not change. once I run an analysis command on my raster, I see the following error:

 elevation file's resolution differs from current region resolution!

please help me to change current region resolution to match my raster map.


texts: I already have found how to import GPX files into PostgreSQL using gdal/org2org. Now I am wondering of dumping into GPX files the results of my modification on dataset stored into PostgreSQL. But I have not found how to do so. Provided I have kept the original structure of tables generated by ogr2org, how can I dump their contents into a GPX file?


texts: Can I create a feature class capable of on-the-fly feature additions?

texts: I am trying to select a subset of roads from a line feature class. The selection is too complicated to use the Select By... option, so I am essentially handpicking several hundred lines to add to a new feature class. This wouldn't take long except one mis-click can clear my selection. To avoid this, I could create selection layers every 100 features or so and then load their data into a master feature class, but it would be much easier to have an ongoing feature class into which I can "drop" or "load" features one at a time, like saving a song to a playlist. Is there an option for this in either ArcPro 2.6 or ArcMap 10.8?


texts: how to project a flat image to spherical mercator EPGS:900013 to use in a Image Layer in openlayers

texts:
I have a flat 400x400 pixels image which represents an area of 400x400 kilometers. I need to draw this image as an overlay on my OpenLayers mapping system. I use OpenStreetMap/Google Maps as base layers. I want to draw the overlay image using an Image Layer and I want to center my overlay image at a specific lat/lon point. The projection of the base layer is EPSG:900913 which is, in fact, the projection used by OpenStreetMap and Google Maps.

What I have understood (please correct me if I'm wrong), is that I can't use a flat, 'unprojected' image as an Image Layer, but I have to do a sort of 'reprojecting' to adjust the image to be compliant with the projection of the map, which is EPSG:900913. Searching here and the net, it seems I should use a library/tool called GDAL to do this, but whereas I am completely new to this I ask you:

Can anyone address me on how to use it?
Could you possibly tell me if what I have just said is correct and the way I described here is valid?
How could I use gdal to translate my flat unprojected 400x400 pixels image to an image projected in EPSG:900913, and suitable to be put as an Image Layer as an overlay?

Thanks.


texts:
Is it possible to show a title "ProjectName" within the map generated by qgis2web? I can put a title in above the map described here and tried to follow the answer here to put it above the zoom buttons. But i cant find the lines of code mentioned in the answer such as var title = new L.Control();
Using QGIS 3.16

Edit. I want to do this using openlayers but want the same effect when exporting with leaflet.


texts: QGIS $area field calculator always produces NULL

texts:
Error occurs under Windows in QGIS 2.8.3 and 2.10.1, answers to similar questions e.g. 57418 have not helped.

The shapefile in question was created as follows:

Vector/GPS/GPS tools/Load GPX file (Routes)
Layer / save as ESRI shapefile [....\epsg4326\routes.shp]

The failing sequence is:

Start a new map/add vector layer [routes.shp]
Open attribute table
Add new column AREA [decimal real width=17 precision=8]
Manually edit both AREA values to verify format.
Select all entries [2 entries]
Field calculator / existing /AREA / geometry / $area

Always get NULL

I have tried countless combinations of width and precision for AREA. I have uninstalled/reinstalled/installed different releases. All to no avail.

What's up Doc?

Beginning to smell a rat in the direction of either the GPS tool or the shapefile. Starting with the shapefile, how may I independently verify it?


texts: Request entire data extent independently from view extent to style GeoServer layer

texts:
I created a modified PointStacker which stacks/clusters the data features according to an indicated attribute of the data set (e.g. state or district) which is a user-indicated parameter in the sld style file.

When I am looking at the entire map extent, it works fine. For instance, I have 20 states and for each of them a cluster point with its number of clustered features is displayed.
However, when I zoom into the map, and a state is not completely within the view anymore, the input data for this state is less than before, thus, the number of clustered features for this state is getting less as well - which is not my intended behaviour for this case. I do not want to cluster only the features within the view, but to cluster the entire data extent - independetly from my zoom level, thus view extent.

Is there any way to accomplish something like a constant data bounding box to request data from independently from zooming? Thanks a lot in advance.

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0"
                       xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
                       xmlns="http://www.opengis.net/sld"
                       xmlns:ogc="http://www.opengis.net/ogc"
                       xmlns:xlink="http://www.w3.org/1999/xlink"
                       xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <NamedLayer>
    <Name>xyz</Name>
    <UserStyle>
      <Title>xyz</Title>
      <FeatureTypeStyle>
        <Transformation>
          <ogc:Function name="gs:PointStackerAttribute">
            <ogc:Function name="parameter">
              <ogc:Literal>data</ogc:Literal>
            </ogc:Function>
            <ogc:Function name="parameter">
              <ogc:Literal>clusterattribute</ogc:Literal>
              <ogc:Literal>state</ogc:Literal>
            </ogc:Function>
            ...
            <ogc:Function name="parameter">
              <ogc:Literal>outputBBOX</ogc:Literal>
              <ogc:Function name="env">
                <ogc:Literal>wms_bbox</ogc:Literal>
              </ogc:Function>
            </ogc:Function>
            <ogc:Function name="parameter">
              <ogc:Literal>outputWidth</ogc:Literal>
              <ogc:Function name="env">
                <ogc:Literal>wms_width</ogc:Literal>
              </ogc:Function>
            </ogc:Function>
            <ogc:Function name="parameter">
              <ogc:Literal>outputHeight</ogc:Literal>
              <ogc:Function name="env">
                <ogc:Literal>wms_height</ogc:Literal>
              </ogc:Function>
            </ogc:Function>
          </ogc:Function>
        </Transformation>

        <Rule>
        ...
        </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>


texts: howto set ogr2ogr imported SRID to 0? Currently defaulting to SRID 1

texts:
I'm currently trying to import a bunch of boundaries into a legacy system with a specific SRID. its a mySQL db.

After import and running some intersect code, i receive this error:

Binary geometry function st_intersects given two geometries of different srids: 1 and 0, which should have been identical.

As it says, there is an srid mismatch of 1 vs 0. The newly imported boundaries are SRID 1, the existing boundaries from another table (cannot be changed) are SRID 0. I have established this with mysql command

select ST_SRID(boundary_polygon) from boundaries;

My question, what EPSG value would i use with ogr2ogr's -a_srs flag. I want to force SRID=0?

Here is the current command i'm using for import:

ogr2ogr -f "MYSQL" MYSQL:"blahdb,host=localhost,user=blah,password=blah,port=3306" boundaries.tab

I'm currently running this command which forces the import to SRID=0 which seems by convention to mean No SRID.

UPDATE boundaries SET boundary_polygon = GeomFromWKB(boundary_polygon, 0);

This then allowed my st_intersects method to work fine and find intersect between two polygons.

However, I still need a way to import as SRID=0 from ogr2ogr instead of SRID=1. So that I don't need this post import SQL step


texts: Bad data beaking IRelationshipclass.GetObjectsRelatedToObjectSet

texts:
I have two workspaces - one a File GDB, the other SDE (SQL Server).  Both workspaces contain the same data with the same relationship classes.  When identifying features in ArcCatalog, the results pane correctly presents the related records.

When using IRelationshipclass.GetObjectsRelatedToObjectSet, I always get an empty result set.  If however I use a loop, calling IRelationshipclass.GetObjectsRelatedToObject, I get the related records.  I see the same result in both workspaces.

I have asked someone with a known good dataset to run my code, and they report that it works perfectly.
So my question is this - what in my workspace could cause this behavior? I need to rectify it.

For the record, I'm using ArcGIS 9.3 SP1.  I've included my test code below, on the chance that there is an issue with my code.

public void test(IWorkspace workspace, IQueryFilter qf, string relClassName, string sourceTableName)
{

    IFeatureWorkspace featureWorkspace = (IFeatureWorkspace)workspace;
    IRelationshipClass relationshipClass = featureWorkspace.OpenRelationshipClass(relClassName);
    ITable originTable = featureWorkspace.OpenTable(sourceTableName);
    var sourceSet = originTable.Search(qf, false);

    var foundSet = new SetClass();
    IRow row;
    while ((row = sourceSet.NextRow()) != null)
    {
        foundSet.Add(row);
    }

    var resultSet1 = GetRelatedRecordsByObject(relationshipClass, foundSet);
    var resultSet2 = GetRelatedRecordsBySet(relationshipClass, foundSet);

    System.Diagnostics.Debug.WriteLine("Result1: " + resultSet1.Count);
    System.Diagnostics.Debug.WriteLine("Result2: " + resultSet2.Count);
}

public ISet GetRelatedRecordsBySet(IRelationshipClass relationship, ISet sourceObjectList)
{
    sourceObjectList.Reset();
    var resultSet = relationship.GetObjectsRelatedToObjectSet(sourceObjectList);

    resultSet.Reset();
    return resultSet;
}

public ISet GetRelatedRecordsByObject(IRelationshipClass relationship, ISet sourceObjectList)
{

    ISet resultSet = new SetClass();

    sourceObjectList.Reset();
    IObject sourceRow;
    while ((sourceRow = (IObject)sourceObjectList.Next()) != null)
    {
        var indSet = relationship.GetObjectsRelatedToObject(sourceRow);

        indSet.Reset();
        IRow indSetRow;
        while ((indSetRow = (IRow)indSet.Next()) != null)
        {
            resultSet.Add(indSetRow);
        }
    }

    resultSet.Reset();
    return resultSet;

}


texts: Updating pgsnapshot database with osmosis - unable to read state from the server

texts:
I set up a local postgis database with osmosis using pgsnapshot schema as described here.

The problem occurs when I try to update the database with latest minutely change sets as described here.

When I run the command osmosis --rri --wxc something.gz (which is supposed to put all the changes between the previous sequence number and the most up-to-date into the file something.gz ) I get the following error:

SEVERE: Thread for task 1-rri failed
org.openstreetmap.osmosis.core.OsmosisRuntimeException: Unable to read the state from the server.
        at org.openstreetmap.osmosis.replication.common.ServerStateReader.getServerState(ServerStateReader.java:116)
        at org.openstreetmap.osmosis.replication.common.ServerStateReader.getServerState(ServerStateReader.java:50)
        at org.openstreetmap.osmosis.replication.v0_6.BaseReplicationDownloader.runImpl(BaseReplicationDownloader.java:290)
        at org.openstreetmap.osmosis.replication.v0_6.BaseReplicationDownloader.run(BaseReplicationDownloader.java:383)
        at java.lang.Thread.run(Thread.java:745)
Caused by: java.io.FileNotFoundException: http://planet.openstreetmap.org/minute-replicate/state.txt
        at sun.net.www.protocol.http.HttpURLConnection.getInputStream(HttpURLConnection.java:1625)
        at org.openstreetmap.osmosis.replication.common.ServerStateReader.getServerState(ServerStateReader.java:97)
        ... 4 more


texts: Loading style from GeoPackage with PyQGIS

texts:
Reading the documentation I do not understand how the the row in layer_styles is specified (and why a table name has to be given)...

layer.loadNamedStyleFromDatabase('/path/to/geopackage/','layer_styles', 'styleName')

... doesn't do the trick.


texts:
I have used pivot table in ArcToolbox of ArcGIS, while I did it, this of my result of summary table, the field is text format, but however it is has numeric value (year value), after using pivot, the value is not shown as it should be, here is my screen shot:

I guess the ArcGIS can't using field with numeric name?, why is this happen?, is this a bug?, Thn_Buat_u is in string format, but why it was changed in pivot result?


texts: ArcGIS Python Field Calculation says failure during processing?

texts:
I have a problem with the field calculator in ArcGIS:

I have:

A complete table of several states of the world with their specific life expectancies.
A world map (not yet joined)

What i have to do:

classify the states in 5 groups by their life expectancies in an extra field in the field calculator.

What i did:

I generated the following code block:

def reclass(a):
    if (a <= 50):
        return 1
    elif (a <=60):
        return 2
    elif (a <=70):
        return 3
    elif (a <=75):
        return 4
    elif (a >75):
        return 5
    else
        return 999

Expression:

reclass(!LIFE_EXP00!)

I added a new field called "Field_1" and i took LIFE_EXP00 from the already given data in the table.

But still ArcGIS tells me after every single time i try it, that there was a failure during processing. What am I doing wrong?


texts: How to overlay Historic map onto Google Earth using QGIS?

texts:
I am very new to QGIS. I have a Plat Map of my home town from the early 1900's that I would like to overlay on Google Earth.

The problem I'm having is the map will not line up. I can get roads to line up in one section but not the whole map. I believe this is due to to a poor scan that is not flat in all areas (not sure).

Is there a way to manipulate the map by referencing multiple areas on both maps? I've read some tutorials on georeferencing, but I'm not really sure if this is the right direction or how to do it.


texts: Deleting records from multiple, similar feature classes by attribute values using SQL statement in ArcGIS Desktop?

texts: I have several hundred geodatabase feature classes which all have the exact same fields. I need to go through each of those feature classes and delete any records which match an SQL query (the same SQL query can be applied to all the feature classes). Alternatively, I could export them all to new feature classes, keeping the records I want. In ArcCatalog I can export them singly using an SQL query to specify which records get exported, but the SQL statement option isn't available when you want to export multiple layers simultaneously (using "Feature Class to Geodatabase (multiple)").


texts: String column in QGIS shapefile that is read as a factor column in R

texts:
I have a shapefile that has 3 columns (ID_1, NAME_1 and Rainfall). The 2nd column contains the names of governorates (type of variable is string).

I have used the sf package in R to read this shapefile and save it in an RDS format, using the code below:

TUN_sf <- st_read("Data/TUN_adm/TUN_adm1.shp")
TUN_spdf <- as(TUN_sf, "Spatial")
TUN_sf <- st_as_sf(TUN_spdf)
simp_sf <- ms_simplify(TUN_sf, keep = 0.01, keep_shapes = TRUE)
saveRDS(simp_sf, "simp_sf.rds")

However, I noticed that the second column is now considered as a column of factors (even though it was originally a column of string).

I would like to convert this column to string but I don't know how to go about it. I tried simple functions like (as.character) but it distorted the table. This is because the file "simp_sf" is not a simple dataframe but a spatial dataframe (it contains coordinates).


texts: GeoServer converting PostGIS timestamp with time zone

texts:
In a PostGIS view i have a field of data type "timestamp with time zone". In GeoServer this field is show with datatype "Timestamp"

My problem is that i want the Timestamp to be show in my application exactly as it is in my Database as UTM. +00.

What i get back from GeoServer via WFS request is:

  2010-08-20T16:56:15.095-04:00

What is stored in the PostGIS database is:

  2010-08-20 20:56:15.095+00

How do i tell GeoServer not to convert the time so that i see the time exactly how it is in the database as UTM? Is it on Apache, System Time, GeoServer config?

Running GeoServer on latest stable Suse linux Ent.


texts: Iterate through features to use selection as input for ExtractByMask

texts:
I need to iterate through each of the features of a feature class successively to use each one as the mask to extract from a raster with the ExtractByMask tool. Basically I'm looking for the "Iterate Feature Selection" function from ModelBuilder for Python scripting.

If I understand correctly, using a search cursor will return only the row of values from the attribute table, not the feature shape itself.

Would iteratively creating a new shape with "Select_Analysis" be a feasible way to go about this instead?

I suppose creating a new shapefile for each feature would take quite long.

Is there another, more efficient way?


texts: Checking SAGA version from a python plugin in QGIS 3.4.0

texts:
I use QGIS 3.4.0 on Windows and I want to verify SAGA version with  processing.algs.saga.SagaUtils.getSagaInstalledVersion(). But I get an error:
module 'processing' has no attribute 'algs'

This was running on 3.2.0

Any other solution?


texts: How to zoom to the extent of a Qgsvectolayer using PyQGIS?

texts: Anyone know how to zoom to the extent of a Vectorlayer in QGIS using python plugin?


texts: Changing symbol size for layer that is symbolized by unique values in ArcMap

texts:
I frequently symbolize layers by unique values (e.g. first rock type encountered in an exploration drillhole).  This typically results in a layer with 15+ unique rock type all symbolized by the same size of simple/circular marker that is colored differently based on the rock type.  The resulting maps that I print are used for interpretation by hand and a common revision that I have to implement is the re-sizing of these layer files.

Is there a method or script that I can use to change the size of all of these symbols in one step instead of manually modifying the size of each rock type's symbol?


texts: How can I use QGIS to convert points with elevations to DXF?

texts:
I'm trying to convert points (with their elevations) into a dxf file for use in AutoCAD LT.  I've tried "save as" to dxf for the layer, which appears to almost work...the dxf has x and y location information, but the elevations associated with the points are not present.

I also tried the "dxfauth" openware tool, but that did the same as the above.  Thanks!


texts:
I know we can create map with a zoom factor like:

map = new ol.Map({
target: 'container',
layers: [
    layer
],
view: new ol.View({
    maxResolution: res,
    extent: ext,
    zoomFactor: factor
}),
controls: [
   new ol.control.Zoom()
]

});

But is there any way to change this zoom factor on the fly? I want to make custom control to set that. Any function like map.getView().setZoomFactor() ?


texts: Removing Polygons Within Polygons After Merge (QGIS 2.18 OS X)

texts:
I've merged together a polygonized raster layer of a certain elevation grade (>25%) and a wetlands buffer. I can't seem to find a good solution to the problem at hand, however. Ultimately I'd like to run a symmetrical difference tool on the merged layer (mentioned in the former) and a parcel layer to find suitable land to build on.

So far I've tried running v.clean (snap, rmdupl, and bpol, all to no avail). Manually cleaning lines is totally out of the question. Dissolving takes a ridiculous amount of time and computing power. Any better alternatives?

Here's what my data looks like:


texts: How to Have Anonymous Users Edit Survey 123 Data in a Form

texts:
I have a Survey 123 form where users enter the status of a project, then they submit it. Then they need to, at the conclusion of the project, head back to the record and edit it. This editing would include changing the drop-down to closed, entering notes etc.
These are anonymous users, not AGOL accounts. If you have an AGOL account then you can log into Survey 123, head to the survey and in the data tab, voilà there you go you have the form when you click on the record. Is there anyway to accomplish for non AGOL users?
So far, my workarounds include creating a link to the feature service of the Survey 123 form in ArcGIS Classic Web Viewer. Doing this allows them to click on one of the entries and then they can edit it, but it's not the prettiest process and it's cumbersome for non-GIS users.
Essentially, I'd like if they could head back to the same Survey 123 form, with it populated with the existing records information and then they click submit and the record is updated.
Is this possible?


texts: Growing a PostGIS geometry by a percentage

texts:
I'm looking for a way to add a buffer around a PostGIS geometry, but the size of the buffer should depend on size of the geometry. That is, I want each geometry to be enlargened by, say, 5%.

The idea is I'm looking for intersecting geometries, but there could be an error of up to 5% associated with each one which I want to take into account.

Anyone know the best way of going about this?

The database has almsot a million rows, so I'd prefer it to be fairly fast.


texts: gdal.VectorTranslate returns empty object

texts:
Having the following GeoJSON:

{"geometry": {"coordinates": [[[0.0, -89.9999999999], [0.0, -87.5], [120.0, -87.5], [120.0, -89.9999999999], [0.0, -89.9999999999]]], "type": "Polygon"}, "properties": {"grid_cell_index": 1}, "type": "Feature"}

I'm trying to load it with GDAL and to perform a VectorTranslate on it. The objective is to guarantee that is compliant to RFC 7946. My code is running on Python 2.7.x with GDAL 2.2.1:

from osgeo import gdal

srcDS = gdal.OpenEx('input.geojson')
ds = gdal.VectorTranslate('output.json', srcDS=srcDS, format = 'GeoJSON', layerCreationOptions = ['RFC7946=YES', 'WRITE_BBOX=YES'])

When I count the objects on the input:

srcDS.GetLayer(0).GetFeatureCount()

Returning 1.

However, when I do the same count after transforming:

ds.GetLayer(0).GetFeatureCount()

It returns me 0, and also the output.json file is empty. I've inverted my point coordinates but it seems that is not a right-hand error. Any ideas on what could be going on?


texts: Clipping Raster with shapefile boundary using QGIS

texts:
I have quastion about Clip Raster with shapefile boundary in QGIS 3.
I would like to receive 'No data' value in places out of the polygon, but as a result of cliping is 0,0,0. I can't defines a value that should be inserted for the nodata values in the output raster, because this value is in the raster.
Any help?


texts:
I'm trying to extract climate variables in QGIS 3.10.2 for sample points with a raster dataset. The raster comes from ClimateNA (https://adaptwest.databasin.org/pages/adaptwest-climatena) and is available as in ASCII grid (.asc) and NetCDF (.nc) formats.

The raster uses the Projection: Lambert Conformal Conic and Datum: D WGS 1984. The best matching CRS seems to be North_America_Lambert_Conformal_Conic, EPSG:102009

I've set the project CRS to this and verified it for all layers. The raster loads with the CRS undefined and must be set manually, and the sample points, which are actually located in California, are displayed at (0,0). The Google base map appears to be correct.

Does anyone know how I can resolve this?

Another user was able to load it correctly, so hopefully its just a simple step I'm doing incorrectly.


texts:
I am using a custom CRS for a project and whenever I load a .shp file, QGIS (3.8.2) creates a new custom CRS entry. As of today, I have about 10000 entries that are all the same and would like to delete them.
Can I do it all at once? Clicking the minus sign is becoming tedious, since there is about a second of delay between each click.


texts: Displaying a symbol label using an expression in MapServer 7

texts:
I want to label the numbered roads at all scales with a symbol instead of text. Also, I would like to label the roads using an expression so I could use the symbol corresponding to the road number.

Right now, I can display a symbol in a label, but the text is still visible and all roads use the same symbol. I can't find a way to use expressions in a LABEL.

So how could I dynamically use a symbol as label, without text, according to a field attribute that stores the road number?

Here is how a part of my LAYER section looks like:

LABELITEM "NomRte"
CLASS
  NAME 'Autoroute'
  EXPRESSION (('[ClsRte]' eq 'Autoroute')and('[CaractRte]' ne 'Bretelle'))

  LABEL
    MINDISTANCE 750
    STYLE    # since to version 6
      symbol 'a10'
      MAXSIZE 24.0
    END # STYLE
  END # end of label

  STYLE
    WIDTH 2
    COLOR 0 0 0
    OUTLINECOLOR 0 0 0
    MAXSCALEDENOM 500000
  END #STYLE


texts:
I read Regular expression substring for labelling

and it got me curious.  I'm trying to do something similar with OSM content in QGIS.

I would like to label features with English names if those names are stashed in the other_tags attribute.  If English names do not exist in the other_tags attribute field, fall back to the name attribute.

I've gotten only as far as finding whether an English value exists:

"other_tags" LIKE '%"name:en"=>"%'

Perhaps someone could explain to me how to put all this in a conditional check within QGIS, calling the name:en value to be rendered as the label, and if the conditional check fails, use the name field.

I am on QGIS 2.14.

UPDATE:

I read through QGIS Multiple CASE WHEN/THEN Statements for Expression Based Labels

and it does not answer my question.  I am not concerned with how to take multiple attributes and stack them in a label.  I am looking at a single attribute field in OSM called other_tags.  The structure of this tag looks like this:

"name:ar"=>"[some arabic text here]","name:ru"=>"[some russian text here]","name:en"=>"Gaziantep"....etc...

What I'm trying to do is take, for example, Gaziantep out of "name:en" out of the overall other_tags field.

I can use CASE expression to roll over to the regular name field if the conditional check fails, but how do I syntactically obtain the "name:en" subvalue?


texts:
I've searched around  but haven't find anything about this topic.
I have several lines that represent the border between adjacent areas, so each area is delimited by several lines. How can I merge the lines to create the polygons representing the areas? Each line has two columns representing the areas that it separates.

Added:

The problem can be divided in two steps:
- Merge lines in a single unified line (ST_LineMerge, ST_Collect, ST_Union)
- Create a polygon from the unified line


texts:
I am trying to code a type of buffer and I tried to find the source code of the buffer function in the ArcGIS and QGIS files in order to understand how it works in depth when it comes to creating buffer around a complex polygons. I would like to know how the vertex are generated and how the arcs close to the edge are created for optimal coverage.

If anyone has in depth knowledge of the source code of a buffer tool without too many external references to other functions, that would greatly help.


texts: Geoserver SLD styling to align labels at the center vertical of roads

texts:
I've used SLD to style road features in Geoserver. As roads include multiple segments the 'group' attribute has been used to control label repeat of roads. The problem is some labels can't be placed on the centerline of roads. In the image below, label 2 is my desire and label 1 is the problematic one.

I applied the 'PerpendicularOffset' to 0 (default) which should do the job, but no success has been made.

<se:TextSymbolizer>
    <se:Label>
      <ogc:PropertyName>name</ogc:PropertyName>
    </se:Label>
    <se:LabelPlacement>
      <se:LinePlacement>
        <se:PerpendicularOffset> 0 </se:PerpendicularOffset>
      </se:LinePlacement>
    </se:LabelPlacement>
    <se:Fill>
      <se:SvgParameter name="fill">#523735</se:SvgParameter>
    </se:Fill>
    <se:VendorOption name="followLine">true</se:VendorOption>
    <se:VendorOption name="maxAngleDelta">90</se:VendorOption>
    <se:VendorOption name="maxDisplacement">200</se:VendorOption>
    <se:VendorOption name="repeat">200</se:VendorOption>
    <se:VendorOption name="group">yes</se:VendorOption>
</se:TextSymbolizer>


texts: Interpreting Firebug results on ESRI Javascript feature layers in the NET tab

texts:
This was cross-posted to the ESRI Javascript API forum yesterday.

I need help in interpreting what I am seeing and wondering why I'm seeing duplicate features.

I'm using this example: http://help.arcgis.com/en/webapi/jav...efinition.html "Feature layer with layer definition applied" on my own server and viewing the results after opening the page in Firefox. I modified the tooltip code so I can see the OBJECTID. Your results may be different from what I describe here. I used the tooltip to find the OBJECTID on screen and then expanded the URLs to find the same OBJECTID. I'm specifically looking at the shaded area next to the I-35 sign on the map.

When I expand the keys under the URLs in the JSON tab on the NET tab I can see where the featurelayer (fields) is being created (the line starts as GET1?). Then there are 3 more lines that when expanded have a key called "features". The first one has 1 feature (0), the second line has 4 features (0,1,2,3) and the third line has 2 features (0,1).

In the second line with 4 features, the 4th feature has an OBJECTID of 5095 and in the third line the 2nd feature has an OBJECTID of 5095. Does that mean that this feature is being downloaded twice? I'm not moving the map around at all.

There are also URLs that when expanded do not show any features. They appear to show the same keys as the ones I described above but the features key cannot be expanded. These have a size of 2.6kb. I'm not sure what is happening here.

In this example the sizes that are being downloaded are not a big deal. In my own data I have some large features and when I run this example against my own internal data and limit the definition expression to one feature I get six URLs returned. I'm assuming the first two are about the creation of the actual feature layer. The next four seem to be pairs of duplicates. The URLs with a size of 865B have a key called features that appears to be empty. The URLs with a size of 4.2KB have a key called features that is listing the feature I defined.

So, I asked myself if in fact it was downloading the feature twice. Am I interpreting these results correctly?


texts:
I am trying to add ECW support into GeoServer 2.12.1 but no luck so far.

Does anyone have a step by step guide to add the ECW support?


texts: How to make customized radar chart in R?

texts:
I am a beginner in R language and trying to make this picture. I tried radar.
I tried using radarchart (radarchart(data)) but I do not know how to modify the picture.

Here is my data:

A       B      C     D     E     F     G      H      I    Group
-0.61 0.13   0.24  -0.30 -0.12 -1.24  0.74   0.55   0.80   Yellow
1.02  -0.40  0.73  0.17  0.68   1.21  -1.35 -0.84  -1.27   Blue


texts: Calculate Sum on One Field Within Several Different Point Feature Classes?

texts:
I have several point feature classes (eg: the location of mall, liquor stores, etc), that I need to get the count of within a specified city district. Basically, I want to great a modelbuilder model that calculates how many liquor stores and malls there are per city district.

I have already created a spatial join between the points and the polygon (which create separate output tables), but need to collapse them so I may do this tabulation. I have tried a bunch of different things but all of them don't allow you to count by polygon.

If you can point me in the right direction, that would be great!


texts:
I'm trying to learn how to use spatialite databases, so have made a few experimental ones. When I delete them, however, they still show up in QGIS's DB Manager. If I try to open them it comes up with a 'database does not exist' type error, but it will not let me delete them, or update itself. Worse still, I created a new spatialite layer in QGIS and named it the same as an old one (but in a different file location). Now DB manager just shows me the old one.

How can I get DB manager to update itself?

I've tried the obvious refresh button, but to no effect.


texts: PyQGIS: standard and advanced linear referencing possibilities

texts:
I'm aware of the existence of the following QGIS plugins:

QChainage
Event layer
linear referencing
LRS

Although, all of these tools offer some linear referencing capabilities, I find that they lack in terms of their ability to be used directly from a standalone python script, are not part of the base QGIS platform, and are not on par with LRS capabilities found in other GIS software. Some of these features include but are not limited to:

Locate features along a road: The ability to extract events from a vector layer into a data table (ex: CSV) of some sort
Overlay events: Take two different event tables (non GIS) that share the same road_id and combine/merge/split/etc them automatically
Make route event layer: combine an existing road layer with a new event table with added bonus to do right vs. left buffer if negative numbers are present.

My question:

What are my options to perform advanced linear referencing using Python QGIS?


texts: Create arc lines in OpenLayers like Kepler.gl?

texts:
I user Openlayers 6. I want to create an arc line in OpenLayers that connects one coordinate to another.
I have a case where when the distance between the points is too short, the line doesn't curve. Even though in this case I want to create an arc line like the kepler.gl site. How can I make an arc line like that?

My Openlayers Code:
var tileLayer = new ol.layer.Tile({
  source: new ol.source.Stamen({
    layer: 'toner',
  }),
});

var map = new ol.Map({
  layers: [tileLayer],
  target: 'map',
  view: new ol.View({
    center: [0, 0],
    zoom: 2,
  }),
});

var style = new ol.style.Style({
  stroke: new ol.style.Stroke({
    color: '#EAE911',
    width: 2,
  }),
});

var flightsSource = new ol.source.Vector({
  wrapX: false,
  attributions:
    'TES',
  loader: function () {
    var url = 'flights.json';
    fetch(url)
      .then(function (response) {
        return response.json();
      })
      .then(function (json) {
        var flightsData = json.flights;
        for (var i = 0; i < flightsData.length; i++) {
          var flight = flightsData[i];
          var from = flight[0];
          var to = flight[1];

          // create an arc circle between the two locations
          var arcGenerator = new arc.GreatCircle(
            { x: from[1], y: from[0] },
            { x: to[1], y: to[0] }
          );

          var arcLine = arcGenerator.Arc(100, { offset: 10 });
          if (arcLine.geometries.length === 1) {
            var line = new ol.geom.LineString(arcLine.geometries[0].coords);
            line.transform('EPSG:4326', 'EPSG:3857');

            var feature = new ol.Feature({
              geometry: line,
              finished: false,
            });
            // add the feature with a delay so that the animation
            // for all features does not start at the same time
            addLater(feature, i * 50);
          }
        }
        tileLayer.on('postrender', animateFlights);
      });
  },
});

var flightsLayer = new ol.layer.Vector({
  source: flightsSource,
  style: function (feature) {
    // if the animation is still active for a feature, do not
    // render the feature with the layer style
    if (feature.get('finished')) {
      return style;
    } else {
      return null;
    }
  },
});

map.addLayer(flightsLayer);

var pointsPerMs = 0.1;
function animateFlights(event) {
  var vectorContext = ol.render.getVectorContext(event);
  var frameState = event.frameState;
  vectorContext.setStyle(style);

  var features = flightsSource.getFeatures();
  for (var i = 0; i < features.length; i++) {
    var feature = features[i];
    if (!feature.get('finished')) {
      // only draw the lines for which the animation has not finished yet
      var coords = feature.getGeometry().getCoordinates();
      var elapsedTime = frameState.time - feature.get('start');
      var elapsedPoints = elapsedTime * pointsPerMs;

      if (elapsedPoints >= coords.length) {
        feature.set('finished', true);
      }

      var maxIndex = Math.min(elapsedPoints, coords.length);
      var currentLine = new ol.geom.LineString(coords.slice(0, maxIndex));

      // directly draw the line with the vector context
      vectorContext.drawGeometry(currentLine);
    }
  }
  // tell OpenLayers to continue the animation
  map.render();
}

function addLater(feature, timeout) {
  window.setTimeout(function () {
    feature.set('start', new Date().getTime());
    flightsSource.addFeature(feature);
  }, timeout);
}


texts: Can anybody give me insight into the grid JSON structure that Mapbox returns for use with overlays on tiles.  What is the significance of !,#,$?