texts:
I was working on a file for a few weeks, saving it on an external drive, in a folder, along with all layers and shapefiles. Today I went to open it and the file is gone. Not in recycle bin or anywhere that I can find. The associated layers are still there. Any idea what could've happened? Luckily this was for school and I've exported most of what I needed into Illustrator, but I'd like to go back in and tweak a few things.
If I got to "open recent", I can see the path of where it was, but it has a little exclamation in a triangle, which I assume means it cannot locate the file.


texts: I am not quite sure how to ask this question...When editing data, is there a way have a File Geodatabase and say, all US states, major cities, and major roads.  When say, New York is chosen as a state, then only the major cities and major roads in New York are available. This just an example, but I guess I am trying to validate data, so when someone is editing they can't select options that aren't in their region.  Using the example above, if New York is chosen, Miami wouldn't be option.


texts: Update table based on closest point between its geometry(Point) and another table's geometry (MultiLinestring)

texts:
I think my question is rather simple, but being new to PostGIS, I can't figure it out myself nor being able to transfer knowledge of other questions to my issue.

The problem: I have two tables (testme, roads) with points as geometries (b.geom2 intestme --> point; a.geom in roads --> multistring)  in my PostgreSQL-database. Now I am trying to find for every point in testme the closest point in roads.

Alter TABLE testme Add Column matched_geom geometry(Point,4326);
Select a.geom, b.geom2,  ST_Closestpoint(b.geom2, a.geom) as the_geom
from roads as a, testme as b

The result of the query should be to append the closest point found (for every point/row in testme according to roads) in a new column called matched_geom in the table testme...

I executed the following query, which doesn't deliver any result.

Alter TABLE testme Add Column matched_geom geometry(Point,4326);
Select a.geom, b.geom2,  ST_Closestpoint(b.geom2, a.geom) as the_geom
from roads as a, testme as b;
WITH closest_points(id, geom) AS
   (SELECT id, ST_ClosestPoint(a.geom, b.geom2) as closest_geom
    FROM roads a, testme b
    WHERE ST_DWithin(a.geom, b.geom2, 100))
UPDATE testme SET matched_geom = closest_geom
FROM closest_points
WHERE testme.id = closest_points.id;

What is going wrong? The table testme contains about 150 rows, can it take that long to determine the closest points or is there still a mistake in the query?


texts: How to export attributes from QGIS to a .txtf file using Python

texts:
I have used this method outlined in:
https://howtoinqgis.wordpress.com/2016/10/31/how-to-export-attributes-from-qgis-to-a-txt-file/

##input_file=vector

from qgis.core import *
import processing
# Load the layer
layer = processing.getObject(input_file)

# Set the path for the output file
output_file = open('C:/.../cities.txt', 'w')

# Get the features and properly rewrite them as lines
for feat in layer.getFeatures():
    msgout = '%s, %s\n' % (feat["name"], feat["pop_max"])
    unicode_message = msgout.encode('utf-8')
    output_file.write(unicode_message)
output_file.close()

It works with QGIS v2.18, but it throws an error in v3 because “processing.getObject” :

AttributeError: module ‘processing’ has no attribute ‘getObject’

Do you know if there is a replacement for GQIS v3?


texts: Calculating average slope of a watershed upstream of many points in ArcGIS (upstream basin slope)

texts:
I am trying to replicate a procedure that the USGS Streamstats performs on points delineated along a stream network. They call it the "upstream basin slope" Basically, for a given point that I have confirmed is placed on a stream network, I want to calculate the mean slope of that point's watershed, including hillslopes and channels.

I have a shapefile of around 100 points, all snapped to a stream network (using snap pour points, then manually correcting errors). Many of the points are close to one another (i.e., many points on the same river). I can use the watershed tool to delineate watersheds for these points, but I can't just as easily use something like zonal statistics to calculate the mean slope in each of those watersheds, since the watersheds don't overlap (i.e., the upstream-most points would be accurate, but the rest would have a slope representative of only the portion of their watershed not taken up by the watershed(s) of the point(s) upstream.

Does anyone know of a way to accomplish this task without manually merging together accurate polygon watersheds for every point? I have a lot of points, and it would be extremely time-consuming to perform a procedure point-by-point.


texts: Geoserver 2.2.5 and Integrated GeoWebcache caching Defaults not saved over restarts

texts:
I am experience an annoying issue with Geoserver 2.2.5 and the integrated GeoWebCache running as a service on Windows. The Caching Defaults Configuration is not retained over GeoServer Restarts or Server Reboots.

This is particularly annoying in a production Environment as a server reboot or GeoServer Service Auto Restart on failure ends up with GeoServer running with no cache at all.

We have our Service set to restart the service on failure, then reboot the computer automatically on the third failure.

The documentation points to a configuration file in the Geoserver data directory named gwc-gs.xml that holds the configuration details. No matter what we do by editing this file, we cannot get the configuration to persist.

Interestingly, I noticed that a file named gwc-gs.xml.tmp was also in the same location. When I delete this .tmp file and then restart Geoserver, it appears again, always containing the default values

So my assumption is that this .tmp file is being created on Geoserver start up and it is being read by the GeoWebCache code, not the actual .xml file.

Question: Has anyone else experienced this or know of this behavior?

Question: Does anyone know of a work around or know if this behavior has been corrected in a later version of Geoserver?

Thanks in advance

Mark


texts: Best method to convert CSV file of lat/lon and depth data to bathymetry polygon?

texts:
I'm trying to create bathymetry polygons based on lat/long and depth data. I have tried several methods and the resulting polygons are clunky.

Make XY Event Layer-->Copy Features--> Project-->
Interpolate with IDW-->Make values Integeter with (Int tool) -->
Raster to Polygon
Dissolve Polygons by depth, so that all areas of equal depth are in same feature in the bathymetry feature class.

This method results as you might imagine is a clunky, blocky polygon.
I have tried several routes, including Contour, Create Filled Contour - I'm sure there must be a simpler way to create a bathymetry polygon?


texts: creating a random raster image using gdal

texts:
I want to create a random image containing say three classes with predefind percentages of each of the classes using gdal in python.

Can anyone provide any sample code?


texts: Is there an architecture for distributed geoprocessing?

texts:
Suppose I have 50 computers on my LAN.  Each computer has a geodatabase for all parcel polygons in a particular state in the US.

I'd like to write a geoprocessing task that finds all parcels valued over x $/acre that are within y feet of another parcel that is valued at less than z $/acre.

I'd like to formulate and run this query without knowing or caring that the data is distributed across 50 computers.  Keep in mind the boundary conditions: I also want the query to return cases where expensive parcels in one state are near inexpensive parcels in another.

Is there an architecture that supports this sort of distributed geoprocessing?

The architecture can be described abstractly, or as an implementation specific to Azure or Amazon Web Services.  Or, preferably, as a typical office where computers sit idle at night with plentiful ArcGIS desktop licenses.


texts: What does the ST prefix stand for, e.g. in ST_Polygonize?

texts: I was just curious, and couldn't find any reference to this anywhere - is it Spatial Types? Spatial Topography? Something else?


texts:
Is there a standard or a best practice way to store certain features about roads:

Road directionality:

One way streets: Roads going in one direction only
Center Turning lanes: roads that exist between two lanes, that either lane can use to turn left or right such as this
Highway interchange Ramps and merge lanes

Road Functionality:

Intersection layout (signalized or unsignalized):

Incorporating the different allowed and disallowed patterns of an intersection into a db, such as:
No-U turns
Turn left on green (or not)
green cycles and times
Actuated versus non-actuated approaches in intersections
Right of way rules: Some intersections might have different Right of Way rules than typical. Take for example a hilly road, where everyone must stop at the stop sign, except for the people who are traveling uphill. They do not have a stop sign, so all traffic must yield to them.


texts:
in the Program "Grapher" from Golden Software, you're writing in a modified script-language, based on BASIC.
I want to run the script I wrote in ArcGIS either using a toolbox or as an arcpy interface.

How do I run the BASIC script in ArcGIS?

I'm working with ArcMap 10.3.1

I found a solution to open "Surfer"(another app from Golden Software), but not Grapher

import win32com.client
def main():
    app = win32com.client.Dispatch("Surfer.Application")
    plot = app.Documents.Add(1)
    app.Visible = True
main()


texts: Using Union on multiple shapefiles to output single result

texts:
Using QGIS 2.2 and Python, I would like to use the Union algorithm on multiple shapefiles (contained in a folder TEST) which would then produce single shapefile containing all the data.

Currently the following code only outputs the union of the first and last file, how can I modify it so that the output of 2 layers is saved and then use that output with another layer etc?

UPDATE: The shapefiles are being processed in terms of the clip_function however I can't load the temp files when running the qgis:union algorithms. Also, an error occurs when it tries to output the Result.

##Model=name
##Select_folder=folder
##Cellsize=number 1000
##Result=output vector

import os
import glob
import tempfile
from qgis import utils
from qgis.core import QgsVectorLayer, QgsMapLayerRegistry, QgsVectorFileWriter, QgsCoordinateReferenceSystem

#   Create grid with defined parameters
outputs_1=processing.runalg("qgis:creategrid", Cellsize, Cellsize, 24108, 18351.157175, 258293.802316, 665638.226408, 1, 'EPSG:7405', None)

def clip_function():

#   Set directory and search for all .shp files
    os.chdir("C:\Users\TEST\DATA\\")
    for fname in glob.glob("*.shp"):

#   Clip .shp files within directory to grid and save files as temp
        outputs_2=processing.runalg("qgis:clip", outputs_1['SAVENAME'], fname, None)
        outputs_3=processing.runalg("qgis:fieldcalculator", outputs_2['OUTPUT'], 'Rank_' + fname, 1, 10, 0, True, 1 , None)

def binary_alg_wrapper(list_of_files,final_file_name):
    run_count=0
    temp_path=tempfile.gettempdir()
    temp_input_name="temp_input_layer"
    temp_output_name="temp_output_layer"
    temp_input_layer_filename=str(os.path.join(temp_path,temp_input_name+".shp"))
    temp_output_layer_filename=str(os.path.join(temp_path,temp_output_name+".shp"))
    for file_name in list_of_files:
        EPSG_Code_for_file=4326
        run_count+=1
        if run_count==1:
            #Copy the first file into a temp location
            temp_output_layer = QgsVectorLayer(file_name, "temp_output_layer", "ogr")
            Add_result = QgsMapLayerRegistry.instance().addMapLayer(temp_output_layer)
            write_error = QgsVectorFileWriter.writeAsVectorFormat(temp_output_layer,temp_output_layer_filename,"system",QgsCoordinateReferenceSystem(EPSG_Code_for_file), "ESRI Shapefile")
        QgsMapLayerRegistry.instance().removeMapLayer(temp_output_layer.id())
    if run_count>1:
        #Copy the result of the last algorithm into a temp location
        old_output_layer = QgsVectorLayer(temp_output_layer_filename, "input_layer",     "ogr")
        Add_result=QgsMapLayerRegistry.instance().addMapLayer(old_output_layer)
        write_error = QgsVectorFileWriter.writeAsVectorFormat(old_output_layer,temp_input_layer_filename, "system", QgsCoordinateReferenceSystem(EPSG_Code_for_file), "ESRI Shapefile")
        QgsMapLayerRegistry.instance().removeMapLayer(old_output_layer.id())

        #Read in a new file and the file stored in the temp location
        input1_layer = QgsVectorLayer(temp_input_layer_filename, "input1_layer", "ogr")
        input2_layer = QgsVectorLayer(file_name, "input2_layer", "ogr")

        #Run the binary algorithm and store the result into a temp location or output location
        if run_count==len(list_of_files):
            processing.runalg("qgis:union", input1_layer, input2_layer, temp_output_layer_filename)
        else:
            processing.runalg("qgis:union", input1_layer, input2_layer, final_file_name)

        QgsMapLayerRegistry.instance().removeMapLayer(input1_layer.id())
        QgsMapLayerRegistry.instance().removeMapLayer(input2_layer.id())

path_res = tempfile.gettempdir()
file_list = [shp for shp in glob.glob(path_res + "*.shp")]
Result=(os.path.join(path_res, "result.shp"))

if Select_folder:
    clip_function()
    binary_alg_wrapper(file_list, Result)


texts: Overlaying elevation data and bathymetry data in R?

texts:
I have two layers that cover approximately, but not exactly, the same area.
One is land elevation of some islands. The other is the depth of the water around the islands. I can resolve the vertical datums between the layers so that is not an issue.

I would like to create a new layer that consists of the land elevation and the depth of the surrounding water. The problems:

The extents aren't exactly the same
The resolutions are very different
The aspect ratios are different

Is there an easy way to create a new layer by looking at the higher resolution layer (land) and if the value is less than 75 say assume it's water then look to the other raster layer as to the depth of water otherwise accept the land elevation as the value in the new layer?

Land

class       : RasterLayer
dimensions  : 488, 1129, 550952  (nrow, ncol, ncell)
resolution  : 6.2e-05, 4.5e-05  (x, y)
extent      : -79.39998, -79.32998, 43.61302, 43.63498  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names       : dem_620483
values      : 74.89935, 87.08463  (min, max)

Water

> to_harbour
class       : RasterLayer
dimensions  : 26, 84, 2184  (nrow, ncol, ncell)
resolution  : 0.0008333333, 0.0008333333  (x, y)
extent      : -79.39958, -79.32958, 43.61292, 43.63458  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0
data source : in memory
names       : ontario_lld
values      : -15.9335, 5.764997  (min, max)

Original Question

I wish to combine elevation data in one raster layer with bathymetric data in another raster layer.

I think using raster::overlay is the correct approach here. However before I get that to work I require both layers to have the same extent.

Here is the code to isolate the part of land I am interested in.

file <- "data/GTA_elevation _data/GTA_DEM/6a7a7e71-f502-4336-bba2-364c7eefd950-SW/dem_620483/dem_620483.flt"

# Convert to class RasterLayer
img_raster <- raster(file)

# Transform RasterLayer to lat/long coordinate system
elevation <- projectRaster(img_raster, crs="+proj=longlat +datum=WGS84")

# Create bounding box
# Longitude (xmin, xmax), Latitude (ymin, ymax)
island_rect <- as(raster::extent(-79.4, -79.33, 43.613, 43.635), "SpatialPolygons")
proj4string(island_rect) <- "+proj=longlat +datum=WGS84"

# Crop just the islands
to_island <- crop(elevation, island_rect)

When I look at the extent it is not exactly what I entered:

> to_island
class       : RasterLayer
dimensions  : 488, 1129, 550952  (nrow, ncol, ncell)
resolution  : 6.2e-05, 4.5e-05  (x, y)
extent      : -79.39998, -79.32998, 43.61302, 43.63498  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names       : dem_620483
values      : 74.89935, 87.08463  (min, max)

This becomes an issue because when I use raster::overlay to combine land elevation and bathymetric data into one raster layer overlay gives an error because the extents are different: Error in compareRaster(x) : different extent

Assume that in the land layer anything lower than 75 is water. I'm not 100% sure how to use the overlay command.

to_island_combo <- raster::overlay(to_island, to_harbour, fun = function(x, y) {
  (x < 75) <- y
  return(x)
})


texts:
Objective : To know shortest path between two points.

Background : Created sample data ,converted into table in PostgreSQL ,in addition created topology,routing using pg routing extension. Data is created in proper format

Now trying to  display data in OL3 but getting an error "source not found". Tried to check error on internet but not getting exact result.

For reference used this link

Software used : Tomcat 8.0 ,PostgreSQL 9.6/PostGIS, Pgrouting 2.6 on Ubuntu 18.04

Please check attached screenshot for details

Code is available on this link


texts: Generating Points Shapefile From Excel File with X Y Values?

texts:
Using ArcGIS 10.2, I have an Excel file (Please review attached zipped Excel file) from which I am able to Generate Points by using the "Add XY Data". The points are generated on the map (using WGS84). I am then able to to export the points to Shapefile from Data --Export Data and then Add the new generated shapefile into the current document.

After this the problem starts:

I can not Select Any Point data from graphic Selection By Location, Clicks on table or Selection by Attributes and ArcMap crashes after any type of Selection happens.
After closing the map I am not able to add the new generated point shapefile into any map! It take several minutes and finally ArcMap display this error!

Can you please help me to fix this issue! I went all through to fix it but didn't work!

Download the Excel File


texts:
Occasionally I'm asked to plot points that relate to sensitive data items. I would quite like to be able to add a varying level of random jitter to points when plotting so that point density is still evident to viewers but without giving away precise locations.

I could do this quite easily by adding/subtracting random values from the co-ordinates in Excel or SQL prior to loading data but this would make it a bit difficult for when I want to do precise analysis based on actual co-ordinates as I would have to swap between real and jittered values (which is just asking for a mistake to happen further down the road).

Is there a way of doing this in QGIS? I've had a look through the documentation and haven't spotted anything so my last hope is that someone knows of a plug-in somewhere.

I know of other techniques to protect sensitive data on maps (binning, assigning to centroid, heatmaps, contours etc.) but none of these are quite right for the data I'm trying to display (individuals who have used a service multiple times, with point size related to number of times they've used that service).


