texts: Legend in map atlas anchored to bottom left or bottom right will not correctly reposition in PDF

texts:
I submitted this as a bug but I'd like to know if anyone here has experienced or has a solution.
In a map atlas, if I anchor a legend to either the bottom left of bottom right and if the legend is set to limit to the features to the current atlas item, it positions itself properly when I traverse the atlas manually, but when I output the atlas to a PDF or individual PDFs, the legend shifts seemingly at random. Sometimes if atlas items have many features in the legend, the legend will fall off the page or position itself higher on the page. Problem occurs with versions 3.14 and earlier.
The only workaround I can see is to traverse the atlas manually and output the PDFs one at a time.
Here are some sample images:


texts: Why do field values change when exporting from QGIS?

texts:
Today I came across a very strange behavior of QGIS. Maybe someone knows what's going on?

I have a table with two fields which are defined as Double, with a Length of 20 and a Precision of 0. Some sample values:

8159080177
8159080742
8159080008
...

The format by the way is shape. Exporting a selection of the features results in a table that has same definitions as stated above but all values for the two columns are

-2147483648
-2147483648
-2147483648
...

I get the same values when I read the features with python using ogr's <feature>.GetFieldAsInteger('<fieldname>'). Reading it as a string returns the correct value.

What is happening here?


texts: ArcGIS Online: possible to build network data set (to be used for distance calculations)

texts:
I am currently using ArcGIS Pro (network analyst, OD cost matrix analysis) to calculate distances between a set of points A and a set of points B. Distances are calculated along my own road map, which I have in a shape-file. To do this, I first build a network data set in ArcMap, which I import in ArcGIS Pro prior to the distances calculations.

I am now trying to move my project from ArcGIS Pro to ArcGIS Online, and I wonder if it is possible to somehow create the network data set (that I need as input in the OD Cost matrix analysis) in ArcGIS Online? Or do I need to have this ready, and publish it to my ArcGIS Online account in order to perform the distance calculations?


texts:
I have a user that recently received a new machine with a fresh install of 10.2. We have a python / arcpy application that worked in 10.1 on the old machine but does not work on the new machine with 10.2. The error message is "ImportError: No module named arcpy".

How do I determine if arcpy is installed on the machine? If it is installed, could  this be a problem with the version?


texts:
When I use both QGis and Geomedia to digitize in the same layer I get some issues in my Oracle table. It looks like Geomedia save my point, lines and polygons as ARRAY, while QGis save it like POINT (to points) and ARRAY (to lines and polygons). The problem cames when a point layer is worked both with Geomedia and QGis --> In my Oracle table all the points digitized with QGis after a point digitized with Geomedia are "snapped" in the Point (coordinate) digitized with Geomedia. Any idea to solve it?
Thanks!


texts: Measuring average, min and max of vertex height of linear shape using QGIS?

texts:
I have a linear shape, all of the vertexes of each segment has a z-value, stored in the geometry.

I need for every record, segment the average, min and max height value aggregated from their vertex heights.

How would I get this information using the field calculator?

I'm using QGIS 3.4 Madeira.


texts: How to save a NetCDF file with proper timestamp in R?


texts: Calculating coastal distances using QGIS?

texts:
I'm trying to work up some historical transport cost maps, across multiple modes (road/rail/sea/river). I have a polygon shapefile of the coast itself, and point data on the cities of interest (some of which are coastal). I'm using QGIS if that is relevant.

Re: coastal routes in particular, I'm not sure what the best way is to generate a coastal distance measure between two port cities A and B.

One approximation is the length of coast between A and B, which I think will be an overestimate given the ruggedness of the coastline. The two ideas I had to improve on this were:

I could use a buffer and then take some kind of smoothed median line parallel to the coastline.
I could simplify the geometries of the coastline to wash out coastal details, and simply use the length of that. (Possibly I could offset the coast by N km and simplify that.)

Does this make sense?


texts: Is Orfeo Toolbox OTB no longer available through OSGeo4W?

texts:
When I run the current OSGeo4W installer on Windows 10, 64 bit, no OTB packages are available:

The OTB binary package status page still lists OSGeo4W as a source for OTB.

~~Does anyone know what is going on here? Is OTB not distributed via OSGeo4W anymore? Or temporarily unavailable?~~

EDIT:

I just discovered the "Hide obsolete packages" checkbox.

So why are OTB packages considered obsolete in OSGeo4W?


texts:
Does anyone know how to get LS-tool mentioned in Zhang et al 2013.
https://www.researchgate.net/publication/256938578_Extension_of_a_GIS_procedure_for_calculating_the_RUSLE_equation_LS_factor

I tried using the Raster calculator in ArcMap 10.3.1 using the following formula by Bizuwerk et al., (2008).

