texts: Is there a QGIS tool or plugin that can identify multiple layers at once?

texts:
For anyone familiar with magic.gov.uk (a UK government web-based GIS) I can shortcut this question by saying that I'm looking for the equivalent of the identify tool. For others, this is a tool which produces a report identifying the value of all the properties of any polygon from all selected layers on a click of the mouse.

So, say I had three layers, one with counties (name, population, area), one with nature reserves (name, area, type of reserve) and one with habitat (type, rarity, conservation status). I would like to be able to click on a point on the map (or draw a polygon), and produce a report which would list all the features from all the layers under that point, together with their attributes.

I'm using QGIS 2.18, but soon 3.


texts: How to downsample a raster using the gdal max method using python?

texts:
I want to downsample a raster image using the max method that is available in GDAL in python.

In other words I would like to recreate the following bash command in python code:

gdalwarp -tr 32 32 -r max input_image.tif resampled_image.tif

I have found the following solution but for some reason this changes the values in the image to a range between approximately 0 and 2 when I load the created tif file.

How can I resample a raster image in python using the max method?


texts: Export each feature to a new geojson individually

texts: I have a shp file in qgis. I would like to export that shp file into geojson. I have seen the save as option. But, my point is to create individual files for each individual shape in the gis project. Is that possible?


texts: Zooming to layer with page definition query set in ArcPy?

texts:
I have data driven pages enabled on a mxd, with a page definition query set on one of the layers to only show features that match the current index feature. This layer looks like this:

After enabling the page definition query, right clicking on the layer and clicking Zoom to Layer yields this:

I would like to perform this function in ArcPy. Setting the data frame's extent to the extent of the layer using the getExtent() method zooms to the extent of the entire layer:

This method does not honour the extent of the page definition query. I cannot use the zoomToSelectedFeatures() method of the data frame because nothing is actually selected; similarly I cannot use the layer's getSelectedExtent() method. The only workaround I see is updating a normal definition query on the layer from ArcPy and then getting the extent.

Is there a way to use the UI's Zoom to Layer method from ArcPy?


texts: Save x,y,z values to a GeoJSON file from QGIS

texts:
I would like to save x,y,z values for polygon vertices as a geoJSON file using Python in QGIS.  I've made several attempts, but all have issues.  I've found that the QgsPoint object must be used as QgsPointXY does not support Z and M values.  Below is a code snippet which includes various tries and their corresponding error messages/issue:
for idx in range(n):
   col = int(round((geom.vertexAt(idx).x() - extent.xMinimum()) / xres))
   row = int(round((extent.yMaximum() - geom.vertexAt(idx).y()) / yres))

   lat = round(self.ur_lat - (row * lat_res), 6)
   lon = round(self.ll_lon + (col * lon_res), 6)
   alt = computeAltFromLatLon
   point=QgsPoint(lon, lat,alt)
   lat_lon_points.append(point)

# Error: type object 'QgsGeometry' has no attribute 'fromPolygonZ'
# lat_lon_poly = QgsGeometry.fromPolygonZ([lat_lon_points])

# Error: type object 'QgsGeometry' has no attribute 'fromPolygon'
#lat_lon_poly = QgsGeometry.fromPolygon([lat_lon_points])

# Error: index 0 has type 'QgsPoint' but 'QgsPointXY' is expected
#    QgsPointXY does not support Z and M values
#lat_lon_poly = QgsGeometry.fromPolygonXY([lat_lon_points])

# File is successfully written here, but only x,y is saved.
lat_lon_poly = QgsGeometry.fromPolyline([lat_lon_points[0],
lat_lon_points[1],lat_lon_points[2], lat_lon_points[3]])

fet.setGeometry(lat_lon_poly)

# print("Write to GeoJSON: ")
writer.addFeature(fet)

Any suggestions for how I can write a geoJSON file with x,y,z values for each coordinate?


texts: Multiple people using one ArcGIS Online Story Map Builder?

texts: Is it possible for multiple people, say a group of students, to collectively work on a story map or to create something with Web AppBuilder?


texts: How to backup my custom defined coordinate systems in Leica Geo Office

texts: How can I use a custom defined cs in Leica Geo Office on another PC (Leica GeoOffice instalation)?


texts: Creating web app using GeoServer and linking ArcGIS Desktop to it?

texts:
I have a database (SQL Server) with a table having lat-long information for sites that keep getting updated daily. I want to create a web application to plot the same on map and add other features to it like buffer creation, shortest distance, heatmap etc.

We have the Standard (editor) license level of ArcGIS Desktop and we have made a geodatabase connection in the same.

How can we proceed to create a web map for the above requirement?

We don't have an ArcGIS Server license.

Is there a way we can publish a service on GeoServer from ArcGIS desktop and then use the same in ArcGIS API for JavaScript considering our data will keep changing on a daily basis?


texts: Attachment from Data Driven Pages - Name Filter

texts:
We are using data driven pages to create a large series of figures.  We used Collector in the field to collect features with data and attachments.

To display the attachments for each feature, we use the "Attachment from Data Driven Pages" and the "Name Filter" option.

The problem we run into is some field devices save attachments as "Attachment1" (Android) or "Photo1" (Apple).  Is there a way to use a wild card or use multiple names in the "Name Filter"?


texts:
Can anyone tell me what all the information on this Geoserver status this means and if there are any issues with it?

I cannot find anything in the documentation that explains what all this is actually telling us.

We have figured out the reload button, but not sure if it is required for EVERY configuration change or not, but we generally press it after changing anything, including styles, GeoWebCache changes, layer changes, etc.

Thanks in advance

Mark


texts:
I have a leaflet map with a circleMarker layer that has a popup. This works as expected.

dotLayer = new L.geoJson(dots, {
  pointToLayer: function (feature, latlng) {
    return L.circleMarker(latlng, style(feature));
    },
  onEachFeature: onEachDot
  }).addTo(map);

The problem is that I have a tile label layer which should be drawn on top for the labels to be legible.

I can put the tile layer to the markerpane and add the layer after the marker layer to show. But like this the popup will no longer work.

var lDarkMatterLabels = L.tileLayer('https://cartodb-basemaps-{s}.global.ssl.fastly.net/dark_only_labels/{z}/{x}/{y}.png',{pane: 'markerPane'});
map.addLayer(lDarkMatterLabels);

I have a jsfiddle here.


texts: AddFeatures to an ArcGIS Online Feature Service?

texts:
Is it possible to add and delete features from an ArcGIS Online-hosted feature service?

I'm able to do this with services hosted on our ArcGIS Server (example here), but don't see the option under "Supported Operations" when looking at the REST endpoint for an AGOL-hosted feature service (example here).

I understand that few people use both ArcGIS Server and ArcGIS Online, but our client is a heavy AGOL user and would like to store all their data in the platform instead of our ArcGIS Server.

If this isn't available, is it part of the roadmap for ArcGIS Online development? It seems strange that the REST endpoint is different from that made available from ArcGIS Server.


texts: Unable to project GPM precipitation data correctly in QGIS?

texts:
After downloading the GPM (Global Precipitation Measurement) data from FTP, I tried opening the image in QGIS and I get a popup message saying CRS was undefined.

As I know the projection is set to WGS84 I added a South Asia boundary with WGS84 projection over the GPM raster data, and the stack of these two layers did not overlay properly.

The small shape on the top left corner of the snapshot above is the South Asia boundary. I have zoomed in to display the small South Asian boundary in the following snapshot.

My boundary data has the correct projection and I have been using it with other remote sensing data and had never faced this problem before?

Both of the layers have the same WGS84 projection, but still they do not overlay correctly, why is it so?


texts: Labeling points with polygon they overlay using QGIS?

texts:
I have a shapefile of a cut-up map of England. I have overlaid points (lat/long) from a CSV as a separate layer.

Is there a way to output what polygon each point is on. So, i'd want to be able to output that point A,B,C are on Polygon 1, D,E,F is polygon 2 etc.

I've counted how many points per polygon but I do not know how to do to swing it back round so it tells me what polygon each point overlays.


texts: lstrip Method not returning desired result

texts:
I can get the desired result for the lstrip method, I want to return the names of my SDE Connections without the full path, but lstrip doesn't remove all before ArcCatalog. What am I missing? Or is there an easier way?

ie: C:\Users\an043\AppData\Roaming\ESRI\Desktop10.2\ArcCatalog\Transport.sde
I want to return: Transport.sde or Transport

This is the code I have been using:

arcpy.env.workspace = r"C:\Users\an043\AppData\Roaming\ESRI\Desktop10.2\ArcCatalog"
workspacesSDE = arcpy.ListWorkspaces("*", "All")
for workspaceSDE in workspacesSDE:
    print workspaceSDE.lstrip("ArcCatalog")


texts: How To Add CSS Class to SimpleMarkerSymbol in GraphicsLayer

texts:
Can you please let me know if there is a way to bind a CSS class .marker to an ArcGIS SimpleMarkerSymbol in a graphic layer?

I have a SimpleMarkerSymbol called schoolMarker which I am styling it using ESRI's options but I need to add more CSS styles like .marker

var schoolMarker = new SimpleMarkerSymbol();
    schoolMarker.setStyle(SimpleMarkerSymbol.STYLE_CIRCLE);
    schoolMarker.setSize(30);
    schoolMarker.setOutline(null);
    schoolMarker.setColor(new Color([255,116,0,0.25]));

Here is the CSS class:

.marker{
   stroke-linecap: round;
   stroke-linejoin: round;
   stroke-opacity: 1.0;
   fill-opacity: 1.0;
   background-color: rgb(0, 89, 190);
   stroke: rgb(0, 89, 190);
   fill: rgb(0, 89, 190);
   fill-opacity: 1;
   stroke-width: 70pt;
   stroke-opacity: 1;
   opacity: 0.12;
 }


texts: QGIS 2.18.4 show '0' value disappear with labels with rules

texts:
I need help to show '0' value with expression.  When I show only label from a row, '0' appear, but I want to add the expression in the color :

CASE
WHEN  "max_44_mol_pharma"  < 100 THEN  color_rgb( 0,0,0)
WHEN  "max_44_mol_pharma"  >=  100 THEN  color_rgb( 255,0,0)
END

and the '0' disappear.  Can you help me to resolve this issue?


texts: What is the recommended ram for an ArcServer 9.3 box?

texts: I've dug around on the ESRI site for a while and can't find a conclusive answer.  Can someone tell me what ESRI's recommenced RAM is for an server running ArcServer 9.3, with hopefully a link to document it?


texts: How to Clear Parameters from Model Builder - ArcGIS

texts:
I'm attempting to clear out the parameters from my model so that the parameter locations don't show up in the tool when launched from ArcToolbox.

ESRI's resource center has a couple of sentences about clearing parameters, but not how to actually do it.

ESRI's attempt at explaining how to clear out the parameters

Does anyone know the steps to clear the parameters of a model built inside Model Builder for ArcMap?


texts: QGIS Expression Builder: returning pair of values (e.g.: coord)

texts:
Seems like this should be simple, meaning I'm probably "trying to do the wrong thing," but: I have a field (specifically, Label Placement's "Offset X,Y") that expects an input of type double coord. I want to do a data-defined override: pull the X and Y values from columns in the attribute table.