texts: Passing GeoJSON API to Leaflet (proj4leaflet) map

texts:
I have created a geographical view of a point using Leaflet based on a coordinate system data. Since I'm using the coordinate system rather than lat and long, I'm using proj4leaflet for the conversion of data.
When I hard code the co ordinates within the JavaScript as a variable as shown in the below code. It is working as expected.
var map = L.map('map').setView([44.97,-93.24], 11);

// MapQuest OSM Tiles

// Attribution (https://gist.github.com/mourner/1804938)
var osmAttr = 'Map data &copy; OpenStreetMap contributors, CC-BY-SA',
    mqTilesAttr = 'Tiles &copy; MapQuest ';

L.tileLayer(
  'http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png',
  {
      attribution: 'Data by OpenStreetMap contributors'
  }
)
.addTo(map);

// GeoJSON layer (UTM15)
proj4.defs('EPSG:26915', '+proj=utm +zone=15 +ellps=GRS80 +datum=NAD83 +units=m +no_defs');

var geojson = {
  'type': 'Feature',
  'geometry': {
    'type': 'Point',
    'coordinates': [481650, 4980105],
  },
  'properties': {
    'name': 'University of Minnesota'
  },
  'crs': {
    'type': 'name',
    'properties': {
        'name': 'urn:ogc:def:crs:EPSG::26915'
      }
    }
  };

L.Proj.geoJson(geojson, {
  'pointToLayer': function(feature, latlng) {
    return L.marker(latlng).bindPopup(feature.properties.name);
  }
}).addTo(map);

But when I try to pass them by an API it is not even hitting my API.
I have tried using the following code
    .then(response => response.json())
    .then(data => {
        L.Proj.geoJson(data, {
            'pointToLayer': function(feature, latlng) {
                return L.marker(latlng).bindPopup(feature.properties.name);
            }
        }).addTo(map);
    })

it throws up an error that Access to fetch at 'https://localhost:44360/api/' from origin 'null' has been blocked by CORS policy: No 'Access-Control-Allow-Origin' header is present on the requested resource. If an opaque response serves your needs, set the request's mode to 'no-cors' to fetch the resource with CORS disabled.
Can someone suggest how to pass a GeoJSON to this Leaflet map by using an API?


texts: Is there a way to create a presence/absence matrix of polygons per cell in a grid using qGIS?

texts:
I'm working in QGIS with species distribution polygons and I want to know how many cells in a grid does the polygon overlap. I've tried creating random points within the polygon and then used vector<Analysis<points in polygon but I end up with gaps.

Can anyone help with this. I'm still learning GIS and I need to do this with over 1000 species.


texts: Do latitude/longitude coordinates ever refer to "real" earth, not idealized sphere or ellipsoid?

texts:
I'm new to the world of geographic information systems. While trying to learn the basics of spatial reference systems, I wonder whether there is any situation in which we can have latitude and longitude coordinates that refer to the "real" earth with its bumps and irregularities, and not an idealized sphere or ellipsoid. I think that such a situation doesn't exist, as all the mathematical methods used to arrive at latitude/longitude values (e.g. in astronomic navigation) need some geometric idealization to be tractable. As far as I understand, the same goes for GPS.

Is it correct that we never have latitude/longitude coordinates referring to the "real" earth, or am I mistaken?


texts: Error when pushing maps to baselayer picker in cesium?

texts:
I want to learn on how to create a baselayer picker and push maps to the baselayer picker in Cesium. The maps which I am going to push into the baselayer picker are accessed directly from the geoserver. I have referred to baselayerpicker sandcastle link, https://cesiumjs.org/Cesium/Apps/Sandcastle/index.html?src=Imagery%20Layers%20Manipulation.html&label=All

I have set up a local imagery server, geo server at port 8090 and I have added maps to the geoserver. I manage to access the map data from the geoserver using WebMapServiceImageryProvider

Here is my code for accessing the map from the geoserver in Cesium

Now I want to add the baselayer picker to allow user to access any of the maps accessed from the geoserver. I refer to the website, https://cesiumjs.org/Cesium/Build/Documentation/BaseLayerPicker.html for reference. After that, I implemented my code. I added two maps to the base layer picker.
my code is found below..

<!DOCTYPE html>

  <meta charset="utf-8">

  <meta http-equiv="X-UA-Compatible" content="IE=edge">

  <meta name="viewport" content="width=device-width, initial-scale=1, maximum-scale=1, minimum-scale=1, user-scalable=no">
  <title>Natural Earth II Bright</title>

the error states

Uncaught TypeError: Cannot read property 'defaultAlpha' of undefined
    at new U (Cesium.js:546)
    at c.addImageryProvider (Cesium.js:547)
    at l.set (Cesium.js:563)
    at l.n [as selectedImagery] (Cesium.js:560)
    at new l (Cesium.js:563)
    at new l (Cesium.js:564)
    at BaseLayerPicker.html:50

Can you explain what is the cause of the error and how to fix it?


texts: How to check if point has three or more lines snapped on to it in ArcGIS?