LS = X/22.1m (0.065+0.045S+0.0065S^2
but I got only 0 value.

I found about LS tool after I researched on other ways to calculate slope length but can't find any place to access LS tool.

I would be really thankful if someone could help me find a way to either access LS tool or any other efficient and precise way to calculate slope length.


texts: Folium - plugins.HeatMapWithTime no longer produces an output

texts:
I have been creating Heat Maps with a Time dimension using both test data and my own datasets using the Folium libraries in Python. I am using Python 3.5 to do this, and the results have been good.

I have been away and one returning I have noticed that these scripts no longer produce a result that is viewable.

The scripts run without error, and the outputs open with the correct basemaps, however, the heatmap and the time sliders are no longer visible.

At first I thought it was the data I was using, but I have a test script taken from Git that I have been using to test some formatting and displays, but this no longer works.

https://github.com/python-visualization/folium/blob/master/tests/plugins/test_heat_map_withtime.py

import folium
import folium.plugins as plugins
import numpy as np
import webbrowser

from datetime import datetime, timedelta

np.random.seed(3141592)
initial_data = (
    np.random.normal(size=(100, 2)) * np.array([[1, 1]]) +
    np.array([[48, 5]])
)

move_data = np.random.normal(size=(100, 2)) * 0.01

data = [(initial_data + move_data * i).tolist() for i in range(100)]

weight = 1  # default value
for time_entry in data:
    for row in time_entry:
        row.append(weight)

print (data)

time_index = [
    (datetime.now() + k * timedelta(1)).strftime('%Y-%m-%d') for
    k in range(len(data))
]

print (len(time_index))

m = folium.Map([48., 5.], tiles='stamentoner', zoom_start=6)

hm = plugins.HeatMapWithTime(
    data,
    index=time_index,
    auto_play=True,
    max_opacity=0.3
)

hm.add_to(m)
m.save (r'c:\temp\Time.html')
webbrowser.open(r'c:\temp\Time.html')

Has something changed with the Folium libraries / anything?
I thought maybe it was a firewall issue at work, but even going offline / using a connection outside of work produces the same result:

A Basemap, with no data to display.


texts:
My organization (a public transit system) will be conducting surveys on our existing infrastructure (bus stops).

My original plan was to send out the surveyors with a  tablet with ArcGIS Collector and Survey123 installed, and a separate GPS unit. The surveyors would select the bus stop in Collector which would open and populate three questions in the Survey123 form, like the unique ID, using a url. The GPS unit will override the tablet's GPS capabilities using a bluetooth connection and the BluetoothGPS app.

The bump I have hit is with Collector.

Collector requires a constant online connection or an enterprise database set-up for offline use, neither of which are available. As of right now, the surveyors could complete the surveys without Collector, but it would be nice to have those pre-populated fields to cut down the rate of human error and save time and effort entering information in those fields for every stop. We have over 5,000 stops that need surveying.

What is an alternative to Collector in this workflow or is there a better workflow?

Note: At this time, I do not have the time or skill set to create a custom data collection app for this project.


texts: Editing and Versioning with ArcGIS Rest API

texts:
I'm developing a web app with ArcGIS REST API.
I created Version A from default version, made some edits on an event layer where I acquired locks for Route A, then deleted Version A. (The API I called in order are: createVersion, locks/acquire, applyEdits, deleteVersion)
Then I created Version B from default version. I try to make some edits on the same event layer, but at this point I couldn't acquire locks for Route A somehow. any ideas why?
If I delete a version, does the locks/edits automatically go away with it, or do I have to release it manually?


texts: Multivalue input from chosen field giving patchy input strings from ArcPy?

texts:
I made a small tool that can be run from a toolbox. It has a validation code to choose multiple values from a specified field from a specified feature layer. It looks quite ordinary:

The validation code is:

def updateParameters(self):
   if self.params[1].altered:
      if self.params[1].value:
         values = set()
         inlayer = self.params[0].value.value
         fieldname = self.params[1].value.value

         cursor = arcpy.SearchCursor(inlayer, None, None, infield)
         for row in cursor:
            values.add(row.getValue(infield))

         self.params[2].filter.list = sorted(values)

It gives me a nice multivalue list of values. But when I check some of them and run a simple script:

import arcpy
# Variables
inlayer = arcpy.GetParameterAsText(0)
infield = arcpy.GetParameterAsText(1)
invalues = arcpy.GetParameterAsText(2)
arcpy.SetParameter(3, invalues)

selected_values = invalues.split(";")

for selected_value in selected_values:
    arcpy.AddWarning(selected_value)

then it returns some unicode strings with an apostrophe and some not. I found out that if there is a gap in the value, it returns an unicode string with an apostrophe. I don´t get it why? There´s a example output of the forementioned script:

In this tool window I can see that these apostrophes already figure as input values!

Why is that, and how do I avoid it?

These inputs with apostrophes troubles me further in my script because

if "withoutgap":
  # returns True but
if "with gap":
  # returns False because the unicode string has apostrophes

running on Windows 8.1, ArcGIS for Desktop 10.1 Basic, Python 2.7


texts: Geomedia Smart Client features, Oracle database

texts:
We have database schema imported in Oracle, but when we want to import it into Geomedia Smart Client (GMSC) certain tables are not visible in GMSC.

On the other hand, VIEW created from these tables has the MDO.SYS Geometry column, is visible in GMSC, but we cannot perform any queries on it.

It seems that the topic:
How does GeoMedia know what feature classes are in a database?

is maybe connected with our problem but we don't now how?
Is there a problem with geomedia's metadata and what? Why can't GMSC see all the features with geometry columns?


texts: How to insert missing vertex of two geometry overlaps area in QGIS

texts:
I need to add the missing vertex of two lines intersection places using QGIS.If any default plugins or tools available for fixing missing vertex in qgis,kindly suggest me.
My Inputs:

In above screenshot one line feature having vertex on highlighted area and another one line has no vertex on that place.I need to place missing vertex on one line.


texts: Is it possible to retrieve the size of a shapefile using python and arcpy?  If so, how?


texts: Moving raster by half cell and estimating values in-between existing cells using ArcMap

texts:
I have a raster at 1-km resolution (GridASCII format with 6 header lines above the space-separated data). I need to produce a raster of the same underlying data at the same 1-km resolution, but with the cell centers offset by 0.5 kilometers. I do not just want to shift the raster (i.e. subtract 0.5 from XLLCORNER and YLLCORNER).
How do I calculate the values at the mid-points using ArcMap 10.6?
The ideal tool I can imagine for this task would be something like Resample (data management), but with specification of XLLCORNER and YLLCORNER permitted.


texts: Setting the extent in cartodb to the sublayer bounds

texts:
How do you set the extent of a sublayer to fit all of the points?

I've been trying to use this bounds function, but it doesn't work:

sql.getBounds('select * from table').done(function(bounds) {
    map.fitBounds(bounds);
});


texts: Raster Attribute Table Where Each Pixel = 1 Row

texts: In ArcGIS, is it possible to create a raster with an associated attribute table where each pixel is represented by a "row" in that table.  I would like to specify unique attributes for each raster pixel, rather than for each "Value" of pixels (the default raster attribute table in ArcGIS seems to be one row per "Value" attribute of pixels).


texts: Get IFeatureLayer from ILayer using ArcGIS Engine?

texts:
I try to select object from layer by name. I think i found sollution here.
First problem i got is that i use GroupLayer. So i made research how to get subLayer from GroupLayer. Sollution is here. Now i have a problem how to get IFeatureLayer if i already have ILayer from where i need to select object.

var mapControl = form.axMapControl2.ActiveView.FocusMap;

var basemapCompositeLayer = mapControl.Layer[0] as ICompositeLayer;
var basemapSubLayer = basemapCompositeLayer.Layer[0];
var featureLayer1 = basemapSubLayer as IFeatureLayer; // always NULL

Or this is because that Layer isn't Featured by default?


texts: How can I show a label for a certain field in an attribute table for a shapefile?

texts:
names = ["Schools"]
layers = arcpy.mapping.ListLayers(mxd, "*", df)

for layer in layers:
  if layer.name == "Schools":
    layer.showLabels = True

My code labels the default field in the attributes table. I want it to label a field in the Schools shapefile called "Name_Type" instead of the default field.


texts: Leaflet.js rendering 'chopped' WMS point symbols for variable sized points

texts:
I am using GeoServer to render WMS for variable sized points based on a data attribute.  This technique is documented in GeoServers SLD Cookbook here.  The SLD that I am creating validates successfully in GeoServer's SLD section of the Admin screen and renders correctly in the OpenLayers preview.  However, when rendering in a Leaflet.js map a specific area of the map consistently returns seemingly incorrect tiles resulting in a 'chopping' like effect.  This may be obvious but I have noticed that the sharp edges are a result of a tile edge where one tile renders a portion of the point and the other is missing the remaining portion of the point.  Since GeoServer is serving complete tiles this seems like Leaflet is retrieving certain tiles that are out of date for that zoom level.

I don't think this is a browser caching issue because the 'chopping' is consistent no matter how many times I recreate the database view and GeoServer services the layer is created from.  One theory I have is that the

Does anyone have any idea why this may be happening and/or how to approach fixing it?


texts: GetFeatureInfo works over HTTP but not HTTPS

texts:
I have GeoServer up and running with ProxyPass allowing access over HTTPS and I have the proxy base URL set to https://maps.tellus-toolkit.com/geoserver/ Everything is working as expected apart from a GetFeatureInfo request for a WMS layer not working via HTTPS but it is working on HTTP.
I can preview the map here:
https://maps.tellus-toolkit.com/geoserver/nhs/wms?service=WMS&version=1.1.0&request=GetMap&layers=nhs%3ASG_NHS_HealthBoards_2018&bbox=5512.998599999584%2C530248.9970999993%2C470332.0037000001%2C1220301.5015999998&width=517&height=768&srs=EPSG%3A27700&styles=&format=application/openlayers
But, a request over HTTPS throws a 403 error
https://maps.tellus-toolkit.com/geoserver/nhs/wms?service=WMS&REQUEST=GetFeatureInfo&SERVICE=WMS&SRS=EPSG%3A4326&STYLES=&TRANSPARENT=true&VERSION=1.1.1&FORMAT=image%2Fpng&BBOX=-12.63427734375%2C54.93345430690937%2C4.63623046875%2C58.955674359706016&HEIGHT=672&WIDTH=1572&LAYERS=nhs%3ASG_NHS_HealthBoards_2018&QUERY_LAYERS=nhs%3ASG_NHS_HealthBoards_2018&INFO_FORMAT=text%2Fhtml&X=833&Y=404
This works though:
http://maps.tellus-toolkit.com:8080/geoserver/nhs/wms?service=WMS&REQUEST=GetFeatureInfo&SERVICE=WMS&SRS=EPSG%3A4326&STYLES=&TRANSPARENT=true&VERSION=1.1.1&FORMAT=image%2Fpng&BBOX=-12.63427734375%2C54.93345430690937%2C4.63623046875%2C58.955674359706016&HEIGHT=672&WIDTH=1572&LAYERS=nhs%3ASG_NHS_HealthBoards_2018&QUERY_LAYERS=nhs%3ASG_NHS_HealthBoards_2018&INFO_FORMAT=text%2Fhtml&X=833&Y=404
Is the problem with Tomcat/Apache reverse proxy setting or is there something in the GeoServer settings I am missing?
Here's the output from Chrome
VM27:1 GET https://maps.tellus-toolkit.com/geoserver/nhs/wms?service=WMS&REQUEST=GetFe…_LAYERS=nhs%3ASG_NHS_HealthBoards_2018&INFO_FORMAT=text%2Fhtml&X=877&Y=242 403
(anonymous) @   VM27:1
send    @   jquery-1.10.1.min.js:6
ajax    @   jquery-1.10.1.min.js:6
getFeatureInfo  @   L.TileLayer.BetterWMS.js:21
fireEvent   @   leaflet.js:6
_fireMouseEvent @   leaflet.js:6
_onMouseClick   @   leaflet.js:6
s   @   leaflet.js:8


texts:
I'm composing a map with aerial footage with mapserver and tilecache. I have tiffs with geo information, then I run gdaltindex to create the shapefiles and then I configure the map file with layers like this:

LAYER
    NAME "viana"
    STATUS default
    TILEINDEX "viana.shp"
    TILEITEM "location"
    TYPE raster
    OFFSITE 0 0 0
     PROJECTION
       "init=epsg:27492"
    END
    PROCESSING "RESAMPLE=BILINEAR"
    PROCESSING "DITHER=YES"
    TRANSPARENCY ALPHA
    MAXSCALEDENOM 54000
    MINSCALEDENOM 423
END

How do I set a default image in mapserver for when I'm have no data for some specific area, like an "Image not available at this zoom"? Like http://12f.p8.sl.pt


texts:
Is it possible to make a copy of a tool from the Processing Toolbox in QGIS and make edits to it?

As an example, I would like to make a copy of the Import geotagged photos tool and then edit the code so that it also brings in extra fields from the EXIF data in the photo.

Note: the goal of this question is figuring out how to edit a processing tool, not to bring in the EXIF data another way.


texts:
I am completely new to the QGIS python api and I am running into an error which I cannot solve when trying to use the dissolve class.
I used the basics of the PyQGIS Cookbook. In my  understanding the problem lies within the "layer". But as I understand it, that is exactly whats assigned to that variable.
Probably it is just a simple error I am running in, as I am not experienced using python.

In essence I just need a script which performs a dissolve on a shapefile with the option "dissolve all" like it is implemented within QGIS.

Thankful for any help.

Andreas

Code:

from qgis.core import *
from qgis.analysis import *

QgsApplication.setPrefixPath("/usr", True)
QgsApplication.initQgis()

data_source = "test.shp"
layer_name = "test"
provider_name = "ogr"

layer = QgsVectorLayer(data_source, layer_name, provider_name)
if not layer.isValid():
    print "Layer failed to load!"

QgsGeometryAnalyzer.dissolve(layer, "test_dissolved.shp", onlySelectedFeatures=False, uniqueIdField=-1, p=None)

Error:

Original exception was:
Traceback (most recent call last):
  File "/home/ubuntu/workspace/testproject/dissolve.py", line 17, in <module>
QgsGeometryAnalyzer.dissolve(layer, "dissolved.shp", onlySelectedFeatures=False, uniqueIdField=-1, p=None)
TypeError: QgsGeometryAnalyzer.dissolve(QgsVectorLayer, QString, bool onlySelectedFeatures=False, int uniqueIdField=-1, QProgressDialog p=None): first argument of unbound method must have type 'QgsGeometryAnalyzer'


texts: GeoDjango: how to calculate distance between two Point()?

texts:
I've noticed that in GeoDjango the distance method behaves differently if it's being used inside an annotation, a query etc.. or if it's used between two points or geometries.

Look at this example:

bss = BusStop.objects.filter(point__distance_gt=(base_point, D(m=1000))).filter(point__distance_lt=(base_point, D(m=1200))).annotate(distance=Distance('point', base_point)).order_by('distance')
In [20]: for bs in bss:
...:     print(bs.amic_code, bs.distance)
...:
1800SBH0781 1002.00849246 m
1800SB30351 1002.56747241 m
1800SB04941 1003.2266972 m
1800SB30331 1008.46770311 m

The distance is correctly expressed in meters, as I expect.

This is instead what happens if I use the distance method between two points:

In [23]: type(base_point)
Out[23]: django.contrib.gis.geos.point.Point

In [24]: bss[0].point.distance(base_point)
Out[24]: 0.014417147670674185

I can't really understand how the distance is being expressed. What I get back is just a float. Having the distance in meters would be nice, but even better if I can have a measurement object where I can get all the other formats: https://docs.djangoproject.com/en/1.10/ref/contrib/gis/measure/

Note: I've read this https://docs.djangoproject.com/en/1.10/ref/contrib/gis/geos/#django.contrib.gis.geos.GEOSGeometry.distance where it says "GEOS distance calculations are linear – in other words, GEOS does not perform a spherical calculation even if the SRID specifies a geographic coordinate system.", but it doesn't explain how to do a spherical calculation.

How can I get the distance in meters like in the first use case?

Update: I've tried to use this library https://pypi.python.org/pypi/geopy/1.11.0 and I've nothing against using it, but... if I try to calculate the distance between one of the points and the base_point I get 1.6 Km that is quite unexpected, given that I have a list of points which have a distance > 1000 meters and < 1200 meters. Something is still wrong.


texts: Create ID for each string attribute in QGIS

texts:
I'm trying to rasterize a vector layer by using an attribute filed stored as string, I got the confirmation through the GDAL mailinglist serv that gdal_rasterize does not work with string field. So my question is: is there an easy way to create an ID column based on a existing string column? Basically that would something like this:

"Forest and Woodland" ----- 1
"Shrubland and Grassland"---2
...

and so on

I only have 7 district strings so I suppose I could write some basic python if loop on the field calculator but what if this value was greater?

I'm using QGIS 2.8 and latest GDAL.


texts: Copy GeoTIFF metadata to TIFF using ArcPy

texts:
I have conducted some image manipulation on GeoTIFFs in a program and the resulting TIFFs now have no spatial information attached.  Their x,y origins are (0,0).  Since these processed images are of the same exact area as the original images, can I simply copy the metadata from the original images to the processed images?

I have seen some gdal solutions, but I'm looking for a solution in Arcpy.


texts: Azimuth, elevation and distance for point-to-point radio links

texts:
After plotting our point-to-point radio sites and links in QGIS using Ordnance Survey map layers I'm wondering whether I can extract some information for this.

For each link I would like to determine the azimuth and elevation of the antennas at either side such that the narrow beam from them will target the opposite site's antenna.

Furthermore, it would be useful to understand the distance of each link. I believe using the length attribute in QGIS gives me the length of a great circle rather than the distance the wave will travel.


texts: Automating Select by location using ArcPy and where_clause?

texts:
I have used the code in Coding python loop to select features of polygon layer than are completely contained in group of 11 layers using ArcPy? to automate select by location analysis using ArcPy.

I have a polygon layer "Provinces.shp" which specify number of cities, and a point layer "Points.shp". The purpose is to find the intersected points in each cities and export these points as separate layers. Thus, I need to first select city names based on "NAME_1" field in the polygon layer and then intersect it with point layer literately.

The code is as below:

import sys, string, os, arcgisscripting
import arcpy

# Script arguments...
A = arcpy.GetParameterAsText(0)
B = arcpy.GetParameterAsText(1)
outfolder = arcpy.GetParameterAsText(2)

# Create a Search Cursor that references the features in layer B
rows = arcpy.da.SearchCursor(B, "*")
row = rows.next()

# Create a numbering system for new output files
filenum = 0

while row:

    # Report a message back to the user about progress
    arcpy.AddMessage("Creating selected set #" + str(filenum) + " from " + A)
    #Create a feature layer to cycle through your selection features
    fc = r'C:/Users/Zarei.f/Desktop/data-module 6/IRN_adm1.shp'
    layerFld = arcpy.AddFieldDelimiters(fc,"NAME_1")

    sql =  """{0} = {1}""".format(layerFld , str(row[6])

    #fidvalue = row.getValue("NAME_1")

    arcpy.MakeFeatureLayer_management(B,"cursor_lyr5"+str(filenum),sql)

    # Process: Select Layer By Location...
    arcpy.SelectLayerByLocation_management(A, "INTERSECT", "cursor_lyr5"+str(filenum), "", "NEW_SELECTION")

    outname = outfolder + "\\" + 'test' + str(filenum)
    # Process: Copy Features...
    arcpy.CopyFeatures_management(A, outname, "", "0", "0", "0")

    #iterate the output file number
    filenum = filenum + 1

    #go to the next row
    row = rows.next()

I got the below error:

ExecuteError: ERROR

    000230: Failed selecting with NAME_1 = Alborz
    Failed to execute (MakeFeatureLayer).

Regarding to this error, the MakeFeatureLayer analysis can't understand the defined sql clause.

I used several ways to define the sql clause but they all raise a same error.


texts: OpenLayers 3 can't draw after loading GeoJson

texts:
I created a function which draws and saves to GeoJson. I got : result from save

Next when I load features it's ok, but when I try drawing I get error:

"layer.getSource().getFeaturesCollection()" == null

But I saved features with type featureCollection. When I use getFeatures() it returns an array of features.

Begin of load layers:

if(map_info != "") {
            var layer_area = new ol.layer.Vector({
                source: new ol.source.Vector({
                    features: (new ol.format.GeoJSON()).readFeatures(map_info['layers'][0])
                }),
            });
            var layer_branch = new ol.layer.Vector({
                source: new ol.source.Vector({ features: (new ol.format.GeoJSON()).readFeatures(map_info['layers'][1]) })
            });
            var layer_region = new ol.layer.Vector({
                source: new ol.source.Vector({ features: (new ol.format.GeoJSON()).readFeatures(map_info['layers'][2]) })
            });
            var layer_ambon = new ol.layer.Vector({
                source: new ol.source.Vector({ features: (new ol.format.GeoJSON()).readFeatures(map_info['layers'][3]) })
            });

            var view=new ol.View({
                center: [{{map.pos_x}}, {{map.pos_y}}],
                zoom: {{map.zoom}}
            });
        }else {
            var layer_area = new ol.layer.Vector({
                source: new ol.source.Vector({ features: new ol.Collection() })
            });
            var layer_branch = new ol.layer.Vector({
                source: new ol.source.Vector({ features: new ol.Collection() })
            });
            var layer_region = new ol.layer.Vector({
                source: new ol.source.Vector({ features: new ol.Collection() })
            });
            var layer_ambon = new ol.layer.Vector({
                source: new ol.source.Vector({ features: new ol.Collection() })
            });
            var view=new ol.View({
                center: [2097096.387839827, 6887410.172738342],
                zoom: 5
            });
        }

Draw functions:

$("#draw_ambon").on("click", function(){
    var layer=manager["ambon"];
    map.removeInteraction(draw);
    select !== null ? map.removeInteraction(select) : 0;
    setInteraction(layer.getSource().getFeaturesCollection(), "Point");
    layer.setStyle(pointStyleFunction);
    layer.setMap(map);

    draw.on("drawend", function(e){
        layer.getSource().addFeatures(e.feature);
        e.feature.setProperties({
            'id': layer.getSource().getFeaturesCollection().getLength(),
            'type': "ambona",
            'name': ''
        })
    })
});


texts:
Are there any CRS that represents south up in a map?
Just like this:


texts:
I have 3 shapefile: pipes (line), junctions (point), customers (point). I want to bring customer demand value from customer to junctions with conditions:

First condition: Nearest road
Second condition: the nearest node

I used "Distances nearest hub" but it only gives the result for second condition.
Is there any way to solve this problem?


texts: Spatial mapping of land use change intensity

texts: I have land cover map for two different time periods 2000 and 2010. I want determine the land use change intensity across space. That is, a spatial distribution/pattern of how the various land use types changes within the two time periods. Which areas/pixels have the highest change. I also want to show the change in ha, i.e. high and low value. Any suggestion on how to go about this using ARCGIS?


texts: How to check spatial autocorrelation correctly in a glm

texts:
I ran a binomial glm as below and from plotting observations i know that there are many observations close to each other and i expect them to be highly correlated. For this matter after selecting the best model, I did the following to check for spatial autocorrelation (SAC) using Moran's I (not sure if u did this correctly too). But after much reading I found you can't use Moran's I to test the SAC of a model with binomial response and should use join count.test.

#best model
best_m <- glm(response ~ predictor 1 + predictor 2 +predictor1*predictor2, family=binomial("logit"), data= data1)

#check Spatial autocorrelation
coords<-as.matrix(cbind(data1[,2],data1[,3])) #x and y geographic coordinates
nb <- knn2nb(knearneigh(coords, 1))
lstw <- nb2listw((knn2nb(knearneigh(coords, k=1))))

# Compute Moran's I using residuals of model and also raw data
moran.test(data1$response, nb2listw(nb, style = "B")) #check for raw data
moran.test(best_m$residuals, lstw) # there is SAC for model residuals


texts: Using 2 different for loops in Python, skips second process

texts:
I am trying to do the following in Python 2.7:

Clip all files in a geodatabase with a single layer (so a batch clip)
Export the resulting clip layers (identified with _clip on the end of their names) to _ready layers, only pulling out features with a Shape_Area under 100,000
Delete the _clip layers to clean up the GDB.

I am doing all of this within 2 File GDBs. I originally created a model in ArcMap 10.6 ModelBuilder to get the basis, then exported to Python and rewrote it so it can be used outside of a toolbox.
I find when I run the below code, it finishes doing the clip successfully (result can be seen in the GDB), then basically skips over the rest with no error or result.
I am very new to Python, so suspect I have missed a core principle. Maybe you cannot have 2 for loops defined like this. If I separate them and run as two different scripts each process works fine, which tells me the construct of each process is correct, but combining them in this way is not.
Can someone guide me as to what I am doing wrong here?
I have done much research on Google to no avail.
I have altered path and file names to generic ones.
#IMPORT MODULES
import os
import arcpy

#SET ENVIRONMENT
arcpy.env.workspace = "E:\\Folder\\GDB_Name.gdb"
arcpy.env.overwriteOutput = True

#SET VARIABLES
inputfc = arcpy.ListFeatureClasses()
clipExtent = "Database Connections\\My Connection.sde\\My Layer"
inputfc2 = arcpy.ListFeatureClasses("*.clip")
outLocation = "E:\\Folder\\Another GDB Name.gdb"
Delete_succeeded = "true"

#PROCESSES
#CLIP
for fc in inputfc:
    arcpy.Clip_analysis(fc, clipExtent,fc + "_clip")

#FC TO FC
for fc in inputfc2:
    arcpy.FeatureClassToFeatureClass_conversion(fc, outLocation, fc + "_ready","Shape_Area > 100000")

#DELETE CLIPS
for fc in inputfc2:
    arcpy.Delete_management(fc, "")


texts: How to create in QGIS shapefile a field that automatically returns the current date

texts: I'm digitizing a shapefile in QGIS and I need to create a field that returns automatically the current date.


texts: QGIS geolocate a csv (insert point for each record manually in map)

texts:
I have a CSV with data (no location, no lat/long). I need to load that CSV into QGIS and create a point for each feature.

I've trying for an hour without success..

I guess there must be a way to create a new empty SHP, upload CSV data into the SHP somehow and after that assign one point in the map for each element..but I can't find a way to do this in QGIS..


texts: How to create a layered PDF document using MapInfo version 12.5

texts:
I am wanting to create an interactive PDF layer using a MapInfo workspace where by layers can be turned on and off in the pdf. I also have a large aerial image I would like included in the pdf .

Any advice/ tips would be greatly apprecaited


texts: Opening netCDF and creating raster file?

texts:
I want to open a netCDF file and create a raster with ArcMap 10.2. The file contains the mean monthly temperature for Switzerland. I found it in meteoswiss web site.

I tried to open it with the multidimension tool and make a netCDF raster layer but when I select the file to open automatically gis put the info in X dimension and Y dimension but I have the error 000237 (One or more dimensions are invalid). Even if I change the variable for the x and y dimensions I have still the same error.

I tried also to drop the file directly into ArcMap and in that case I have just a point (and a message saying that there is no coordinate system so I added the coordinate system, WGS 1984, but nothing changes).

I found this question that I found helpful:  Importing NetCDF file to ArcGIS Desktop as point feature layer?
but my problem is that I can't create the file with the coordinate.

How can I solve my problem?


texts: Producing offset after georeferencing Soviet topographical military maps from printed control point coordinates in QGIS

texts:
I have 'acquired' a collection of scans of (previously top-secret) Soviet-era military topographical maps at a set of 1:25.000. I want to georeference a set of these scans.

I know the coordinate system and projection (Pulkovo 1942, Gauss-Kruger Zone 8, EPSG:28408) and am selecting this when pulling the scans into QGIS's Georeferencer Plugin.

To do the georeferencing, I'm using the coordinates printed at each corner of the map, which has been standardised at 7.5' divisions of latitude and 5' divisions of longitude.

I'm using WGS84 as the output coordinate system. However, when I then load the resulting GeoTIFF, it produces a very small latitude offset (approx 2") when compared to all other data I can find, including my own GPS traces and GCPs.

Within the area I'm interested in, the offset seems to be the same for every map, but I'm wondering if the issue is with my georeferencing process or if the maps themselves were all produced with this error?


texts: Gdalwarp producing odd results around 0 degree longitude

texts:
I am warping ECMWF grib data (full globe, 0.25 degree resolution) to match the SEVIRI geostationary satellite projection. I use this command:

gdalwarp -overwrite -ot Float32 -t_srs "+proj=geos +a=6378169 +b=6356583.8 +lon_0=0.0 +h=35785831" -te -5570248.832537 5570244.832537 5567248.429179 -5567252.429179 -ts 3712 3712 -wo NUM_THREADS=6 -wo SOURCE_EXTRA=1000 -q -r cubicspline --config CENTER_LONG 180 /data/ECMWF_T2_20181223T00_000.grib /data/cubicspline.tiff

The output projection looks good, and if I try this with the ECMWF land-sea mask then I can see that it is correct for most of the view.

However, along a vertical centred on zero degrees longitude I see this artefact.

Any ideas what is causing that, and how to resolve it? It artefact appears for all resampling methods except nearest neighbour.

(edit) In response to the comment, here's the original data - visualised using QGIS. Obviously there is some problem with how QGIS views ECMWF files, which start at 0 degrees longitude and go to 360 degrees.
Viewing the file in non-GDAL based viewers (such as Panoply) shows no problem.


texts: Converting Raster to Vector, Google Earth Engine,

texts:
By running the code, one obtains:

FeatureCollection (Error) Image.reduceToVectors: First band ('x') of image must be integral.

I don´t understand this error. My goal is to convert the pixel values into vector data, e.g. Points.

var coorsfeat = mean.reduceToVectors({
  reducer : null,
  geometry : geometry,
  crs :mean.projection().getInfo().crs,
  crsTransform :mean.projection().getInfo().transform,
  geometryType : 'centroid',
  labelProperty : 'Precip'
})

Link to GEE


texts: Create a color scale/legend for choropleth map using geotools (or other open source Java library)?

texts:
I'm experimenting with choropleth maps using geotools.

Is there an easy way to add a color scale to the map (using geotools itself or some other os library)? By a color scale I mean a bar that shows how the feature attribute values map to color, with labeled ticks either automatically generated, or at break points of my choosing.

Ideally, I'd like to generate a separate color scale that could be rendered independently of the choropleth map, so that I can use a single color with layouts containing multiple maps.

I'm not absolutely wedded to geotools, and would consider other open source libraries, though I'd much prefer to stick with JVM languages (Clojure, Groovy, Java, ...).


texts: Addig title to marker or shape using Leaflet Draw?

texts:
I am trying to create a participatory mapping site where users can draw pts, lines and polys using Leaflet draw.

Is it possible to create a pop window where users can add a title or description of the feature they draw -or- after they draw, click on the feature and add a label or description?


texts: Apply .qml style to a folder of tif images

texts:
I am trying to apply a .qml style from QGIS to each individual .tif image, contained within a folder, using python, but I can't get it to work. Is there a way to do this without loading the images into QGIS?

Here is my code which isn't working:

import glob
for layers in glob.glob('E:/Clipped/*.tif'):
    layers.loadNamedStyle('E:/style.qml')


texts: I asked a question before about whether it was possible to scale labels and I was correctly pointed to the IAnnotateLayerTransformationProperties.ReferenceScale property. What would like to know now is: does anyone have pointers on the correct algorithm to use to calculate that number (the reference scale) and the font size in order to best match the bounds of your layer and map? (I want to automate the labeling process as much as possible without asking the user questions that they probably can't answer (like me now.))


texts: Create a Shapefile datastore in GeoServer using REST API

texts:
I am trying to create a datastore for shapefile in GeoServer using REST API but I am unable to create a datastore for shapefile.

http://localhost:8080/geoserver/rest/workspaces/cite/datastores

<dataStore>
 <name>Test</name>
  <type>Shapefile</type>
  <enabled>true</enabled>
    <workspace>
    <name>cite</name>
  </workspace>
  <url>file:///home/ubuntu/mapfiles/test/shapefile.shp</url>
</dataStore>

ERROR: 500 Server Error

When I follow same procedure for GeoTIFF files everything works fine.


texts: How to create layer group in geoserver with gsconfig?

texts: How can I create layer groups using python gsconfig ?


texts: Exporting location coordinates as .csv file

texts: Is there a way to export the coordinates of regions that have been drawn on a map in google earth engine to a .csv file?


texts:
In the Python API of GEE, is it possible to delete a known asset ?
I can crawl in but don't know how to delete one :
assets = ee.data.listAssets({'parent': folder})

for asset in assets:
    # something like
    if asset.description = 'toto':
        asset.delete()


texts: Call an input layer in field calculator for another layer

texts:
My model start by asking two input layer LAYER1 and LAYER2.
Then, in field calculator for LAYER1, i want to join an attribute value from LAYER2. Using this code :
attribute( get_feature( 'LAYER2','ID1',"ID1"),'FIELD1')

The problem is that I'm not able to know the name of the input layer LAYER2
So I try to use the layer_property in order to see if the modeler recognize the name of yhe valeur but he won't:
layer_property( @LAYER2 ,'layerName')

But it always return a NULL value. Is there is something that you have to done to return the name of the layer?


texts: Using Zonal Histogram for each pixel in a raster by creating a buffer?

texts:
I was looking to use the Zonal Histogram for each pixel in a raster by creating a buffer (let's say of 500 meters) for a LULC map of 100 meter resolution. I'm using ArcGIS 10.6., Model builder.

One approach I tried was to convert raster to point, then  created buffer map using the point shapefile. Then I use the iterate feature selection on buffer map and gave it as input to Zonal Histogram along with original LULC map. The problem starts here, I don't get either list of tables. I tried to merge them in model itself but still it doesn't work.


texts: Is it possible to serve vector tiles directly from S3?

texts:
I have seen that it is possible to use S3 as a storage backend for vector tile servers such as Tegola. These work by fetching the tile from S3, then serving it through an application server.

Is there any possible mechanism by which tiles would be fetched directly from S3 in the browser? Having a massively scalable infrastructure like S3 piped through a tiny server seems like a bottleneck to be avoided.


texts: Sophisticated generation and calculation of Vector entities from Raster Data?

texts:
The first image below is shown a raster wish represents the altitude value. I created this raster by applying this process, bu using this process we can manually select an upper and lower area (so 2 vector entities) under some constraints like:

1 The difference of altitude between the two entities is upper than 200 m
2 The maximum distance between the tow entities is lower than 1000 m
3 The area of each entities is :   1.5 ha <area< 2.5

The second image below is shown a manuell detection of the upper and lower area wish is too hard to perfom.

My purpose is to automatically generate vector entities where every proposition (two entities) is saved on one vector layer with the ID, area and the distance between two entities so if it's generate 1000 propositions (1000 upper and 1000 lower) my result will be 1000 shp.
I can control the form of area using the Miller index using this formula:

  C=(4π*A)/P^2

  where A is the area and P is the perimeter

  C is equal to 1 if the form is circular and 0.75 if the form is square

  So I can define a threshold to delete all proposition wish have a Miller index mower than 0.6 for example

Is there any way to solve this problem using GRASS commands or Python?


texts: How to compute latitude and longitude of a point at a distance from another point in some direction?

texts:
I have the coordinates of a point on Earth and I want to know the coordinates of a point that is at a distance in meters from the first point in an angle in degrees from the north but I haven't found the way to do it.

Does anyone know how to do it?


texts: Clipping points with complex polygon geometry using GeoPandas

texts:
I have a huge point dataset with almost 300.000 points and I want to intersect them with a country geometry that has a complex geometry (multiple islands and borders).
What I've done so far is to construct a GeoPandas dataframe and then clip it, however the process takes almost 50 minutes whereas if I perform the same task with ArcGIS the entire process takes only 43 seconds. I've been searching for a way to improve the process without success.
gdf = gpd.GeoDataFrame(
    merged_data,
    crs='EPSG:4326',
    geometry=[Point(xy) for xy in zip(merged_data.longitude, merged_data.latitude)])
poly = country_boundary.geometry.unary_union
points_clip = gdf[gdf.geometry.intersects(poly)]


texts:
I'm trying to show a label if a field + another field = a certain value. Not sure the proper way to do this. I've tried cint([field1]) + cint([field2]) = x  but that doesn't work either.  It says the expression is valid but the feature that meets the criteria is not labelled.

Function FindLabel ( [SQFT_CONC] , [DEPTH_CONC] , [SQFT_ASPHT] , [DEPTH_ASHT]  )
  Dim strInput
  if [DEPTH_CONC] + [DEPTH_ASHT] = 10 then
     strInput = "IT WORKS"
  end if

  FindLabel = strInput
End Function


texts:
Background:
I am dynamically populating a feature layer. On this feature layer I have set the symbol (as you can see) to a nice little .png image. Everything seems to work fine around the layer except there is a blue cross box around my selection. I have tried setting the feature layers selection symbol to no avail. After reading through the api docs I don't see what else it could be.

Can someone please tell me:

what this symbol is
where it is coming from
and possibly how to remove it

Here is my original code that isn't working as expected:

    var featureLayer = new esri.layers.FeatureLayer(featureCollection, {
    mode: esri.layers.FeatureLayer.MODE_ONDEMAND,
    infoTemplate: template,
    outFields: ["Name"]
});

var infoSymbol = new esri.symbol.PictureMarkerSymbol({
    "angle": 0,
    "xoffset": 0,
    "yoffset": 12,
    "type": "esriPMS",
    "url": "/Images/map-marker.png",
    "contentType": "image/png",
    "width": 24,
    "height": 24
});

var renderer = new esri.renderer.SimpleRenderer(infoSymbol);
featureLayer.setRenderer(renderer);

var selectionSymbol = new SimpleMarkerSymbol(SimpleMarkerSymbol.STYLE_SQUARE,
    10,
    new SimpleLineSymbol(SimpleLineSymbol.STYLE_DASHDOT, new Color([255, 0, 0]), 2),
    new Color([255, 255, 0, 0.5])
    );
featureLayer.setSelectionSymbol(selectionSymbol);


texts:
How can I export a raster file out of QGIS (version 2.4) into a PostGIS database?

I am looking for a module or python script.


texts: gdal2tiles generates small picture inside 256x256 tile

texts:
I've used gdal2tiles for generating tiles in number of different polar projections and always experienced the same problem — strange small picture inside normal 256x256 png with the rest of space transparent.
I'm using -p raster option.

That's how one particular tile looks like http://grab.by/uNHw

And here is the whole map view.

Looking for suggestions. Gladly will switch to other utility/program that can generate tiles properly in non-mercator projections.

Source geotiff I use are arctic and antarctic maps:

http://www.ngdc.noaa.gov/mgg/bathymetry/arctic/ibcaoversion3.html
http://www.ibcso.org/data.html

Sample command I use: gdal2tiles.py -z 3 -p raster ibcso_background_hq.tif


texts: Does a geotagged image contain information about its accuracy?

texts: When I view my current location on Google Maps on my Android it shows me a circle of uncertainty. If I take a geotagged picture it will off course contain information about its Lat or Lon. Does it also contain information concerning accuracy of this geotagging?


texts:
Well, I'm a Windows user and a newbie with Web mapping technology, but I've had success with qgis server and the "getmap" and "getcapabilities" work well.

Now I'm trying to have the same success with Qgis Web Client, but appeared a little problem: The metadata, legend and layers appear, but the map is not showing up (instead appears a red broken layer). What could be the problem?

I don't know, but maybe I could have done something wrong in files like qgis-web-client.conf, or globaloptions.js, so i'll show you screenshots of it.

Qgis-web-client.conf:

GlobalOptions.Js:

http://jsfiddle.net/nYfLZ/

Or, do I need to edit GISProjectListing.js ?

http://jsfiddle.net/rWV7v/

GISProjectListing.js Screenshot:

Or it could be something else?


texts:
How do I add in multiple fields in the attribute table?  In the attribute table I have crop types such as corn, cotton, sorghum....etc (as shown image below), in square meters. I want to generate new fields so I can calculate into acres.  I am using python that would have add fields, but I would still have to provide field names, one by one (and there are a lot of field names to enter).  Can I just have GIS to generate fields that are based on the existing fields without me listing them?

Image of table is attached.


texts:
I am mapping the municipalities 911 numbers on our road network and built a simple python expression to find the lowest value between 2 fields.  Then it shows that value on the line followed by an arrow indicating direction of numbers ascending.

Spot checks revealed random segments that did not present the correct value.

If anyone can spot the problem, please let me know.

def FindLabel ( [L_FIRST], [R_FIRST] ):
    if [L_FIRST] < [R_FIRST]:
        return [L_FIRST] + '>'
    elif  [R_FIRST] < [L_FIRST]:
        return [R_FIRST] + '>'


texts:
My task is to upload geojson, add data from corresponding csv, drop some columns, then save it back to geojson.

However, after I did everything, added csv data (gj.join(df) ), etc,
I can't save it to geojson getting this Keyerror:

KeyError Traceback (most recent call last)
in ()
----> 1 x = newGj.to_json()
2 x

/usr/local/lib/python2.7/site-packages/geopandas/geodataframe.pyc in to_json(self, na, *kwargs)
249 return json.dumps(
250 {'type': 'FeatureCollection',
--> 251 'features': [feature(i, row) for i, row in self.iterrows()]},
252 *kwargs )
253

/usr/local/lib/python2.7/site-packages/geopandas/geodataframe.pyc in feature(i, row)
245 'properties':
246 dict((k, v) for k, v in iteritems(row) if k != self._geometry_column_name),
--> 247 'geometry': mapping(row[self._geometry_column_name]) }
248
249 return json.dumps(

KeyError: None


texts:
I seem to be having a problem previewing a layer group in GeoServer.  I have three shapefiles that I converted from a DXF file and I can preview each successfully, however when I try to preview the layer group nothing appears apart from the OpenLayers border and nav controls.

Any ideas?


texts: Does SVG offer better performance than plain HTML to render text labels in Leaflet?

texts:
I'm using Leaflet built-in L.DivIcon to add text labels to my map. Text labels are rendered in simple HTML divs and can be rotated and scaled using CSS based on values obtained from geojson data properties.

I would like to know if using SVG (with for example d3.js) to render my labels would allow faster panning and zooming than my actual HTML/div implementation ?

Thanks for any suggestions.


texts: Setting height and width to TileWMS layer in OpenLayers 3

texts:
I'm trying to set up a TileWMS layer in OpenLayers 3 as:

var baseMap = new ol.layer.Tile({
    title: "MayBase Map",
    source: new ol.source.TileWMS({
        url: mapUrl,
        params: {
            'LAYERS': 'GIS_Demo',
            'TILED':     true,
            'SRS':'EPSG:900913',
            'WIDTH':256,
            'HEIGHT':256},
        serverType: 'geoserver'
     })
 });

But when I actually request I get the URL like:

http://localhost:8080/geoserver/gwc/service/wms?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&FORMAT=image/png&TRANSPARENT=true&LAYERS=GIS_Demo&TILED=true&SRS=EPSG:900913&WIDTH=205&HEIGHT=205&CRS=EPSG:3857&STYLES=&FORMAT_OPTIONS=dpi:81&BBOX=8624542.775473006,1452915.0336446315,8629434.745283257,1457807.003454883

Here even though I am setting the height and width of 256 each its somehow overriding the default values of 205 (for more view docs).

This is basically causing the GWC error as : 400: The requested tile dimensions 205x205 do not match those of the grid set (256x256)

Why is this happening? is this a bug? or is there any other way of specifying the width and height?


texts: Select by location (multiple shapefiles in ArcGIS ModelBuilder)

texts:
I have looked around the forums and seen some helpful questions and answers around model building, but with ModelBuilder being relatively new to me, I have struggled to fully understand everything that has been said.

Where I work we have a process of using OS Codepoint data to form our boundary area. So for example the North East of England.
Usually this is a manual approach we do a few times a year, but I decided it would be a lot easier if we just ran a model to do this.

So we have a polygon feature called 'NE_Area' which contains codepoint polygons for the North East area.
I have created a model that uses the NE_Area shapefile, selects all records and deletes them.

Where I have got stuck is when I need to use a selection by location tool to select multiple codepoint polygons (all of which make up the NE area, TS, DL etc.), which intersect a polygon shapefile we have which covers the North East  Postcode area. The query shown below is an example for the Stockton Borough polygon we create:

Once all the relevant code point polygons are selected I copy and paste them into the NE_Area shapefile.

I was hoping in the model I could configure it to use all of the individual codepoint polygons, select all which intersect with the North East postcode area polygon, and finally copy and paste the selected records into the North East polygon.
Can I do this as a single tool in the model? I have been playing around in model builder and I don't know if I can... the only work I can think of doing it is if I use a select tool on each individual postcode area, export it out to a shapefile, then eventually merge them all into one.


texts:
I want to add a digital elevation model (DEM) layer to my leaflet map and access the height data through javascript, i.e. I cannot use .png files or qgis2webplugin but rather .tif or some compressed webmaptile file.

Any efficient idea or plugin recommendation?


texts: Scale dependent Pop-up, Is Arcade the answer?

texts:
Is there a way (possibly through an Arcade Expression) to control ArcGIS Online Pop-up visibility based on scale?
For example, I only want the Pop-up to be available when the user is zoomed in to a certain extent.
I can't create a copy of the data and change visibility because it will impact the filters I have set up.


texts: Matching X and Y axis coordinates of two GeoPandas geometries to overlay them

texts:
I have two GeoPandas GeoDataFrame objects with their corresponding geometry columns.

GeoDataFrame 1: counties --> has polygons and multipolygons of
Dominican Republic counties. The source data was obtained from the Dominican Republic 2010 census website and consist of shape files. It was apparently made using ArcGis.
GeoDataFrame 2: hospitals --> contains
longitude and latitude coordinates (Points) for medical centers in
the Dominican Republic. The source data was obtained from the Dominican Health Ministry website and consist of a Microsoft Excel file. The coordinates appears to align well with Google Maps.

I am trying to overlay the hospitals (Points) on top of the counties (Polygons). The problem I have is that the two geometries don't match. The hospitals are in longitude/latitude coordinates, but I'm not sure about the counties coordinates.

How can I make the two coordinates (X and Y axis) match so that I can overlay them?

The easiest way to understand my problem is to take a look at the plots:

https://drive.google.com/open?id=14xYmNJh958q1XANdq7dz98wn1Nrxymyi

Some sample code below. Note the total_bounds attribute.

import pandas as pd
import geopandas as gpd

# Import counties data
counties = gpd.read_file('ShapeFilesCenso2010/MUNCenso2010.shp')
counties.crs = {'init':'epsg:3395'}

# Import hospital coordinates
hospitals = pd.read_excel('hospitals.xlsx')

# After processing the excel file I turn it into the geometry below:
print(hospitals.geometry.head())

0    POINT (-69.89062 18.49381)
1    POINT (-69.88906 18.47107)
2    POINT (-69.89158 18.49421)
3    POINT (-69.92358 18.45247)
4    POINT (-69.86830 18.50938)
Name: geometry, dtype: geometry

print(hospitals.crs)

{'init': 'epsg:3395'}

print(hospitals.geometry.total_bounds)

[-72.      ,  17.      , -68.4513  ,  19.884642]

# The output from the counties geometry
print(counties.geometry.head())

0    POLYGON ((397122.651 2050791.921, 397142.013 2...
1    MULTIPOLYGON (((298172.939 2034259.341, 298374...
2    POLYGON ((333950.281 2046005.250, 334156.656 2...
3    MULTIPOLYGON (((300647.656 2050014.250, 300638...
4    POLYGON ((303931.923 2102011.319, 304050.231 2...
Name: geometry, dtype: geometry

print(counties.crs)

{'init': 'epsg:3395'}

print(counties.geometry.total_bounds)

[182215.7656, 1933511.9638,  571429.3273, 2205216.25]

ADDENDUM:

Checking an XML file that came with the counties (census) data, I see the following attribute and wonder if it has anything to do with the problem:

<Process Name="DefineProjection_2" ToolSource="C:\Program Files\ArcGIS\ArcToolbox\Toolboxes\Data Management Tools.tbx\DefineProjection" Date="20080220" Time="081209">
DefineProjection D:\DIVTER_ENERO2008\dm_1207_v1.shp PROJCS['WGS_1984_UTM_Zone_19N',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-69.0],PARAMETER['Scale_Factor',0.9996],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]] D:\DIVTER_ENERO2008\dm_1207_v1.shp
</Process>


texts: Where can I find a histogram of the grade of the Earth's surface?

texts:
I've found a histogram of the elevation of the Earth's surface on Wikipedia:

However this doesn't give any information about the distribution of grade. For example, the entire surface could be made entirely of little hills and have a high grade everywhere, or the surface could be composed entirely of perfect plateaus, putting the average grade at 0°. Obviously both these scenarios are untrue, but it illustrates how this information cannot be determined from the elevation histogram alone. Does anyone know where I can find a similar histogram for grade?


texts: Creating file geodatabases in script tool after specifying path and names as GetParameterAsText

texts:
I am trying to create a script tool in ArcGIS Pro that uses 2 clipping feature classes to clip multiple other feature classes and save them to two new geodatabases I create in the script.
I use 6 parameters as GetParameterAsText to specify the conditions:

sourceWorkspace (0)
ClipFC_North (1)
ClipFC_South (2)
out_folder_path (3)
out_gdb_N (4)
out_gdb_S (5)

Below that I try to create the new geodatabases:
arcpy.CreateFileGDB_management(out_folder_path, out_gdb_N)
arcpy.CreateFileGDB_management(out_folder_path, out_gdb_S)

Thats where my script fails and doesn't continue to clip the feature classes in the source geodatabase.
What am I doing wrong?
Here is my code up to the point it fails as a script tool.
import arcpy

out_folder_path = arcpy.GetParameterAsText(0)
out_gdb_N = arcpy.GetParameterAsText(1)  # North.gdb
out_gdb_S = arcpy.GetParameterAsText(2)  # South.gdb
sourceWorkspace = arcpy.GetParameterAsText(3)
clipFC_N = arcpy.GetParameterAsText(4)
clipFC_S = arcpy.GetParameterAsText(5)

print(out_folder_path, out_gdb_N, out_gdb_S)

# Execute CreateFileGDBs
arcpy.CreateFileGDB_management(out_folder_path, out_gdb_N)
arcpy.CreateFileGDB_management(out_folder_path, out_gdb_S)

Here is the error message I get:
Start Time: Friday, December 4, 2020 11:28:44 PM
Failed script ClipScript...Traceback (most recent call last):
  File "D:\Dan_GIS\Programming\DAN_Scripts\SH_36.py", line 21, in <module>
    arcpy.CreateFileGDB_management(out_folder_path, out_gdb_N)
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\management.py", line 21978, in CreateFileGDB
    raise e
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\management.py", line 21975, in CreateFileGDB
    retval = convertArcObjectToPythonObject(gp.CreateFileGDB_management(*gp_fixargs((out_folder_path, out_name, out_version), True)))
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\geoprocessing\_base.py", line 511, in <lambda>
    return lambda *args: val(*gp_fixargs(args, True))
arcgisscripting.ExecuteError: ERROR 999999: Something unexpected caused the tool to fail. Contact Esri Technical Support (http://esriurl.com/support) to Report a Bug, and refer to the error help for potential solutions or workarounds.
Failed to execute (CreateFileGDB).
Failed to execute (JasonClipScript).
Failed at Friday, December 4, 2020 11:28:45 PM (Elapsed Time: 0.69 seconds)


texts: What happens if you set QGIS as the default program to open a *.SHP file?

texts:
In the screenshots of shapefiles in this question, Accessing downloaded shapefiles/attribute data?, you can see that the .shp files have the QGIS logo as its icon.

I think this means that on that particular computer, someone tried to open a .shp file directly, Windows asked them what program to use, and they chose QGIS as the default program to open that type of file.

Now, that user is having some unusual difficulties getting shapefiles to properly load into QGIS. Polygon shapefiles are loading into QGIS as point layers with no data in them. I wonder if these problems be related to having QGIS set as the default program to open a .shp file.

I don't want to take the risk of testing this out myself, in case it causes the same problems on my computer. Hopefully someone out there has seen this done before.

So my question is: If you've ever set QGIS as the default program to open a *.SHP file, did it cause any problems?

Eg, shapefile corruption, causing QGIS to misinterpret the data type, etc.


texts: Update multiple leaflet locations if value

texts:
i have a vehicle tracking api that I am attempting to consume in leaflet. My idea is that is the vehicle driving value is true, then update the location. Ideally, it would be if the location is different, but starting small.

Below is what I have that works for one vehicle, but am unsure how to apply this to 150 or so vehicles

var ism;

function update_position() {
    $.getJSON('http://localhost:5000/live_events', function(data) {
        var latitude = myjson._items[0]['latitude'];
        var longitude = myjson._items[0]['longitude'];
    if (!ism) {
        ism = L.marker([latitude,longitude]).bindPopup("I am driving").addTo(map);
    }
    ism.setLatLng([latitude,longitude]).update();
    setTimeout(update_position, 1000);
});
}
if  (myjson._items[0]['isDriving'] = "True") {
  update_position();
}

So this would be updating every 5, 10, 15 seconds, whatever and loading under 200 vehicles. Does anyone know an optimized way to accomplish this and have it work for all vehicles?


texts: Referencing another layer in Field Calculator inside Graphical Modeler

texts:
I want to use the Field Calculator inside the Graphical Modeler to create a new column in an input vector layer (LAYER1) using data in another layer (LAYER2). LAYER2 is the output of another algorithm inside the model.

The formula is as follows:

attribute(get_feature("LAYER2", 'WKT', geom_to_wkt(start_point($geometry))) , 'NODEID')

Actually, I am trying to implement an algorithm similar to what's described in Getting points-touching-a-line attributes back to lines layer using QGIS?.

When used outside the model, it works perfectly.

When used inside the model, the column is always created with NULL values.

Any ideas on how to reference LAYER2 inside the Field Calculator for this to work?


texts: Field names change when selecting by attribute/location and then copying feature in ModelBuilder?

texts:
I am trying to create a tool to make my workflow faster for a routine job task.  Basically, I have a USGS 7.5 min quad layer and a feature (species) layer, and I want to ultimately create a shapefile (must be a shapefile, because I am sending these out to people who use them in a variety of GIS programs) of the species found in a certain quad and that contain the associated attribute data.  I am using ModelBuilder, and I have gotten it to go through the process of selecting the quad of interest (with Select by Attribute), use this quad selection to select the species that intersect that quad (with Select by Location), and export as a shapefile (with Copy Features).  The problem is that the resulting shapefile, while it contains all the attribute data, has lost its field names.  The field names become 'Request1', 'Request2', and so on (instead of 'Species Name', 'Common Name', etc. respectively).

I have tried using the 'Make Feature Layer' and 'Feature to Feature' tools, but without any change.  I have a Basic license, and so do not have access to the 'Alter Field Names' tool.

Any thoughts on what I'm missing or on a workaround?


texts: Linking to external files on AGOL in Windows 10

texts: I have a link field in a feature layer on AGOL. I want to link to files that reside on my desktop. When I click on the link nothing happens. If I copy the link to IE or Chrome, I get a search page list of options, not the file I want. This is a sample link that I am trying to use: "file:\aiotxtem3fp1\geodata$\elevation\Texas_Lidar_1m\Rasters\Texas_Lidar\USGS_NED_ArcherJack_1m_2014\USGS_NED_one_meter_x50y370_TX_Archer_Jack_2014_IMG_2016". The word "file" preceding the link is recommended by ESRI. I set this configuration up last July and this link worked perfectly: "\aiotxtem3fp1\geodata$\elevation\Texas_Lidar_1m\Rasters\Texas_Lidar\USGS_NED_ArcherJack_1m_2014\USGS_NED_one_meter_x50y370_TX_Archer_Jack_2014_IMG_2016" before the newest upgrade in AGOL and Windows 10.  Any ideas or suggestions as to the problem or a solution.  I don't have access to a server to put my files on.


texts: Why does function raster.Zscore in package spatialEco not work in R after installing?

texts:
I want to calculate the z-scores for a rasterstack. I downloaded the package spatialEco from https://cran.r-project.org/web/packages/spatialEco/index.html
Now it is not possible to write and run the function raster.Zscore.
Do somebody know this problem and can help?
That is the official description of the code:

Calculates the modified z-score for all cells in a raster
Usage
raster.Zscore(x, p.value = FALSE, file.name = NULL, ...)
Arguments x A raster class object
p.value
Return p-value rather than z-score raster (FALSE/TRUE) file.nameName
of raster written to disk...Additional arguments passed to writeRaster

If I run my code
raster.Zscore(rasterstack, p.value = FALSE)

I get following answer:
Fehler in raster.Zscore(rasterstack, p.value = FALSE) :
  konnte Funktion "raster.Zscore" nicht finden
means that--> function not found


texts: extracting the outline of raster imagery using GDAL

texts: Salams I want some help regarding developing a method using GDAL to extract and represent the footprint of any raster file (or group of raster data) in a vector format... any ideas please ??


texts:
Does anyone know how I can link to a location in Bing Maps and have it open in 3D mode?

A link formatted like this:

http://www.bing.com/maps/default.aspx?lvl=15&cp=37.9609257054894~-87.3786394334287

will open that location using their standard 2D map.  Is there anything I can add in to have it open in 3D instead?

Thanks,
Seth


texts:
(1) How to create a multiband raster from a set of single band rasters (except when somebody talks about doing a r.composite).

I attempted creating a raster3D and it turned out that it is not possible to display it with the GUI.

(2) How to subset a raster dataset (at least I do not see a direct command for doing so)

Edit (following the comment to my humble question): If my observations are wrong, please let me know the answer to these questions.


texts: ArcPy.da.UpdateCursor error The requested operation is invalid on a closed state

texts:
I am getting an error when trying to use the arcpy.da.UpdateCursor().
The error I am receiving is as followed

  The requested operation is invalid on a closed state

I am in the python windows and trying to work with the current document

upsql = "OBJECTID = 3456"
with arcpy.da.UpdateCursor("GravityMain",["FromManHole"],upsql) as upcurs:
     for updat in upcurs:
         if updat:
             updat[0] = '222'
             upcurs.updateRow(updat)
del upcurs

How can I use the UpdateCursor with the current document mxd = arcpy.mapping.MapDocument("CURRENT")

I looked at arcpy.da.Editor but it wants the full path to the database.

Also, I am using a MSSQL database.

If I input the following I still get the error

fc = "C:\Users\name\AppData\Roaming\ESRI\Desktop10.3\ArcCatalog\Sewer.sde"
fcc = "C:\Users\name\AppData\Roaming\ESRI\Desktop10.3\ArcCatalog\Sewer.sde\EC_Sewer.DBO.Sewer_Network\EC_Sewer.DBO.ssGravityMain"
workspace = os.path.dirname(fc)
edit = arcpy.da.Editor(workspace)
edit.startEditing(False, True)
edit.startOperation()
upsql = "OBJECTID = 345"
with arcpy.da.UpdateCursor(fcc,["FromManHole"],upsql) as upcurs:
   for updat in upcurs:
      if updat:
            updat[0] = '222'
            upcurs.updateRow(updat)
edit.stopOperation()
edit.stopEditing(True)
del upcurs


texts: Set geodatabase as outfolder in ArcGIS 10.7 when creating tool in Python toolbox

texts:
I am trying to specify a file location for an output file when creating a python toolbox tool, but the tool continues to output the file to the default geodatabase. I have tried to set the env.workspace within the tool, but clearly this is not the right approach. Here is a code snippet for the parameters I have created to specify the outfolder and the outfile name:
def getParameterInfo(self):
    """Define parameter definitions"""
    # First parameter. Input file
    inFile = arcpy.Parameter(
        displayName="Input File",
        name="inFile",
        datatype="GPFeatureLayer",
        parameterType="Required",
        direction="Input")
    inFile.filter.list = ["Point"] ## restrict to point files

    # Second parameter. outfolder
    outFolder = arcpy.Parameter(
        displayName="Output Folder",
        name="outFolder",
        datatype="DEWorkspace",
        parameterType="Required",
        direction="Input")

    # Input File Parameter
    outFile = arcpy.Parameter(
        displayName="Output File",
        name="outFile",
        datatype="DEFile",
        parameterType="Required",
        direction="Output")

Unfortunately, the tool continues to write to the default database.
How do I change this behavior?


texts: What is the easiest way to detect geometry type from text in PostGIS?

texts:
I have some records with different types of 'geometries (poly,line,points)' stored as string. I'm trying to use ST_GeomFromText as I want to do some intersection with them. What is the easiest way to detect whether I should put ST_GeomFromText(Polygon or LineString or MultiPolygon or Point in the SQL besides splitting the string and counting the coords pair?

Example data:

71.1776585052917 42.3902909739571

-71.1776585052917 42.3902909739571,-71.1776820268866 42.3903701743239,
-71.1776063012595 42.3903825660754,-71.1775826583081 42.3903033653531,-71.1776585052917 42.3902909739571

-71.1776585052917 42.3902909739571,-71.1776820268866 42.3903701743239,
-71.1776063012595 42.3903825660754,-71.1775826583081 42.3903033653531

So the first one should be point, the second a polygon, the third is a linestring.
But programmatically i don't know how to differentiate them, unless I split the string up count them up and if it's 1 pair it's point, etc...


texts: Solutions for uploading maps to GPS (forestry maps: roads, trails, hydro, polygons)?

texts:
I'm using a QGIS 2.2.0 Valmiera on a Windows 8.1 system Desktop. I'm a forester and I would like to have a cheap solution to load maps into the handheld Garmin 60, 62 or 64 series GPS. Free would be even better, if possible.

There is a software that costs around $350 that creates an image (.img) that is easily loaded into GPS, but it's expensive AND I don't think it works with QGIS. The best solution I've found so far would be Mapwel, which costs $45 or $90. You just take the shapefiles you need and merge them, then load them to your GPS.

Does anyone know of a cheaper solution?


texts: How to turn a point layer into classified proportional circles in QGIS 2.12?

texts:
I'm trying to style a point layer representing cities around the Mediterranean following these two rules :

The points should be represented as proportional circles, following a classification in 5 classes ;
The circles should be coloured according to the five classes, on a graduated scale from yellow to red.

I tried different things, including mixing classification and data defined styling, but didn't succeed yet.


texts: Merge all layers become a single layer in WMS

texts:
I am using ArcGIS Server now to create my WMS, but my WMS has many layers and I use another program (The Leica Spider Business), it can add the WMS in the basemap, but you just can choose only one layer for the WMS.
Is there any method to merge all layers become a single layer for WMS?


texts:
I am using Landsat data to determine NDVI for a specific region of interest (roi). However, the roi has buildings and other artificial objects that I wish to mask.
I have followed through the example at Masking NetCDF time series data from shapefile using Python.
However, this method masks all data outside of the shapefile, and retains the data from inside of the shapefile.
Has anybody any solutions for doing the reverse of this, and thus masking the data INSIDE the shapefile?


texts: How to make an arcpy tool to take multiple inputs csv/text lists, instead of a single input?

texts:
It's been only few days that i've been working in python environment.

Based on this code, i've been trying to batch few process for our QAQC task's. At certain times, like when i use Describe.(path_in_a_csv) i am always getting this error "Input value is not valid type". Even here, only with the help of Philip i was able to make the code run successfully. I tried the same thing arcpy.ExportMosaicDatasetPaths_management(path_in_a_csv,"Output"). The error was "not able to execute the tool".

Though i've been trying to understand this for the past couple of days, i couldn't able to crack it. It'll be of great help if anyone shed a bit of light on how to go from here.

import csv
import arcpy
from os import path
csvpath = r"D:\QAQC\Chk_Rep_File_list.csv"
with open(csvpath, "r") as csvfile:
   fcs = [r[0] for r in csv.reader(csvfile)]
for fc in fcs:
   desc = arcpy.Describe(fcs)

RuntimeError: Object: Describe input value is not valid type


texts: Filter image collection to one specific date (Earth Engine)

texts:
I'm using GEE and I have many codes were I filter my imagecollection by using filterdate (and then I have date range) and then I work with the first image in the imagecollection, when I know from the beginning which date I want .
I want to change my code so I can choose one specific date.

for example, if I want to get image from march 15th, my current code does this:

  var dataset = ee.ImageCollection('COPERNICUS/S2')
                    .filterDate('2020-03-14','2020-03-16')
                    // Pre-filter to get less cloudy granules.
                    .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
                    .select('B1','B2','B3','B4','B8','QA60')
                    .filterBounds(geometry)
                    .map(maskS2clouds);

  });

instead of filterDate with two dates, I would love to have something like that:

  var dataset = ee.ImageCollection('COPERNICUS/S2')
                    .filterDate('2020-03-15')
                    // Pre-filter to get less cloudy granules.
                    .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
                    .select('B1','B2','B3','B4','B8','QA60')
                    .filterBounds(geometry)
                    .map(maskS2clouds);

  });

But this doesn't work. I remembered in the past I mange to do this but can't find the script and also in the internet I find only date-range filter.

My end goal: to be able to insert one specific date.


texts:
I've been using GIS softwares in the past few years by I'm new using AutoCAD.

Well, I have contour lines in AutoCAD and I want to export it into QGIS (I also have access to ArcMAP).

I saved as .dxf (AutoCAD 2007/LT 2007 DXF) in AutoCAD. It opened fine in QGIS. However, it seems that didn't export the contour values (Take a look in the attribute table).

Am I missing some step?


texts: Extracting Cities, Towns and Villages boundary coordinates

texts:
How can I extract the polygon coordinates of cities,towns and villages from OpenStreetMap?

Are there any other free databases available?


texts: Translating geometry from EPSG:23030 to EPSG:21818, then to lat long, using PostGIS?

texts:
I'm from Colombia, South America. I'm a computer programming student and I have to develop a prototype for a GIS system for my small home town using Postgis >= 2.0 as database back-end.

The local public works office gave me the old Arcgis shape files: one for urban roads, one for rural roads, one for mountain streams, one for lakes, and so on. All shape files contain the *.shp, *.shx and *.dbf files, except the *.proj file.

This is the problem: The local goverment made these shape files long time ago using the EPSG:23030 reference system which locates in Western Europe (A very bad mistake with they still work and are not interested to change due to time and budget resources).

However, Colombia is located within the EPSG:21818 Colombia Mainland. Using http //epsg.io/, curiously, I figure that coordinates x,y (in meters since EPSG:23030 is UTM) of the roads, streams and lakes correspond to locations in Nigeria, Africa.

You can find a point located on the central park of my home town http //epsg.io/ using these coordinates: (-10153166.72, 1702714.02).

So my question is: How can I use Postgis, for example, using the st_transform function, to translate the imported shape files data from EPSG:23030 to EPSG:21818? How can I "convert" the EPSG:21818 to lat, lon coordinates later? Is it possible to do all this with Postgis?

I uploaded the shapefiles (no *.proj available), sql files (with EPSG:23030) and some pictures showing a sample point in UTM: 1106289.573,1131273.240 which I've taken as an example (it's an upper left corner on the central park).

Since, the original shape files didn't have the projection file, I set it in Postgis with this command:

SELECT UpdateGeometrySRID('viasurb02','geom',23030);

When using QGis we can see some rural roads distort due to low  bad quality when capturing the data, so I considered the urban roads layer (viasurb02.shp) for taking the sample point. Huh, That point is located in Benin, Africa.

www.mediafire.com/?s1s6ae8bh8i3kev

I've tried the suggestion about importing the shape files by setting the EPSG:23030 SRID:

 shp2pgsql -s 23030 -W LATIN1 -D -I viasurb02.shp viasurb02_b | psql -d mycity -U me

Then updated the SRID of the table:

ALTER TABLE viasurb02
ALTER COLUMN geom
TYPE Geometry(MultiLinestring, 21818)
USING ST_Transform(geom, 21818);

However, the traslated geometry is impossible for EPSG:4326 because it is inf/inf

Then, I've tried the suggestion assuming that the shapefiles already has the EPSG:21818:

shp2pgsql -s 23030 -W LATIN1 -D -I viasurb02.shp viasurb02_b | psql -d mycity -U me

And now, the same sample point is located in North Venezuela i58.tinypic.com/2jdny2x.jpg.

I don't know if the local public works office created in the late 1980s, using a bad projection or using Arcgis as it just it were Autocad.

Are there additional steps I can take to make another traslation geometry?


texts: Vector to Raster using GDAL (gdal.RasterizeLayer) : Error in Output Raster(generating output raster with NAN values)

texts:
For rasterizing a vector layer, I have tried  gdal.RasterizeLayer method for converting vector(shapefile) to raster(tiff). But its giving raster output with NAN values(complete black image). I need some help in which where I am doing wrong. And I want to print one attribute value from shapefile to raster & remaining values should be zero in outRaster.

My code ( using Gdal )is

    NoDataVal = -9999

    # Open the data source and read in the extent
    inPolygonShp = r"E:\polygons.shp"
    outputRaster=r"E:\OutRaster.tif"
    inGridSize=float(2)/110575      # nearly I am converting to 2 meters gridsize since its in GCS coordinates (vector layer)
    shpDS = ogr.Open(inPolygonShp)
    shpLayer = shpDS.GetLayer()

    # Create the destination data source
    xMin, xMax, yMin, yMax = shpLayer.GetExtent()
    xRes = int((xMax - xMin) / inGridSize)
    yRes = int((yMax - yMin) / inGridSize)
    rasterDS = gdal.GetDriverByName('GTiff').Create(outputRaster, xRes, yRes, 1,    gdal.GDT_Byte)

    # Define spatial reference
    rasterDS.SetProjection(shpLayer.GetSpatialRef().ExportToWkt())
    rasterDS.SetGeoTransform((xMin, inGridSize, 0, yMax, 0, -inGridSize))
    rBand = rasterDS.GetRasterBand(1)
    rBand.SetNoDataValue(NoDataVal)
    rBand.Fill(NoDataVal)

    # Rasterize
    err = gdal.RasterizeLayer(rasterDS, [1], shpLayer, burn_values=[200], options = ["ALL_TOUCHED=TRUE"])

    # for rasterizing with a attribute value of polygon
    # err = gdal.RasterizeLayer(rasterDS, [1], shpLayer, burn_values=[0], options = ["ATTRIBUTE= Height"])

Here I feel its generating a raster file but its failed to overwrite with shapefile values.
I need some help in this regard.


texts:
I seem to have put too much detailed information in my previous question, confusing the masses.  My apologies.  Let me try again from a different angle.

I want to rotate labels interactively/visually, using edit tools.  I don't want to use the Data defined Rotation method, found under Layer Properties > Placement.

Is there a way to interactively rotate labels in QGIS using any of the edit tools?

While Data defined Rotation works well, it is tedious to set the text rotations for a feature, having to manually enter and "guesstimate" the text angle value into the attribute fields.  Especially for point and polygon features.


texts: OpenLayers: how to save added point in server side

texts:
I'm new in OpenLayers.
I have simple application, which use a picture as map,
on this picture users can add, move, or delte point, and then save.
I have implemented adding, moving, deleting functionalities, only
except save.
codes are in below :

codes here

anyone who knows how to save change on server side? save to database,
file?


texts: Accessing the child of a group with owslib in QGIS

texts:
How can I iterate a layer group with owslib in QGIS in a script? If is not possible, are there alternatives?
Here is a sample code of what I'm trying to achieve:
2,5,8 has 1 as parent.
from owslib.wms import WebMapService
wms = WebMapService('https://neo.sci.gsfc.nasa.gov/wms/wms', version='1.1.1')
list(wms.contents)) #prints wms layers,not the groups also [['MOD14A1_M_FIRE', 'MOD14A1_E_FIRE'...]
# the above is listing only the children not the parents
# I want to iterate also the parent with 0->1->[2,5,8] within a script

img = wms.getmap(layers=['1'], #this is not working bcs a parent obviously
                 srs='EPSG:6343',
                 bbox=(708542, 460329, 731172, 493289),
                 size=(100, 100),
                 format='image/tiff')


texts: Migrating several shapefiles from one server to another?

texts:
my organization is a large ESRI user and we are under going some hardware changes. My question is:

Is there a way to move several shapefiles scattered across many folders (in one server) to another, completely new server, in an efficient manner without compromising the integrity of associated MXD files?

I'm aware that I can cut and paste the folders over but that would take ages as I'd have to reset each MXD file's layers  data source to the correct path. This is the main concern.

Perhaps there might be python solution however my scripting ability is weak.

Any ideas?


texts: Calculate total area of umnasked pixels of a raster image in Google Earth Engine

texts:
I have a partially masked image, and my goal it to compute an area of all unmasked pixels. The way how they describe process of area calculation in this example is by using a vector mask.

So do I have to vectorise all those unmasked pixels and then calculate an area of the raster within that vector, or is there an easier way?

var alos_dem = ee.Image("JAXA/ALOS/AW3D30_V1_1"),
    imageVisParam = {"opacity":1,"bands":["AVE_MSK"],"max":8,"gamma":1};

// Get the target layer
var alos_ave = alos_dem.select('AVE_MSK');
// Select only pixels with value 1
var alos_mask_cond = alos_ave.eq(1);
// Mask out all other values
var alos_mask = alos_mask_cond.updateMask(alos_mask_cond);

// calculate area in square meters here


texts: How to get landmass polygons for bounding box in Overpass API?

texts:
Can anyone show me how to get area polygons (just any landmass) with bounding box using Overpass?

I just want to differentiate between water and land - if you think it's better to assume everywhere is dry land and query for polygons of water to specifically render them instead, then please demonstrate this if possible.

EDIT: To clarify, I can't do massive queries (I'm splitting them up into bounding boxes), so I'm really pleased with the relation"type"="boundary"; solution, but I just need to tailor that query to filter out or ignore bodies of water. There must be a query that just returns useful "land only" polygons for a specific bounding box.


texts: Displaying points on OpenLayers with correct CRS

texts:
I have previously used QGIS but not for a while.

I just installed 2.14 (Essen) and the OpenLayers plugin.

I want to display a basic csv of thousands of points (mostly in USA) each with a latitude and longitude (mostly near the East Coast, e.g.

40.749473   -73.994262

I set project to Pseudo Mercator.
I import the csv as a layer (delimited), set it correctly and identify longitude and latitude as X/Y (correctly).
I turn on OpenLayers OSM (tried with Google too)
All with the same CRs (Pseudo Mercator - EPSG 3587), on the fly transformations

But... my points still all end up in the sea off Africa (see second image in link below)... but they are admittedly distributed in the right places, just as if the USA was really tiny and floating in the sea (the big cluster in the first image is NYC, then one/some in Europe and a few in West USA).

http://imgur.com/a/4P0bj

I feel like I must be doing something obviously wrong here, any tips?


texts:
I am trying to apply basic authentication in GeoServer but somehow it is not reflecting while consuming the WMS service through OpenLayers.

I am accessing it thorough below urls:

http://admin:geoserver@localhost:8080/geoserver/localhost/wms?.. // psw as plain text (GetMap)
http://admin:Z2Vvc2VydmVy@localhost:8080/geoserver/localhost/wms?.. // psw as base64 encrypted. (GetMap)
http://admin:Z2Vvc2VydmVy@localhost:8080/geoserver/wms?request=GetCapabilities //Check for capabilities

Either way it is not working. Below are the screen shots from GeoServer :

Filter Chain :

default chain config (with basic authentication as well)

I created a Child role and created a child user with this role, Also I created service access rule as below

Please suggest if I am missing anything.


texts:
I want to capture requests from QGIS WMS. So I set up Tinyproxy and configured QGIS: Settings\Options\Network and checked Use proxy for web access, set up localhost for Host, 8888 for Port and HttpProxy for Proxy type.

Now I am able to trace logs

sudo tail -f /var/log/tinyproxy/tinyproxy.log

That's what I got in the log when loading new WMS layer taken from example:

CONNECT   May 23 16:26:21 [9237]: Connect (file descriptor 6): localhost [127.0.0.1]
CONNECT   May 23 16:26:21 [9237]: Request (file descriptor 6): GET http://mrdata.usgs.gov/services/ca?version=1.1.1&SERVICE=WMS&REQUEST=GetCapabilities HTTP/1.1
INFO      May 23 16:26:21 [9237]: No upstream proxy for mrdata.usgs.gov
CONNECT   May 23 16:26:21 [9237]: Established connection to host "mrdata.usgs.gov" using file descriptor 7.
INFO      May 23 16:26:22 [9237]: Closed connection between local client (fd:6) and remote client (fd:7)
CONNECT   May 23 16:26:22 [8787]: Connect (file descriptor 6): localhost [127.0.0.1]
CONNECT   May 23 16:26:22 [8787]: Request (file descriptor 6): CONNECT mrdata.usgs.gov:443 HTTP/1.1
INFO      May 23 16:26:22 [8787]: No upstream proxy for mrdata.usgs.gov
CONNECT   May 23 16:26:22 [8787]: Established connection to host "mrdata.usgs.gov" using file descriptor 7.
INFO      May 23 16:26:22 [8787]: Not sending client headers to remote machine
CONNECT   May 23 16:26:29 [8788]: Connect (file descriptor 6): localhost [127.0.0.1]
CONNECT   May 23 16:26:29 [8788]: Request (file descriptor 6): GET http://mrdata.usgs.gov/services/ca?version=1.1.1&SERVICE=WMS&REQUEST=GetCapabilities HTTP/1.1
INFO      May 23 16:26:29 [8788]: No upstream proxy for mrdata.usgs.gov
CONNECT   May 23 16:26:30 [8788]: Established connection to host "mrdata.usgs.gov" using file descriptor 7.
INFO      May 23 16:26:30 [8788]: Closed connection between local client (fd:6) and remote client (fd:7)
CONNECT   May 23 16:26:31 [8785]: Connect (file descriptor 6): localhost [127.0.0.1]
CONNECT   May 23 16:26:31 [8785]: Request (file descriptor 6): CONNECT mrdata.usgs.gov:443 HTTP/1.1
INFO      May 23 16:26:31 [8785]: No upstream proxy for mrdata.usgs.gov
CONNECT   May 23 16:26:32 [8785]: Established connection to host "mrdata.usgs.gov" using file descriptor 7.
INFO      May 23 16:26:32 [8785]: Not sending client headers to remote machine
CONNECT   May 23 16:26:32 [8793]: Connect (file descriptor 6): localhost [127.0.0.1]
CONNECT   May 23 16:26:32 [8793]: Request (file descriptor 6): CONNECT mrdata.usgs.gov:443 HTTP/1.1
INFO      May 23 16:26:32 [8793]: No upstream proxy for mrdata.usgs.gov
CONNECT   May 23 16:26:32 [8793]: Established connection to host "mrdata.usgs.gov" using file descriptor 7.
INFO      May 23 16:26:32 [8793]: Not sending client headers to remote machine
INFO      May 23 16:27:02 [8785]: Closed connection between local client (fd:6) and remote client (fd:7)

The layer is loaded and displayed. When I begin to move the map or zoom in the image is successfully reloaded, but that's all I have in the log:

CONNECT   May 23 16:29:32 [8789]: Connect (file descriptor 6): localhost [127.0.0.1]
CONNECT   May 23 16:29:32 [8789]: Request (file descriptor 6): CONNECT mrdata.usgs.gov:443 HTTP/1.1
INFO      May 23 16:29:32 [8789]: No upstream proxy for mrdata.usgs.gov
CONNECT   May 23 16:29:32 [8789]: Established connection to host "mrdata.usgs.gov" using file descriptor 7.
INFO      May 23 16:29:32 [8789]: Not sending client headers to remote machine

Why are there no GetMap requests in tinyproxy log? How can I trace all the requests?

QGIS 2.18.4


texts: Understanding ST_RasterToWorldCoord() in PostGIS?

texts:
@Sorin helped me putting raster into PostGIS at Changing table name of command using raster2pgsql?

Now, I want to output raster tile, then view it on QGIS. But the problem is the output of ST_RasterToWorldCoord(raster rast, integer x, integer y) makes me confused. I have 2 questions here:

According to the docs, the function returns (long,lat) of the upper left corner of raster. So what is the use of the last 2 params (x,y)?
I run the sql, and got the result of a raster (long,lat) = (563985, 2346615). I don't understand the result here. I think the long, lat values should be with in (-90 -> 90, -180 -> 180).
I output the raster to file, and view it on a map. But the coordinate of the raster seems to be wrong, so it cannot be displayed.

Here my query:
    SELECT st_rastertoworldcoord(rast,1,1)
    FROM landsat

I have a raster before putting in postGIS. I can display it in QGIS, and it interscts/overlap some polygons, which come from a SHP file.
However, after cutting it into tiles and putting it into postGIS then run ST_Intersects(rast, geom). The function return 'false' in all rows.

My raster was imported with SRID 4326

(long,lat) = (-180 -> 180, -90 -> 90).


texts: Is there any interoperability between GRASS and PCRaster?

texts: After having looked at PCraster river network it seems that GRASS makes better work. But PCRaster is the best for building models. So is it really possible to use python-grass scripts for delineation and them convert them to .map and work through PCRaster?


texts: arcmap envelope coordinates from webbrowser tool

texts:
So I have my own php file with textfields which I run in my add in. Now I want to get the variables from the textfields and automatically give the value of the textfields to my env.PutCoords(x,y,x,y)
How do I do this?


texts: Filter by attribute names from another layer in QGIS

texts:
I have some difficulties with a task I have.

Data: layer 1 (tmc_segments), layer 2(regional_roads(earth)). Layer one contains simplified geometry of the road network (just simple lines with data about movement). Layer two contains some roads from Layer one, but these roads have more detailed geometry. Objective: to get from layer 1 only those roads that are present on layer 2. Layer one contains two distinctive attributes:road_numbe and road name. Layer two contains one attribute Name, wich contains both number and name.

As you can see, road geometry is very diferent, so Select by location not working that well.


texts: Including Less Boundary Detail in Zoomed out Map of ArcGIS Desktop?

texts:
I have a large map in ArcGis with a great deal of boundary detail which I don't need.

My issue is when my map is zoomed out the outline of these extremely detailed boundaries block other important features on my map.

What is the best way to either remove this extra detail or make it so the boundaries aren't drawn with such detail?


texts: Selecting multiple fields in QGIS Modeler

texts:
I am trying to make a model in the Processing Modeler of QGIS.
I need a step where a layer with polygons is dissolved based on two fields.
At the moment I do this manually. When I open the dissolve algorithm i can select the two fields.
In the modeler version I get a drop down menu or the option to enter the fields. If I try one field, all works well, If I try two fields, the algorithm dissolves my polygon layer as if I didn't enter any field to dissolve on.
I am almost sure there is something wrong in the notation of how I give in the two field. It is not mentioned in the help document. I must have tried any combination I can think of, but not the correct one.
Does anybody have the answer?


texts: Questions about SRTM , .hgt format and the sea level/coastline

texts:
What a nice place here ! The questions/answers system is very elegant.
I have read and get some useful informations allowing me to write a little .hgt reader...
I am working on C.G.I and create 3d images and visual FX. Then, writing this script bring me to be aware of my huge lack of knowledge about geographic science.

During my tests, I noted a difference between the terrain made from a specific hgt tile and the actual satellite photo ( googlemap ref) of the same area.

You can watch my terrain generation result

tile :N37E028

The blue part is MY water level ( 0 meter in the hgt file).
But if you watch the same zone in googlemap, if the overall look is good, you note big difference in the coastline.

Why that difference? Is the sea level above 0.0?


texts: Choosing specific ArcSDE version to create replica geodatabase from?

texts:
I'm trying to automate the process of exporting a replica geodatabase using the geoprocessing tools and eventually a Python script.

I can't figure out how to create the replica from a specific version.

Is this possible?


texts: Reshaping automatic polygon edge based on reference layer using ArcGIS for Desktop?

texts:
i have two polygon layer, one for Region administrative"Reference layer", and the second is administrative border is governorates.

Governorates border is not identical with Region border. So there is any idea to perform reshape automatic Governorates border with Region border

some_text http://www91.zippyshare.com/scaled/yuf4eBnj/file.html


texts: How to automate exporting map for many fields in the same layer?

texts:
I need to make a few ~75 frame time-lapse animations where each frame is a map of site status in a given week. Each week's statuses are a different field in a single shapefile. I export a frame of each week to jpg files and put them together in a video editor. I've done this manually when there were fewer frames, but now it's taking too long to make each frame individually. I'm using ArcGIS 10.2.

Is there a way to automate the creation and export of these time series maps? I need to do the following many times:

Set the display symbology for one field at a time, one field for each week
(e.g. cell value of 0 = blue symbol, 1 = red symbol, 2 = no symbol)
Export the map to a jpg, with the file name for that image
corresponding to which week is displayed (e.g. 2013-01.jpg, 2013-02.jpg).
Preferably be able to adjust this to add more weeks on the end for when the shapefile gets updated every few months with new data


texts: how to create map package programatically using C#?

texts:
I'm using .NET SDK to create ArcMap toolbar with arcgis 10.1. how do i create map package .mpk programatically?

http://resources.arcgis.com/en/help/main/10.1/index.html#//006600000403000000


texts: OMI data projection not lining up with WGS84

texts:
I am using data from Aura's OMI of global NO2. I need to compare values of NO2 measurements from OMI with measurements from ground based stations, for which the pixel location needs to be matching with the ground based stations location. However, my data layers do not line up (see image).

The coordinate system of the shapefile (world) is WGS 84 with datum World Geodetic System 1984. The coordinate system of the .tif file is GCS_WGS_1984 with datum D_World_Geodetic_System_1984. Trying to reproject a layer results in error 000365 "Invalid geographic transformation".

Does anyone know how I can solve this?