But, how do I return the two values as a single coord pair from the expression builder? Trying [x,y], (x,y), {x,y}, etc. all give a syntax error...


texts: How to configure Geoserver WMS to make it printable from QGIS/ArcGIS?

texts:
I have installed geoserver 2.3.1 - I am using this to create a WMS of some data held in PostGIS.
I can view the WMS in QGIS and ArcGIS, but it won't print from ArcGIS (blank sheet) and in QGIS I get the message:

Some WMS servers (e.g. UMN mapserver) have a limit for the WIDTH and HEIGHT parameter. Printing layers from such servers may exceed this limit. If this is the case, the WMS layer will not be printed

So I assume I have to set something in Geoserver or on the server to allow larger print sizes - where can I do this?


texts: With Modelbuilder, how to iterate field value so final output table is empty if input is empty?

texts:
I have created a model which will output two Excel tables (the green circles at the end) per state. I am using "Iterate Field Values" to go through each state to create the tables. While this gives me each table I need, the problem is that when a table is empty, the model STILL outputs a table for that state - with rows for every state (all input).

Goal:

1) to have NO table output if there is no data in the table

or

2) to output a table that has empty fields

I would prefer the latter option. For Iterate Field Values I tried both checking 'Skip Null Values' and leaving it empty. Neither resolved the issue.


texts: Trying to add a WMS layer to Openlayers 3 but not showing anything

texts:
I have started with OpenLayers 3 recently and I need to add a layer from an own WMS. I have been trying this:

map = new ol.Map({
    target:'map',
    renderer:'canvas',
    layers: layer,
    view: new ol.View({
        projection: 'EPSG:4326',
        center: [4, 40],
        zoom:2
    })
});
var layer = new ol.layer.Tile({
    source: new ol.source.TileWMS({
    preload: Infinity,
    url: 'http://ide.cime.es/menorca/wms/base_referencia/',
    params:{
        'LAYERS':"RE007RUS_municipis", 'TILED':true
    }
})
});

But this is showing nothing, I don't know what to do. With OpenLayers 2 I had the map working, but I need it with OpenLayers 3.


texts:
Is it possible to spatially pre-clip layers (say with a IPolygon object) before they are loaded into ArcMap?

This would be a bit like a query definition but spatial and for raster data.


texts: Shapefile of European Parliament constituencies?

texts: Is there any shapefile (or similar vector data set) of European Parliament electoral constituencies? Some countries (e.g. Germany) are basically one constituency, but others (e.g. UK) are not. Is there any shapefile of all the constituencies in one file?