texts:
I inherited a shapefile with lines and points. I need to make selection of all points that have three or more lines snapped at same point. The problem is that the network is pretty large and it is impossible for me to do it manualy.

Is there a way to accomplish this in ArcMap? To make selection of points with 3 or more lines that are conneted to it.


texts: Getting attribute value from vector layer after mouse click in PyQGIS?

texts:
I am trying to write a piece of PyQGIS code that runs after a mouse click to collect the Map Coordinate position, and the value of a polygon vector layer within the project under the cursor position.

The vector layer is called 'Parish Boundaries'. The attribute I would like the value of is called 'name'. The first bit in the canvasReleaseEvent makes sure that the correct layer is selected, regardless of how its named.

I then have got very confused as to how to get the attribute value.

So far, I have the following:

class FaultTool(QgsMapTool):

def __init__(self, canvas):
    QgsMapTool.__init__(self, canvas)
    self.canvas = canvas

def canvasReleaseEvent(self, event):

    parish_layer = None
    layers = self.canvas.layers()
    for each_layer in layers:
        if each_layer.source() == "C:/MapData/QGIS/Base Maps/Other Layers/Parish Boundaries.shp":
            parish_layer = each_layer.name()

    x = event.pos().x()
    y = event.pos().y()
    parishName = # This is the bit I can't work out!
                 # something like....event.pos().attributes()[0]...?

    point = self.canvas.getCoordinateTransform().toMapCoordinates(x, y)


texts: How to zoom to selected feature with QGIS Python Console

texts:
I want to select a feature by ID and zoom to it. I tried the following but the canvas doesn't change. Any ideas?

canvas = qgis.utils.iface.mapCanvas()
cLayer = canvas.currentLayer()
cLayer.select(51)
canvas.zoomToSelected()


texts: Select Issues with Converting Feature Class to Shapefile with only certain fields

texts:
This is not a copy of previous questions like Selecting certain fields from feature class to create new feature class in ArcPy? but builds on them.

I have two approaches and two different issues.

Project: Take numerous feature classes in several SDE's and create a shapefile of each FC but only containing a selection of fields.

Method 1:

I copied the entire Feature Class to a new GDB, made a feature layer from that feature class and iterate through the layer to create a list of "delete fields" by excluding all of the already established "keep fields." I then create a new shapefile from the feature layer that deletes all fields I don't want. This works fine except when subtypes are discovered.

I wrote the following to handle subtypes:

subtypes = arcpy.da.ListSubtypes(tempLayer)
sKeys = subtypes.keys()
for keys in sKeys:
    if keys != 0:
        print("Subtype found at code {}.".format(keys))
        sys.stdout.flush()
        arcpy.RemoveSubtype_management(tempLayer,keys)

However, I still get errors saying fields cannot be deleted due to subtypes. When I look into/print the subtypes found in arcpy.da.ListSubtypes(tempLayer) ...no subtypes are found.
    So, with Method 1: What am I missing in identifying and deleting subtypes?


texts: Using ArcGIS ModelBuilder to reclassify by raster properties (mean)?

texts:
I'd like to reclassify raster using ModelBuilder based on the raster properties.

For example assign 0 to all values less than mean+2, and assign 1 to other values

I could do it manually by calculating the proper value, but for sure some of you know a better way to do it.


texts: Do you know a source where I can get a table for geocoding the numbers from German "Amtlicher Gemeindeschlüssel" (http://de.wikipedia.org/wiki/Amtlicher_Gemeindeschl%C3%BCssel)?


texts: Symbology based on related table ArcMap10.5

texts:
I have a parcels layer and related table of water meters. I did labeling of water meters IDs and I pulled it from the related table through the script.

Can I symbolize polygons based on related table records?

So far I was able to do that through joining tables, but it is not ideal for the end user they do not want those tables together.

I read it is possible in ArcGIS Online, but I am working in ArcMap 10.5.


texts: Which route will provide the best education and return on investment for a GIS career?

texts:
I currently have a B.S. in biology, and am interested in a career change (I am currently a high school science teacher.) I am trying to decide which route will provide the best education and return on investment: doing a post-bacc in GIS, getting a GIS certificate, or getting a master's in GIS since I already have a bachelor's degree. I'd love to hear any and all opinions!

Edit: To be honest, I wasn't more specific because I'm not really sure what my options are. My background would point me more in the direction of something environmentally/conservation related, but as I don't know much about the field, I would hesitate to pigeon-hole myself. I am trying to figure out the best way to find out as much as I can to not only know where my best fit would be, but to also make myself marketable.


texts: Leaflet markerClusterGroup 'undefined is not a function'

texts:
I have to create a map with more than 5000 points from a geojson. So I would like to use markercluster plugin with Leaflet.

I wrote the following code with help from the documentation and from other posts (I'm new to webmapping and to programming languages) and the chrome console returns that "L is not defined" for the leaflet.markercluster.js file and that "undefined is not a function" for my line of code with L.markerClusterGroup().

Here is the code :

var clusterGroup = new L.markerClusterGroup();

var points = L.geoJson(gares, {
    pointToLayer: function (feature, latlng) {
      return L.marker(latlng)
      .bindPopup(feature.properties.NOM);
    }
});

clusterGroup.addLayer(points);

map.addLayer(clusterGroup);

Is there someone to help me understanding why it doesn't work ?


texts: Plot image on map based on corner coordinates

texts: I have a color image and I know geographic coordinates of its 4 corners. Now I'd like to place this image on basemap in GMT using any projection. Any idea how to do that?


texts: Obtaining weight of evidence plugin for ArcMap

texts:
I read  this article:
An integrated user-friendly ArcMAP tool for bivariate statistical modelling in geoscience applications
but I did not find this tool "bivariate statistical modelling , BSM ".
I sent an email to its authors but they did not respond,
Is there any other way that I can obtain this tool?


texts: Finding the "perfect" track without using a road twice using QGIS

texts:
As a part of our work we have to scrutinize the network of roads in our forest for their timeliness. We were given the data (roads as lines and nodes) in the form of vector files. The problem with these files is that they are possibly outdated and due to that we have to countercheck them.

My idea was to create a track with some kind of analysis to lead us through each tile without having to drive some of the roads multiple times. Is there a a tool or some kind of extension to help me calculate the most efficient path?

So far I have only seen approaches to find the shortest track between two or multiple points, but this would not work for me since my prime intention is not the shortest path.

To make my problem clearer since I can not comment on my former posted question due to my lack of reputation (my registered account is fjaellraeven, Pascal was just a guest account):

At the moment I am using QGIS 2.14 Essen.

The factor I measure efficiency with is time. The situation is as follows:

The roads in the local forest are not up-to-date anymore and that is why we have to stocktake them. For this purpose I want to generate a track which is time efficient. From my point of view the most efficient path would be the one where I have to drive each road only once. The "Travelling Salesman" deals with this problem quite well.

We do not have to file the tracks we are driving so it is not about accuracy in the first place. If we encounter a missing/deviating road or one which is not accessible anymore we will digitize it later. The main intention is to check the existing roads for their timeliness.

I hope my "problem" is a little bit clearer by now. Keep in mind that the forest we have to drive along in has an extension of multiple thousand square kilometres making it kind of a huge deal if we drive inefficient and without a plan.


texts: Is it possible to highlight water steam sources on Landsat 7 image?

texts:
I am trying whether it is possible to view hot steam sources like power plant cooling towers and resulting cloud using satellite images. I was advised to try Landsat (needed free data for initial research)

I downloaded two Landsat 7 images from various regions with two big power plants (about 4 cooling towers each) and tried to use common band combinations (3-2-1, 4-3-2, 4-5-3, 7-4-2, 5-3-1) in SAGA GIS to create composite images. I though it could show the vapor cloud at least near the power plant but I am not able to identify it. It is possible to somehow filter the images and remove most of the background to highlight what I need?

thanks

PS: the goal is to get an approximate result whether it has sense to start more detailed research in this area or not. I have QGIS and SAGA on all of my computers.


texts: Calculate distance between buffer and geographic points for grouped subsets (and get tally)

texts:
The data I used can be found here.

After importing my data, I select a subset of points (Nest2017 in column type) to create 30 m buffers around:

library(dplyr)

nests<-df %>%
    filter(type %in% "Nest2017")

library(maptools)
library(sp)
library(rgeos)

nests$lon <- as.numeric(as.character(nests$lon))
nests$lat <- as.numeric(as.character(nests$lat))
coordinates(nests) <- c("lon", "lat")

proj4string(nests) <- CRS("+init=epsg:3578") #setting the projection of the points

spTransform(nests,CRS("+init=epsg:3578")) #tranforming the coordinates

gBuffer(nests, width=30, byid=TRUE) #puts a circular buffer around each individual point

I know gBuffer()requires the points in UTM, but my points go across two zones (7 and 8), so I used the Yukon Albers projection instead.

Now I am trying to figure out how to calculate the distance from the buffer around each Nest2017 (in column type) to each point (Foray) within the same NatalMidden group.

A subset of my data looks as follows:

lat                  lon                NatalMidden   squirrelID    type
60.9577819984406    -138.0347849708050  -27           NA            Nest2017
60.9574120212346    -138.0345689691600  -27           NA            NatalMidden
60.9578209742904    -138.0346520338210  -27           23054         Foray
60.9575380012393    -138.0348329991100  -27           23054         Foray
60.9576250053942    -138.0339069664480  -27           23054         Foray
60.957643026486     -138.0338829942050  -27           23054         Foray
60.9575670026243    -138.0348739866170  -27           23054         Foray
60.9600000176579    -138.032592013478   -515          22780         Foray
60.9600180387497    -138.032631995156   -515          22780         Foray
60.9599519893527    -138.032342987135   -515          NA            NatalMidden
60.959974033758     -138.032317003235   -515          NA            Nest2017

So, for example, squirrelID 23054 was located (Foray) multiple times (type column) and I have a corresponding latitude (lat) and longitude (lon) for each Foray. I am trying to calculate the distance between each Foray (type column) and the 30 m buffer (assuming I correctly created the buffers) around its corresponding Nest2017 (type column) for each individual (squirrelID) separately. The common field that links squirrelID and type is the NatalMidden column.

  Is there a way I could work within the dplyr framework to group_by(squirrelID) and then calculate the distances between (and tally) each Foray and its corresponding Nest2017 30 meter buffer (the common field being NatalMidden for both the Foray and Nest2017)? If not, how might I be able to do this?

My ultimate goal is to create new columns for:

distance between each Foray and its corresponding 30 meter buffer (centered on Nest2017) for each squirrelID
number of Forays inside corresponding 30 meter buffer (centered on Nest2017) for each squirrelID
number of Forays outside corresponding 30 meter buffer (centered on Nest2017) for each squirrelID


texts: How to solve bitwise error in Google Earth Engine?

texts:
I've wrote the below code for chlorophyll detection using MODIS ocean reflectance data but I've got this error:

Image.bitwiseXor: Bitwise operands must be integer only.

code link: https://code.earthengine.google.com/31dc6a60f330ee27c8e2a3e75eddee3c
Map.centerObject(table);
Map.addLayer(table);

var chl = modis
.filterDate('2020-05-01','2020-06-01')
.filterBounds(table)
.map(function(img){
  var bands = img.multiply(0.0001).clip(table);

  var R = bands.expression('log(b9 / b12)',
  {'b9': bands.select('sur_refl_b09'),'b12': bands.select('sur_refl_b12')});

  var index = bands.expression('10.0 ^ (0.283 + (-2.753 * R) + (1.457 * R^2) + (0.659 * R^3) + (-1.403 * R^4))'
  ,{'R':R});

  return index.copyProperties(img,['system:time_start','system:time_end']);
}).mean();

print(chl);
Map.addLayer(chl);


texts: Slim down large project with many layers to prevent crash in QGIS

texts:
I am afraid my project might stop opening at some point because of too many layers and some costly rendering options.
Is there a way to still keep the layers in the layer tree without loading (or even "searching" for) them when opening the project? I am looking for a way to open a project as if it only had, say, 10 layers and some kind of pseudo-layer-tree with some kind of phantom layers that QGIS will only care about when I select them for working on. Would it mayhaps be better to convert my project file from .qgz to .qgs?
The difference to this approach in contrast to the browser in QGIS is that layer settings and options are preserved. The difference in contrast to the current layer tree is that nothing is done with the project layers, not even their file path is checked. This approach is something between the browser state and the current layer tree.
Is my question understandable and if yes, what are the obvious disadvantages to this technique in your opinion?


texts: My SHP file does not appear on my base map

texts:
I am using QGIS 2.18.3 with a base map named Natural_Earth_quick_start_for_QGIS

I use Australian Albers EPSG: 3577

When creating a bunch of polygons on my base map, my QGIS suddenly stopped working. So I had to re-start QGIS and added vector layers that I was working on again. After retarting and changing the default coordinate to Australian Albers EPSG: 3577, almost all of the layers that I have just created previously do not show up, except for the layer named Cambodia Chamic1.

What seems to be the issue here? how to make them show up.

This is a list of the files associated with my created polygons

=========================================================================
=========================================================================

@ Joseph

I did check the Coordinate reference system of my layers at Properties > General > Coordinate reference system. All of my layers have CRS as WGS 84, EPSG: 4326 which is not the same as my chosen projection Australian Albers EPSG: 3577.

But the layer named Cambodia Chamic1 still shows up on my base map, whereas other layers having the same setup do not show up.

When I change the Coordinate reference system of all of my layers from WGS 84, EPSG: 4326 to Australian Albers EPSG: 3577, the layers that have not shown up still do not show up

and the only layer that's showing up (Cambodia Chamic1) suddenly disappears.

I don't know what's going on right now?


texts: How to transform camera calibration parameters to other coordinate system?

texts:
I have all camera calibration parameters (X0, Y0, focal length, R0, P0, P1, K0, K1, coarse rotation, rotation, translation) and position of the image in WGS84 (XYZ in meters) and other parameters like fi (azimuth), omega, kappa.

I have laser points in the coordinates system different than WGS84, so each time when I want to check if the point is in the image I need to transform my point into the WGS84 coordinate system.
To avoid unnecessary calculations I want to transform my position of the image to coordinate system like in my laser data.

The problem is that I also need to transform camera calibration parameters and fi, omega and kappa. I don't know which exactly values need to be transformed and how I can achieve it.

I'm not asking about the soft, but the method to do that because I need to implement this solution to my program. Also, I can't transform image position before calculate camera calibration parameters, because I will only have final values.

The question is: what parameters I need to transform and which mathematical formulas I need to use?

I've tried to figure it out myself, but unfortunately, I don't have enough knowledge about it.


texts: rasterize produces a raster full of NA values

texts:
I am having trouble converting a polygon to a raster. The desired raster is created, however all its values are equal to NA.

require(rgdal)
require(raster)

#download a shapefile with 2010 census data
tmp_dl <- tempfile()
download.file("http://files.hawaii.gov/dbedt/op/gis/data/blkgrp10.shp.zip", tmp_dl)
unzip(tmp_dl, exdir=tempdir())
HIshp <- readOGR(tempdir(), "blkgrp10")

#add a new attribute/field with population density of each polygon
p_areas <- sapply(HIshp@polygons, function(i)i@area)
p_pops <- HIshp@data$POP10
HIshp@data[,3] = p_pops/p_areas
names(HIshp@data)[3] <- "POPDENS"

#here is the raster format that I require for the result
new_ras <- raster(nrow = 1520, ncol = 2288)
extent(new_ras) <- c(-159.816, -154.668, 18.849, 22.269)
crs(new_ras) <- "+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs"

rasterize(HIshp, new_ras, field="POPDENS")

And here is the sad result:

> rasterize(HIshp, new_ras)
class       : RasterLayer
dimensions  : 1520, 2288, 3477760  (nrow, ncol, ncell)
resolution  : 0.00225, 0.00225  (x, y)
extent      : -159.816, -154.668, 18.849, 22.269  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=4 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
data source : in memory
names       : layer
values      : NA, NA  (min, max)

How can I make the values convert properly?


texts: Removing addresses / points within certain distance of each other and have the most points remaining?

texts:
I have a data set of nearly 2000 points. I would like to find out which points remain and the most possible if no points can be within x kilometers of another point. So, if one point was randomly chosen, then the next closest to be chosen would need to 1) be greater than x km away and 2) the next closest.