texts: Can GlobalMapper open Personal Grodatabase (MDB( without ArcGID license?

texts:
I am wondering if GlobalMapper can read an .mdb geodatabase without an ESRI license?  I have encountered some apps that can not read a .gdb without the license but can read an .mdb, so wondering if GlobalMapper is one of them.

I don't have the software to test it out, and if I download it I have an active ESRI license so can't test it out, and have not been able to find any info on the GLobalMapper support pages or user forum.


texts: With a routing model, how do I break up the input locations into start/end pairs?

texts:
I have the following model:

How do I break up the feature set of points (Input Locations) into start/end pairs, and iterate through each pair


texts:
I try to generate a UUID (edit widget) for a vector shapefile in qgis 2 but the field value is always NULL. Does the field requirea specific datatype? I used integer. In the manual I found the info "what uuid-generator is" but not how to use it.

Thanks for any help,
Patrick


texts: Loop through several arcpy.MakeQueryLayer_management?

texts:
I have several make query and feature class to feature class arcpy tools that create the necessary feature classes from my SQL queries. Is there a way to incorporate this into a for loop?

arcpy.MakeQueryLayer_management(input_database=db,
                                out_layer_name ='test1lyr',
                                query = test1,
                                oid_fields = "test1id")

arcpy.FeatureClassToFeatureClass_conversion('test1lyr', scratchGDB, 'test1')

arcpy.MakeQueryLayer_management(input_database=db,
                                out_layer_name ='test2lyr',
                                query = test2,
                                oid_fields = "test2id")

arcpy.FeatureClassToFeatureClass_conversion('test2lyr', scratchGDB, 'test2')

arcpy.MakeQueryLayer_management(input_database=db,
                                out_layer_name ='test3lyr',
                                query = test3,
                                oid_fields = "test3id")

arcpy.FeatureClassToFeatureClass_conversion('test3lyr', scratchGDB, 'test3')

arcpy.MakeQueryLayer_management(input_database=db,
                                out_layer_name ='test4lyr',
                                query = test4,
                                oid_fields = "test4id")

arcpy.FeatureClassToFeatureClass_conversion('test4lyr', scratchGDB, 'test4')

arcpy.MakeQueryLayer_management(input_database=db,
                                out_layer_name ='test5lyr',
                                query = test5,
                                oid_fields = "test5id")

arcpy.FeatureClassToFeatureClass_conversion('test5lyr', scratchGDB, 'test5')


texts: Using GRASS command line tools within QGIS?

texts:
How do I use the GRASS command line tools within QGIS?

I use the QRASS plugin to open the mapset, and then I open the GRASS tools window. I click on 'GRASS shell' under modules and it brings up a bash terminal, which says:

bash: /sw/bin/init.sh: No such file or directory
bash-3.2$

I assume init.sh was supposed to open some GRASS command line program, but it failed to find it. Does anyone have any suggestions for how to procede?

Some background, any advice will be gratefully received: I'm trying to teach myself GRASS using the third edition of "Open Source GIS: A GRASS Approach" by Neteler and Mitasova. Unfortunately it appears that the book is somewhat out of date, designed for GRASS 6-3, while I have GRASS 6-4. This means that some of the GUI stuff doesn't work (particularly d.mon, I believe it has to do with a new Python WxGUI interface) so I can't work through following the examples exactly. I'm hoping that the book will still be mostly relevant in terms of the modules available and the basic GRASS command syntax and use, but I can't follow the instructions that relate to the GUI. For this reason I've decide to ignore the GRASS GUI completely and to just use the QGIS GRASS plug in.


texts: Calculating centroids in decimal degrees using ArcPy without reprojecting vector data?

texts:
I have vector data not in geographical coordinate system (i.e. not in WGS 1984).

How to calculate object's centroids using arcpy without reprojecting the data itself?

I have tried to do this but it doesn't work (centroids are in meters):

arcpy.CalculateField_management(feature, fieldName, "!SHAPE.CENTROID.X!","PYTHON_9.3")
arcpy.CalculateField_management(feature, fieldName, "!SHAPE.CENTROID.Y!","PYTHON_9.3")

arcpy.CalculateField_management(feature, fieldName, "!SHAPE.CENTROID@DECIMALDEGREES!.split()[0]", "PYTHON")
arcpy.CalculateField_management(feature, fieldName, "!SHAPE.CENTROID@DECIMALDEGREES!.split()[1]", "PYTHON")


texts:
The Dissolve tool is not working for me when I use ArcGIS Desktop 10 to dissolve a shapefile.
It says:

ERROR 999999: Error executing function.

It works perfectly when I use version 9.3.
Any help?


texts: Adding a text input field as coded value in domain?

texts: I want to see this in my final table -- I click on a box in a field and i get a list of text values in the dropdown box but the last value should be 'other'. By selecting this 'other' I should be able to enter a text by myself.


texts: Skill requirements and resources to run a private GIS consultancy/ GIS service providing firm?

texts: I'm planning to make a private GIS service providing firm in the coming years. I just wanted to know what skills and resources are required to set up such a firm? I'm aware that the software, hardware and data collecting equipment would cost me. Can you give me some idea of the things i will be needing? For you information, i am a Geography hons student, and presently doing a diploma course in Geo-Informatics. I also have some knowledge in programming.


texts: How to create random points inside polygons in QGIS?

texts: Has anyone created random points within polygons in QGIS 3.0? I created multiple polygons and would like to create multiple random points within each polygon. I am trying to have the number of points within each polygon to be determined by the area of the polygon (I would like it to cover 20% of the polygon). In older versions of QGIS there is the option to select point density and to select the percentage of points with relation to the polygon. However in QGIS 3.0 there is no option to select the percentage of the polygon to cover with points. Is there a way to create these random points within polygons having them cover 20% of the polygon?


texts: Interpolating parts of raster layer in QGIS?

texts:
I have the following raster layer

The part in black have a value of 255, while the yellow parts have other values that goes from 0 to 255 depending on the color. The black part are the one that i want to interpolate since they have a value of 255 because no value were present, so it is like  they are null values.

Is there a way to interpolate/assign a value to the black parts and obtain a raster layer without them?


texts: Measuring basin length in ArcGIS Desktop?

texts:
I am using ArcGIS Desktop 10 to compute Morphometric parameters for 2 catchments. I have already calculated the Linear parameters, (stream order,stream length, Bifurcation ratio) I have also calculated some of the Aerial parameters ( drainage density, drainage texture, stream frequency ) I am currently calculating the Form factor and Enlongation ratio. I have calculated the Area of the basin. The next thing to do is calculate the Basin length. That is where i am stuck.

Is there a way to get it in ArcMap?

Or is basin length same as flow length?

What method or measure can i use to determine the length of the basin?


texts:
I have this table:

And I want to create a field called "COSN2" which after the creation I want to fill it with three characters from the field "COS207".
Than I wrote the following script:

import arcpy, os
from arcpy import env, UpdateCursor
from arcpy.management import AddField

arcpy.env.overwriteOutput=True
arcpy.env.workspace="C:\Users\Thomas\Desktop\Trabalho AE\Trabalho.gdb"
fc="COS2007"

arcpy.AddField_management(fc,"COSN2","TEXT","","","","NULLABLE")
arcpy.CalculateField_management(fc,"COSN2","LEFT([CONS2007],3)","VB")

cursor=arcpy.SearchCursor(fc)
for row in cursor:
    print row.COS2007+" ; "+row.COSN2

Then the eclipse returns me an error:

Traceback (most recent call last):
File "F:\MSIGOT\AE\AnaliseEsp\TrabalhoAE\COSN2.py", line 10, in <module>

arcpy.CalculateField_management(fc,"COSN2","LEFT([CONS2007],3)","VB")
File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\management.py", line 3360, in CalculateField
raise e
arcgisscripting.ExecuteError: ERROR 999999: Error executing function.
Failed to execute (CalculateField).

I've tried to change "VB" (above) to "PYTHON" and "PYTHON_9.3" but the error persists.

Can anyone explain this?


texts: Blurring raster data across tile boundaries?

texts:
I have a raster dataset that is made up of 280 tiles, 51Gb total. I need to smooth/blur the data, retaining the same resolution. I would like to blur across tile boundaries, where possible, so that points on the edge of each raster contain information from the neighbouring raster(s). I would like the output to be in the same format or smaller tiles.

I there an easy way to do this? I see that it's possible to blur using a KernelFilteredSource in VRT, but it's not clear how that can be made to work across tile boundaries.

I guess the other option is to write a script to load each tile, plus the margins from the surrounding 8 tiles, blur the result, trim, and re-save. I'm just wondering if there's a more efficient way to do it before I try that.


texts: I'm looking for a way to continue a line feature in QGIS. I have checked  the Digitizing Tools plugin (both regular and advanced) for a Prolong Line tool, but I haven't had luck. In ArcGIS there's a Continue Feature button when editing but I haven't found one so far in QGIS. I'm looking for a way to do this while editing/digitizing rather than through code. Can this be done?


texts: I want to create composite band image through batch processing using ArcPy

texts:
I have multiple Sentinel 2a raster images in separated folders at same location, which means every image is in separate folder. Now I want to create composite band of each raster. I have written a code and its given an error, its says
File "C:\Python27\ArcGIS10.3\lib\genericpath.py", line 91, in _splitext
sepIndex = p.rfind(sep)
AttributeError: 'list' object has no attribute 'rfind'
How do I specify arguments in composite band?
Here is the code
import arcpy, os
from arcpy.sa import *

arcpy.env.workspace = r"L:\Arcpy\Data\Satellite_image\T44QPJ_20200405"
arcpy.env.overwriteOutput = True

input_raster = arcpy.ListRasters()
outputDir = r"L:\Arcpy\Data\Layer_stack"
print(input_raster)

blue = r"L:\Arcpy\Data\Satellite_image\T44QPJ_20200405_B02_10m.jp2"
green = r"L:\Arcpy\Data\Satellite_image\T44QPJ_20200405T045651_B03_10m.jp2"
red = r"L:\Arcpy\Data\Satellite_image\T44QPJ_20200405T045651_B04_10m.jp2"
nir = r"L:\Arcpy\Data\Satellite_image\T44QPJ_20200405T045651_B08_10m.jp2"
swir_1 = r"L:\Arcpy\Data\Satellite_image\T44QPJ_20200405T045651_B011_20m.jp2"
swir_2 = r"L:\Arcpy\Data\Satellite_image\T44QPJ_20200405T045651_B012_20m.jp2"

# the first part of the file name before extension (if present) then add *.img
output = os.path.join(outputDir,os.path.splitext(input_raster)[0].split('_')[0] + '.img')
arcpy.CompositeBands_management(";".join([blue,green,red,nir,swir_1,swir_2]), output)

I am stuck in CompositeBands_managemnet tool. How can I solve this problem?


texts: Looking for a pythonic way to calculate the length of a WKT linestring

texts:
I was quite unsatisfied with Calculating Length of Linestrings in WGS84 in Miles. It kept me wondering if there is a more convenient, Pythonic way to calculate the length of a WKT linestring according to a given SRID.

I have in mind something like:

srid="WGS84"
line="LINESTRING(3.0 4.0, 3.1 4.1)"
print length(line, srid)

I'm looking for an accurate answer, not sin\cos approximations.

Any ideas?


texts:
I would like to project my shapefile information on the maps. I converted my shapefile into a JSON file. How can use the Json file to create a featurecollection object in the ARCGIS Javascript API. The Json file contains 50 thousand features. Have given a sample set of my feature collection below. I do not want to embed the entire JSON file within the javascript file.

{ "type": "FeatureCollection", "features": [ { "type": "Feature", "properties": { "ROADNAME": "AAA", "LINK_INDEX": 1}, "geometry": { "type": "LineString", "coordinates": [ [ 163.554342, -77.702882 ], [
163.553160, -77.702327 ] ] } } , { "type": "Feature", "properties": { "ROADNAME": "KIAN TECK WAY","LINK_INDEX": 2}, "geometry": { "type": "LineString", "coordinates": [ [ 162.527145, -77.686591 ], [
162.522273, -77.686598 ] ] } } }


texts: Constructing SQL statement for ArcPy SearchCursor

texts:
I'm struggling with a real simple cursor problem, and I figure it's probably a syntax issue. But I cannot see what it is.
tmpHeader = ['OID@', 'SHAPE@', 'ProjectID','Status','ProjectName','Description','ProjectStartDate','ProjectEndDate']

for row in arcpy.da.SearchCursor( "SourceFeatures", str(tmpHeaderHolder)):
     ...
     ...

And then I get this weird error...What puzzles me is the syntax of the resulting SQL statement in the error message
Traceback (most recent call last):
  File "loader.py", line 77, in <module>
    for row in arcpy.da.SearchCursor( "SourceFeatures", str(tmpHeader)):
RuntimeError: An invalid SQL statement was used. [SELECT ['SHAPE@', 'ProjectID', 'Status', 'ProjectName', 'Description', 'ProjectStartDate', 'ProjectEndDate'], OBJECTID FROM MyFeatureClass]


texts: What does it mean by scale in reduceRegion function in Google Earth Engine

texts:
I want to calculate the mean over an area in Google Earth Engine. I am unable to exactly understand the meaning of scale in the function reduceRegion (https://developers.google.com/earth-engine/apidocs/ee-image-reduceregion) I checked https://developers.google.com/earth-engine/guides/scale but am still not clear. In other examples on stack exchange people have shown examples with landsat and have shown them with scale as 30. As I am trying to understand this, it should be equal to the spatial resolution of the data. And as the landsat data is at 30m spatial resolution, the scale being used is 30?
I am specifically trying to use the data https://developers.google.com/earth-engine/datasets/catalog/MODIS_006_MOD09GA
Is the following code right
region = ee.Geometry.Polygon([

  [[74, 14], [87, 14],[87, 28], [74, 28], [74, 14]]

  ]);
collection_ndvi = ee.ImageCollection('MODIS/MYD09GA_006_NDVI').select('NDVI').filterDate(dates,datee ).mean()
 data = collection_ndvi.reduceRegion(ee.Reducer.mean(), monsoon_region, scale=500).get('NDVI')

I don't know if I am doing this right by using a scale of 500. But the page https://developers.google.com/earth-engine/apidocs/ee-image-reduceregion says that
scale   Float, default: null
A nominal scale in meters of the projection to work in.

and the data (https://developers.google.com/earth-engine/datasets/catalog/MODIS_006_MOD09GA) is at 500m


texts: GeoServer error on styling line width relative to zoom level

texts:
I've applied MBstyle (gs-mbstyle-2.13.2, gt-mbstyle-19.3) extension for GeoServer (2.13.2) to style layers in mvt format. It works fine for simple styling but the problem occurs when I want to set the line-width of features based on the zoom level of the map display. I've used the following snippet but Goserver throws me this error:

  org.geotools.mbstyle.parse.MBFormatException: "line-width" number from
  JSONArray not supported

{
  "version": 8,
  "name": "cyan line",
  "layers": [
    {
      "id": "road",
      "type": "line",
      "source-layer":"road",
      "paint": {
      "line-color": "#0099cc",
     "line-width": [
         "interpolate",
         ['exponential', 2],
         ['zoom'],
         5,1,
         20,8
     ]
    }
    }
  ]
}

Is it something wrong with MBstyle version or something else causes the problem?


texts: Why does a raster in Gauß-Krüger has a little offset to a shape in UTM?

texts:
I have raster data (TIFF files) where I think it's spatial reference is in Gauß-Krüger ("DHDN_3_Degree_Gauss_Zone_3") and some shape files in UTM ("ETRS_1989_UTM_Zone_32N_8stellen"). As you can see in the screenshot the raster data (black drawings on top) has a slight offset to the UTM shape (light blue in background) of ~200 meters to north and ~90 meters east.

Example coordinates from one of the tfw files:

0.26051843980496
0.00000000000000
0.00000000000000
-0.26051843980496
3648545.83601748060000
6038270.01816578210000

I am wondering, what the reason is. I think of two possibleties for the raster:

I suggested a wrong spatial reference. But how could I find the right one?
The data itself is inaccurate. Then I need a way to define an offset for this raster. I tried adding the offset to the FalseEasting and FalseNorthing in IProjectedCoordinateSystem5, but it didn't help.

Does anyone have a clue?

Update/Solution

Thanks to Mintx I found the solution. The data itself is accurate but needs to be transformed to fit to the UTM project. Here is how I set the on the fly transformation on the map:

IMapGeographicTransformations transf = _map as IMapGeographicTransformations;
if (transf != null && transf.GeographicTransformations != null)
{
  ISpatialReferenceFactory2 srFactory = new SpatialReferenceEnvironmentClass();
  IGeoTransformation geoTransformation = (IGeoTransformation)srFactory.CreateGeoTransformation((int)esriSRGeoTransformationType.esriSRGeoTransformation_DHDN_To_ETRS_1989_8_NTv2);

  transf.GeographicTransformations.Set(esriTransformDirection.esriTransformForward, geoTransformation);
  transf.GeographicTransformations.Set(esriTransformDirection.esriTransformReverse, geoTransformation);
}


texts: Concatenate and sequentially order unique groups of lines using SQL

texts:
I have a segmented line (road) dataset that potentially runs through multiple area polygons (localities). I want the road dataset that is grouped by the unique field (RoadMntnc) to have a concatenated list of the localities it passes through in the direction of flow of the road.

The following partially works

select RoadMntnc, lower(Locality) from clip_DSC_Roads group by RoadMntnc, lower(Locality) order by RoadMntnc

But it lists every locality each segment is in and doesn't follow the direction of travel of the road.

I need the result as

Miallo_Bamboo_Creek_Rd | Miallo to Whyanbeel to Bamboo
Kingston_Rd | Whyanbeel

So if it is multiple localities (pink below) it should be in the direction that the road is.

To clarify a bit -there is a field with locality that changes in the dataset when a segment of road is inside different localities. It doesn't come from an external dataset -however if it has to then this won't be a problem either. This is why there is no 'intersect' in the SQL.

A QGIS or Mapinfo sql process would be fine as well.

UPDATE

Based on Concatenate or group by values using spatial join in QGIS?

If I import the data into Spatialite and run

select RoadMntnc, group_concat(Locality) from clip_DSC_Roads group by RoadMntnc

I get

But how do get just the unique values and also check if it's in the right order?

I have tried

select RoadMntnc, distinct group_concat(Locality) from clip_DSC_Roads group by RoadMntnc

but get a syntax error.


texts: Does a feature need to be in the Spatial Index (QGIS) in order for NearestNeighbor to find its neighbours?

texts:
To be more specific: I have a feature class for which some features have an entry for a variable ('KnownBj') and some not (have a 9999 entry). I want for each feature to find the 5 nearest neighbors that have an entry in 'KnownBj'.
What I currently do is to prepare two dictionaries - A has all the features and B only the features with a 'KnownBj' .I add B to a SpatialIndex, and loop through A searching for neighbours with .NearestNeighbor:

canvas = qgis.utils.iface.mapCanvas()
layer = canvas.currentLayer()
iter = layer.selectedFeatures()
feature_dict = {f.id(): f for f in iter}
featDictKnownBj = {k: v for k, v in feature_dict.iteritems() if v[KnownBj] != '9999' }

#Prepare Spatial Index
SpatIndex = QgsSpatialIndex()
map(SpatIndex.insertFeature, featDictKnownBj.values())

#SEARCH NEIGHBOURS
NNDict = {}
for f in feature_dict.values():
    NNs = SpatIndex.nearestNeighbor(f.geometry().centroid().asPoint(),5)

Now this seems to work, my question is - am I missing sth? Should I add the current feature to the Spatial Index, then search, then remove it and so for all, or should this be perfectly ok?


texts: My office does a lot of field work where we log points via the Collector app. Those points are then always exported to a geodatabase feature class for further editing and analysis. If edits and changes are made to the feature class, we then have to go back and re-publish the file to the server so that everything matches and the file is up to date in both places. Is there a way to just link the two so that when changes are made to one, changes are automatically made to the other?


texts: Difference between "heap size" and "SOC"

texts:
About ArcGIS Server Administrator Directory and concerning the Server Machine Properties, what are the differences between :

App server maximum heap size (in MB)
SOC maximum heap size (in MB)

Current Version: 10.5.0

Platform: Windows Server 2012 R2-amd64-6.3


texts:
I am looking to run

https://polar.ncep.noaa.gov/waves/examples/usingpython.shtml#example_2

example on wave model data.
It refers to downloading grib files from

https://nomads.ncep.noaa.gov/pub/data/nccf/com/wave/prod/

When I browse into a current date subfolder, what is the naming convention?

I see

multi_1.ak_10m.t00z.f000.grib2

I assume it must be something like

multi_1.<LOCATION>_??m.t??z.f<FORECAST_IN_HOURS>.grib2

LOCATION I see ak - Alaska? , ep - East Pacific? , glo - Global?

What are the others?


texts:
I know PostGIS Topology supports splitting road edges while creating a network.

I was wondering if PostgIS topology can do the opposite and merge two lines?

Does it support merging two consecutive road segments that touch each other at a common (degree-2) node?


texts: How to load a raster into the Map Canvas using QGIS Action

texts:
Is it possible to write an Action in QGIS that will load a raster into the current map canvas based on a url/path of the raster as specified in the attribute table of a shapefile?

The shapefile is a index-tile.


texts: Cloud masking of Sentinel-2 product using SNAP

texts:
I am using Sentinel-2 images in tropical region where almost all image contains clouds which needs to be masked out.

What is the best way to remove the clouds from the image in all the bands, is SNAP is the right tool?
After masking out the clouds how do I interpolate the empty cells of the raster?


texts: Flipping coordinates in projection results in distortion. A simple mistake

texts:
Edit: I made an embarrassing rookie mistake late at night (excuses, excuses) and flipped my lats and longs. As a result, people searching for an answer to my original question will not be helped by this. To resolve this, I've changed the title to this post to reflect the mistake caught in the accepted answer.

I'm attempting to project some data to a planar coordinate system and the results look very wrong. Here's a reproducible example with a map of New England:

library(rnaturalearth)
library(sf)
library(dplyr)
library(ggplot2)

states_sf <- ne_states(country = "united states of america", returnclass = "sf")
states_names <- c("Maine", "New Hampshire", "Vermont", "Connecticut",
                  "Rhode Island", "Massachusetts")
states_sf <- filter(states_sf, name %in% states_names)

ggplot(states_sf) +
  geom_sf() +
  theme_bw()

But when I project and replot the data, it doesn't look as expected.

# Lambert equal-area projection
laea_centered <- "+proj=laea +lat_0=-70.30744 +lon_0=43.15268"
coord_shift <- "+x_0=12019341.4 +y_0=1730136"
proj_ref <- " +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
laea_proj4 <- paste(laea_centered, coord_shift, proj_ref, sep = " ")

states_projected_sf <- st_transform(states_sf, crs = laea_proj4)

ggplot(states_projected_sf) +
  geom_sf() +
  coord_sf(datum = st_crs(states_projected_sf)) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

Just for grins and to further illustrate how absurdly I've managed to mangle this, here's the projected map with geographic coordinates and graticule.

ggplot(states_projected_sf) +
  geom_sf() +
  coord_sf() +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

Can someone explain where I went wrong and how to do this better?


texts: Filter on a array of value in OpenLayers

texts:
I'm running into an issue.

I have an array of value coming from a JSON POST.
After retrieving these values and storing them into a var, I need to filter an OpenLayers layer with this values.

I'm using this :

mapFilterbyStation = new OpenLayers.Filter.Logical({
        type: OpenLayers.Filter.Logical.AND,
        filters: [
            new OpenLayers.Filter.Comparison({
                type: OpenLayers.Filter.Comparison.EQUAL_TO,
                property: "id4geoserv",
                value: stationsParse
            })
        ]
    });

where stationsParse is an array of values.

Problem, is that the value to filter can only contains : number or string, according to the OpenLayers docs : http://dev.openlayers.org/releases/OpenLayers-2.10/doc/apidocs/files/OpenLayers/Filter/Comparison-js.html

Is there a way to achieve what I want using filters ?

I was thinking about iterate on each value stored in my array, and create a filter for each of them, but this solution looks ugly for me.

I'm surprised there is nothing to filter on multiple values from the same parameters in OpenLayers.

So if anyone has a tips ;) !


texts: Labeling on Symbology in Table of Contents of ArcMap?

texts:
I have to recreate a map from one I have found on the internet as close as possible for an assignment and can anyone tell me if there is some way to change the labeling in the TofC?

I'm doing a map of Africa and the climate. I found one I like online but the climate is broken down if different labels.

For instance the online map shows under the heading of Climate: Tropical Wet, Tropical wet & Dry, Arid, Semiarid... but the map layer in AcMap is showing Equatorial, Intrazonal, Subequatoral and so on.

I don't know how to add my screen shot here sorry. But what I wanted to do was leave the layer name the same: Climate but to change the sub names under Climate to match the map that I had found online. I can change the colors so they are almost alike but need to have the sub name headings the same


texts: How to create a consecutive alphanumeric column in a QGIS attribute table?

texts: I am trying to create a column in the attribute table such that each cell lable is a prefix (M) follwed by a unique number ie M1, M2,.. etc


texts: QGIS distance matrix in meter instead of decimal degrees

texts:
I'd be thankful for your help.
i created a point layers of  store locations and customers in Florida by given coordinates (WGS 84 EPSG: 4326). Now i need a distance matrix between stores and customers but the result is in decimal degrees i think. I understand i have to reproject the 2 layers so that i can work with meter but how dows this work in QGIS? Which projection should i use and which setting should i use?


texts: Error Creating Project in QGISRed Plugin

texts:
I am trying to create a new project in QGISRed Plugins using QGIS 3.10 A Coruna, and it gives me an error message:

CRITICAL    Error : The specified directory service attribute or value
does not exist.

I tried creating a new project in QGISRed before and it worked normally. I googled some articles about this issue. Maybe it needs read/list permission in the directory; but I have no clue about QGIS.
I am a bit of a beginner to all of this.


texts: How to display the current polygon area when editing in QGIS 3?

texts:
I use QGIS 3.10 to create new polygons manually.

I want to see the area of an object interactively during its creation.

I know that after creating a polygon, you can automatically calculate its area and display this value in the attribute table, but this option does not suit me.

How to do that?


texts: Setting options for individual layers within ArcGIS Server Map Service?

texts:
I've been working on an internet mapping site to include Esri basemap layers, and in-house map services broken out by layer.

I am also still fairly new to javascript so this might be a javascript question.

Here is the code for the mapservice and the layers I am accessing.

var patLayer = new esri.layers.ArcGISDynamicMapServiceLayer("http://gis.co.frederick.va.us/ArcGIS/rest/services/Planning/Planning_PAT/MapServer",{id:'PAT',opacity:.6});
    patLayer.setVisibleLayers([7,6,1,2,3,4]);
    legendLayers.push({layer:patLayer,title:'Planning Layers'});

What I am not sure of is how to set the options only for specific layers and not just the entire map service e.g.opacity?


texts: Importing .CSV into QGIS 3.10 missing some data points

texts:
I am trying to display a delimited text layer from a .csv (utf-8) in QGIS 3.10.  I have done the same thing successfully in ArcMap 10.7.1 with no problem, but when I add the delimited text layer in QGIS over 700 of the 1900+ data points do not show up.  I tried to resolve this by creating a shapefile in ArcMap and bringing that into QGIS, but exact same problem occurred.  When I open the attribute table in QGIS the data is all there, it simply isn't being displayed. I am posting some images.  The first on is the result in Arc, then Q and then I selected the missing data in Arc so you can get a glimpse of what the data look like.

The data were collected using WGS_1984 and the Maine/NH towns is projected in NAD83_StatePlane for NH in Feet. Any help solving this problem would be greatly appreciated.  Is it a CRS issue?  Are the data in the .csv messed up?  If so, why does it work in Arc and not Q?
Here is the text of the data the first line of data show up fine, the second does not show up.

id,username,lat,lng,hotspot,observed,country,county,state,municipality,postcode,observation_notes,Sample Number,Sample Location,Was the sample filtered?,Type of Filtration System,Water filtration description,Beryllium (Be) (ug/L),Chromium (Cr) (ug/L),Manganese (Mn) (ug/L),Iron (Fe) (ug/L),Nickel (Ni) (ug/L),Copper (Cu) (ug/L),Arsenic (As) (ug/L),Selenium (Se) (ug/L),Cadmium (Cd) (ug/L),Antimony (Sb) (ug/L),Barium (Ba) (ug/L),Thallium (Tl) (ug/L),Lead (Pb) (ug/L),Uranium (U) (ug/L),Lab data entered on,SEPA School Name,City,State,ZIP code,Well Type,Previous arsenic test,Year,Do you treat your water?,Image 1
103070,anonymous,43.181041,-71.999461,,2019-12-09 23:35:00+0000,United States of America,Merrimack County,New Hampshire,Bradford,3221,,2019-1259,Kitchen,No,No filter,,0,0,1.41,9.3,0.64,487.15,0.06,0.11,0.07,0.08,19.55,0,5.73,0.01,1/23/2020,Kearsarge High School,Bradford,New Hampshire,3221,I don't know,No,2019,,
102564,anonymous,44.821972,-67.20856,,2019-11-23 21:00:00+0000,United States of America,Washington County,Maine,Whiting,4691,,2019-815,Other,No,No filter,,0,0.3,1.64,8.13,0.03,0.18,0.67,0.05,0,0.03,3.11,0,0.02,0.04,12/20/2019,Machias Memorial High School,Whiting,Maine,4691,Dug well,No,2019,,

Not sure if this is the best way to share this.  Thank you all for the help.


texts:
I wanted to find out mean of the selected points using QGIS. I tried out selecting using attribute table but could not find an option for mean in field calculator.

I have attached a screenshot with this.


texts: How to connect nodes via shortest path along network lines with GRASS v.net modules?

texts:
My goal:

I would like to connect several nodes via the shortest path along a network --> The red nodes in the next image should be connected along the shortest path of the yellow lines.

My approach:

I have tried different v.net modules, such as v.net.steiner directly on the nodes and the network file --> This doesn't seem like the shortest way along the lines to connect the nodes.

My tried solutions:

v.clean (snap) to close possible small gaps in the topology
add length per line segment via field calculator and use this as a 'weight' in v.net modules
use other v.net modules (e.g., v.allpairs, v.net.path)

My impression:

I think all modules consider the shortest path from line intersection to line intersection in the general direction to the next node. As a consequence, nodes are connected via short line segments that might not result in the overall shortest path (see 3rd image).

My question:

Which module/approach would connect the nodes via the overall shortest path?


texts: How to identify parameters sent from ArcMap to geoprocessing service on ArcGIS Server?

texts:
I have a question regarding the use of ESRI geoprocessing services and maybe some of you can give me a hint for the following "problem"...

I have set up a geoprocessing service which should solve a Vehicle Routing Problem on a given network. If I call the service via Rest URL with the expected parameters the service does not calculate a route, but does also not print any error messages. The service itself works correct. I have tested the functionality by calling the service via ArcMap by entering the expected parameters via the GUI. The parameters are send to the service on the server, the results are calculated there and finally visualized in ArcMap. Thus I think the service works correct.

Now the question: Is there any software that can show me how the request is build up / how the parameters are send from ArcMap to the server? My idea is to copy the syntax and to use it in the Rest URL. I have already tried Charles (http://www.charlesproxy.com/) but it doesn´t show any paramters. Another idea was to look at the log files on the server, but also here I cannot identify any parameters.


texts: Minimum number of samples for kriging interpolation

texts:
I get some data with number of specimens with a request to interpolate it using kriging method.
After some investigation it appeared that kriging results (performed in ArcGIS Geostatistical Analyst with default parameters) are not satisfactory. Interpolated values are far different from measurements (especially top ones) and the surface doesn't look reliable. Here's the picture:

I suppose that the main problem is insufficient number of samples.

How many points should we use to get reliable results?
Or maybe kriging method is not appropriate for such diversed values?


texts: Tool for smoothing values based on neighboor analysis in QGIS?

texts:
I have a raster with areas of different values. Lets say area 1 has values that range from 1 to 10 while area 2 has values that range from 30 to 40. Now I also have one value of 35 in area 1.

Is there a tool that can detect those irregular values and revalue them all based on a neighboor analysis?


texts:
I haven't found a way of setting the selection method (Create new, Add to current, Remove from current, Select from current selection) in arcpy. Is it possible to do this and - if yes - how?

I am switching between these methods regularly and having to press a modifier key is somewhat error-prone, so I'd like to create a python add-in toolbar that allows me to select the method directly.

Also, I want to select interactively because the features are more or less picked by gut feeling.


texts:
I was looking in Fiona to get for each feature its extent but I didn't find how.

I have try to far to do something like below

import fiona

with fiona.open('countries/ne_10m_admin_0_countries.shp', 'r') as source:
    for f in source:
        geom = f['geometry']
        print geom

I was expecting to get a method for this at the f variable level. After some documentation reading, I've seen that f is a pure python record. So at the end, how with Fiona can I get the extent (or bounding box) of each feature geometry?

PS: I already know the pure GDAL/OGR python solution so I expect a Fiona solution please


texts: Inverse of Euclidean distance in IDRISI?

texts:
I am trying to do the reverse / inverse of a euclidean distance in IDRISI (TerrSet). For example, I am doing distance to roads analysis. The output from the distance raster shows roads as 0 and areas closest to roads as 0.1, 0.2, etc and the areas furthest from roads as a very high value. However, I would like the raster to be the reverse of this - the areas closest to the roads should have the highest values and the areas furthest from the roads should have the lowest values.

How do I create a raster that reverses the number in IDRISI?

I would use reclass but this would be difficult to use as the values are continuous and not integers.

Is there a way to do this in reclass or raster calculator that I am not aware of for continuous values without overgeneralizing them?

Or is there a different tool that I can use?


texts: difference between st_accum and st_collect

texts:
Can anyone please tell me the difference between ST_Accum and ST_Collect of PostGIS.
https://postgis.net/docs/ST_Accum.html

https://postgis.net/docs/ST_Collect.html


texts: QGIS point label expression - using a value from another layer

texts:
I have a point data with depth values in CSV format as well as a polygon layer with a 'Critical_Depth' attribute. I'm trying to use the label rule expression to change the colour of the text depending on whether the point 'Depth' value is above or below the 'Critical_Depth' for that area.
I can achieve this by saving the CSV as a geopackage and adding a field 'Critical_Depth' that is calculated using a spatial join, but i'm wondering if i can simplify the process.
Above - Depths CSV data
Above - Polygons with Critical_Depth value
Above - Raw xyz data
Above - What I'm trying to achieve


texts:
Here is the next function to select and modify the selected polygons. I have an error in var "vectorLayer = new ol.layer.Vector": missing: after property id. that someone are correcting my code and giving me a solution of how polygons can be modified

Function ModifySnappingZone(MapIdentifier){

  var selectf = new ol.interaction.Select({
      var vectorLayer = new ol.layer.Vector(

    layers:[vectorLayer],
     });
   map.addInteraction(selectf);

  selectf.on('select', function(evt) {

//here I want to modify the polygone (deform)
});
}

//function draw(works already):
function DrawSnappingZone(MapIdentifier){
        var coordinatearray = [];
        var snappingpoints = prompt("Wieviele Punkte sollen gesetzt werden?");
        MapIdentifier.on('click', function(evt)
        {
            var single_coordinate = ol.proj.toLonLat(evt.coordinate);
            coordinatearray.push(single_coordinate);
            //console.log(coordinatearray);
            console.log(coordinatearray.length);
            //DrawSnappingZone(coordinatearray, 5);
            if (coordinatearray.length == snappingpoints)
            {
                coordinatearray.push(coordinatearray[0]);
                RouteAnzeigenAusKoordinaten(coordinatearray, KarteGlobal_o, 'red');
            }
        });


texts: Starting a QGIS project for a fantasy world map from raster images

texts:
I want to use QGIS to georeference two fantasy (D&D) raster maps, one is the map for the entire world (Toril), one is the map for one of the continents (Faerun).

This is the post with the world map.

This is the website with Faerun map.

I have already aligned the continent map to the world map in gimp, and marked the corners on the world map where the continent map should go (in order to be able to easily find them on the georeferencer plugin). I have also recentered the whole map to have the equator and the prime meridian at the center of the image (I don't know if this was needed).

This is the result.

Now I would like to import both of them in a QGIS project and (sorry for improper terminology) place them on a sphere/world.

I need to specify a projection, but I have no idea what kind of projection I am looking at in the world map, and I also think I need to create a .pngw file to tell the georeferencer plugin about the measures in this map, but it's not clear how should I calculate the needed values. For the moment we can assume that Toril has the same size of Earth.

I have already seen this question Get correct measurements for fantasy world map but I am starting with only two raster images with no georeference data of any kind.

How should I proceed?


texts: select files in folder by geographic location

texts:
In a single folder I have a large set of aereal images (.tif), all georeferenced in EPSG:5677. They cover a "large" total area.

I need to load a small set of them, covering only a fraction of the total area in QGIS.

As the filenames are just unsorted numbers, I wonder wheather there is an option to select files by geographic location while they are not loaded in a project. I work in QGIS, a R-solution would also help.


texts:
I have GeoJSON files that contains altitude values. I have used these GeoJSON files to be shown on leaflet.js web map.
I need to show values of highest and lowest point in two divs when I clicking on each line.
Here is my code that to show GeoJSON lines on the map:

$.getJSON('data/routes/New/Lwn.geojson', function (collection) {
    Lwn = L.geoJson(collection, {
      'style': function () {
        return {
          'color': '#0000FF',
          'weight': 3,
          'opacity': 0
        }
      }, onEachFeature: function(feature, layer) {
        layer.on('click', function(e) {
          showElevation(e.target);
        })
        layer.on('mousedown', function() {
          layer.setStyle({color: '#FFFF00', opacity: 1});
        })
        layer.on('mouseup', function() {
          Lwn.resetStyle(layer);
        })
      }
    }).on('click', function (e) {
      $('#overlay').slideDown();
      $('#Name').html("Name: "+''+e.layer.feature.properties.name+'');
      $('#Distance').html("Distance: "+''+e.layer.feature.properties.Lwn_length+''+" Meters");
      $('#Surface').html("Surface: "+''+e.layer.feature.properties.surface+'');
      $('#Route-Symbol').show();
      document.getElementById('Route-Symbol').src = `https://hiking.waymarkedtrails.org/api/symbols?osmc:symbol=${e.layer.feature.properties.imgSource}`
    });

});


texts: QGIS 3.14 Python error when executing Clip raster by extent tool

texts:
I am not proficient in Python, therefore I am not able to parse the error message or understand it. I am trying to clip a DEM (raster, .asc format) by the extent of town boundaries (vector). Fortunately, I can do this with FME without a problem.
Input parameters:
{ 'DATA_TYPE' : 0, 'EXTRA' : '', 'INPUT' : 'W:/xxx', 'NODATA' : None, 'OPTIONS' : '', 'OUTPUT' : 'TEMPORARY_OUTPUT', 'PROJWIN' : '977760.940000000,986582.310000000,6262610.720000000,6273061.720000000 [EPSG:2154]' }

GDAL command:
gdal_translate -projwin 977760.94 6273061.72 986582.31 6262610.72 -of GTiff W:/xxx/Dpt_83_asc.asc C:/Users/xxx/AppData/Local/Temp/processing_gHDFAW/c7c4be4bb5c34603a7ed682f43cbae6c/OUTPUT.tif
GDAL command output:
Traceback (most recent call last):
File "C:/PROGRA~1/QGIS 3.14/apps/qgis/./python/plugins\processing\algs\gdal\GdalAlgorithm.py", line 138, in processAlgorithm
GdalUtils.runGdal(commands, feedback)
File "C:/PROGRA~1/QGIS 3.14/apps/qgis/./python/plugins\processing\algs\gdal\GdalUtils.py", line 108, in runGdal
for line in proc.stdout:
File "C:\PROGRA~1\QGIS 3.14\apps\Python37\lib\codecs.py", line 322, in decode
(result, consumed) = self._buffer_decode(data, self.errors, final)
UnicodeDecodeError: 'utf-8' codec can't decode byte 0x82 in position 68: invalid start byte

Execution failed after 0.81 seconds


texts: Define line height for features in Autocad Map

texts:
What is the procedure to assign a specific lineheight to a feature in Autocad Map. I understood the style definition but I can only assign a thickness between 0 and 5 instead I want a specific height (1mm).

Thanks a lot


texts: Auto adjustment of table size automatically in ATLAS from QGIS

texts: I would like to know if anyone knows any way to configure my ATLAS to automatically resize the attribute table if, for example, one field is having a long attribute than usual. I know in the layout of ArcGIS is possible but i am not finding the way in QGIS.


texts: How to make a polygon on QGIS based on a PDF image?

texts: I'm trying to take an image out of a PDF and make a shapefile on QGIS. How can I do that?


texts: Getting max feature value with GDAL Python RasterizeLayer

texts:
I'd like to know how to use the Python/GDAL function RasterizeLayer to create a raster layer, where each pixel give the maximum value of a integer attribute. According to docs and forums, if I use :

gdal.RasterizeLayer(out_raster_ds,  # output to our new dataset
                             [1],  # output to our new dataset's first band
                             layer,  # rasterize this layer
                             None, None,  # don't worry about transformations since we're in same projection
                             [0],  # burn value 0
                             ['ALL_TOUCHED=TRUE',  # rasterize all pixels touched by polygons
                              'ATTRIBUTE=id']  # put raster values according to the 'id' field values
                             )

This will use the value of "id". But in my case, each raster cell intersects many spatial data (from PostGIS), and I want to keep the highest value.

Any idea ?


texts:
I am using SAGA --> Grid Calculus --> Raster Calculator on QGIS version 2.6.  The formula I am using is:

ifelse(gt(a,0), (a), b)

where a is my layer with values in a greater extent than that of b.  The resulting layer doesn't seem to abide by the else condition of my formula contains only 0s.

Dug into error log and found a bunch of these type errors:

  ERROR 1: Can't load requested DLL: C:\PROGRA~1\QGISBR~1\bin\gdalplugins\gdal_ECW_x


texts: Changing lakebed symbology to relative depth in ArcGIS 10.5 (TIN or raster file)

texts:
I have several TIN files representing empty lakebeds at altitudes ranging from 5000m to 3000m. I have calculated the TINs using ArcGIS hydrology tools and an ASTER v3 DEM projected to WGS_1984_UTM_Zone_43N.
The DEM data as well as a folder containing the TIN data for several of the lakebed depressions can be found here: https://drive.google.com/open?id=1ERFdsqDGLH1a_FbxwawE_gPm0Au0Q9vT

I want to present the depth of the lakebeds in an appealing way.
I can change the symbology type to elevation but then the same color ramp is used for every lakebed, regardless if it is 10m or 100m deep - the deepest part in every lake has the same color. Therefore, at first they all seem to have the same depth:

I would rather have the parts of all lakes with a depth of 5m, 10m, 15m etc. to have the same color coding so that I would be able to tell the depth difference between the individual lakes at a glance. Is there a possibility to tweak the appearance of a TIN to display not the absolute altitude but its individual depth?


texts: Can't add a new field to an attribute table

texts:
I haven't been able to find a question that matches my situation. I have a .csv file that I made, and I loaded it into arc to join it to another layer. Due to another issue for another question, I have to add a field. Problem is the add field button is greyed out, and I'm not in an edit session.

I can't just add a new field into the .csv, there's another issue blocking me from doing that. Is there a reason why I can't add a new field?

I'm pretty new to using arc, and I haven't been able to find any answer to this online.


texts: Georeferencing vs Georegistration of an Image, Vocabulary Clarification

texts:
Question
This is a question about the meanings and possible differences between georeferencing and georegistration. Language can evolve over time so I will ask this question in two parts.

Is there technically or originally a difference between these terms?
If they are technically different, do you find those in the field usually conflate the two (as in a shift in meaning, evolution of language)?

My Understanding
Say you have an aerial photo (below left) with no spatial reference and a basemap (below right) with roads and other features that does have a spatial reference.

If you had coordinates for road intersections or other known control points, you could click on the aerial photo, enter the coordinates to say this pixel is at these coordinates and that would be georeferencing. However, more common it seems is to click on a feature in aerial photo, and then click on the same feature in a basemap that already has a spatial reference. In this case I would call this georegistration or just registration, but it seems like most others I talk to would call this georeferencing.
I know I am being terribly pedantic, using one term more generally is fine, but for my own sanity is there technically a difference?


texts: Joining attribute data from multiple polylines using ArcGIS for Desktop?

texts:
I've got two polyline shapefiles; one showing the location of cross sections taken across a river and the other containing the river centre line that has been split by the cross section layer to create multiple features. Each cross section contains a unique ID along with some additional information in the attribute table (most of which isn't too important - mainly the ID). Each river centre line feature also contains a unique ID.

I'd like to spatially join the river centre line features to both the upstream and downstream cross section features - the result being that each centre line feature now contains information on the upstream and downstream cross section IDs,  The output shapefile needs to contain a single feature for each river channel with the connected upstream and downstream section names in the attributes.

Using the detailed MS Paint example image below, channel 1 should show it is connected to A and B, channel 2 connected to B and C, channel 3 connected to C and D etc.

I have access to ArcMap 10.2.

I had previously attempted to use the spatial join tool as suggested by @mmoore, however this tool creates a duplicate of each river channel polyline in the output. Each feature then contains either the upstream or downstream section name. Is there a way to either produce the required output using the join tool, or maybe by dissolving the output somehow?


texts: How to find all features in one table that intersect features in another table in PostGIS?

texts:
I have two tables, boundary polygon table and location point table. I need to find all boundary polygons that contain at least 1 location point. What is the proper way to do this in PostGIS?

I tried:

SELECT *
FROM boundary, location
WHERE ST_Intersects(boundary.geom, location.geom)

But this returns each combination of boundary and location that intersect, when all I want is the boundary returned once. For example, I got 76 results, but only have 10 boundary entries. Do you know a better way to do this?


texts:
I want to create a 3d map and have found the osm2world converter.
But i do not understand whether I can use it to convert the whole world map from 2d to 3d, for example on 17th and 18th zoom levels, so that when the zoom is below 17 or 18 I see a 2d map, but it switches to 3d at zoom level 17.  Alternatively it could switch to a 3d map using a button.

I am using osm, leaflet, extjs4 on the client.  On the server I am using pyramid framework.

Can anybody help me find out how to create such a mechanism?


texts:
I have imported base and rover data to Leica Geo Office and did a selection to the reference and rover station in the GPS-Processing Tab.

Why are the Process button and Data Processing parameters still inactive?


texts: Check Personal Geodatabase attribute against external MySQL database?

texts:
I have a Personal Geodatabase to which editors enter values from some field inputs.

I would like some of the entries, eg Project No ( whose name is unique and has some some related attributes like contractor etc ) to be checked against an MySQL database and pop up properties of that Project number from External DB to be available for comparison during data entry..

This would help to avoid wrong Project name entry , double check whether its the Same project which is currently being entered, based on related properties.

I had a conversation with the Developer of the External Web application from where we want to fetch the data - Its a PHP application having a MySQL database at backend. So basically I would need to have a secure Login from the ArcGIS editing Window/Session.

Does this warrant a different approach?


texts: Is it possible to import movie into QGIS?

texts: I would like to 'geotag' a movie in a QGIS layout. Is there a plugin that can accomplish this or have I missed something in the new release of QGIS?


texts: Invoke FME batch process/Process using c#

texts:
I want to invoke FME batch process and general FME processes using c#, I have no idea about FME, can anyone please let me know if anything like this is possible?or is it possible in python, any help on this is most welcome.
I'm thinking of making a system process file using c# which can automatically initiate FME process or FME batch at midnight.
Thanks


texts: Fixing "The SDK seems invalid" error when setting up PyCharm and PyQGIS?

texts:
I wanted to set up PyCharm the way this blog post explains: Setting up PyCharm for PyQGIS and Qt.

I did everything as instructed but when it comes to adding project interpreters I got this error message:

Does anyone have a solution to this?


texts: How to duplicate features based on attribute value pyshp?

texts:
I have a shapefile in which there are a series of polygons that I have to attribute to the name of a project which are recorded in 6 attribute fields.
For each polygon there might be one to 6 projects associated.
I need for each instance of a project name to have one polygon (up to 6 identical but for the name of the project) associated so that I can later dissolve the polygons according to the projects names.

I've used the pyshp module but I encounter a problem with the writing of the new polygon shapes based on the existing polygons. I read the existing using shapes[i].parts but it returs [0] to most of the shapes. Could you please give me a hand solving this issue ?

Here is my code :

import shapefile

#### Set input and target feature
inputShp = r"Pre-treatment\FP_BVC_Polygons.shp"
targShp= "Pre-treatment\FP_BVC_Polygons_test"
sf= shapefile.Reader(inputShp)
shapes= sf.shapes()
records = sf.records()

fields=sf.fields
w=shapefile.Writer(shapefile.POLYGON)

for i in range(0,len(records)):
    for j in range(12,19):
        if records[i][j] != "0":
            listGeom.append(list(shapes[i].points))
            listRecord.append(records[i])

del i,j

for i in range(0, len(listGeom)):
    w.poly(points=listGeom[i])

for i in range(0,len(fields)):
    w.field(fields[i][0],fields[i][1],fields[i][2],fields[i][3])

for rec in listRecord:
    w.record(rec[0],rec[1],rec[2],rec[3],rec[4],rec[5],rec[6],rec[7],rec[8],
         rec[9],rec[10],rec[11],rec[12],rec[13],rec[14],rec[15],rec[16],rec[17],rec[18])

w.save(targShp)


texts:
I am trying to clip a raster to the shape of a polygon as shown in image 1.

But the result that I get is as shown in image 2.

Why is this happening?


texts: How to speed up raster to polygon conversion in R?

texts:
I've decided to process my Landsat data in R instead of ArcGIS - due to my missing knowledge of python and because of (assumed) high computation capacities of R. I want to :

import r1 raster to R,
import shp1 convert raster r1 to shp r.to.poly (dissolve = TRUE)
intersect converter raster r.to.poly with my polygon shp1
calculate area of every created polygon of intersected shp

Thus:
# read shp
shp <-readOGR(dsn = "C://...",
    layer = "m")

#read raster
r1<-raster("r1.tif")

# convert raster to polygon, dissolved neighboring same values
r.to.poly<-rasterToPolygons(r1, dissolve = T)

# define the same projection
proj4string(shp) <- proj4string(r.to.poly)

# use intersection from raster package
int.r <-raster::intersect(r.to.poly,shp)

# calculate area per polygon
int.r$area <-gArea(int.r, byid = T)

# export shapefile
writeOGR(int.r, dsn = "C:/...",
          layer = "...", driver="ESRI Shapefile", overwrite = TRUE)

So far, so good, but it takes about an hour to run the single conversion! moreover, when I tried FOR loop, my R on Windows crashed twice... It runs on mac, for the moment. Where the problem could be and how can I increase computation speed? Am I running out of R memory? The raster size on my disk is only 779 580 byte, size of shp is 1 729 532 bytes, thus are small. Also, make the same task in ArcGIS takes only couple seconds.
I've found some related discussion here: Increasing speed of crop, mask, & extract raster by many polygons in R? but as I have only about 10 rasters to process I don't want to start with parallel processing...


texts: Esri cache white tiles levels L10 L11 L12 with Geowebcache

texts:
I'm trying to publish a WMTS using the basemap from Esri cached tile. The cache is called by geowebcache 1.7.2 and displayed via a openlayers view.

The problem I deal with is that the levels L10/L11/L12 of my basemap are displayed in white whereas all of the other levels of tiles are working well.
Nothing appears in the logs of geowebcache or in firebug (F12).
Do you have any idea about this problem?


texts: Google Earth Engine: checking if ee.List() object includes a specific element

texts:
I want to include ee.List.contains() in a conditional if statement, perform something similar to the following operation, but am getting some unexpected results. In the following code, test.contains(1) returns false, however if(test.contains(1)) seems to be returning true. Can anyone explain why this is occurring and recommend an alternative approach?

var a=5
var b= 3
var c=2
var test=ee.List([a,b,c])
print(test.contains(1))
if(test.contains(1)){print('yes')}
else {print('no')}


texts: How do partition tables improve query performance in PostGIS?

texts:
I have static road, country and state datasets in my PostGIS database. Road networks are partitioned by their containing country. Countries, States and other areas reside in their own table and are partitioned by 'version' or 'release'. All spatial data is fully indexed with GiST on all geometries.

The data is used only for look ups (read-only), with updates being performed periodically as a full database refresh.

All nearest-neighbour queries require a minimum of two spatial queries. One to derive the country where the query is located - which allows me to select the correct partition table. Then a second spatial query performs the actual NN search.

Assuming that spatial queries like nearest-neighbour are O(n log n) then partition tables should result in improved performance. That should be true if the spatial query which finds the correct partition table doesn't take longer than the query on an unpartitioned table.

I have been told that partition tables always the way to go when working with global road networks, but to me I think surely it can be just as, or even more efficient to lump all the world's roads into a single PostGIS table.

In what situations should partition tables be used for spatial data?


texts:
How can I get this ArcMap SQL Trim Def query code correct?

I have a attribute field that looks like:
LOC: R1, THEN A BIG LONG BUNCH OF WORDS
or
LOC: L10, THEN A BIG LONG BUNCH OF WORDS

What I need to display on the label for the asset is just R1 or L10 (bearing in mind these answers are different per asset e.g. R1, R2, R3, R4, R5, R6, R7, R8, R9, R10 and then the same for the L's.

All my searches online keep finding regular SQL answers which ArcMap doesn't like.

Basically I don't want to see LOC: or everything after L10 (for example) till it runs out of words.

Here's where I stopped and tried more web searching for a better answer

LTRIM(COLUMN_NAME, 1, 5) AND RTRIM(COLUMN_NAME, 9, 250)

What can I do?


texts: Compare multiple fields using SearchCursor dictionary and UpdateCursor

texts:
I have two feature classes. I need to transfer attributes from a point feature class to a line feature class based on a combination of fields. They do not have matching field names. Here is an example of what I want to compare between my SearchCursor dictionary and my UpdateCursor table:

if row[0] in intdict AND row[2] in intdict != row[2] AND row[3] in intdict == row[3] AND row[4] in intdict = row[4] THEN row[6] = row[1] in intdict

intdict = {row[0]:row[1:] for row in arcpy.da.SearchCursor(interSectPts,['FID_ADDRSTAND','STREET', 'POINT_X','POINT_Y'])}
with arcpy.da.UpdateCursor(addrLyr,['FID_ADDRSTAND','STREET','STARTX_2','STARTY_2', 'ENDX_2', 'ENDY_2','FROM_','TO_']) as cursor:
     for row in cursor:
         if (row[0] in intdict):
            update row with value
            update row with value
            update row with value
            etc.....

To add context, I am changing the values of a road segment in the TO_ and FROM_ columns to show the street names of intersecting roads. For example a segment named 'Main St' could intersect '1st st' on the west end and '2nd st' on the east end...thus the attribute table would look like:

STREET   |    FROM_     |     TO_
Main St  |   1st St     |    2nd St

I basically an comparing the x,y coordinates of the intersection point feature class and its street names with the x,y of the starting and ending coordinates of street segments. When an intersection occurs, I am grabbing the values of the street name. My primary question is how to structure my 'if' in the update cursor...the rest, I feel like I could do....


texts: Finding cumulative area in ArcGIS for selected polygons

texts:
I have 50 polygons of 5 different classes viz. A,B,C,D,E. and a field in the attribute table for area of each polygons. I want to add a field which give the cumulative area of the polygons of a particular class. e.g. if I query the attribute table for class A polygons, the cumulative area field give the cumulative area of the class A polygons. Likewise for class B polygons and so on. How can I do this in ArcGIS? The field calculator in ArcGIS gives the cumulative area of all the polygons, but not of particular class after query?

Thanks.


texts: What does the ArcGIS10.2 sde_state_id suffix 'L' mean in delta tables, when read using pyodbc in Python?

texts:
Enterprise geodatabase sits in MS SQL Server database. I'm using Python's pyodbc SQL driver to read SDE delta tables. I'm running the following query:

select
        delete_table.SDE_STATE_ID,
        delete_table.SDE_DELETES_ROW_ID,
        delete_table.DELETED_AT,
        add_table.sde_state_id,
        add_table.OBJECTID

        from db.schema.a123 add_table
        full outer join
        db.schema.d123 delete_table
        on add_table.sde_state_id = delete_table.DELETED_AT
        and add_table.OBJECTID = delete_table.SDE_DELETES_ROW_ID

        where  add_table.OBJECTID = 9005656
        or delete_table.SDE_DELETES_ROW_ID = 9005656
        order by add_table.sde_state_id desc

from which I get an SDE_STATE_ID from the add table that has an L appended to it.
'sde_state_id': 3646872L

What does the L mean?


texts: Time enabled service has different time extent in ArcGIS Online?

texts:
I have a time enabled layer that I have published to ArcGIS for Server 10.3

There are times ranging from 5am to 10am (all on one day).

This flows through to the ArcGIS Server service as expected:

However, when I then add this same ArcGIS Server service to the embedded map viewer within my ArcGIS Online account, it comes up with the wrong time range by default.  If I change the time slider to the correct time window and choose to increment it in 5 minute intervals (which is how I am successfully doing it within ArcGIS Desktop) it fails to show any data.

How can I configure this so that ArcGIS Online has the right time extent to match the ArcGIS Server map service,


texts: Using Rule-based/Geometry generator style in QGIS

texts:
I have layer for comments in qgis, symbology style is  rule based with Geometry generator.

Problem is with displaying name in layers tab because the line before text is too long

How do I solve the problem so that the line was as short as on all other layers?


texts: Incorrect extent when importing shapefile to PostGIS

texts:
Using PostgreSQL 9.6, pgAdmin 4 v3.5 and QGIS

I’ve found quite a bit of info on this issue but no clear answer. I’m attempting to import a shapefile that was created using ArcGIS (don’t know what version) by importing it using the loader for PostGIS. The loader appears to import it correctly and the multipolygon displays correctly in the pgAdmin’s geometry viewer but without having a background map. I have read that when importing Esri shapefiles into PostGIS this is a normal phenomenon but not what to do about it (see SQL after import via loader: below).

I’ve tried adding a SRID by using the 2 following lines and after which appears correctly in the geometry_columns as having SRID 4326.

UPDATE "public"."400_map_grid" SET geom = ST_SetSRID(geom, 4326);
ALTER TABLE "public"."400_map_grid" ALTER COLUMN geom TYPE geometry(MultiPolygon, 4326);

But now the multiploygon no longer appears in the viewer, but the background map now does.

When I import it into QGIS it show to have WGS 84 (4326) but the extent appears to be incorrect (see attached image1):

2977708.89153, 9976225.51161 - 3236908.88834, 10200225.50756

Other shapefiles I have that I didn’t receive from ArcGIS shows the WGS 84 (4326) with extent:

-97.97695, 29.92911 - -97.52100, 30.75366 (see attached image2).

How do I correct the extent to match the latter?

SQL after import via loader:

-- Table: public."400_map_grid"
-- DROP TABLE public."400_map_grid";

CREATE TABLE public."400_map_grid"
(
gid integer NOT NULL DEFAULT nextval('"400_map_grid_gid_seq"'::regclass),
objectid bigint,
perimeter numeric,
grid_id character varying(20) COLLATE pg_catalog."default",
geom geometry(MultiPolygon),
CONSTRAINT "400_map_grid_pkey" PRIMARY KEY (gid)
)
WITH (
OIDS = FALSE
)
TABLESPACE pg_default;

ALTER TABLE public."400_map_grid"
OWNER to postgres;

-- Index: 400_map_grid_geom_idx
-- DROP INDEX public."400_map_grid_geom_idx";

CREATE INDEX "400_map_grid_geom_idx"
ON public."400_map_grid" USING gist
(geom)
 TABLESPACE pg_default;


texts: Value of longitude of central meridian for sinusoidal projection

texts:
I am trying to project a PDS 3 dataset into sinusoidal projection using ENVI. How do I calculate the longitude of central meridian value?

It is required for projecting the dataset. I have four corner reticules selenographic coordinates of the dataset.


texts: Problems building a world map showing population size in QGIS 2.10 using Cartogram plugin

texts:
I am using the Cartogram plugin to build a world map showing population size. I have all required data in the attribute table.

My source layer  looks like this:

The plugin starts processing, but the output layer looks like nothing happened:

What could the problem be?
I  used this dataset for polygon layer and this  downloaded as a CSV file, with CRS 4326.


texts:
I am making a choropleth map of population estimate by counties in United States. I am almost done but the problem is a few counties have population count in seven figures and they are the ones that are not being filled. Can somebody help me fix that? I think something's wrong with this part of the code.

var color = d3.scaleThreshold()
.domain(d3.range(2,10).map(function(d) { return d * 1000; }))
.range(d3.schemeReds[9]);

Here's a link to my Gist on Github
https://gist.github.com/sabeehtaqi/3c680626970a66ad63058b08f2a042a4


texts: Raster to vector (two attributes per feature)

texts: I have a raster map of land-use which I want to convert to shapefile. I've used the raster to vector tool from QGis, and I got the polygons for each land-use classes; in the attribute table I got one land-use per polygon and I want to know if there is a way to have two attributes per polygon, the inner and outer land-use (or left and right land-use).


texts: Dynamic leaflet scale / zoom / bounds in R Shiny Applications

texts:
I've dug into a few older stack exchanges on this topic, but am at a standstill.

I am attempting to let leaflet dynamically resize in a large shiny app when the application is resized. I'm doing that through an event listener tied to the leaflet's bounds, which are handily stretched when the ui is resized.

It's doing some resizing, but waits for certainly cut points, then the map resizes in a large choppy jump. With chrome developer console on and watching, the svg gets a new translate3d with a scale parameter that only catches and runs at 0.5 and 2, e.g. the translate 3d works pixel by pixel, then the scale pops in at some threshold at runs either 0.5 or 2, then disappears again.

<svg class="leaflet-zoom-animated" width="768" height="313" viewBox="-22 -6 768 313" style="transform: translate3d(-22px, -6px, 0px);">....</svg>

I've tried working the leaflet options as best I can (see code below, see here: http://leafletjs.com/reference-1.2.0.html#map-option), but can't seem to get around this large jump in zoom.

Below is the code excerpt I have so far. Anyone have experience in this?

output$leaflet_test <- renderLeaflet({
bb = bbox(NC_counties); c(bb)
#bb = bb+c(+1,+1,-1,-1)*.8; c(bb)
#bb = c(-84, 34, -77, 36.0)
# E W
# S N
leaflet(NC_counties
        ,
        options=leafletOptions(attributionControl=F, zoomControl=F, # turn off controls
                               zoomSnap=0.01, #zoomDelta=0.01,
                               trackResize=F, #doesn't seem to work
                               boxZoom = F, doubleClickZoom = F, dragging=F,
                               transform3DLimit = 2^5,
                               keyboard=F, scrollWheelZoom=F,
                               zoomAnimation = F  #?

                              # zoom to limits? https://rdrr.io/cran/leaflet/man/mapOptions.html
                                )
        ) %>%
  #setView(-80.793457, 35.782169, 6.5) %>%
  fitBounds(lat1 = bb[2], lng1 = bb[1], lat2 = bb[4], lng2 = bb[3]) %>%
  addPolygons(fillColor = ~pal(ALAND),
              weight=1, opacity = 1, color = "white", fillOpacity = 0.9,
              highlight = highlightOptions(weight = 5, color = "#666", dashArray = "",
                                           fillOpacity = 0.7, bringToFront = TRUE),
              label = ~lapply(labels, HTML)#,
              #labelOptions = labelOptions(style = list("font-weight" = "normal", padding = "3px 8px"), textsize = "15px", direction = "auto")
  )
})

# https://github.com/rstudio/leaflet/issues/448 or https://groups.google.com/forum/#!topic/shiny-discuss/P80waHG_hhU
# See observers? observe() or leafletProxy()
# or here https://stackoverflow.com/questions/34082452/shiny-send-an-resize-event-when-ui-element-gets-resized
# or the leaflet documentation on proxy: https://rstudio.github.io/leaflet/shiny.html
#observe({leafletProxy("leaflet_test") %>% fitBounds(lat1 = box[2], lng1 =     box[1], lat2 = box[4], lng2 = box[3])})

observeEvent(input$leaflet_test_bounds, {
  bb = bbox(NC_counties) # consistent for all leaflets. Resize them all here.
  #bb = bb+c(+1,+1,-1,-1)*.8; c(bb) #bb = c(-84, 34, -77, 36.0)
  leafletProxy("leaflet_test") %>% fitBounds(lat1 = bb[2], lng1 = bb[1], lat2 = bb[4], lng2 = bb[3])
})


texts: I have difference between Oracle sdo_geom.SDO_Distance and Google Maps

texts:
SELECT sdo_geom.SDO_Distance (
SDO_Geometry (2001, 4326, SDO_Point_Type ( -34.59593955425856, -58.43453386049867,null), null, null),  --  Parque Centenario
SDO_Geometry (2001, 4326, SDO_Point_Type ( -34.60643742373489, -58.376673444758225,null), null, null)  --  Obelisco
,1)
FROM dual;

-- Google Maps show: Distancia total: 5.10 km (3.17 mi)
-- Query returns: 6473,92625918709 m
Does anyone know why there is this difference?


texts: Create time series chart with linear trend using NDVI values at specific locations (9000) stored in a feature collection in Google Earth Engine

texts:
I try to find a correlation between a linear trend of two NDVI time series before and after a specific day in time. And that I want to repeat (iterate/map?) for a lot of points globally distributed. I created a feature collection (asset in GEE) with all the points (coordinates), the date of a conflict (for which I want to see if there is a change in NDVI), and the NDVI value of that point. I used the GIMMS NDVI coarse resolution dataset (found in GEE) for the NDVI values.

Function to create the feature collection worked with following code.

  .map(function(f) {
    var GIMMSndvi = ee.ImageCollection('NASA/GIMMS/3GV0')
                      .filterDate('1989-01-01', '2013-12-16')
                      //.filterBounds(Conflicts_ALL) //*********** this doesn't wrok since the entire GIMMS is one image, you should use clip()
                      .map(function(img) {return img.clip(Conf25)})
                      .select(['ndvi'])
      .map(function(img) {
        var d = img.reduceRegion({
                    reducer: ee.Reducer.first(),
          geometry: f.geometry(),
          scale: 8000
        });
        return f.set(d).set('GIMMSdate', img.date().format('yyyyMMdd')); //Do I need this operation still? Return of date?  /////************sure, otherwise how do you know the time of your value
      });

    return GIMMSndvi; //global variable, use later in script

  }).flatten()

Now I'm stuck in the next step. I need to create a function which goes through the FC (598 rows per point), creates an NDVI trend for all values till the day of conflict and safe that slope. Then the same for the second half of the table. Then I have to find the difference between these two and safe it. And this calculation should be repeated for all the points (around 9000). And in the end, I want to have one overall trend estimate.

I have been thinking and trying a lot the past weeks, and wonder if this task is doable in GEE. Should I use plots add the trend lines and save the slopes, e.g. I could plot all time-series and average them to an overall trend (the problem is I need always two trend lines per point). Or should I do it rather with the linear trend function (ee.Reducer.linearFit())?

Could someone could give me a hint or two on how to do it?


texts: Extract street right of way width from ROW polygon using street centerline

texts:
In the image below, I want to extract the width of each street's ROW, which would be twice the distance of each streets centerline to the nearest parcel boundary.

Its slightly difficult to see, but I have a sample measurement drawn on Division ST of about 60ft.

Eventually, I want to have a field in the street centerline called row_width which would be this measurement of 60ft, etc. Any shortcuts to handling this short of measuring each street's actual row?


texts: Display and aggregate time series based on map selection

texts:
My data is geospatial and consists in records of concentrations that are measured every day on fixed locations. I would like to display the time series of these records when I select one point of the map. If possible, I would like also to do an aggregation of the measures if I select many points.

I have been scrolling through QGIS documentation and look at plugins (TimeManager is great but does not do the job), widget, actions, tip preview, but I have not found a solution yet.

I precise that I run QGIS 3.2 and that my data is stored in PostGIS and every record has: lon, lat, date, measure, name where lon, lat are fixed for a given location name.

This can be done with Tableau, but I struggled to reproduce something similar with QGIS.

Example with Tableau: locations are selected on the right and time series of the selected locations are aggregated and display below


texts:
I'm running into a situation (using R on MacOS) where sf() is giving different levels of precision on the CRS. This causes problems for geospatial operations, as indicated below with st_crop().

Can anyone explain what's going on here? I'm admittedly puzzled.

Reproducible example:

# Download data
download.file("http://sedsh127.sedsh.gov.uk/Atom_data/ScotGov/ZippedShapefiles/SG_UrbanRural_2016.zip",
              destfile = "data/SG_UrbanRural_2016.zip")
unzip("data/SG_UrbanRural_2016.zip", exdir = "data")

urbanrural_sf <- st_read("data/SG_UrbanRural_2016.shp") %>% st_transform(27700)

download.file("https://borders.ukdataservice.ac.uk/ukborders/easy_download/prebuilt/shape/Scotland_ca_2010.zip",
              destfile = "data/Scotland_ca_2010.zip")
unzip("data/Scotland_ca_2010.zip", exdir = "data")

admin_lev1_sf <- st_read("data/scotland_ca_2010.shp") %>% st_transform(27700)

centralbelt_region <- st_bbox(c(xmin = 224479.2, xmax = 642963.5,
                      ymin = 347475.0, ymax = 711014.5),
                    crs = st_crs("+init=epsg:27700")) %>%
  st_as_sfc()

Running st_crs shows how an extra decimal of precision is present on the centralbelt_region polygon (scroll to the right to see). Otherwise CRS are exactly the same:

> st_crs(centralbelt_region)
Coordinate Reference System:
  EPSG: 27700
  proj4string: "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.1502,0.247,0.8421,-20.4894 +units=m +no_defs"
> st_crs(admin_lev1_sf)
Coordinate Reference System:
  EPSG: 27700
  proj4string: "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs"
>

This causes problems when I try to run:

> urbanrural_sf_simplified_centralbelt <- st_crop(urbanrural_sf_simplified, centralbelt_region)
Error in geos_op2_geom("intersection", x, y) :
  st_crs(x) == st_crs(y) is not TRUE

Would be really glad for any insights folk have on this.


texts: Creating buffer rings of equal population size using ArcMap

texts:
I am trying to create a shapefile of equal population size buffer rings. Similar to the procedure of creating multi-ring buffers of distance around a point, find the intersects of the buffers with a polygonal layer, and then dissolve with the required statistics (e.g. https://support.esri.com/en/technical-article/000010582).
I have a point layer (Atlanta) and a polygon layer (Atlanta metropolitan area census block groups). I want to create a new polygon layer of equal population size (say 5% of total population) rings around the point layer, using population attribute from the census block groups. So a total of 20 rings around the city, each equal of 5% of population. Distance and area of each ring obviously differs.
So far I have tried to use ArcMap Grouping Analysis and modifying the process in the link above.
Any suggestions on how to implement this?
I can use Python for this if needed.
I expect to have an output similar to this:


texts: Why is my SSF file opening in Pathfinder Office in read only mode?

texts:
I have a new SSF file that I just copied from my Trimble GeoExplorer 2008 series GPS unit to my desktop.  I also copied the shapefiles from the unit.
When I double-click on the SSF file, I receive the following error.
Changes cannot be made to this file as it was created with GPSCorrect.
The file has been opened in read only mode.

It never actually opens a file though.  I need to be able to run a correction on it and then apply the correction to the shapefile that the data was collected for.


texts: Openlayers popup - no scroll bars on ios or android

texts:
I have an openlayers app with a FramedCloud popup. The popup has a lot of lines of data, so scrolling is needed to get to the bottom. All works well on my pc with Firefox. However, out in the field with an iPhone or Android (running Chrome), I get no side scroll bar, so cannot get down to the end of the form.

Here's some popup code that gives the problem:

function onFeatureSelect(event) {

    var feature = event.feature;

    var content = "" + feature.attributes.art.value +
         " (Id: " + feature.attributes.trädnummer.value +
         ")
Risk: " + feature.attributes.risk.value +
         "
" + feature.attributes.anmärkningar.value +
         "
Åtgärd: " +
         feature.attributes.rekommendationer.value;
    content += "

A string";
    content += "

A string";
    content += "

A string";
    content += "

A string";
    content += "

A string";
    content += "

A string";
    content += "

A string";
    content += "

A string";
    content += "

A string";
    content +=  '
<input type="submit" value="Spara"">';

    popup = new OpenLayers.Popup.FramedCloud(
        "chicken",
        feature.geometry.getBounds().getCenterLonLat(),
        null,
        content,
        null,
        true,
        onPopupClose);

    feature.popup = popup;

    map.addPopup(popup);
}

Does anyone have a solution to this?

Thanks

Mini


texts: How to create a text file containing layer names in QGIS?

texts:
My QGIS ToC contains layer groups and some layers that do not belong to any group.

I would like to create a text file containing the name of those layers that do not belong to any group.

Is there an automated way to do it?


texts: Merging layers with different objectid types

texts:
I have a problem that's stumped me. I'm merging flood hazard layers from FEMA. For most states, the objectid is a real number. California (we're always the weird ones) has an objectid that's 'integer64'. I don't know enough about objectids to know which one is correct or if there's a correct option. Merging the layers into a single layer runs into problems with the difference in object ids. Can anyone suggest a solution to this problem? Is there way to convert one or the other objectid fields to match or a setting to select to enable the different types to work together?

I'm using QGIS.  I pulled the data from FEMA's Esri Geodatabase. I then clip these layers to my study area and save 'em as shapefiles. I'm working with those, trying to merge them into a single layer


texts:
I'm very new to GIS and to pgRouting. We're implementing a Routing system in our app and we need to route between points that we'll define ourself.
So far, we have a working pgRouting on Postgis (on Postgres) and a Spring boot application to call pgRouting's methods. We managed somehow to download from Geofabrik and import OSM files into the running database using osmttools to reduce it's size, but there still are way too much points and it takes too much time to proceed a pgr_dijkstra call.
Now, during the import phase, I need a way to remove all points and leave only those we're.