I would like to have the optimal result but it could also be good to be able to have the first point chosen at random and run that a few times to get an average result.


texts: Find Direction/Location of point in polygon using Python ArcGIS10.1

texts: I have a point located inside a polygon, I want to determine its location or direction i.e. if it is in the northwest corner or southwest. I have the x,y of the point. Any idea on how to do this? I have the shape file of the point and polygon. Any help on this is most welcome.


texts: Find most contrast pixel on raster in certain range of pixels around it?

texts:
How to assign to each pixel of a raster value of the biggest difference between values of all pixels around it to a mask on 1 picture (red pixels and everything that after them is not included, black pixels and everything that inside of them are included).

Here is an example of how it works on 2 picture. I wrote such pseudo-code with Python assuming that as input i get list with all closest pixels using as a Central pixel lower left a corner 2 pictures of the. I do not understand how this can be transformed into a script for QGIS or Raster Calculator.

I tried in R the neighbour's module, but every time I set a radius greater than or less than 3 it throws an error. Also, you can take pixels not only for this mask but also for a circular radius of 128 meters. (For images of Landsat with a spatial resolution of 30x30). With what should I implement this?

near_pixels_list = [4,
                    1, 3, 3, 5,
                    8, 4, 1, 3,
                    0, 0, 7, 7,
                    4, 3, 0, 2, 1]
center_pixel = 4
output = max(list(map(lambda x: max(center_pixel, x) - min(center_pixel, x), near_pixels_list)))


texts: Creating Annotations from point features using ArcPy?

texts:
Because in ArcGIS 10.0 there is still no simple toolbox function to convert labels to Annotation, I try to go a substitude way with the callable Toolbox function [TiledLabelsToAnnotation][1] on a polygon shape, which covers all of my data area (= 1 polygon).

The script is now running, because of the simplifications for the posting!

import arcpy

# Settings
scale = "30000"
folder = "D:\Temp\LBL.gdb"
arcpy.env.workspace = folder

# Mapdocument
map = r"D:\Temp\MAP-LBL.mxd"

# Local variables:
dataFrame = "Layer"
grid = "Grid"
group = "Group"
suffix = "Anno"
fid = "tile"
link = "STANDARD"  # not linked = "STANDARD" / linked = "FEATURE_LINKED"
unplace = "GENERATE_UNPLACED_ANNOTATION"

# Process: Tiled Labels To Annotation
arcpy.TiledLabelsToAnnotation_cartography(map, dataFrame, grid, folder, group, suffix, scale, "", fid, "", "", link, unplace)


texts: How to execute st_read function including multibyte character?

texts:
how to execute st_read function including multibyte character?

library(sf)

nc <- st_read(system.file("shape/nc.shp", package="sf"))
st_write(obj = nc, dsn = enc2utf8("データ"),layer = "nc", driver = "ESRI Shapefile")
st_read(dsn = enc2utf8("データ"), layer = "nc")

Cannot open data source C:\Users*****\nc\データ
Error in CPL_read_ogr(dsn, layer, as.character(options), quiet, type,  :
  Open failed.


texts: How to create a single depth contour line from a bathymetric raster?

texts: I have bathymetric data (raster form) that I'm trying to tease out specific depth contours.  For example, I would like to "create" the 150-m depth contour from the bathymetric data in the Gulf of Alaska.  By using the Raster Calculator, I was able to create a raster of just the 150 m depth cells.  However, I would like these to be connected as a "smooth" contour line.  I'm having a hard time thinking of how to go about this.  I hope this makes sense!


texts: How to Change Precision of Export / Add Geometry Columns field

texts:
I would like to display the coordinates of a number of points on a map.  I can use Vector > Export / Add Geometry Columns which adds coordinates to 6 dp.  However as I am using a UTM-type projection this level of (false) precision is not necessary.

I can copy the values to another predefined column in Field Calc, but would need to do this every time I updated or added a point.  Is there a way of changing the default precision value for Export / Add...  ?

Edit - I am using QGIS 2.0.1 Dufour


texts:
I am doing the classification of mangrove trees according to their species, height, and etc by using QGIS. I have done the species classification part. I would like to know how can I differentiate these trees by their height?

Let's say I have a few categories:

0-5m
6-10m
11-15m

How can I classify them using QGIS?

Is remote sensing able to classify them automatically or should I go to the site and collect the data and input as attributed data (e.g. CSV)?


texts: Arcpy/REST: how to save raster dataset from URL

texts:
I have a script which downloads features from an ArcGIS MapServer via the REST API. Getting vector features is relatively trivial since the features can be retrieved as JSON. However, when retrieving a map as an image, I'm getting stuck at the point where I have retrieved a response.

The relevant portion of my code, which gives me a valid image URL, is this:

# Export the map to JSON.
urlGetImage = '%s/export?bbox=%s&bboxSR=%s&format=png&f=json' % (
    theBasicUrl, bbox, bboxSR)
urlGetImageResponse = urllib2.urlopen(urlGetImage)
jsonImage = json.loads(urlGetImageResponse.read())

# Read the image URL from the JSON and save the image.
urllib.urlretrieve(jsonImage['href'], os.path.join(tempFolder, 'test.png'))

Note the parameter where I retrieve the response as JSON. This is not essential. I don't care what I retrieve it as, so long as I can retrieve the generated raster somehow.

A sample response is this:

{"width": 400, "href": "https://[SERVER NAME]:[PORT]/arcgis/rest/directories/arcgisoutput/imagery/DEM_MapServer/_ags_map26932bd6681a4546a951125dd7c60e88.png", "scale": 20326865.193372156, "extent": {"xmin": 2996004.971208999, "ymin": -3232131.5223958646, "ymax": -1080867.3202355732, "xmax": 5147269.17336929, "spatialReference": {"wkid": 102100, "latestWkid": 3857}}, "height": 400}

This gives me a path to download an image from, as well as the relevant spatial reference WKID. I then save it to an image.

My question is this: how do I go from here (an image on the drive) to a (referenced) raster dataset? Or have I already gone wrong in getting to this point?

EDIT:

To clarify, I need to work with a MapService, not an ImageService, as that's how the rasters will be published. When I export the result above as a KMZ (rather than JSON), and then load it in Google Earth, it appears in the correct location. So in theory, getting a raster image from the map server would seem possible. However, trying to convert this into a raster layer via arcpy.KMLToLayer_conversion gives:

  "ERROR 000401: No features were found for processing."

... and even if I extract the image information in the KML, download the image, and replace the path in the KML (which by default points to the generated image on the server), it still complains with this error.


texts: Creating colormap where colors diverge at value 0 in QGIS?

texts:
I calculated a raster as difference of two other rasters. Now I want to create a colormap where one colour represents negative values and one colour positive ones. The two colours should become lighter as they approach zero and around zero the colour should be white.

Here the colour scale I would like to use:

My problem in using the plugin "1 Band Raster Colour Table" is that it creates palettes where the colours diverge around some value above zero. In my case my values go from -6880 to 5539 with a mean of 104.


texts: How to load a GeoJSON (in a javascript variable ... ) in MapBox GL JS?

texts:
I need to load in MapBox GL JS a layer based on a GeoJSON and I started form this example

https://www.mapbox.com/mapbox-gl-js/example/geojson-line/

I've my data in a Javascript variable

   var theRouteJsonText = "{ \"type\": \"geojson\", \"data\": { \"type\": \"Feature\", \"properties\": {}, \"geometry\":";
   theRouteJsonText = theRouteJsonText + JSON.stringify(theCoords);
   theRouteJsonText = theRouteJsonText + "} }";

....  and its value is, for example, ...

{
    "type": "geojson",
    "data": {
        "type": "Feature",
        "properties": {},
        "geometry": {
            "coordinates": [
                [9.123777, 39.213546],
                [9.124838, 39.214753],
                [9.125036, 39.219553],
                [9.123548, 39.220278],
                [9.122217, 39.219429],
                [9.121313, 39.219682],
                [9.120233, 39.221339],
                [9.119572, 39.221669],
                [9.118786, 39.221757],
                [9.118277, 39.223214],
                [9.118174, 39.223181],
                [9.116678, 39.225676],
                [9.116099, 39.226343],
                [9.111633, 39.228077],
                [9.11111, 39.228516],
                [9.110992, 39.229156],
                [9.111139, 39.231757]
            ],
            "type": "LineString"
        }
    }
}

I'm trying to add the GeoJSON layer to the map in this manner ....

map.addSource("route", theRouteJsonText);

map.addLayer({
    "id": "route",
    "type": "line",
    "source": "route",
    "layout": {
        "line-join": "round",
        "line-cap": "round"
    },
    "paint": {
        "line-color": "#888",
        "line-width": 8
    }
});

When execute my Chrome console says

Suggestions / examples / alternatives?


texts: Conditional Question on Survey123 Connect

texts:
Anyone here familiar with ArcGIS Survey123? I want my next question to appear based on to the previous question answer. I tried the expression ${previous_question} = 'true' however I don't seem to get it right. When I published the survey by applying the expression, the question seems to be missing from the whole survey.
The question: 32. How often do you use the bus in Mena City?

Choice of answer:
B)During the weekdays
C)During weekends
D)Everyday of the week
E)Only when my regular transport is not available
F)When it is a faster mode of transportation
G)Almost never

The expression I used:  ${_33_how_often_do_you_use_the_bu} = 'Almost never'
I want the next question to be shown based on the answer from B-F
Appreciate if someone could point me to the right direction.


texts:
I am trying to get some polygon features drawn on my map. point features are working fine. even the coordinates of my polygon are correct, i draw them as point features on my map to check that.

i am using the latest version of openlayers (3.2 i guess).

here the code i am using for populating the polygons:

mapView = new ol.View({
      center: [15.438445,47.071028],
      projection: 'EPSG:4326',
      zoom: 17
    });

            // MAP ERSTELLEN
            map = new ol.Map({
                target: 'map',
      layers: [
        new ol.layer.Group({
          'title': 'Basis',
          layers: [
            new ol.layer.Tile({
                        title: 'Bing',
                        type: 'base',
                        source: new ol.source.BingMaps({
                        key:
        'my key here ..',
                        imagerySet: 'Road'
                        }),
                    })
          ]
        }),
        window.overlayGroup
      ],
      view: mapView
            });
    sehenswuerdigkeitLayer = new ol.layer.Vector({
            title: 'Sehenswürdigkeiten',
            type: 'overlay',
      projection: 'EPSG:4326',
            source: new ol.source.Vector()
        });

....

  var polygonGeometry = new ol.geom.Polygon();
          var coordinates = [];

          console.log(geoJson.coordinates);
          $.each(geoJson.coordinates, function(i, polyCoordsArr) {

            $.each(polyCoordsArr, function(i2, polyCoords) {
              if (polyCoords.length == 0) {
                return true;
              }

              var geschaeftFeature = new ol.Feature({
                geometry: new ol.geom.Point(
                  [polyCoords[1],polyCoords[0]]
                )
              });
              sehenswuerdigkeitLayer.getSource().addFeature(geschaeftFeature);

              coordinates.push([polyCoords[1],polyCoords[0]]);
            });
          });

          console.log(coordinates);

          polygonGeometry.setCoordinates(coordinates);
          sehenswuerdigkeitFeature.setGeometry(polygonGeometry);
sehenswuerdigkeitLayer.getSource().addFeature(sehenswuerdigkeitFeature);

the loops are neccessary because i get those coordinates from my database, also i have to switch those coordinates since openlayers use [longitude,latitude] and not [latitude, longitude] as point.

the result of my code gives following object output:

here you can see that the points / coordinates of my polygon are valid (visible as blue circles, ignore the green round circles, thats just other features)


texts:
I want to create an archive of mail I receive and send. Part of my goal is to provide public statistics about this archive, like a world-map with a pin for each location I sent mail to or received mail from. In order to protect the identities of those I receive mail from and send mail to, I would like to anonymize the addresses. This is, instead of pointing to the exact address mail went to, the pin should be placed over the general region the address is in. This region should follow administrative divisions so people have an idea where the mail went to.

The word “region” is intentionally fuzzy as I am not sure what these regions should be. Equating regions to countries is coarse for the statistics I would like to provide. I thought about using ISO 3166-2 codes for the regions, which has the advantage of a standardized and short textual representation, but I haven't found a dataset that resolves ISO 3166-2 codes to coordinates yet.

For instance, when I send mail to Deutscher Bundestag, Platz der Republik 1, 11011 Berlin, Germany, an anonymized address could look like this, from fine to coarse:

11011 Berlin, Germany
Mitte, Berlin, Germany
Berlin, Germany / DE-BE
Germany

Is there a free set of data points that I could use for these anonymized addresses? The dataset should contain for each region:

The English or local name of the region
A pair of coordinates in the center of the region
optionally a short, unique textual representation of the region

It should be obvious for a third party what region is meant when provided with the name.


texts:
I'm creating table in PostgreSQL using these commands:

 shp2pgsql -s 32632 watershed.shp watershed > watershed.sql
 psql -f watershed.sql cambridge
 psql -f watershed.sql cambridge

But I had this FATAL

  password authentication failed for user "ASUS"

How to resolve this?


texts: Create vector extents for PDF files exported within QGIS

texts:
I'm interested to create a vector point layer representing the centroid of some 500-1000 PDF files exported from QGIS to find PDFs geographically rather than file explorer.

The algorithm 'Tile Index' seems most appropriate as the directory file path is retained. Tile Index will leave me with a polygon extent of the PDF from which I will extract the centroid. The issue is that by selecting >1 PDF as the input, no output is created. Below is the result of the algorithm run twice for each PDF shown. They are labelled with a file path for demonstrative purposes.

Is there a bug or PyQGIS workaround to allow more than one PDF input?


texts:
I have some shapefiles that I am reading and then loading into tables in the database. The shapefiles with no zoom (eg. "au_postcode") have no shapes with type null. When I tried to load a shapefile with zoom (eg. "au_postcode_100"), I saw some null shapes.

It would be great if someone can please explain to me how the null shapes work and why would a shape file be null. It is confusing because to me a shape is at least a point and if there is a physical location, there would be a shape for it.

I am new to GIS and if something in this question does not make sense or I've used the wrong terminology, please let me know.


texts:
I'm having trouble adding users to a map in QGIS Cloud Pro. Is there any alternative way anyone knows of to add users to a map other than using the 'edit Map' form?

The 'edit Map' form seems very unresponsive trying to add users from what appears to be a drop-down list of all users on the platform. I just can't add anyone. After typing in a letter into the 'select a user to add' box, the web page just hangs.


texts: Changing decimal digits of extent in raster layer using Python in QGIS

texts:
I want to change the decimal digits of the raster layer extent in QGIS 3.4.3 with Python, so it ends with 0.5.

My code is as follows:

input = vlayer
output = "Directory"
extent = vlayer.extent()
xmin = extent.xMinimum()
xmax = extent.xMaximum()
ymin = extent.yMinimum()
ymax = extent.yMaximum()
coords = "%f,%f,%f,%f" %(xmin, xmax, ymin, ymax)
out_raster = processing.run("gdal:rasterize", {'INPUT': input,'FIELD':'ID','BURN':0,'UNITS':1,'WIDTH':1,'HEIGHT':1,'EXTENT':coords,'NODATA':0,'OPTIONS':'TFW=yes','DATA_TYPE':5,'INIT':0,'INVERT':False,'OUTPUT':output+"vlayer.tif"})

When I ask for the coordinates of the extent with extent I get '-7081.840000,124283.820000,144044.740000,217761.860000'
but I want them to end with 0.5. So in the end the extent is a little bit larger. The values of the extent may then be -7100.5,125000.5,145000.5,220000.5
Round does not work.

Is there a way in Python to change these values?


texts: Selecting an amount of features that together sum up to a certain value?

texts:
I am looking for a way in QGIS to randomly select features within my shapefile. The selected features together need to sum up to more or less  2000m2.
So, it is not selecting 10% of the features, but 10% of the sum of one attribute of the features (20 000 m2 in this case)


texts:
I'm trying to upgrade to the newest version of qGIS but get the following error message when I open the installer:

QGIS requires the GEOS framework version 3.4 or newer.

I've had no luck maintaining my version of geos via the geos download page, and have instead installed it via homebrew.  Is it possible to tell qGIS where the homebrew maintained version of GEOS lives?

I would maintain osgeo via http://trac.osgeo.org/geos/ but keep getting error messages when I run the newest installer (3.4.2) and all advice I've found elsewhere on debugging the installer suggest that I shouldn't bother messing around w/ all the dependencies and should instead use homebrew.

So this brings me back to where I started: is it possible to have homebrew maintain a version of osgeo that qgis can continue to locate?  Or, am I forced to keep trying to debug the geos installer?


texts: How to change the condition of an interaction

texts: How can I change the condition of an exiting interaction? Or do I have to remove the interaction and create a new one with the desired condition? For example I would like to switch the condition of the DragZoom interaction back and forth from ol.events.condition.shiftKeyOnly to having no condition.


texts: Distance of shape points from given lat-long

texts:
I have a point shapefile, and in its attribute table I have some coordinates (DISP_LON, DISP_LAT) which are slightly different from the actual location of the points.

I want to calculate the distance between actual point location and attribute table given co-ordinates, and want to insert in attribute table; preferably in Python, because I don't have ArcObjects SDK.

Initially I was trying to create line from the attribute coordinate and the Shapefile point coordinate. and then calculate line length.

To do some rough work I calculated the point lat-log in two new field name RP_LAT and RP_LONG and tried below code to create line by calculating distance between these two sets of coordinates:

import arcpy

in_rows = arcpy.SearchCursor("C:\Temp\PA.shp")
point = arcpy.Point()
array = arcpy.Array()
featureList = []
cursor = arcpy.InsertCursor("C:\Temp\Line.shp",["SHAPE@"])
feat = cursor.newRow()
for in_row in in_rows:
  # Set X and Y for start and end points
  point.X = in_row.RP_LON
  point.Y = in_row.RP_LAT
  array.add(point)
  point.X = float(in_row.DSPL)
  point.Y = float(in_row.DSPLT)
  array.add(point)
  # Create a Polyline object based on the array of points
  polyline = arcpy.Polyline(array)
  # Clear the array for future use
  array.removeAll()
  # Append to the list of Polyline objects
  featureList.append(polyline)
  # Insert the feature
  feat.shape = polyline
  cursor.insertRow(feat)
del feat
del cursor

But the result is unexpected, it created strange long lines and empty records in a line shapefile.

My ultimate goal is find to distance of point coordinates and attribute given coordinates.


texts: Loading non-spatial data via CSV using PyQGIS for standalone script?

texts:
I am trying to load a .shp file, and a non-spatial data table .csv file.  Then join the table to the shapefile.  Both files have a primary key.

Here is my code:

countyLayer = QgsVectorLayer(countySHP, "County_Layer", "ogr")
QgsMapLayerRegistry.instance().addMapLayer(countyLayer)

Loading the shapefile works fine:

table = QgsVectorLayer(inputFile, "table", "delimitedtext")
QgsMapLayerRegistry.instance().addMapLayer(table)

Loading the table returns "File cannot be opened or delimiter parameters are not valid" in the delimitedtext log in QGIS.  I've tried adding the "?delimiter=;" to the end of the filename and what not still get same errors.  These errors don't stop the code from running however.

processing.runandload('qgis:joinattributestable', countyLayer, table, 'AFFGEOID', 'GEO_ID','memory:layer')
joinedLayer = QgsMapLayerRegistry.instance().mapLayersByName('memory:layer')[0]

This next part actually works... even with the errors above.  It creates and loads a new shapefile WITH the added attributes like I want it too.  However, the table never appears in the layers panel, and when I populate a list of all layers in the map layer registry, only the two shapefiles are present and no table.

I'm fairly new to PyQGIS, so perhaps I am completely overlooking something.

Does the table get stored somewhere else and not the layers registry?

When I first encountered the problem I thought it was something to do with the 'processing.runandload', but when I comment that out, the table still never shows up in the layers panel.  This is a stand alone script so I can't use any iface stuff.  I just don't understand where this .csv is going.  The cookbook has no information on loading delimited text without geographic data, so maybe there is another way to open it.


texts: Burning zeros instead of nodata in GDAL rasterize

texts:
When rasterizing a vector layer I would like to achieve an output that contains 1 for areas covered by the shape and 0 for areas not covered. So far GDAL always converts areas that are not covered to nodata. I tried to change that with the following GDAL command but it does not work properly. I guess he burns nodata as zero but also the metainformation that zero = nodata so I cannot use the layer for raster algebra afterwards (all zells with nodata are excluded, however i'd like them to be zero in the algebra operation).
This is my gdal command:
gdal_rasterize -l GMW_clipped_2010 -a pxlval -tr 0.000269959 0.000269959 -init -999999999.0 -a_nodata 0.0 -te -94.9017777777779 10.918977291452409 -78.1278774106337 23.2768888888896 -ot Byte -of GTiff /home/GMW_clipped_2010.shp /tmp/processing_46ba07b600a54c5581ad63a22017470a/b2dc5cbedbf746a2b9dbbdb406bd35e8/OUTPUT.tif

any help would be much appreciated.


texts:
I'm trying to create a point geometry from a lat, lon pair read from a json, and then write to a database table in postgreSQL that has already been created with the appropriate variable types. I don't think there's anything wrong with my code but it will not write to database.  When I use the command line in pgadmin and put numbers in, it works fine.  What am I doing wrong? Relevant code below:

        points=item['points']
        trip_id = int(trip_ID)

        for point in points:
            lat = float(point['lat'])
            lon = float(point['lng'])
            dist = float(point['dis'])
            elev = float(point['ele'])

            try:
                conn = psycopg2.connect(host=PG_HOST, port=PG_PORT, dbname=PG_DBNAME, user=PG_USER, password=PG_PASS)
                cur = conn.cursor()
                sql ='INSERT INTO mmf_2d (trip_id, dist, elev, geom) values(%s, %s, %s, ST_SetSRID(ST_MakePoint(%s, %s), 4326))'
                cur.execute(sql, (trip_id, dist, elev, lon, lat ))
                conn.commit()
            except:
                print('Save to DB failed!')
                break

Resolution: make parameters in cur.execute string and no spaces after %s, thanks everyone for the help!


texts: Exporting data to new table using ArcPy?

texts:
This is my script and I need to create a new attribute table which contain the information on the print line. Could someone help me?

import arcpy
from arcpy import env
import datetime

env.workspace= "C:/temp/temp.gdb"
arcpy.env.overwriteOutput = True

fc_pnt = "C:/temp/temp.gdb/point"
fc_buf = "C:/temp/temp.gdb/poly"

max_dist = 300
final_bff = "C:/temp/temp.gdb/final_bff"
final_ins = "C:/temp/temp.gdb/final_ins"
fields_buff = ["SHAPE@","id_local"]
fields = ["OID@", "SHAPE@AREA"]
fields_final = ["Shape_Area", "SC_11"]

with arcpy.da.SearchCursor(fc_pnt,fields_buff) as curs:
    for row in curs:
        pnt_cc = row[0]
        id_local = row[1]
        arcpy.Buffer_analysis(pnt_cc, final_bff, str(max_dist))

        for i in fc_buf:
            arcpy.Intersect_analysis([fc_buf,final_bff],final_ins,"ALL","","INPUT")

            with arcpy.da.SearchCursor(final_ins, fields_final) as cursor:
                for row in cursor:
# i need to pass this print to a new table
                    print(str(id_local) + " " + "{0}, {1},".format(row[0], row[1] ))


texts:
I'm trying to use the SAGA tool ''split Lines by points''.
I have created the  line layer using '' Polygon to lines''

My goal is to split the lines. I created points using Extract points

I want to have new lines based on those points. I tried to use ''split lines by points'' but then I get something like this

I have changed the epsilon but it still has the same problem.

My end goal is to have ''independent'' lines that I can divide into segments of 12.5 meters.


texts:
I've been trying to set up a few timed WMS layers in GeoServer. They're simple enough, of the form

ID(PRIMARY KEY) REGION DATE VALUE

so for every region and date you have some value, e.g. the precipitation in that day.

When I set up the layer without time dimensions, I get a working layer where I can see all my regions (although each of them is present multiple times, once for every corresponding date).

Then, I enable the time dimension, with attribute = DATE and presentation = LIST. When I call GetCapabilities, it really contains something like

<Dimension name="time" default="current" units="ISO8601">
  2003-12-31T00:00:00.000Z,2004-12-31T00:00:00.000Z,2005-12-31T00:00:00.000Z
</Dimension>

However, when I request GetMap, results vary.

with MS SQLServer backend, I get blank image if the type of DATE column is date, but I get the correct map image if it is datetime
with PostGis backend, I always get blank image, whether I use date or timestamp

( one day granularity is enough for my list, so to represent 2003-12-31 I'm ok with using any timestamp within that day )

Note that GeoServer considers all of these datatypes compatible with the time dimension and allows the user to set them as the time attribute.

Has anyone encountered such problems? Am I misunderstanding something about the WMS-T?

EDIT: I've recently found out that timestamptz is the type that works with PostGis. Still, why...?


texts: How to get hypsometric curve from elevation data?

texts: I have my elevation data but cannot  make an attribute table from it and so I'm not sure how to get any of that data into a table in order to make a hypsometric curve.


texts: OpenLayers 6 how to add line breaks to text labels?

texts:
I'm trying to add simple line breaks to my text labels in ol 6 using \n but somehow both text line drift far apart as soon as I zoom out. It's weird. I created a simple fiddle, just click on the button "set style".
The relevant part:
var createTextStyle = function (feature, resolution) {

  return new ol.style.Text({
    text: 'first line' + '\n' + 'second line'
  });
};

function polygonStyleFunction(feature, resolution) {
  return new ol.style.Style({
    stroke: new ol.style.Stroke({
      color: 'blue',
      width: 1,
    }),
    text: createTextStyle(feature, resolution),
  });
}

The first line should be above the second text line but it's not.
https://jsfiddle.net/a7cbgtd9/
The \n should not be a problem here because it works in this example:
https://openlayers.org/en/latest/examples/vector-labels.html
Im even using the same ol version.


texts: WFS request returns FeatureCollection with no feature - OpenLayers

texts:
Created the WFS feature request as following:

this.map.getViewport().addEventListener("click", (evt) => {

  //Point transformation to EPSG:25830
  var proyeccion = getProjection('EPSG:25830');
  var p = new Point(this.map.getEventCoordinate(evt)).transform(this.map.getView().getProjection(),proyeccion);

  //WFS Feature Request
  var featureRequest = new WFS().writeGetFeature({
    srsName: 'EPSG:25830',
    featurePrefix: 'CP',
    featureTypes: ['CadastralParcel'],
    outputFormat: 'application/json',
    filter: containsFilter('MultiPolygon',p,'EPSG:25830')
  });

  fetch('https://inspire.navarra.es/services/CP/wfs', {
    method: 'POST',
    body: new XMLSerializer().serializeToString(featureRequest)
  }).then((response) => {
    return response.json();
  }).then((json) => {
    var features = new GeoJSON().readFeatures(json);
    this.vectorSource.addFeatures(features);
  });
});

This makes the getFeature request correctly using POST, with the following request payload:

<GetFeature xmlns="http://www.opengis.net/wfs" service="WFS" version="1.1.0" outputFormat="application/json" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd">
    <Query typeName="CP:CadastralParcel" srsName="EPSG:25830">
        <Filter xmlns="http://www.opengis.net/ogc">
            <Contains>
                <PropertyName>MultiPolygon</PropertyName>
                <Point xmlns="http://www.opengis.net/gml" srsName="EPSG:25830">
                    <pos srsDimension="2">586784.8369967575 4722127.891344394</pos>
                </Point>
            </Contains>
        </Filter>
    </Query>
</GetFeature>

But the request above returns an empty FeatureCollection when clicking a point where there is a feature. The type of Features that I was supposed to get were Polygon.

{"type":"FeatureCollection","totalFeatures":"unknown","features":[],"crs":null}

Guess it could be because of the projections, but don't know where the mistake is.

EDIT

I'm using projection EPSG:25830. Defined like this:

proj4.defs("EPSG:25830","+proj=utm +zone=30 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs");
register(proj4);

Added to the map's view and layers as following:

const proyeccion = getProjection('EPSG:25830');

view: new View({
    projection: proyeccion,
    center: fromLonLat([-1.629950,42.63],proyeccion),
    zoom: 9
})

//WMS
source: new TileWMS({
    url: 'http://localhost:3000/idena',
    params: {'LAYERS': 'catastro'},
    projection: proyeccion
})


texts: Plotting scatterplots with subset of rasters in R according to land cover

texts:
I have three rasters:

> raster1
class      : RasterLayer
dimensions : 2803, 5303, 14864309  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : 60.85, 105.0417, 15.95833, 39.31667  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
source     : memory
names      : layer
values     : 0, 31.72  (min, max)

> landcoverraster
class      : RasterLayer
dimensions : 2803, 5303, 14864309  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : 60.85, 105.0417, 15.95833, 39.31667  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
source     : C:/Users/XX.tif
names      : landusemaskedme
values     : 1, 12  (min, max)
attributes :
       ID
 from:  1
  to : 12

I would like to plot scatterplots between raster1 and raster2 according to landcover types in landcoverraster. But first I would like to subset a random subset of raster1 and raster2 (as they are very large rasters and when I plot scatterplots it gives a big chuncky scatterplot) according to landcover type. How can I subset a random raster according to land cover type and plot a scatterplot between raster1 and raster2? I am looking a desired output similar to this

What I initally have done (without subsetting), and came up with this plot:

plot(raster1[landcoverrraster==1],raster2[landcoverraster==1])

Reproducible rasters:

library(raster)

ras1 <- raster(matrix(c(1,1,1,2,2,2)))
ras2 <- raster(matrix(c(1,1,1,2,2,2)))

#Generating landcover example data
raster2 <- raster(matrix(c(1,1,1,2,2,2,3,3,3),ncol =3))
raster2 <- as.factor(raster2)

rat <- levels(raster2 )[[1]]
rat[["landcover"]] <- c("land","ocean/lake", "rivers")
levels(raster2 ) <- rat


texts: Geoserver Extension access restriction with proxy nginx

texts:
I have installed Geoserver on my server. And I have proxy server(installed nginx) to access Geoserver.

location /geoserver {
    proxy_pass http://127.0.0.1/geoserver/;
}

http://mydomain.com/geoserver ==> to ==>    http://localhost:8080/geoserver

I installed geoserver extensions (Shape, Excel, DXF) that users can download geographic data.

http://localhost:8080/geoserver/wms?request=GetFeature&typeName=Polygons&BBOX=-103.929,44.375,-103.633,44.500&outputFormat=dxf

But someone can download all geographic data in a city of BBOX . I want to restrict access download. How can I do this with proxy server or other solutions?


texts:
I would like to extract a named subset of the text records from a shapefile using Python/pyshp. Is there a clean/efficient way to do this?

Ideally I'd like to use dictionary-like access such as:

import shapefile
sf = shapefile.Reader("MyFile.shp")
data = sf.record['desired_record_name']

But of course this doesn't work. So far my solution is the brute force approach of looping over all entries:

import shapefile
sf = shapefile.Reader("MyFile.shp")
for r in sf.iterRecords():
    if r[0]=='desired_record_name':
        data = r
        break

But this leaves a lot to be desired, especially for large/many files.

Note: Although I like the simplicity of pyshp, I am not tied to it, so other free & open-source solutions using python are also of interest.


texts: Reading label settings with PyQGIS in QGIS 3?

texts: I am struggling to read (with PyQGIS) even the most basic settings from a labeled shapefile in QGIS 3. I know this was pretty easy in QGIS 2 with .readFromLayer but it does not seem to work in QGIS 3. Is there another way that works in QGIS 3?


texts: Measuring distance from points to polygon: won't measure in metres!

texts: I'm using ArcGIS 10.2. I have some point data (sampling points) and I want to know how far each of these is from the nearest protected area (polygon data). I have tried doing this with both the Join function as well as the near function, it gives me values but they are in decimal degrees, I want it in metres. I have tried changing the projection which does not help, the data is from Namibia so I tried projections: Africa equidistant and WGS 1984 Zone 33S. Any suggestions would be great


texts: Getting ArcGIS for SharePoint up and running?

texts: I am having some issues getting ArcGIS for sharepoint up and running.  I'm not overly fluent with sharepoint, but I don't know what's going on.  I have it installed on my sharepoint along with ESRI maps for sharepoint, and I cannot get either to load when I add them to the site.  At first, the ArcGIS had a map background show up, but now the box shows up with nothing in it.  I can edit the actual properties that pop up, but I cannot use any of the edits on the top of the screen such as add data, change basemap, or anything really.  Not sure if something happened but hoping someone can assist.


texts: ArcPy Not Able to Open Current Map From Stand Alone Code While Map Is Open

texts:
Running an Arcpy stand alone script on PyScripter I followed everything in this answer to ArcPy gives error Object: CreateObject cannot open map document?

As:

Start up the Streets.mxd in ArcMap
Run the mxd=arcpy.mapping.MapDocument("CURRENT") inside the ArcMap (without no error)

3- Add the same mxd=arcpy.mapping.MapDocument("CURRENT") to my code in PyScripter IDE as

import arcpy
####mxd = arcpy.mapping.MapDocument(r"C:\Temp\Street.mxd")
mxd=arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
targetGroupLayer = arcpy.mapping.ListLayers(mxd, "Main", df)[0]
addLayer = arcpy.mapping.Layer(r"C:\Data\NetMap.gdb\Streets\OneWay")
arcpy.mapping.AddLayerToGroup(df, targetGroupLayer, addLayer, "BOTTOM")
arcpy.RefreshActiveView()
arcpy.RefreshTOC()

but I am getting this runtime error message:

  Traceback (most recent call last):   File
  "C:\Users\bhosseini\Documents\ArcPy\Trace With Graphics\module3.py",
  line 3, in
      mxd=arcpy.mapping.MapDocument("CURRENT")   File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\arcobjects\mixins.py", line 625,
  in init
      super(MapDocumentMethods, self).init(mxd)   File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\arcobjects_base.py", line
  47, in init
      for arg in args)) RuntimeError: Object: CreateObject cannot open map document


texts:
When making a sum between two raster layers, for every cell, if one of the two cells have no data in that cell, the opeartion is not done in that cell, example:

In QGIS raster calculator I do:

Plus (+)

And I get:

The red circle in the middle is like cut in the result, QGIS it's not adding every pixel on each raster, it makes a kind of intersecction...

To try to avoid this, I need to fill every no-data cell in both raster layers with a 0 value. How can i do this in QGIS?


texts: Enumerating messages in arcpy.GetMessages()?

texts:
Messages are present in GetMessages() but it's a formatted string. I want to add some of my own newlines and alter the content on a per-message basis.

How do I access the individual messages in arcpy.GetMessages()?


texts: Join table field with shapefile programatically (via PyQgis)?

texts:
The following is currently unsuccessful at fulfilling a table field join to a point shapefile in Qgis.
All instances of the join object will print the proper values, and no error is flagged when processing - but the target layer does not exhibit any effect.

Something missing?

    Heads = QgsVectorLayer(headCSV, 'Heads', 'ogr')
    QgsMapLayerRegistry.instance().addMapLayer(Heads)
    self.iface.legendInterface().setLayerVisible(Heads,False)
    grid = ftools_utils.getMapLayerByName(unicode('Grid'))
    joinObject = QgsVectorJoinInfo()
    joinObject.targetFieldName = str('ID')
    joinObject.joinLayerID = Heads.id()
    joinObject.joinFieldName = str('cellID')
    grid.addJoin(joinObject)
    grid.reload()

UPDATE:

The following returns 'None'.  I expect this is actually supposed to return the field names to be joined?  (e.g. - a clue to the mistake)

    print joinObject.joinFieldNamesSubset()

FURTHER UPDATE:

Added the following, and now the previous update command returns the accurate fields to be joined - yet the destination layer does not show 'joined' fields...

    joinObject.setJoinFieldNamesSubset(['someFieldName'])


texts: Transformation parameters for Everest_Bangladesh and Nepal_Nagarkot to WGS84

texts: What are the parameters for transformation of Everest_Bangladesh and Nepal_Nagarkot Geographical Coordinate System datums to WGS84 datum?


texts: How to find a least cost path through several points?

texts:
I want to do a least cost path analysis and the least cost path has to go through several different points instead of having just a start and end destination. Is this possible?

I am using ArcGIS 10.2.2


texts: Snapping line segments together to remove gaps between them

texts: I have a series of line segments that were created based on the location of an existing line.  There are small gaps of around 0.1 feet between the different line segments.  I need to somehow snap the 2 ends of the line segments together to remove these small gaps.  Manually snapping them together is not a option as there are literally 100s of these gaps in my dataset.  If arcgis has a tool that will do this that is great but I'm willing to try python or VB whatever will get the job done.  I have included a screencap that illustrates my problem.  The green line is the original and the black segments are what I'm trying to stretch.  I tried geometric networking but it is moving the whole line segment around.  I just need to stretch the line.


texts: Coordinate transformation from ITRF2008 to ITRF96


texts: Re-running Calculate Field tool in ArcGIS Pro?

texts:
Calculated a field with a little ArcPy function yesterday and now I want to change it. Just realized that the function I defined, once run and the geoprocessing pane is closed, disappears.

Is there a way to add the function to the gdb toolbox so I can edit and re-run at will?


texts: OWSLib vs shp2img output of MapServer Layer

texts:
I'm working on a map layout engine using Python and ReportLab and have come across a problem and was wondering if anyone can give me some insight.

Scenario:

I'm trying to generate an image using OWSLib (see example code).  The layer that I'm trying to render has an opacity of 40%.  When I generate the image using the sample code the background is transparent but the features are not but if I generate the image using shp2img -m /www/cgi-bin/mapfiles/40e6b3a8.map -o test.png the features are 40% opaque.  Any ideas how I can make this work using OWSLib or if I'm doing something wrong?

OWSLib example code

from owslib.wms import WebMapService
wms = WebMapService('http://localhost/cgi-bin/mapserv?MAP=/www/cgi-bin/mapfiles/40e6b3a8.map&', version='1.1.1')
img = wms.getmap(   layers=['vlyr_oqiJw4Ib4U'],
    srs='EPSG:26910',
    bbox=(685400,5616300,689500,5619500),
    size=(800, 600),
    format='image/png',
    transparent=True
)
print(img.geturl())
out = open('c:/temp/test2.png', 'wb')
out.write(img.read())
out.close()


texts: Which registry key holds the license type in ArcGIS 10.x?

texts:
I'm unable to change my license type (ArcView, ArcEditor or ArcInfo) in ArcGIS 10.0 using the Desktop Administrator due to permissions restrictions.

I'm trying to use the back-door method of changing the registry key in Regedit. In older versions of ArcGIS this was found under > HKEY_LOCAL_MACHINE > Software > ESRI >
License but I'm not seeing this in ArcGIS 10.0

The post Opening ArcGIS on a specific license level mentions a similar problem on Citrix, however I'm unable to set a system variable due to the same permissions issues.

Which registry key holds the license type in ArcGIS 10.0, and is it still possible to change the license type using this method?


texts: Seeking ArcPy equivalent to IPixelBlock in ArcObjects?

texts:
Trying to figure out a way of doing this via ArcPy, but don't see anything remotely related to this.  My raster dataset is a GRID.

    'pPixBlock is an IPixelBlock

    'sets raster to be one column and one row in size
    pSizePoint.SetCoords(1, 1)

    pPixBlock = praster.CreatePixelBlock(pSizePoint)


texts: Setting fill style in polygon in OpenLayers?

texts:
I have created a polygon using Vector layer in OpenLayers 3 now I have to fill that polygon with red color. But I cannot do this using style option in OpenLayers 3.

Any one have any idea how to do this?


texts: When using feature-linked annotation in ArcGIS 10, can I create multiple annotations from one feature?

texts: I have a road that I want to create multiple annotations for. Is it possible to attach more than one annotation to the same feature?


texts: "LineStrings must have at least 2 coordinate tuples" error in QGIS polygonize

texts:
I am using QGIS version 2.18.9. I wanted to convert DXF into polygon. First, I used the tool called "v.in.dxf" under the processing tool. Using that tool, line shapefile was created. Then I ran the "v.clean.advanced" tool under the processing tool. Then output of that tool was used as the input to the polygonize tool. While running that "polygonize tool, the following error appeared:

  LineStrings must have at least 2 coordinate tuples See log for more details

Can any one explain this?


texts: What is this PostGIS query doing to show great circle connections?

texts:
I'm trying to draw a flight map in CartoDB. My starting point is a single lat long point. We'll pretend it is this:

ST_GeomFromText('POINT(-71.060316 48.432044)', 4326)

And I want draw lines between that point and the points in my CartoDB table. I had a hard time even getting started, but someone pointed me to this map:

http://bl.ocks.org/andrewxhill/raw/8406976/

Which is fancier than what I need (I don't need to redraw everytime you click someplace) but includes the basic query:

SELECT ST_Transform(ST_Segmentize(ST_MakeLine(
   ST_Transform(the_geom, 953027),
   ST_Transform(
          (SELECT ST_PointOnSurface(the_geom)
          FROM lower_48_zips
          WHERE cartodb_id = "+data.cartodb_id+"), 953027)), 100000), 4326)
   the_geom FROM us_ski_areas WHERE mindist < 1000

Someone also pointed me towards an example which uses a very similar query:

UPDATE experimental.airroutes
SET the_geom =
(SELECT ST_Transform(ST_Segmentize(ST_MakeLine(
       ST_Transform(a.the_geom, 953027),
       ST_Transform(b.the_geom, 953027)
     ), 100000 ), 4326 )
 FROM experimental.airports a, experimental.airports b
 WHERE a.id = airroutes.source_id
   AND b.id = airroutes.dest_id
);

But I'm having trouble pulling either of these apart. I'm decently conversant in MySQL but totally new to PostGIS and I would love help taking apart what is happening here and replacing source_id with a specific point.

I thought I could do something like

UPDATE my_table
SET the_geom =
(SELECT ST_Transform(ST_Segmentize(ST_MakeLine(
       ST_Transform(my_table.the_geom, 953027),
       ST_Transform(ST_GeomFromText('POINT(-71.060316 48.432044)', 4326), 953027)
     ), 100000 ), 4326 )
 FROM my_table
);

But I know I need a "where" clause in there and I'm not even sure what it should be.


texts: Please any body help me for creating lines from points in electric network for instance electric poles are there in point feature with "polename" field like "P01", "P02", "P03" etc. I need a tool that create lines on Target (Feature Class) from P01 to P02 then P02 to P03 in continue.


texts: Python OGR: Filter Shapefiles by Lat Long and Export To CSV

texts:
Using osgeo's Python API for OGR, I'd like to implement the following process:

Import ESRI Shapefiles from a given directory (call it input_dir)
Import a CSV file that contains lat-longs of interest (call it poicsv)
Join the CSV file long, lats with the Shapefile layer and extract the feature attribute values for that long, lat

I am currently up to the step where I have a single Shapefile loaded, and the CSV file loaded with long/lats in a list of tuples (ie. [(-123.76073, 41.24946), (-123.84451, 41.24946), ...]). Per long lat, I am trying to use SetSpatialFilter on the layer by the following code:

for ll in longlat:
    point = ogr.Geometry(ogr.wkbPoint)
    point.AddPoint_2D(ll[0], ll[1])
    layer.SetSpatialFilter(point)
    for feat in layer:
        feat.items()

However, upon sanity checking I've noticed there is no API for removing the Spatial Filter in my for loop, and I am also concerned about the memory leak created by creating a wkbPoint through the ogr.Geometry function in every iteration. Will the point object be cleaned up every iteration, and how could I reset the Spatial Filter every iteration?


texts: How do I extract only the visible placemarks from a KML?

texts:
I've been given a KML file containing over 1000 placemarks, a small subset of which were selected by a user in Google Earth and made visible (by checking the checkbox next to the feature name). I've looked in the KML using a text editor and placemark visibility is set by the <visibility> KML tag. Is there a way to save only those features in the KML that are visible, or export them to some other format?

I'd like to generate a table of the features containing the names (lat/longs would be nice too) of only those features that are visible (turned on) in the KML, or at least a table that shows me the Visibility of each placemark so I can filter out the ones that are not visible. I've tried converting the KML using Xtools, ArcGIS' built-in KML conversion, and QGIS and I get a layer containing all the features in the KML, even the ones that are turned off. I have Google Earth Pro but haven't found a tool in it to do this.

ETA: One more problem: the feature names are not contained in the placemark's <name> tag, they're in a SimpleData element.


texts: Geoserver administration - defining different privileges

texts:
I need to create separate admin accounts for Geoserver. I'd like to have one superadmin account that would be able to read/edit/delete any style, layergroup, workspace, store or layer - that is what the default admin is already able to.

I also need to have several accounts that would be able to create workspaces, layers and stores, would not be able to edit styles and would be able to edit and delete only the workspaces, layers and stores they created (e.g. if superadmin creates a workspace, he is the only one to delete it).

Is this possible to achieve with Geoserver?


texts: L.esri.Cluster.featureLayer does not call style function

texts:
Using the L.esri.Cluster.featureLayer plugin for Leaflet to cluster features from Esri does not call the style method despite the documentation stating that this is the way to style vector layers.

Docs: https://esri.github.io/esri-leaflet/api-reference/layers/cluster-feature-layer.html

Options has a style method like the one for vanilla Leaflet feature layers. L.esri.featureLayer does use the style option properly, but the cluster version does not appear to be doing this - or I am doing something wrong.

Sample code:

var bikePaths = L.esri.Cluster.featureLayer({
    url: 'https://services.arcgis.com/uCXeTVveQzP4IIcx/ArcGIS/rest/services/Bike_Routes/FeatureServer/0',
    onEachFeature: function(feature, layer)
    {
        // this works - called for each feature
    },
    pointToLayer: function (feature, latlng)
    {
        // this works - called for each point feature
    },
    style: function (feature, layer)
    {
        // never called - despite having linear and polygon features
        // the linear and polygon features do show up fine.
        // they don't cluster, but that's to be expected (they are vectors)
        // wth documentation?
    },
  }).addTo(map)

Anybody have a solution for this?


texts:
I have to write a script in Arcgis which I would like to use the SpatialJoin tool (arcpy) but I do not have the appropriate extension.

It may be a tool (or example of scripts) for doing the same thing but using the basic license?


texts: Landsat oli mosaicking for LULC classification

texts: I have to classify land use land cover data for western ghats- 1,07,770 sq km area.The landsat images for the area cover 15 tiles . I have to use minimum cloud cover images . But for the same the images available at minimal cloud cover are from different dates . Can i mosaic these images of different dates (say dec1 2014 , 30 oct 2014 ) to achieve YEARLY lulc without any problem ?


texts: Downloading images on Planet Explorer using free student account?

texts:
I am logged in to @Planet Explorer with a free student account. I performed a search using an uploaded shapefile of my AOI. I selected a 3-month mosaic to reduce the amount of cloud cover and I am satisfied with the result.

How can I now download the actual images?

There is no download button on the interface!


texts: How to search and move to feature in QGIS?

texts:
After solving how to search text string I need to move to the selected feature.
Is it possible to do that?


texts: Does anyone know how to perform PCA using PL/pgSQL in a PostgreSQL function, considering correlation matrixes of potentially high order? Maybe a specific library for this problem.


texts:
I have read that Australia is updating its local reference frame.

If a small territory has a complete CORS GNSS network aligned to ITRF2000, does it need to update the CORS network to the latest ITRF2014?


texts: What is the proper way to backup an ESRI File Geodatabase that is published to ArcGIS Server?

texts:
I have an ESRI file geodatabase (v10) that is published to an arcgis server map service. When the service is running, the fGDB is locked. Do I need to stop the service to get a clean backup? Or, is there a way to perform a backup through an arcpy script or through Catalog? I'm currently using Windows robocopy to transfer the fGDB to a backup drive. Here's the output showing the locked files:

 New File           0 Bikepaths.CFP0026.4968.5140.sr.lock
 New File           0 BuildingFootprints.CFP0026.4968.5140.sr.lock

etc, etc...


texts: Add Google Maps Layer to OpenLayers Page

texts:
<script src="openlayers.js"></script>
<script src="http://maps.google.com/maps/api/js?v=3&amp;sensor=false"></script>
    <script type="text/javascript">
....
  var google = new OpenLayers.Layer.Google(
            "Google Streets", // the default
            {numZoomLevels: 20}
        );
....
   map.addLayers([
      main_layer,
      google
   ]);
....
</script>

Before adding the Google Maps Layer like this, everything worked. Now I get an TypeError: google.maps is undefined.  What is wrong?

Please tell me if I totally misunderstood the concept of adding a google maps layer to openlayers


texts: Delete overlapping points between two shapefiles QGIS

texts:
I have a point shapefile of all human settlements in South Africa (ranging from tiny hamlets to major cities). In another shapefile, I have point features as well, representing settlements located inside an urban area.

I'd like to join the cities falling inside the urban extent to the shapefile of all settlements. However, if I do so I will get repeated features (cities in the urban extent also exist in the human settlement shapefile).
Is there a way to avoid having this happen?
The first picture represents all the human settlements, and the second only the settlements inside an urban extent.


texts:
I'm ArcGIS beginner.

At the moment I have 2 layers:

roads for whole country
zone of the city

I would like to delete the roads that are outside the city. Any idea how can I do it?