texts: Listing all of contents of geodatabase with ArcPy?

texts:
This seems like a simple problem, but I have not found a way to accomplish it. I am writing a script that needs to record the paths of everything inside of a geodatabase to a text file, and later copy those files to a new location and zip them up.

I have been able to use arcpy.ListFeatureClasses to grab feature classes, and arcpy.ListRasters, and arcpy.ListDatasets, but how do I obtain a list of, say, Address Locators in my gdb? I haven't found a way to do it with arcpy. There should be a simple function that lists all of the contents of a geodatabase. Is there one I'm missing? Or is there any other way to accomplish this in python?


texts: Splitting polygon with line using ArcPy?

texts:
In looking at Buffering with physical barrier using ArcGIS for Desktop?, it occurred to me that I am not sure how one would go about using geoprocessing tools in ArcGIS to split a polygon with a line programatically.

Manually, you would use the Cut Polygons tool or the Split Polygons tool on the Topology toolbar, but how would you accomplish the same task using modelbuilder or python groprocessing scripting tools?

Right off the bat I think of all of the tools in the Analysis toobox like Union, Identity, etc, but those are all Polygon-Polygon tools, NOT Polygon-Line tools.  Even the Split tool is Polygon-Polygon.

Any ideas?


texts:
I'm running some process for about 48 hours and I'm not sure about what is going on or maybe better - how should I interpret the log.

Please have a look at attached part of the log.
Just after successful ? clipping (at least is something I can read from log ?) the FME executes "performing low-level intersection" and repeat that part at least about 200 times in my case. Is it correct ?
Should I patiently wait to finish that part of the process or it somehow stuck in that place ?

Thanks in advance for help !

Finished clipping 303150 / 525807 clippees against all clippers
Performing low-level intersection at phase #1... 100% done
Setting up lines for phase# 2 ... 100% done
Setting up lines for phase# 2 ... 100% done
Performing low-level intersection at phase #1... 100% done
Setting up lines for phase# 2 ... 100% done
Setting up lines for phase# 2 ... 100% done
Performing low-level intersection at phase #1... 100% done
... Last line repeated 2 times ...
Setting up lines for phase# 2 ... 0.48% done
Setting up lines for phase# 2 ... 100% done
Setting up lines for phase# 2 ... 100% done
Performing low-level intersection at phase #1... 100% done
... Last line repeated 16 times ...
Setting up lines for phase# 2 ... 0.48% done
Setting up lines for phase# 2 ... 100% done
Setting up lines for phase# 2 ... 100% done
Performing low-level intersection at phase #1... 100% done
Setting up lines for phase# 2 ... 0.48% done
Setting up lines for phase# 2 ... 100% done
Setting up lines for phase# 2 ... 100% done

Thanks for remarks. Two data sources: Geotiffs (about 13k of tiles) and water areas features (about 600k) in shapefile. The bounding boxes of geotiff's in the first part of the workflow should clip the WA features. According to info from Dale process all DEM tiles at once raster feature's don't hold any data... so I believe I'm not processing DEM data at all and in the process only them bounding boxes are used.
Then it should be simply use of clipper with  clippers and clippees both being vector features. This is also a place of the process that I believe (interpreting the log) the process stuck so I'm not describing the further part of the workbench.
If any won't find the extra information enough I'm able to share with the workbench which will be more self-explanatory for the rest of the process (can't find any button to attach it here).


texts: Google geocoder doesn't return U.S. county information

texts:
I'm geocoding a bunch of U.S. addresses with various levels of completion. In some cases, Google's geocoding API (https://maps.googleapis.com/maps/api/geocode/json?address=) doesn't return a county (administrative_area_level_2) although the data appear specific enough to determine a county.

For example, the following sample data (in quotation marks) don't return a county although, to the best of my knowledge, they can be assigned to a single unambiguous county:

Street-level resolution data

"Las Vegas Blvd S. 89109 Las Vegas" (Google recognises "South Las Vegas Boulevard")

"155 West Dewey Drive WI 53149 Mukwonago" (Google doesn't recognise "West Dewey Drive")

Postal code-level resolution data (Google recognises the postal codes)

"83714 Boise", "11702 Babylon"

Is this a bug? Is there a way to get Google to return a county? Or is there another way to do it? In order to process at least the subset of single-county postal codes, is there a postal code-to-county geocoder out there?


texts: Request JSONP from Geoserver using AJAX request in Openlayers 3

texts:
I am creating a web application using OpenLayers 3 and I am looking to import my vector layers into this. I currently have my data in Geoserver and I am trying to use the following code, but when I make the below request using jQuery:

    var vectorLoader = function(extent, resolution, projection) {
    var url = 'http://XXX185:8080/geoserver/wfs?service=WFS&' +
        'version=1.1.0&request=GetFeature&typename=YYY:QGIS&' +
        'outputFormat=text/javascript&format_options=callback:loadFeatures' +
        '&srsname=EPSG:3857&bbox=' + extent.join(',') + ',EPSG:3857';
    $.ajax({
      url: url,
      dataType: 'jsonp'
    });
  };

  var loadFeatures = function(response) {
    var features = vectorSource.readFeatures(response);
    vectorSource.addFeatures(features);
  };

  var vectorSource = new ol.source.ServerVector({
    format: new ol.format.GeoJSON(),
    loader: vectorLoader,
    strategy: ol.loadingstrategy.createTile(new ol.tilegrid.XYZ({
      maxZoom: 19
    }))
  });

  var serverVector = new ol.layer.Vector({
    source: vectorSource,
    style: vectorStyle
  });

I get this error:

Uncaught SyntaxError: Unexpected token <
jquery-latest.min.js:4 Resource interpreted as Script but transferred with MIME type application/xml: "http://XXXpc185:8080/geoserver/wfs?service=WFS&version=1.1.0&request=GetFea…SG:3857&callback=jQuery111109427054924890399_1428089051875&_=1428089051876".

I imagine it is something to do with the way the vectorLoader is requesting the WFS?


texts: Is there a way to clip the Google Satellite Layer from OpenLayers plugin?

texts: I have imported the Google satellite imagery though the OpenLayers plugin. I would like to clip it to a polygon shapefile of New York State. Is it possible to do this? I tired the GDAL Tools and it wont let me select the imagery layer.


texts: Copying attributes of one row to another row in polygon shapefile in ArcGIS?

texts:
I have a single polygon shapefile in ArcMap with 25*25m grids, and I want to assign the attribute values of a single grid/polygon to another one close by i.e. copy attributes from a row to another row. Is there a way to do this other than Merge or manually as I want to retain the 25*25m grid size of polygons?

The light red colour square polygon in the middle has to be dark red as it is lying in between a bunch of dark red ones. There are numerous such spots in my data that I have to change so I am looking for an easier and faster way to change the attributes and trying to avoid doing manually. I have to retain the 25*25m polygon size as well. Any solutions?


texts: Wrong size when use Mapproxy to serve Mapnik tiles as ESPG: 900913

texts:
In MapProxy Demo TMS page, EPSG:4326 looks fine, except the map is tilted and I don't want that. Tried to make EPSG:900913 work, However, no matter what I try. the map is too small, at zoom level 18 I can not see street, the whole city remains a little dot on screen for all zoom levels. Or sometime I got nothing.

Background:

Street data stored in Spatialite in Lat Lng format, and all geometry
has srid 4269
Uses MapProxy 1.5 and Mapnik 2.2.0
tried to append 900913 to mapnik\share\proj\epsg file (seems mapnik don't support 900913 natively). randomly give value such as

<900913> +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs <>
<900913> +proj=longlat +datum=WGS84 +no_defs  <>

randomly change all srs string in mapnik xml config to srs strings above
Map looks fine in TileMap with SRS "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs", map looks small and tilted with SRS "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0.0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs +over"
I thought Lat Lng data is absolute format, any projection should works fine with my data.
Clear Cache and restart server every time.


texts: Is there a way to identify empty feature classes in ArcMap without deleting them?

texts: Is there a tool or script in ArcMap (currently using 10.5.1 basic) to identify empty feature classes without deleting them? I have to keep all feature classes in the geodatabase, even if they are empty, but I would like to know which ones actually have data without going through every feature class. There are up to 522 feature classes per geodatabase, with 100+ databases.


texts: How to designate a column of polygons as a georeference

texts:
I uploaded a dataset with a column of polygons, but the map is blank. It's been a few years since someone has asked this question and looks like the dashboard has changed since then. How do I designate that column of polygons as the georeference for a map?

I read there was a 'georeference' button near the 'export' button on the previous dashboard. The polygons are pandas objects like this:
POLYGON ((1006076.28338623 262138.1090087891
Which is the same as I see here in this map


texts: Loop using SearchCursor and extract by attribute tool

texts:
I have multiple rasters with 14 different categories. I want to use SearchCursor to loop through each raster's attribute table, obtain "Value" and use it in ExtractbyAttributes tool. This is the code I've so far, but it is giving me an invalid SQL statement error which I don't understand as I'm disregarding the optional SQL statement. I'm using 10.4.

inRas = "C:/ABCD_10000_P6_Mos"

rows = arcpy.SearchCursor(inRas, "", "", "", "")
for row in rows:
    luClass = row.getValue("Value")
    print (luClass)
    rasExtract = ExtractByAttributes(inRas, luClass)
    rasExtract.save (inRas + "_" + str(luClass) + "_P6_Split")

    ERROR:
    Runtime error  Traceback (most recent call last):
    File "<string>", line 7, in <module>
    File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\sa\Functions.py", line 1324, in ExtractByAttributes     where_clause)
    File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\sa\Utils.py", line 53, in swapper     result = wrapper(*args, **kwargs)
    File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\sa\Functions.py", line 1320, in Wrapper     out_raster)
    File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\geoprocessing\_base.py", line 506, in <lambda>     return lambda *args: val(*gp_fixargs(args, True))
    ExecuteError: ERROR 999999: Error executing function.
    An invalid SQL statement was used. An invalid SQL statement was used. [VAT_ABCD_10000_P6_Mos] An invalid SQL statement was used.
    [SELECT * FROM VAT_ABCD_10000_P6_Mos WHERE 1] Failed to execute (ExtractByAttributes).

Updated code:

inRas = "C:/ABCD_10000_P6_Mos"

rows = arcpy.SearchCursor(inRas, "", "", "", "")
for row in rows:
    luClass = row.getValue("Value")
    sqlQuery = "Value = luClass"
    print (luC)
    rasExtract = ExtractByAttributes(inRas, sqlQuery)
    rasExtract.save (inRas + "_" + str(luClass) + "_P6_Split")


texts: Visual connection between two linked frames in QGis Print composer

texts:
I am still a beginner.

I created a printing preview in Print composer, with a general map on the background and several small frames on top, containing zoomed-in portions of the general map.

With the "Overviews" tab, in the "Item properties", I linked the small frames to the relative area in the general map, which now is also surrounded by a frame. I would now like to know if there is a way to create a visual connection between these two frames (e.g. an arrow connecting the frame on the general map to the small frame with the zoomed-in portion) without having to draw it manually. I can not find a function for that.


texts: Understanding regression function output in ArcGIS Kriging?

texts:
Below is "regression function" output in ArcGIS Kriging.

I'm assuming the first value is my r2, but does the 2nd value represent error rate, residual, or something else?


texts: How to remove a layer and its legend in a standalone pyqgis application?

texts: In continuity in my question How to add a legend to a canvas in a standalone PyQGIS application? > In every addition of a new layer in my standalone application, the legend lengthens and becomes overloaded. So I would like to know how to remove one or several layers with his(their) legend(s). Sorry for my english.


texts: Conversion of tilestache cache into mbtiles

texts:
I have read and I am looking for someone with documentation on how to achieve converting a cache generated by tilecache into mbtiles using mainly mbtiles. I see the mbtiles has an import function but it requires a metadata.json which tilecache does not have and also it seems the --scheme option has been depricated.

I am trying to use mb-util -i tiles world.mbtiles and hopefully I would like to open the mbtiles in qgis


texts:
Just found out that QGis 3.x has a metadata editor. That is great! Is there a tool to convert the metadata saved by QGIS (.qmd file) to a ISO-xml standard metadata file?

Why am I asking this: I though of using QGis to create the xml metadata file that we will use to populate a Metadata Catalog


texts: Zooming to feature not in center using OpenLayers 3?

texts:
inside my webgis application with Ol3 I've to zoom to a single feature.

My problem is that I've a panel at the bottom of my screen over the Map div. So the zoomed features often finished half hidden by this panel. The result is like this:

Obviusly I want to see the whole feature when I zoom.
I tried this dirty workaround:

map.getView().fit(feature.getGeometry(), map.getSize());
var geom_centre = map.getView().getCenter();
var lon = geom_centre[1] - 50;
var north_centre = [geom_centre[0],lon];
map.getView().setCenter(north_centre);
map.getView().setZoom(map.getView().getZoom() - 1);

Essentially I made the "classical" zoom and then I move to north the centre of the view and go back of one zoom level.
It doesn't work well because the features have different dimensions and I move them to north of arbitrary length (50) and the smallest features are pushed outside view to north.

How can I zoom to a feature and move it to north without push it outside the view?


texts: Unable to install QGIS 2.2 successfully?

texts:
I'm having problem in installing QGIS 2.2. (64 Bit).

After downloading & installing, when I click the OSGeo4WShell, I get a message:

  C:\***.exe not a valid win32 application

and ends with:

  Access denied

What should I do?


texts: adding field in multiple shapefiles and adding oid

texts:
I am new to Python and programming in general. However, I was trying to find a solution to iterate on multiple (>700) shapefiles in one folder. Each shapefile needs a new field with an ID.
I was able to insert the new field, but not perform the sequentially numbering.

I am working in the QGIS console (QGIS 2.2.0-Valmiera).

from PyQt4.QtCore import QVariant
import glob

def fillfields(list, idx):
    i =0
    for item in list:
        new_values = {idx:item }
        provider.changeAttributeValues({i: new_values})
        i += 1

directory = r'D:\mypath'

files = glob.glob(directory+r'\*.shp')

for layer in files:
    layerLyr = QgsVectorLayer(layer,'test','ogr')
    layerLyr.dataProvider().addAttributes( [ QgsField("oid", QVariant.Int) ] )
    layerLyr.commitChanges()

    provider = layerLyr.dataProvider()
    idx = layerLyr.fieldNameIndex('oid')
    n = layerLyr.featureCount()
    field1 = range(n)
    fillfields(field1, idx)
    layerLyr.updateFields()


texts: I need a tool which displays the sum of selected lines as I go in QGIS. Something like the Measure line-tool. I´m new to PyQGIS. The logic seems quite simple, I just don´t know where to start. I guess this could be a good starting point https://gis.stackexchange.com/a/211162/161290?


texts: Creating polyline feature layer for min. dist. between two point feature classes ArcGIS Pro

texts:
I have been struggling to figure out how to build a Python script with the following utility using ArcPy for ArcGIS Pro: Taking two user-supplied point feature classes (A and B) and creating a third, polyline feature class C, that contains the lines connecting each point in feature class A to the nearest point in feature class B, while also recording the distance between the points in the layers’ properties table.

I am new to Python programming in this context.


texts: How to apply new WFS-Filter features like intersect

texts:
I was pleasantly surprised to find very much enlarged filter support in QGIS-Dev (1.9.0_Master).

One application would be a WFS-Filter which combines $geometry and Intersect in order to load only the data from WFS which would be visible in the actual extent.

Something like this:
intersects( $geometry, geomFromWKT( 'LINESTRING(765145.88 6274561.22, 776031.67 6284189.52)' ) )

Found that one in
http://www.3liz.com/blog/rldhont/index.php?post/2012/11/06/Spatial-operators-in-QGIS-expression

But the above mentioned statement is not a valid expression in the "Expression string builder".

What is the correct way to phrase it?

Thanks


texts: I have my openlayers-based website, and we are planning to add a timezone map. We want to show a clock with the current time on top of each area. Is that even possible with Openlayers and Javascript?


texts:
I use some django project with PostGIS for database and leaflet map.
I have add simple leaflet popups like this :

and here the snippet code to call popups :

onEachFeature:function (feature,layer) {
    layer.bindPopup(feature.properties.kalcode.toString())

}

but because I want to show more details from my fields I want to add one more html popups left down of map  popups like this image :

Any idea what change I need about this tasks ?


texts: How to speed up a query for nearest linestring to point data

texts:
I have processed a Strahler stream order layer with a lot of river linestrings (type: Multilinestring) and their corresponding strahler number (see information here). Besides I have a point layer.

I want to know which linestring is the nearest to the distinct points from the point layer to extract the strahler number for every point.

My query here is very slow. For around 500 points and 2000 linestrings it takes over 6 minutes:

SELECT DISTINCT ON
    (lname)
    lnumber,
    lname,
    strahler,
    min(ST_Length(geom)) AS distance
FROM (
    SELECT
        ST_MakeLine(ST_ClosestPoint(strahler_streams.geom, landslide.geom), landslide.geom) AS geom,
        strahler,
        lname,
        lnumber
FROM landslide, strahler_streams) AS foo
    GROUP BY strahler, lname, lnumber, geom
    ORDER BY lname, distance;

I know there is the Indexed Nearest Neighbour Search in PostGIS.

But when I try this...

SELECT
    landslide.lnumber,
    strahler_streams.strahler
FROM landslide a, strahler_streams b
    WHERE lnumber=2114
        ORDER BY b.geom <#> a.geom LIMIT 1

... it only works when I select one single point (WHERE lnumber=2114)

How can I speed up the query for the nearest neighbor search (linestring, point)?


texts: mapbox - how to generate a random coordinate inside a polygon

texts: I've seen the leaflet plugin for mapbox which finds the polygon within which a specific point / coordinate lies, but how would I do the opposite - generate a random point / coordinate within a specific polygon?


texts:
I have a lot of historical orthophotos, like the year 1950 of a region

I have also of the same region the orthophotos of more recent years such as 2012

The question is how to compare and quantifying the geometric quality of both orthophotos?
For knowing what is the geometric quality of the historical orthophotos


texts:
I'm hoping to find out the methodology for mapping increases and decreases in point based occurrences from multiple tables.

Eventual output: A map displaying point density/kernel density increases and decreases of event/location based burglary calls for Police service over three years. (2012, 2013, 2014).

More basically, in a given city, where did burglary calls for service increase or decrease over the 3 years?

The data I have are the three separate years with XY coordinates for the burglary occurrences in excel .csv format.

Field calculating the attributes between two tables is straight forward enough but I wasn't sure how to account for or calculate for an additional table.


texts: How can I modify DEM elevation data within QGIS using a polygon?

texts:
How can I modify an existing DEM using a polygon?  I have a DEM raster with terrain varying between 355 and 365 meters.  I have a polygon I've drawn over the area.  Ideally, I would like to uniformly reduce the DEM by 1 meter from its current height within the polygon (so new values will range from 354 to 364 meters).

I found this previous question and answer, however the workflow seems to have changed as I'm having trouble following the answer in QGIS 3.4.  Also, that answer addresses how to alter the DEM to have a new absolute height based on the elevation values of the polygon, but now how to create new height values relative to the old height values.  The reason I want to do this is because the terrain in question is on a slope, so I cannot simply assign a single elevation value to the entire area bounded by the polygon.


texts: The ordering of layers with CartoDB/Leaflet

texts:
Just want to get a bit of clarity about this query from before:

Layer order in CartoDB

I have a map with a Basemap and I want to draw points between features based on layer featureclicks, The good news is I've got that far.

However i'd like my newly created layer to sit between the marker layer and the tile layer. Is this even possible?

Heres a link to my code on github


texts:
I have a script that tries to read from the arc environment using

import arcpy
tolerance = arcpy.env.XYTolerance

as recommended here:  http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//002z0000000w000000

Alas tolerance always ends up as None.  Even when I set it explicity in the Environments... part of the script dialog box.

Is this a bug?


texts: By default zoom to, setview functions show the center point of the map in the center point of div measured by pixel of browser. I want to show the center at slightly right when left panel is visible. In http://2gis.ru/ site, they used the similar function. How can I do this?


texts: What is CID field created by Create Random points?

texts: In ArcGIS the result of Create Random Points tool when using a feature class to constrain the extents has an attribute called CID. What do the values in this field represent?


texts: Find borders between all unique country pairs with R

texts:
I would like to obtain a border as a line between all unique country pairs that are contiguous. Here is how to do it:

library(spData)
library(sf)
library(rnaturalearth)
library("rmapshaper")

#load countries dataset:
countries_sp <- ne_countries(scale = 10, continent= "Europe")

# find borders
border_sf <- st_as_sf(countries_sp) %>%
  ms_innerlines() %>%
  as_tibble() %>%
  st_as_sf()

# create ids to visualise lines found by ms_innerlines
border_sf$ids <- rep(1:length(border_sf$geometry),1)

# plot
tm_shape(countries_sp, bbox = c(c(-8.750803, 36.14111), c(40.15954, 70.07531))) +
  tm_fill() +
  tm_shape(border_sf) +
  tm_lines(col = "ids")

Example shows how to find have one line for the border between France and Germany, one for France and Switzerland, France and Italy for each contiguous country pair in the sample.


texts:
How do arcpy searchcursors sort records if there is no 'sort_field' specified?

It seems, on the surface, to sort ascending by the primary key. However, I don't want to base my script around this idea, and then after 6 months of using it find out that its not true for some reason.


texts: Kartograph - converted svg not the same as before

texts:
I'm using the following script to convert a map from NaturalEarthData. But the result is very sobering. It seems like something went wrong with the scaling, but I got no errors when executing the script. What am I doing wrong??

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from kartograph import Kartograph
K = Kartograph()

config = {
    "layers": {
        "mylayer": {
            "src": "ne_110m_admin_0_countries.shp",
            "charset": "utf-8",
            "attributes": "all"
        }
    },
}

K.generate(config, outfile='world.svg')


texts: Batching model output problem in Arcgis 10.2!

texts:
I built a simple model for projecting raster from one coordinate system to another and it works fine but when i try to batch it, so i can convert more then one raster it generates only the first raster, in tiff, and on the same directory rather than the one i specify in the model using values and path.

And when i use iterate raster it continues forever and generate only the raster one, the second and then that's it, it continue endlessly and it keeps deleting the first raster.

This is the other model

and the iterate parameters

Project raster parameters and output parameter


texts:
I have a vector layer in qgis connected to PostgreSQL database table.

Suppose I encrypt the table column using pgp_sysm_encrypt or some other encryption in PostgreSQL

Is there any way i can say qgis to decrypt every column using my key


texts: cannot get rid of sr.lock files while FOR LOOPING dbfTables

texts:
I am new to Python and to this community as well, hope someone could help me.
I am running  script to extract and merge toghether the tables of some Raster. Each Attribute Table presents only one row and the same field names.
My script works fine until the very end of it, when I try to merge all the tables. It seems the tables are locked (I can see an "sr.lock" file for each Table I want to merge in my folder) and I cannot access them.
Strange thing: I do access them if I simply use a sample (for instance) of two of these so created tables in the Python window within arcMap and type their names instead of using a FOR LOOPING... Guess it might have to do with the
Here is my script:
    #import relevant modules, create geoprocessing dispatch object
    import win32com.client, sys, string, os, arcpy

gp = win32com.client.Dispatch("esriGeoprocessing.gpDispatch.1")

# Remember to change this to wherever your files are stored
gp.workspace = "G:\\Chile_2012\\MODIS10A2.V005\\provvisoria\\ATTRIBUTE_TABLE"
out_dir = "G:\\Chile_2012\\MODIS10A2.V005\\provvisoria\\ATTRIBUTE_TABLE"

try:
    rasters = gp.ListRasters("*", "all")
    rasters.reset()
    rst = rasters.Next()

    while rst:
        # Create the new field
        gp.AddField_management (rst, "FILENAME", "text", "", "", "50")

        # Apply the filename to all entries
        gp.CalculateField_management (rst, "FILENAME", '"' + rst + '"')
        rst = rasters.Next()

except:
    print gp.GetMessages ()

listRaster = arcpy.ListRasters ()

for raster in listRaster:
    arcpy.TableToTable_conversion (raster, out_dir, raster+ '.dbf')

listTable = arcpy.ListTables ()

for table in listTable:
    arcpy.Merge_management (table, 'Cl_Rio_Aconcagua.dbf')

The code is taken from others, I simply adapted it for my needs. Anyway, the error I am receiving so far doesn't have to do explicitly with LOCK topics, but since I have those files in my folder I need help to find out how to get rid of them.

ERROR 000732: Input Datasets: Dataset 'Cl_Rio_Aconcagua_2007121.img.dbf' does not exist or is not supported
Failed to execute (Merge)


texts: Postgresql ST_Point problem (Can't see point symbol on map)

texts:
i am using Arcgis for Server and Arcgis Portal webmaps with Postgresql db.
Created a Feature Class with ArcCatalog and shared as Feature Service, so i am trying to add data to this Feature Class with a Sql Query. (Geometry column name "shape")

i made success adding data with these queries;

INSERT INTO sde.data_table (objectid, name, globalid, shape)
VALUES (21,'dataname',sde.next_globalid(),sde.ST_Geometry('point (622000 4122000)', 32635)
);

or

INSERT INTO sde.data_table (objectid, name, globalid, shape)
VALUES (21,'dataname',sde.next_globalid(),sde.ST_Geometry('point (622000 4122000)', 32635)
);

or

INSERT INTO sde.data_table (objectid, name, globalid, shape)
VALUES (sde.next_rowid('sde', 'data_table'),'dataname',sde.next_globalid(),sde.ST_Geometry('point (622000 4122000)', 32635)
);

My problem is; I can see data with Arcmap's and Webmap created with Portal Attribute Table, but i cannot see point symbol on the map, it is there but no symbol, also i cannot identify with a mouse click. Just like a hidden object.

PostgreSQL 9.4.5 /
Arcmap 10.4 /
Arcgis For Server /
Arcgis Portal


texts:
If I am looking to build a network based on road features and attributes, is Network Analyst my best bet?

Is it my only bet?

The world of FOSS (Free and Open Source) is new and exciting to me, but I find it difficult to navigate as a newbie.

My intended uses would include tracing the network, finding intersecting features at given locations along the network, and creating network decision rules for route creation.

Answers so far have been great, but could anyone share their personal experiences and preferences as far as routing interfaces are concerned?


texts: ArcGIS Raster Calculator won't work in main model

texts:
I've created a sub-model in ArcGIS, which stretches my raster data to new values.

The strange thing is, that the model works, if I run it itself - as soon as I implement it in the main model, I get this error message:

ERROR 000539: Error running expression: rcexec()
Traceback (most recent call last):
  File "", line 1, in
  File "", line 5, in rcexec
  File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\sa\Functions.py", line 244, in Con
    where_clause)
  File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\sa\Utils.py", line 47, in swapper
    result = wrapper(*args, **kwargs)
  File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\sa\Functions.py", line 238, in Wrapper
    where_clause)
  File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\geoprocessing_base.py", line 498, in
    return lambda *args: val(*gp_fixargs(args, True))
RuntimeError: ERROR 010240: Could not save raster dataset to C:\Users\Mattias\Documents\01_UNI\Diplomarbeit\GIS\DATA\data_out\scratch_IV.gdb\plus_ras with output format FGDBR.
Failed to execute (Raster Calculator).

In my opinion, there's nothing wrong with my syntax - else it would not work in the sub-model directly...

Anyway, my map algebra expression is:

Con("%to_be_stretched%","%min_r_raster%",Con("%to_be_stretched%","%max_r_raster%",((("%to_be_stretched%" * 1.0) - "%min_in_raster%") / ("%max_in_raster%" - "%min_in_raster%")) * ("%max_r_raster%" - "%min_r_raster%") + "%min_r_raster%","VALUE >= %max_in%"),"VALUE  <= %min_in%")


texts: edit node coordinate in qgis, moving a point without having to make another layer

texts:
so I wanted to make a script, but I am not good using the qgis api...

To sum up the issue I have this :

I want that the 3 points of the yellow polygon are snapped by the closest point of other layers (example an orange point will snap to a point of the yellow layer).

Right now in my script I have the list of all the points coordinate, list of the selected figure, I found out the closest distance. But I am stucked on trying to move the points of the other layers to make them snap to the nodes of the yellow layer... And I don't want to make new layers to do this and I don't want to do it manually, I want to edit the existing coordinates of the points of the other layers...

Here is my script right now :
https://pastebin.com/KYjFdEyn

I found out on the qgis api that there is method setgeometry but I have no idea how to use, even when I read the doc... Unless I am looking at the bad thing, I am not even sure it will be able to do what I want.

Does anyone have an idea on how to edit the coordinate of points of a layer ?


texts: Calculate NDVI by avoiding use some unusual values of surface reflectance of NIR and Red bands included in an image collection

texts:
I want to use a imagecollection to calculate NDVI, the following code shows how I did. However, I know some values of surface reflectance of Near-infrared band and Red band are less than 0, or larger than 10000, or is Null, so how to calculate the NDVI by avoiding use the Near-infrared band and Red band with those unusual values?

var collection = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')
    .filterDate('2019-01-01', '2019-12-31') // the period of interest
    .filterBounds(fc) // the study area
// B4 is near-infrared band and B3 is red band
var NDVI1 = function(image) {
  return image.expression('float(b("B4") - b("B3")) / (b("B4") + b("B3"))');
};
var NDVI = collection.map(NDVI1)
Map.addLayer(NDVI, {}, 'ndvi')
Map.centerObject(fc)

this is the code link:
https://code.earthengine.google.com/202776c3d784fed6c7a04947f6b769f6


texts: Error installing QGIS Open Layers plugin on 2.14 and 2.18

texts:
When attempting to install the plugin Open Layers I get the following error

"Plugin installation failed:Failed to unzip the plugin package. Probably it's broken or missing from the repository. You may also want to make sure that you have write permission to the plugin directory: C:/Users/username/.qgis2/python/plugins"

I have confirmed that the user has write permissions to the directory and have also tried to manually install the plugin to the folder. This is on a windows 10 computer.The plugin will successfully install on another machine that is using windows 8.

Has anyone else experienced this issue? If so, what is the solution?


texts: Split polygons from a shapefile according to other polygons in another shapefile keeping the overlapping polygons separated in Arcgis

texts:
I have two shapefiles: the first with a series of polygons that overlap; the second with a series of polygons that don't overlap. I need to separate the first ones according to the line of the second, keeping the overlapping polygons separated.

Graphically:

This is the initial situation: the shapefile contains several polygons that overlap and each polygon can be crossed by a line of a polygon that is in another shapefile (the red rectangles in the picture)

What I need to obtain at the end are the old polygons splitter according to the red one. In the example each polygon will be divided in two parts:

1a will be the yellow on the left + the green part;
1b will be the yellow on the right + the orange part;
2a will be the grey on the left + the green part;
2b will be the grey on the right + the orange part.

What I tried is to use the union tool and then the dissolve tool: with the first one I obtain six different polygons (one different polygon for each line that cross, no matter if it comes from the red polygon or the overlapping one).

Then I tried to dissolve indicating as "dissolve fields" the FID of the original two layers (FID_shp1 and FID_shp2) but basically nothing happened.

Any suggestion? should I use something different? I mean, is there a way to avoid that the overlapping lines became a separator and only the "red lines" are taken into consideration?

EDIT 1: Just to make the comprehension easier. This is the initial situation:

Each of this polygons can be crossed by other polygons that are in another shapefile

I need to arrive at the end at this situation

i.e. four different polygons that basically are the initial polygons divided by the red line of the polygons of the other shapefile.

The problem with the union tool are two-fold:

The union tool take all the lines as a separator, included the ones of the overlapping polygons;
Once I unified, the attributes of the second file are not copied but it will be transformed in something different (for example, all the FID of the second shapefile became -1. That's why maybe the dissolve tool doesn't work properly).


texts: Creating a new raster layer based on mean of other raster layers, ignoring NA values in R

texts:
I am trying to create a new raster layer that is made up of the mean values from 7 other raster layers. I get results using overlay function mean:

mean <- overlay(2000, 2001, 2002, 2003, 2004, 2005, fun=mean)

The results I get however are skewed by a large amount 0 values in many of the cells. As the  data within the rasters changes each year some cells could have 0 values one year and a value of 56 the following, for example any particular cell may contain the following:

+----------+---------+
|    Year  |  Value  |
+----------+---------+
| 2000     | 0       |
| 2001     | 54      |
| 2002     | 78      |
| 2003     | 102     |
| 2004     | 0       |
| 2005     | 0       |
+----------+---------+

So while the mean is:39, I am trying to achieve the mean ignoring the 0 values so the result should be: 78.

I have tried changing the 0 values to NA:

2000[2000 == 0]<-NA
2001[2001 == 0]<-NA
2002[2002 == 0]<-NA
2003[2003 == 0]<-NA
2004[2004 == 0]<-NA
2005[2005 == 0]<-NA

Which works fine, however the bit that I can't get to work is the calculating of the mean as before but with ignoring these NA values. I have tried

mean <- overlay(2000, 2001, 2002, 2003, 2004, 2005, fun=mean,  na.rm=TRUE)

but the results this returns are incorrect. So I am just looking for a way to find the mean whilst ignoring zero values or NA values.


texts: How do I compare multiple years worth of geomorph changes in a river while drawing one profile line?

texts:
I have many years worth of DEMs of a river I'm studying.

I'm very, very new to ArcGIS (think one intro class), and learned how to use the 3D analyst tool.

What I'd like to be able to do, is to draw one profile, and from that profile line, show graphs of the DEM changes for each year.

What was previously done, was a line for each year was drawn as close as possible to original line of year previous, but it's not directly on top of previous year.

Is it possible to stack these profiles and compare multiple years in the same graph?


texts: Using Virtual Raster File Viewer of QGIS?

texts:
I have created a virtual file using Build Virtual File functionality of QGIS.

Is the file linked strictly to where the raster files are stored or can it be viewed on another computer?


texts: Are there Plugin Compatibility Issues between QGIS 1.6.0 and 1.7.4?

texts:
I've been working on a QGIS plugin on Quantum GIS 1.6.0 Copiapo with Python 2.5 that calculates additional values for eachf feature of a vector layer and adds them to the .dbf. We've just upgraded to QGIS 1.7.4 Wroclaw with Python 2.7. But since the upgrade the plugin doesn't work anymore in 1.7.4.
While on 1.6.0 all functionalities work fine, running the plugin in 1.7.4 actions on a number of QPushButtons and some other functionality within the ui do not trigger any actions. I've tried both the old and the new-style connections, like:

QObject.connect(browseCalcButton, SIGNAL('clicked ()'), self.browseCalcFunction)

browseCalcButton.clicked.connect(self.browseCalcFunction)

This seems to be true for all QPushButtons that DON'T either call...

... a handler:
attributesButtons[i].clicked.connect(self.createHandler_assign_attribute(i, key))
... an existing pyqt method:
eA_Fertigstellen.clicked.connect(self.eAui.hide)

... but rather a function within the script:

browseCalcButton.clicked.connect(self.browseCalc)

Are there compatibility issues between QGIS 1.6.0 and 1.7.4?
Are there any adjustments neccessary?

Has anyone had similar problems?


texts: Using azimuthal equidistant projection with buffers in QGIS?

texts:
I am a journalist.

My goal: A globe that shows the range ok North Koreas ICBMs - the azimuthal equidistant projection centered on NK.

My steps so far:

Created Standard QGIS project (EPSG:4326)
Imported a world shapefile (EPSG:4326). The center of the shapefile is in Greenwhich.
Imported a CSV with the lat/long of Pyongyang
Used MMQGIS to create a buffer (13,000 km)
Created a Custom CRS with these details:

+proj=ortho +lat_0=39.03385 +lon_0=125.755504 +x_0=0 +y_0=0 +a=6371000
+b=6371000 +units=m +no_defs

Didn't look so nice (and since the US was not visible it was useless for me)

Created a different Custom CRS with these details:

+proj=aeqd +lat_0=39.03385 +lon_0=125.755504 +x_0=0 +y_0=0 +a=6371000
+b=6371000 +units=m +no_defs +proj=aeqd  +R=6371000 +lat_0=39 +lon_0=125

My problem:

Being not a GIS pro, I would be quite proud of myself - if there wasn't a part of the buffer (=range) missing - exactly at the line where my original shapefile "ended."

For the record, I checked posts like:
 - Map Range of Airplanes
 - Manipulating Azimuthal Equidistant Projections in QGIS
 - How do I project a Azimuthal Equidistant global map into a circle in QGIS?

I also tried to first create a new shapefile centered on the pacific using QGIS display world country shape files centered on pacific ocean using Robinson, Miller Cylindrical or other projection but that also did not work.


texts: Leaflet markercluster - cluster/uncluster toggle

texts:
Is it possible to toggle between clustered and unclustered with a leaflet MarkerClusterGroup?

I know I can disable clustering at max zoom levels using the disableClusteringAtZoom param below

markerHolder = new L.MarkerClusterGroup({
    disableClusteringAtZoom: 18,
    maxClusterRadius: 80,
    spiderfyDistanceMultiplier: 1,
})

I would like to add a button that clusters/unclusters on click


texts: Viewing Domain values as dropdown list in data grid in arcgis javascript api

texts:
I'm trying to view layer data in a data grid, everything is fine, but when the data contains domain values, its place stays empty, and does not view any data.

Can anyone can help how to view it or am I supposed to build the table manually?


texts: How to disable spatial index creation by ArcSDE for Oracle Spatial

texts:
Working with ArcSDE 10.1 and Oracle spatial 11g.

I set the feature class in load_only_io before loading shapefiles. I am using geoprocessing tool (append) in ArcGIS to load shapefile to geodatabase. After loading ArcGIS is creating spatial index on sdo_geometry.

How to disable this automatic index creation?

I want to load shapefiles and later I want to create index myself.


texts: Showing two layers at same time in QGIS?

texts:
I want to see two QGIS layers at the same time.

one is lines which are boundary lines and the other is polygon which are district.

when i was surfing about this question, it is impossible to show two layers with different type.


texts: QGIS GeoServer response "Server replied: Forbidden"

texts:
I have the following problem. I have created a WCS with Mapserver (MS4W) and would like to integrate it into QGIS. The URL request in the browser via KVP works perfectly. However, I get the error message with QGIS: "Answer not understood: (...) Server replied: Forbidden"

What could be the Problem?

The Error Logs in Apache, Mapserver und QGIS have no clues.


texts: ArcPy Network Analyst solve function, shows 'ExecuteError: Failed to execute. Parameters are not valid. ERROR 000824: The tool is not licensed.'

texts:
I want to find the shortest path between two points. I first make a route layer from an analysis network and then add two point locations.
While python shows the following errors when I use the solve function. Can anyone tell me the reason?
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000824: The tool is not licensed.
Failed to execute (Solve).
Code works when I use Python window in ArcMap and I actually have all ArcGIS licenses
codes:
'''
import arcpy
import os

arcpy.env.workspace = r"C:\Users\hyan1\Desktop\try\network\network_for_wageningen.mdb"
network = r'C:\Users\hyan1\Desktop\try\network\network_for_wageningen.mdb\network2_ND'
outNALayer = arcpy.na.MakeRouteLayer(network, "WorkRoute", "length")
outNALayer = outNALayer.getOutput(0)
subLayerNames = arcpy.na.GetNAClassNames(outNALayer)
points_layername = subLayerNames["Stops"]
routes_layername = subLayerNames["Routes"]

### add location, solve and export
input_point = r'C:\\Users\hyan1\Desktop\Route_2\T11_1_11_pro_lyr_1.shp'
output_name = os.path.split(input_point)[1]
arcpy.na.AddLocations(outNALayer, points_layername, input_point, "", "", "DateTimeS")
arcpy.na.Solve(outNALayer)
arcpy.FeatureClassToFeatureClass_conversion("WorkRoute\Routes", r"C:\Users\hyan1\Desktop\Route_2\shortest_route", output_name)

'''


texts: Plotting xarray over background map using cartopy

texts:
I am trying to plot an xarray DataArray over a background map using cartopy.
When including the background map, the DataArray is not visible, although the map has an appropriate extent.
I am adding the background map using cimgt and geo_axes.add_image().
This method has worked for me before when plotting points with plt.scatter().
The data array has dimensions x (100 values) and y (140 values), and contains np.nans and floats.
Here is my code:
# DataArray creation
# ...

# create figure
fig = plt.figure()

# create geo axes
projection = ccrs.epsg(32630)
geo_axes = plt.subplot(projection=projection)

# add open street map background
# when commenting the two following lines, the data array is plotted correctly
osm_background = cimgt.OSM()
geo_axes.add_image(osm_background, 10)

# plot dataset
xr.plot.imshow(
    darray=dataset.velocity.mean(dim="time"),
    x="x",
    y="y",
    ax=geo_axes,
    transform=projection,
)

# show plot
plt.show()

I can't seem to find what is causing the issue. Has anyone encountered this before, or can anyone see what I am doing wrong?


texts: How to access pgRouting Layers via script

texts:
I would like to provide input to the pgRouting Layers through a script.

The scenario is as follows:

I have a bunch of points whereof every point is representing a starting point and I have another bunch of points representing destinations.

Each of the start points need to be routed to every destination point (shortest path) and I want to retrieve the length of the route (distance), and if possible the travel time.

I can iterate over the start, and destination points, but how do I get it done that I feed the x/y data to pgRouting Layers and extract the required information from the return values?


texts: Checking ArcSDE Compatibility with SQL Server

texts:
I am trying to find if I can use ESRI sde 10.5.1 with a SQL Server 2016 SP1.

I am not able to find this information. What I've found so far:

On the Microsoft SQL Server database requirements for ArcGIS 10.5.x page, it states that the 10.5.x version of ArcGIS is compatible with SQL Server 2016.
The ArcGIS Server 10.5.x system requirements page sends the user to the aforementioned page to get the info.
The Supported Environment Policy says that they usually support the main environments

Where or how can I find out if these two are compatible ?

(Note: The preferred answer would give the place where I can find this information, so that other users with different versions could use it too)


texts: Recentering a tiled raster table in PostGIS?

texts:
Does anyone know a nice and easy way to recenter a tiled raster table in PostGIS?

By recentering I mean that I have a raster that is centered on the Pacific, but I need to recenter it on Greenwhich/Europe. In a sense the effect would be to offset pixels while "rolling" pixels that fall off the right edge back around to the left edge. It would be easier if the entire raster data was just a single big raster, but due to large memory size I have chopped up the raster into several subrasters/tiles with the raster2pgsql -t option (which is the PostGIS approach to raster data), which means the operation would have to be done on a tile by tile basis making the solution slightly more challenging.


texts: Select by attributes to create feature classes

texts:
I am using OSM point data and am trying to pull out schools, hospitals and place of worship. I have created my feature layer and am using a "for loop" to make from each amenity its own shapefile.

What I currently have will create the 3 shapefiles all with hospitals populating every shapefile. Do I have to create a new query for each amenity and then make a new feature layer from that?

amenities = ['school','hospital','place_of_worship']


texts:
There are lots of comments on different forums about different types on online GIS data/map sharing tools.

What I need to know is what ArcGIS online CAN'T do compared to the functionality of the open source options?

I have big geodatabases that contain dynamic data that is constantly being updated and shared across multiple offices.  The geodatabases contain raster data with information that needs querying - Lidar/bathymetry data.  We also have vector features that could be as many as 2000 features in any one feature layer.

I have heard that it is a struggle to get raster layers into AGOL and they will just be an image with no data behind them

I have also heard that there is a limit on the number of features AGOL can display at any one time?

It seems a total minefield as to which one to use but I would prefer the free option if the functionality is the same or better in opensource.


texts: How do I get zipcodes from longitude and latitude on python?

texts:
I have a data frame of latitude and longitude coordinates on this CSV file:Longlat.
I used this code to try to get the zipcodes:
import copy

def get_zipcode(df, geolocator, lat_field, lon_field):
    location = geolocator.reverse((df[lat_field], df[lon_field]))
    return location.raw['address']['postcode']

geolocator = geopy.Nominatim(user_agent='myusername') #My OpenMap username

zipcodes = longlat.apply(get_zipcode, axis=1, geolocator=geolocator, lat_field=longlat['LATITUDE_X'], lon_field=longlat['LONGITUDE_X'])

I got the error:
KeyError
"None of [Float64Index([39.0962320000896, 39.1462010000896, 39.1347670000896,\n 39.1076250000897, 39.0928490000897, 39.1648900000896,\n 39.1846440000895, 39.0970790000897, 39.1491220000896,\n 39.1145560000896,\n ...\n 39.1039560000896,

How do I fix it?


texts:
I want to offset my label appearing in the interactive map generated by QGIS2web plugin.
Unfortunately the variables offsetX and offsetY don't work...
 var size = 0;
 var placement = 'point';

 var style_Area_5H_Schematic_7 = function(feature, resolution){
 var context = {
    feature: feature,
    variables: {}
 };
 var value = ""
 var labelText = "";
 size = 0;
 var labelFont = "20px \'MS Shell Dlg 2\', sans-serif";
 var labelFill = "#000000";
 var bufferColor = "";
 var bufferWidth = 0;
 var textAlign = "top";
 var offsetX = 60;
 var offsetY = 75;
 var placement = 'point';
 if (feature.get("Name") !== null) {
    labelText = String(feature.get("Name"));
 }
 var style = [ new ol.style.Style({
    image: new ol.style.RegularShape({radius: 15.6 + size, points: 8,
        radius2: 7.8, stroke: new ol.style.Stroke({color: 'rgba(9,4,2,1.0)', lineDash: null, lineCap:
  'butt', lineJoin: 'miter', width: 1}), fill: new ol.style.Fill({color: 'rgba(9,4,2,1.0)'})}),
    text: createTextStyle(feature, resolution, labelText, labelFont,
                          labelFill, placement, bufferColor,
                          bufferWidth)
   })];

   return style;
  };

I've changed these values to:
 var offsetX = "6";
 var offsetY = "8";

but still the same. No result at all.
Does anyone know how to solve it?

UPDATE:
Labels are not editable for the following code:
 var size = 0;
 var placement = 'point';

 var style_Area5chamber_14 = function(feature, resolution){
 var context = {
    feature: feature,
    variables: {}
 };
 var value = feature.get("Chamber Type");
 var labelText = "";
 size = 1;
 var labelFont = "15px, sans-serif";
 var labelFill = "#787878";
 var bufferColor = "";
 var bufferWidth = 0;
 var textAlign = "left";
 var offsetX = 0;
 var offsetY = 0;
 var placement = 'point';
  if (feature.get("Chamber No") !== null) {
    labelText = String(feature.get("Chamber No"));
}

  var style = categories_Area5chamber_14(feature, value, size, resolution,
   labelText,
                      labelFont, labelFill, bufferColor,
                      bufferWidth, placement);

    function categories_Area5chamber_14(feature, value, size, resolution,
  labelText,
                   labelFont, labelFill, bufferColor, bufferWidth,
                   placement) {
             switch(value.toString()) {case 'BT':
                return [ new ol.style.Style({
         stroke: new ol.style.Stroke({color: 'rgba(35,35,35,1.0)', lineDash: null, lineCap: 'butt', lineJoin: 'miter', width: 0}),fill: new ol.style.Fill({color: 'rgba(0,0,0,1.0)'}),
    text: createTextStyle(feature, resolution, labelText, labelFont,
                          labelFill, placement, bufferColor,
                          bufferWidth)
})];
                break;
      case 'Voneus':
                return [ new ol.style.Style({
    stroke: new ol.style.Stroke({color: 'rgba(35,35,35,1.0)', lineDash: null,
     lineCap: 'butt', lineJoin: 'miter', width: 0}),fill: new
    ol.style.Fill({color: 'rgba(252,0,0,1.0)'}),
    text: createTextStyle(feature, resolution, labelText, labelFont,
                          labelFill, placement, bufferColor,
                          bufferWidth)
  })];
                break;}};

   return style;

   };

where I have 2 cases.
The label comes out only for one of them, but apart from the font color remains not editable.

is there any way to sort it?


texts: Adding custom names to a GeoPandas legend

texts:
I have a shapefile which has an attribute table with a column I would like to make a map/plot of. The attribute values are numerical (integer). I have made two dicts to map the colors and names I want to these integers.
Palette = {0: 'black',
       20: '#FFBB22',
       30: '#FFFF4C',
       40: '#F096FF',
       80: '#0032C8',
       90: '#0096A0',
       112: '#009900',
       114: '#00CC00',
       116: '#007800',
       124: '#A0DC00',
       126:'#648C00'}

names  = {0: 'NAN',
       20: 'Shrubs',
       30: 'Herbaceous',
       40: 'Cultivated',
       80: 'Permanent Water',
       90: 'Herbaceous Wetland',
       112: 'Closed Forest: Evergreen',
       114: 'Closed Forest: Deciduous broad leaf',
       116: 'Closed forest: Other',
       124: 'Open forest: Deciduous broad leaf',
       126:'Open forest: Other'}

However, while I can map the colors to the right values, I cannot get the legend to show the right names. The legend comes up empty and I get a message that No handles with labels found to put in legend.
My code is:
fig, ax = plt.subplots(figsize=(5, 5))

# Loop through each attribute value and assign each
# with the correct color & width specified in the dictionary
for ctype, data in map_df.groupby('landcovermode'):
    color = Palette[ctype]
    label = names[ctype]
    data.plot(color=color,
      ax=ax,
      label=label,legend=True)

# Place legend in the lower right hand corner of the plot
ax.legend(loc='lower right',
  fontsize=15,
  frameon=True)

ax.set_axis_off()

How do I get the legend to read my labels from the dict?


texts: How to get city data only and not county?

texts:
Most of the data I've come across only comes in state/county/block group/block format. However, I'd like to pare down the shapefile to only observe the block groups within a city, not the county of San Diego. I have the Census ID of the block groups but don't know how to limit it to city data only.

How do I know which census tracts belong within the 2000 Census definitions of San Diego city? Suggestions?


texts: What spatial functions are available to shapefiles as "virtual layers" in QGIS DB Manager?

texts:
Can anyone provide a list of the spatial functions available in the QGIS DB Manager plugin when using shapefiles as Virtual Layers?

I've been trying to use things like ST_DWithin and even ST_Buffer but the queries fail - though ST_Intersects and other functions work just fine.


texts: What does "You do not have necessary license to run tool" mean when ArcGIS extension installed and activated?

texts:
I just purchased ArcGIS for Home Use 10.2 and installed and activated the extensions with the codes I received from Esri. The ArcGIS Administrator says Yes for Installed and Authorized, and even has an expiration date.

However, when I try to run the tool it states "You do not have the necessary license to run the tool" and that there is no license available or installed.

Am I missing something ?


texts: QGIS 3x rotation of the label, when another expression has been set

texts:
Following this query:
QGIS 3x conditional rotation of the label
I would like to rotate one label, which already has been assigned to one expression.
The previous expression looks like this:
 CASE
 WHEN "fid" = 18 then '5,-5'
 End

and applies to the offset for a whole rule-based label.
I tried to make something like this:
  CASE
  WHEN "fid" = 18 then '5,-5'
  End

  CASE WHEN "fid" = 17 then 60
  end

but the expression isn't valid.
I modified it a bit:
  CASE
  WHEN "fid" = 18 then '5,-5'
  WHEN "fid" = 17 then 60
  End

but it returns no reaction at all.
How can I change the rotation for this one label particularly?
I tried also via geometry generator:

but all my labels belonging to the rule disappeared.
The last way was here:

but it looks like I need some auxiliary columns in my data attribute table.
Can anyone help me to solve it?


texts:
I’m trying to import OpenStreetMap (OSM) data into postGIS it doesn’t appear to be working in that features do not appear to be in the correct place..

Using the QGIS  OSM add in  and a coastal shape file imported into PostGIS setting the SRS to 4326 on import (all my data uses WGS84 Longitude & Latitude) the two data sets are in the same space.

Importing  the OSM data set  into postgis with the following code & adding it as a QGIS layer I would expect lines points etc of the OSM file & PostGIS layers to be adjacent if not on top of each other but they are not,
I can’t actually see the coastal shapefile & the OSM PostGIS layer together without turning one layer off and selecting zoom to layer on the other.

BTW I have attempted to use the -E 4326  switch with out -m but it failed
So what am I doing wrong?

createDB osm
createlang plpgsql osm –U username
psql –U username –d osm  -f /path/to/postgis1.5/postgis.sql
psql –U username –d osm  -f /path/to/postgis1.5/spatial_ref_sys.sql
psql –U username –d osm  -f /path/to/osm2pgsql/900913.sql

osm2pgsql –S /path/to/osm2pgsql/default.style –U username –d osm  -l –c –m /path/to/ cambridgeshire.osm.bz2

I then tried following the instructions on http://spatialmounty.blogspot.co.uk/2010/12/transform-osm-data-in-postgis-to-wgs84.html to re-project the data to WGS84 Long Lat it failed at The first line with invalid SRID

SELECT AddGeometryColumn( 'planet_osm_line', 'geom_wgs84', 4326, 'LINESTRING', 2);

My database  has the following tables

Geography_columns  view
Geometry_columns table
Planet_osm_line table
Planet_osm_point table
Planet_osm_polygon table
Planet_osm_roads table
Spatial_ref_sys table

Can somebody give me the steps / syntax i need to follow to accurately import OSM data into PostGIS so that it uses ESPG:4326?

Any one know of any ground truth data I can use to test the process?


texts:
I would like to split a line with a smaller line on QGIS in order to have the rest. I used split line with line but it doesn't work.

Sorry I have difficultes to explain because I'm not an English speaker. I have a big line and a smaller line which are on the space. I would to have the subtraction of them. I would like to split lines in a layer with lines in another layer using processing tools. When I try I only have the full segment as a result and not the subtraction of both

This is the bigger line :

This is the smaller line :

I would like the subtraction (the red circle) of them  (in yellow it's still the bigger line and in green it's still the smaller line)


texts: QGIS Raster Calculator, joining two raster

texts:
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...

If i select in the QGIS layer selector the two layers to be added (using blending effects = multiply) i get visually pretty much what I spected to see when making the raster calculation adding:

How could I make a real adding calculation between two rasters?


texts: Automatically deleting lines within generated polygon using QGIS?

texts:
How do I automatically delete lines within a generated polygon?

what I have: is a polygon with lines within. The lines I want to be deleted automatically -- like a auto clean up.

This is what I did:

opened a postcode shapefile from https://www.suche-postleitzahl.org
Selected features (Berlin postcodes)
Merged selected features
Struggling to find a way on how to clean up the polygon in the layer.

Example


texts: R - How to get raster values across a latitude band

texts:
Let's say I have a raster like this:

library(raster)
r <- raster(ncols=15, nrows=15)
values(r) <- runif(ncell(r)) * 1:ncell(r)
plot(r)

I need to:

find the latitude point closest to a given value, say -50. In the
example above, that latitude point would be -48;
extract the values for all the longitude points across that latitude band (i.e. the red line in the figure below).

The desired output would be a data frame with the longitude points and the extracted values.

Any hints on how to achieve that?


texts: Spatial join list lat,lon and multipolygon Python

texts:
I have a dataframe with lat, lon columns in WGS84.

I also have a multipolygon layer (GeoJSON, I can also convert it to a SHP) with all countries' boundaries, where the name of the country is in the attribute table.

I would like to add to the dataframe a column where for each lat, lon I'll have the country name.
What would be the efficient way to do that (assuming I have 2000 lat lon pairs)?
UPDATE:
Thanks to @BERA I've created this function:
import geopandas as gpd
from shapely.geometry import Point

df_countries = gpd.read_file(r"C:\countries.geojson")

def get_countries(df, lat_col, lon_col,df_countries):
    df_latlon = df[[lat_col,lon_col]].copy()
    df_latlon['Coordinates'] = list(zip(df_latlon[lon_col], df_latlon[lat_col]))
    df_latlon['Coordinates'] = df_latlon['Coordinates'].apply(Point)
    df_latlon = gpd.GeoDataFrame(df_latlon, geometry='Coordinates')
    df_latlon = df_latlon.set_crs(epsg=4326)
    df_latlon = gpd.sjoin(df_latlon, df_countries[['CNTRY_NAME','geometry']], how='left')
    return df_latlon


texts: How to filter a vector layer in openlayers?

texts:
I'm using openlayers (v3.20.0) in an angular 2 project.
I have a vector layer which has its source bound to a geoserver.

My task is to filter the features displayed in this layer.  I have a collection Ids (which are in features' attributes) representing the features needed to be shown.
But I can't find how to filter the layer.  There seems to have been an easy way to do this in openlayers v2, but not in v3.

What am I missing ?


texts:
I have a shapefile of a set of points that I want to space out. I DO NOT want to do this using line interpolation as I run the risk of placing a point where there wasn't one before.

I can load this into Postgis or QGIS and would like to make sure that no point is closer than 100 metres to its neighbour.

I guess I need to calculate the KNN value for each point to its neighbours as a starting point but am not sure what strategy to use to remove points to achieve the desired result.

The graphic shows a sample of the data, I am trying to overlay it onto the map to show cycle routes as points on top of the road network. At the moment the points are too dense and as such are obscuring the road below. Snapping to roads IS NOT an option as the tracks cover all sorts of different tracks/roads/bridleways etc..


texts:
I have a question relating to a comment in a closed discussion (Does a value in srtm dem file refer to the center or corner of a cell?).

The comment was following:

"When talking about any DEM, it's resolution is important. SRTM has a resolution of 90 m for non US areas.
Theoretically speaking, you can't really say anything about the height variation within that pixel. So just as PostGIS rasters (or rather All Rasters), it is the value for the entire area of 90m*90m, and not for just some specific point within that square.
However, many interpolating algorithms, assume the height is for the center of the pixel, and will interpolate the heights for points in between the pixel centers."

...................

Does this mean that it is not possible to estimate correctly the height/slope of the landscape from 90m DEM? E.g. if I have calculated slope based on 90m DEM. What is meant by interpolating algorithms, should one interpolate 90m SRTM to get correct height?


texts: Tabulate Intersect with Overlapping Polygons

texts:
I am trying to create a table that summarizes soil information for the given buffer areas of 125m and 500m radii. Both the soil data and buffer data are in vector form. Within each buffer, there can be several different soil polygons varying in total area covered. My survey areas are circular buffers saved as a polygon shapefiles. The problem with this is that at the 500m buffer spatial extent my buffers have a high amount of overlap. Likewise, my 125-meter spatial extents have no overlap. When I run the Tabulate Intersect (Spatial Analyst) tool, the values for overlapping polygons totals to more, in some cases substantially more, than 100%.
My goal is to tabulate the intersect of the 500m buffer extent with the soil layer to get the percent coverage of each soil type within. I'm looking for a more simple way to do this without creating several shapefiles with the same buffer, but only ones that don't overlap.


texts: QGIS 3.4.11: intersecting vector data/list all land cover classen

texts:
I have a question regarding QGIS 3.4.11.
For research purposes I need to intersect two vector data - one polygonized raster (with land cover classes)/mulitiple polygons shp and one point shp (buffers are set at certain distances).
I want to know which land cover classes intersect with the buffers. But when i open the attribut table only one class is listed even if more classes intersect.
Fig 1 shows the intersection with the fid 1117.
Fig 1:

Fig 2 shows the atttibut table. In column DN the land cover classes are noted. As you can see in Fig 1 multiple land cover classes intersect with the buffer yet just one class is listed.
Fig 2:

How can I solve this Problem?


texts:
Is there a JavaScript parser to obtain GeoJSON from GML?

I have found only Python and OGR solutions.


texts:
--The problem--
I've successfully added mbtiles from Mapbox free online source to OpenLayers, but cannot figure out how to get the mbtiles I've created/stored locally to show up in OpenLayers.

--The process--
Style data in Tilemill, export mbtiles, add mbtiles to Tilestream running on localhost, then add to openlayers using the localhost:8888 url.

--The code--

This works:

// Add mapbox mbtiles from web
map.addLayer(new OpenLayers.Layer.TMS('geography-class',
'http://a.tiles.mapbox.com/mapbox/', {
maxResolution: 156543.0339,
type: 'png',
layername: 'geography-class',
isBaseLayer: true
}));

This doesn't work:

// Add tilestream mbtiles local data layer
map.addLayer(new OpenLayers.Layer.TMS('overlayWorld',
'http://localhost:8888/2.0.0/overlay_62ee43/{z}/{x}/{y}.png', {
maxResolution: 156543.0339,
type: 'png',
layername: 'overlayWorld',
isBaseLayer: false
}));

That localhost url comes from the 'TILE URL' field in Tilestream. I've tried multiple different edits to the localhost url but cannot get it right.

Any help is appreciated.


texts: Cannot visualize a WMS in Geoserver from a JPG file

texts:
I'm using a JPG file (8bit unsigned), together with a .jgw and the .jpg.aux.xml having the projection. The image displays well in ArcGIS and QGIS.

Now I want to create the image in Geoserver 2.11.0 as a 'WorldImage' store, and then create a Web Map Service around it. But the image doesn't work in 'Layer preview'. You see nothing.

When I try to consume the WMS in ArcGIS, this is the error I get.

Any ideas?


texts: Reducing the number of mosaic elements using gdal_merge

texts:
I have 4 thousand mosaic elements. So i would like to reduce their count. (for example to 1 thousand, merge every 4 elements)
How can i do this? gdal_merge give me one big raster


texts: I am interested to know why ‘Flow Direction Raster’, derived from ‘Filled DEM’, generates values ranges from 1-128, with pixel depth of 8 bit unsigned integer which is readily capable to store values 0-255, As such, in the Tool’s Help, it is clearly stated that “the output of the Flow Direction tool is an integer raster whose values range from 1 to 255”. My ‘filled DEM’ carries 16bit signed integer pixel depth already, which has capacity to store values range from -32768 to 32767!


texts:
In layout view of ArcMap it is easy to add a scale bar.

In ArcScene there is no layout view and no insert a scale bar tab.

How do I add a scale bar in ArcScene?


texts: How to use the field calculator's centroid() Geometry function in qgis 2.0?

texts: There is an option in the field calculator for centroids under geometry which when you double click gives you "centroid( ".  What information is required here in order to use this function?


texts: Is it possible to query/download satellite images with Python?

texts:
I was wondering if there exists a Python library/API for downloading satellite images. Ideally, I would want to enter the latitude and longitude of a bounding box and get a high resolution image covering that box. Whether it is from Google Maps or other services doesn't matter.

Any ideas?


texts: App using ArcGis Engine throws error at customers

texts:
I developed a program using the ArcGIS Engine 10.0. On my development computer it runs nicely, but when I deploy the binary to the customer computer the program throws an exception when trying to initialize the license.

The exception: "ArcGIS version not specified. You must call RuntimeManager.Bind before creating any ArcGIS components."

This error does show after I called RuntimeManager.Bind, and I'm quite sure I did not accidentially instanciate any other ArcGIS objects before calling RuntimeManager.Bind.

How can I find the source of this problem? Why does it only show up at the customer?

The code I'm using to initialise:

ESRI::ArcGIS::RuntimeManager::Bind(ESRI::ArcGIS::ProductCode::Engine);
ESRI::ArcGIS::esriSystem::AoInitialize^  aoi = gcnew ESRI::ArcGIS::esriSystem::AoInitializeClass();

ESRI::ArcGIS::esriSystem::esriLicenseProductCode productCode = ESRI::ArcGIS::esriSystem::esriLicenseProductCode::esriLicenseProductCodeEngine;
if (aoi->IsProductCodeAvailable(productCode) == ESRI::ArcGIS::esriSystem::esriLicenseStatus::esriLicenseAvailable)
{
  aoi->Initialize(productCode);
}


texts:
I have created two different plugins in QGIS using python. I would like to combine both in one plugin.

How can I do that?


texts: arcpy cannot extract a subdataset from a .3D file

texts:
I am working with data which represents the atmosphere as rasters at 29 different altitudes. This data is stored in a ".3D" file. When I drag and drop into ArcMap it asks me to select the layer I want. When I use arcpy's "extractsubdataset" tool it says "ERROR 000840 The value is not a raster layer". I'm pretty sure it is because the tool only accepts ".HDF" and ".NITF". Is there any other way to automatically extract the data?

I need to extract over 100 subdatasets from over 100 ".3D" files so I can't do it manually. If I do multiple subdatasets at the same time from drag and drop they all get the same name and it is hard to distinguish them apart.

The data is the National Centers for Environmental Prediction North American Regional Reanalysis

Not sure if I am allowed to have links, but here is a link to the data: https://rda.ucar.edu/datasets/ds608.0/index.html#cgi-bin/transform?xml=/datasets/ds608.0/metadata/wfmd/3HRLY%252017%25NARR3D_201708_1012.tar.GrML

For syntax I have the input and output as strings to where the file is located and the subset index as an int.
arcpy.ExtractSubDataset_management(filepath,outputfolder+"test.tif",1)

Does GDAL work? If it does I can get IT to install it on my workstation (I've never used GDAL).


texts: Connecting QGIS 2.14.3 to SQL 2012 both on LocalHost or Remote

texts:
I am working on establishing MSSQL (2012) spatial connection from my QGIS 2.14.3, its giving me an error like:

  (Error opening connection: [Microsoft][SQL Server Native Client 11.0]Named Pipes Provider: Could not open a connection to SQL Server [53]. [Microsoft][ODBC Driver Manager] Driver's SQLSetConnectAttr failed [Microsoft][SQL Server Native Client 11.0]A network-related or instance-specific error has occurred while establishing a connection to SQL Server. Server is not found or not accessible. Check if instance name is correct and if SQL Server is configured to allow remote connections. For more information see SQL Server Books Online. [Microsoft][SQL Server Native Client 11.0]Login timeout expired QODBC3: Unable to connect).

Please anyone with solution.


texts: Installing PostGIS gives error: could not find libpq?

texts:
After giving the command ./configure, it is showing that:

  pQserverversion in -lpq not found. configure: error: could not find libpq.

What does it mean?

I am installing PostGIS 1.5 using Red hat Linux (32 bit) version 5.3. PostgreSQL version 8.3.


texts: libiconv-2.dll missing when running shp2pgsql

texts:
My system is Windows 7 with PostgreSQL 9.4 and PostGIS 2.2.

I'm running shp2pgsql program from command line (cmd window, not PSQL interface).

shp2pgsql "D:\file.shp" schema.table > file.sql

I'm getting the following error message (sorry this is in french but it says System error: the programme cannot start for libiconv-2.dll file is missing on the computer. Try to reinstall the program in order to fix the bug):

I have then completely reinstalled PostgreSQL 9.4 and PostGIS 2.2. The libiconv-2.dll is properly installed inside the bin directory (C:\Program Files\PostgreSQL\9.4\bin). I don't understand why this file is not found on the computer and I don't know what else I could do...


texts: I'm a journalist based in Paris, and for a project I'm working on, I would like to display how public debt rose dramatically across Europe during the crisis. My plan is to show each country's ratio of debt/GDP between the years 2000 and 2010 inclusive. Any country with a ratio under 50% will be shown in green, anything from 50-60% in lighter green or yellow, etc. I'm just not sure of the most effective way to portray this, using Cartodb. I'm really new to the software and to geographic visuals in general. Thanks for your help!


texts: Copernicus ERA5 crs/projecition for NetCDF

texts:
i'm trying to work from netCDF from this link:

https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels?tab=form

but i can't find the correct EPSG.

In documentation page i found that :

  If you request the data in NetCDF format, it is automatically converted and interpolated from the above grids to a regular lat/lon grid.

So i try to open it as WGS84 Lat/Lon but when i try to view under another layer i obtain something like this


texts: Convert Planet Imagery Surface Reflectance product DN to surface reflectance values

texts:
I have a PlanetScope Ortho Scene Product (3B) which includes an AnalyticMS.tif and AnalyticMS_SR.tif files. Based on what I have read so far, the AnalyticMS_SR.tif file is scaled surface reflectance in Digital Number. How do I convert the DN to surface reflectance value?

I have read that for AnalyticMS.tif, I can convert the DNs into TOA reflectance or radiance by multiplying them with the given reflectance coefficient or radiometric scale factor. Is there something like that for the the Surface Reflectance product?


texts: How to determine the location that has the most intersections with other polygons in a certain distance/buffer?

texts:
I have a very scattered map with a lot of very differently shaped polygons. Now I want to determine the location where a circle polygon of a certain size intersects the most surface of these scattered polygons.

It would be even better to be able to determine not only one best location but maybe even some more locations that are the most intersecting to get something like a heatmap.

Here is an exemplary picture:


texts: Move moved long and lat to their own position

texts:
I have a dataset with point that have fields containing latitude and longitude.
I have displaced some points (up to 2 mts) and saved the modification. The field containing the original lat long in the data is not changed but the geometry of the point is changed.
How can i move the point back to their location defined in the fields long and lat that haven't been changed ?


texts: Topographic Position Index unexpected result in QGIS

texts:
I am using Windows 7 enterprise and QGIS version 2.12.2 Lyon.  Have found an unexpected result when I use QGIS to generate a topographic position index (TPI) raster from a 30m by 30m DEM.  I find a pattern of what I believe to be aberrantly valued pixels forming a square grid pattern that are approximately 900 m on each side.  I expected the resulting raster to not have this grided appearance.

UPDATE: I have looked at my hillshade, slope, aspect, terrain ruggedness index and roughness rasters generated using the same tool and they all have this same "grid" (although it is much less apparent than in the TPI raster).

Here's a screenshot of the original DEM I used:

Here's a screenshot of the resulting TPI.  Notice the grid pattern formed over what I suspect to be an accurate TPI raster layer.

Here's how I created the TPI Raster: Raster – analysis – DEM:


texts: Avoiding ArcGIS Online limitation on number of records/features per shapefile?

texts:
I know that to publish a shapefile to ArcGIS Online you must zip the shapefile and it must contain a maximum of 250 records.

I have a shapefile with 32,000 records.

How do I divide this into 250 records per shapefile?

Original Question (prior to Google translation and touch up above)

So che per pubblicare uno shape su arcgis online si deve zippare lo shape ed esso deve contenere 250 record massimo...
Io ho uno shape con 32000 record...come faccio a suddividere tutti questi record a 250 per volta?


texts: Any way to get smaller chunk of raster grid from big `netCDF` file in R?

texts:
I have climate data in big netCDF file which covers 60 years time span climate data observation (please take a look into this link). However, I want to subset this climate data for only recent 20 years while discarding the time period that I don't want to use. I used raster package to read netCDF file in R, but to get a smaller chunk of raster grid from original netCDF file is unknown to me.

goal:

I want to subset or get a smaller chunk of raster grid (only 20 years' period climate data observation) from big netCDF file (it covers 60 year's climate data).

Is there any efficient fast way to get this done in R?


texts:
I am looking for latitude and longitude coordinates of french subdivision.

I am not used to GIS but I saw that ISO 3166-2 contains default code of such subdivisions, their name etc but where can I download the complete data in a readable format (with python preferred) with name, code and coordinates ?


texts: Portal for ArcGIS WMTS layers are optimised & un-cacheable

texts:
When a WMTS in Portal for ArcGIS 10.3.1 is opened (via Portal) in an mxd, the WMTS layer properties are not standard. The; Source, Layers, and Advanced tabs are missing. When the portal WMTS layer used as an input into arcpy.ExportMapServerCache_server() it fails without any errors.

I've discovered that if you copy and paste the layer into the same data frame, the layer regains all of it's normal properties. I think the portal extracted layer is actually put in a hidden optimized base map layer.

I need to use ArcPy to cache this service locally as a part of a greater user workflow. However beofre I proceed I need a method to handle this "optimized portal WMTS layer".. Has anyone experienced this and then figured out a way around it?


texts: How to calculate average slope and compass point of polygons in QGIS (3.18)?

texts:
I am using QGIS to rank about 800 PV areas according to their potential. I have a vector (.shp) layer with certain areas in polygons (inc. attribute table, every polygon has it's own ID).
Now I would like to find out the avg. slope and orientation (compass or in degree 0-360) of every polygon. I already have a raster (.TIFF) with the slope (resolution 25m) over the full area.
Is it possible to calculate slope and orientation somehow automatic and import them as columns in the attribute table?
If there is no possibiltiy to calculate with a function in the attribute table, what is the best way to calculate manual for every single polygon?


texts: What expressions for QGIS field-calculator?

texts:
I do not understand what language I have to use in field-calculator and labels on QGIS:

SQL or Python?

And where I can find the possible script/expressions to use?

in the image you can see the options for expressions in QGIS 1.8 italian version

Thanks


texts: Eliminating gaps between polygons using ArcMap?

texts:
I work in ArcMap 10.4 with an Advanced level license. I'm searching for a fast and effective method for eliminating gaps between polygons, so that one polygons follows boundary of the other. In other words - changing shape according to the other polygon (as seen in the example picture - purple polygon's boundary should follow more detailed blue polygon's boundary).

Is there a tool to do it?


texts: UpdateCursor not updating when run from PyCharm?

texts:
I have a script in ArcPy that iterates through a set of intersection points and determines if all intersecting streets have the same attribute value in the CROSSLTS field, and if they do, sets the the intersection point's MarkForProcessing field value to 1. if the intersecting segments do not have identical CROSSLTS values, it sets the MarkForProcessing field value to 0.

When I run this script from the Python window of ArcMap, it takes 2 hours and I get an "out of memory" message before it completes. when i run it from PyCharm, it takes roughly 10 minutes, but the MarkForProcessing field is never updated.

Can someone could take a look?

import arcpy
workspace = 'C:\SRC\Bikestress\LTS classification.gdb'
arcpy.env.overwriteOutput = True

##check if all elements in a list are identical, returns true if so.
def checkEqual(lst):
    return lst.count(lst[0]) == len(lst)

arcpy.MakeFeatureLayer_management(r'C:\SRC\Bikestress\LTS classification.gdb\Intersections','INTERSECTIONSLYR','',workspace)
arcpy.MakeFeatureLayer_management(r'C:\SRC\Bikestress\LTS classification.gdb\roadmidsplit2','STREETSLYR','',workspace)

Intersectionfields = ["CONTROLTYPE","GLOBALID","MarkForProcessing","SHAPE@"]
StreetsFields = ["CROSS_LTS"]

count = 0
with arcpy.da.UpdateCursor('INTERSECTIONSLYR', Intersectionfields) as IntersectionCursor:
    for Intersection in IntersectionCursor:
        count += 1
        if Intersection[0] not in ["4STOP", "PED", "TRAFFICSIG"]:

            # select the intersection of the current row, and all adjacent segments
            arcpy.SelectLayerByLocation_management('STREETSLYR', "INTERSECT", Intersection[3], "", "NEW_SELECTION")

            CROSSLIST = []

            #make sure selection is not empty
            desc = arcpy.Describe('STREETSLYR')

            if desc.FIDSet == '':

                with arcpy.da.SearchCursor('STREETSLYR', StreetsFields) as StreetCursor:
                    for Street in StreetCursor:
                        CROSSLIST.append(Street[0])

                if len(CROSSLIST) != 0:
                    if not checkEqual(CROSSLIST):
                        Intersection[2] = 1
                    else:
                        Intersection[2] = 0

            IntersectionCursor.updateRow(Intersection)
        if str(count)[-2:] == "00":
            print(str(count) + " records processed")
print('writing intersections layer...')
arcpy.CopyFeatures_management("INTERSECTIONSLYR", r'C:\SRC\Bikestress\LTS classification.gdb\IntersectionsUpdated')
print("done")


texts: How to convert a shapefile (shp) to raster via GEE

texts:
I would like to convert a shp to raster via GEE. My SHP has several polygons, which have a value between 1 and 2 in the attribute table.
I would like the output raster to assume values 1 and 2 according to the value of the polygons.
I know there is a "toReduce" command, but I didn't understand how to do that, even with the tutorial provided by Google.
Context
I would like to add a mask for an image rating. In this shp, the value 1 would be ignored by the classifier.


texts: How to install PostgreSQL 9.2 with PostGIS 2.0 on Ubuntu 11.10 (or higher)?

texts:
I have read this question: https://dba.stackexchange.com/questions/25033/install-postgis-v2-for-postgresql-v9-2, but it didn't resolve my problem.

I added the PostgreSQL backports PPA and installed the postgresql-9.2 package, but there is no package for PostGIS.


texts: Sort and create attribute field to the shapefile

texts:
I am creating a shapefile layer and adding an 'Area' attribute field using the following script:

 # Create shapefile layer
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource(shp_output_path)
srs = osr.SpatialReference(wkt=prj)

dst_layer = dst_ds.CreateLayer(dst_layername, srs=srs)
area_field = ogr.FieldDefn('AREA_ha', ogr.OFTReal)
dst_layer.CreateField(area_field)

for feature in dst_layer:
    geom = feature.GetGeometryRef()
    area = geom.GetArea()
    feature.SetField("AREA", area)
    dst_layer.SetFeature(feature)

dst_ds = None

How can I select the 10 largest area features and delete the rest?


texts: What is the most natural way to save a graphic view of a map generated by an Overpass Turbo query?

texts:
As a preface, I'm not very familiar with Overpass Turbo, so I may be giving too little or too much info. Please ask for clarification as needed.

Brief version of the question:

Someone built a query that draws on a map, and I want to save that drawing for future reference. What would be the most convenient way to do that? I don't want to have to re-run the query as it is pretty slow. Can I get a hyperlink to this data, or should I export it in some form, that could either be imported back to the website, or to another graphical tool? (A screenshot is not detailed enough for this situation.)

Here's more background if you're curious:

The game Pokemon Go, developed by Niantic, heavily uses data from OSM to determine where certain things will happen in the game. One example is using areas with certain tags, such as parks, and making them a "nest" in the game, where a specific Pokemon will appear frequently.

Certain individuals have managed to determine the exact timestamps of the data that Niantic uses for certain things, as well as the tags that are used. One such individual built a couple queries and posted them on Reddit here. It is a map drawn by this query, over my hometown, that I want to refer to.


texts:
I am having an issue, which has essentially seemed to have split into two different issues. I was initially trying to work with TIMED satellite data, starting off with GUVI 1C data. Rummaging around http://www.timed.jhuapl.edu and other  sites, they made it clear the the data format was netCDF - which is unfamiliar to me. Regardless, it was downloaded, and I got the following (for Win 10 OS, 64-bit):

[![enter image description here][1]][1]

So initially speaking, where is the actual application to use netCDF?

Unsure, I downloaded Panolpy to help out.

Then I proceeded to download some TIMED GUVI L1 files:

And now I can't seem to open the L1C file - although according to the site it was supposed to be a netCDF file.


texts: Quick way to isolate layer and hide others in QGIS?

texts:
Similar to how CAD works, is there a way to quickly hide all other layers than the one you are working on?

Working on QGIS 2.18 on mac


texts: Determining percentages of shape layer in each grid cell using ArcGIS Pro

texts:
We have five different shape layers (water, forest, ...) for a certain area (which is supposed to be a square about ~50 sq km). What we need is to put something like a grid over our shape layers (each cell = 1 sqm). Then, we would like to know the percentage of each shape layer for every grid cell.
We are using ArcGIS Pro.
It would be also great to put the results into a PMG file or similar so that we have a grid as a layer where each cell has a color that represents the percentage of a layer (for example: a grid cell with five different colors for the cells, presenting the percentage of water).


texts:
I am trying to add the WFS layer to my map to make it clickable via L.control, unfortunately, I am not using this properly as when I click on the other layers in control, the popup feature is still there on the map while WMS layer disappears.

I guess the solution is putting WFS as a layer in control, but I am not sure how.

My Code:

L.control.layers({
        "osm": osm.addTo(map),
    {
         "points_Mohsen": // WMS Layer
           L.tileLayer.wms("http://localhost:9090/geoserver/localhost/wms", {
            layers: 'localhost:DP_View',
            format: 'image/png',
            transparent: true,
            version: '1.1.0'
        })
        ,
        "polygon_Mohsen": L.tileLayer.wms("http://localhost:9090/geoserver/localhost/wms", {
            layers: 'localhost:TEST_JSON',
            format: 'image/png',
            transparent: true,
            version: '1.1.0'
        })
        , 'drawlayer': drawnItems
    }
     ,   { position: 'topright', collapsed: false }).addTo(map);

    map.addControl(new L.Control.Draw({
        edit: {
            featureGroup: drawnItems,
            poly: {
                allowIntersection: false
            }
        },
        draw: {
            polygon: {
                allowIntersection: false,
                showArea: true
            }
        }
    }));

  map.on('overlayadd', onOverlayAdd);

    function onOverlayAdd(e) {

        //  alert(e.layer.wmsParams.layers);

        if (e.name == "points_Mohsen") {
            if (map.hasLayer(WFSLayer))
                map.removeLayer(WFSLayer);

            var defaultParameters = {
                service: 'WFS',
                version: '2.0.0',
                request: 'GetFeature',
                typeName: 'localhost:DP_View',
                outputFormat: 'text/javascript',
                format_options: 'callback:getJson'

            };

            var parameters = L.Util.extend(defaultParameters);
            var URL = owsrootUrl + L.Util.getParamString(parameters);
            ////layerVisibility[L.Util.extend(defaultParameters)] = true;

            ajax = $.ajax({
                url: URL,
                dataType: 'jsonp',
                jsonpCallback: 'getJson',
                success: function (response) {
                    WFSLayer = L.geoJson(response, {
                        style: function (feature) {
                            return {
                                stroke: false,
                                fillColor: 'FFFFFF',
                                fillOpacity: 0
                            };
                        },
                        onEachFeature: function (feature, layer) {
                            popupOptions = { maxWidth: 200 };
                            layer.bindPopup("Popup text, access attributes with ParcelID:" + feature.properties.P_id + "</br>"
                                , popupOptions);
                        }
                    }).addTo(map);
                }
            });
        }
 };


texts: SQL query builder in QGIS 3 doesn't display columns names

texts:
I'm creating sql queries using QGIS 3.8 DB Manager.

The problem i'm facing is that when I write a query, I don't have any option to display the columns names :

Is there any way to display the columns names like in QGIS 2?


texts: How to enable snapping when moving multiple features at once?

texts: Is there a way to enable snapping when moving two or more features at once?


texts: Connecting QGIS 3.4 to a MS Access database

texts:
I try to connect to MS Access and i have already created an ODBC connection.
Now, I have the the problem that I cannot open the database in QGIS (but the connection was succesful).


texts: Dynamic Tiling Isn't working(create tiles on demand isn't working)

texts:
I have a map where users are going to primarily being viewing a particular area, so I cached this area and checked off the dynamic tiling option inside of the service properties this worked fine for a while every time the map would load and a user went outside my tiled area the system behaved as expected it started the tiling process tiled that area they went to and stopped. Then the dynamic tiling quit operating, it won't start at all no matter where I go on my map I of course have checked the web and this site but have found no answers every time  my end user goes goes to the area outside my tiled area they get a massive blank spot all the way around the tiled area. Any suggestions as what to try will be greatly appreciated.

EDIT: I'm using arcgis server and my client is a webpage that we built to use the map, it uses the javascript api.


texts: Significance test for Moran's I using Monte Carlo Simulation?

texts:
I calculated the residuals for my model and measured Moran's I for the residuals of the model.

How can I pragmatically determine the significance of Moran's I using Monte Carlo simulation?


texts: Where can I find free shapefiles of oil/natural gas pipelines in Mexico?

texts: I'm looking for free shapefiles of oil or natural gas pipelines pertaining specifically to Mexico or even trans-national pipelines between the U.S and Mexico.


texts: How to make area from specific node or relation, using it's element id number in Overpass QL?

texts:
I can make an area in Overpass QL, using name and relation attributes (tags) filter of element, something like this:
area["name:en"="Moscow Oblast"]["boundary"="administrative"]["admin_level"="4"]->.region;
(node(area.region)[place="city"];node(area.region)[place="town"];);
out;

But searching via name is not so stable, and will broke if element name is changed. So much better (for performance too) is to make an area from element id. Here is example of searching relation via it's element id:
rel(51490);
out;

And can't understand, how to pass element type and id to area query? Something like this:
area(rel(51490))->.region;
(node(area.region)[place="city"];node(area.region)[place="town"];);
out;

or other variants like area["id"=51490]->.region; or area["type"="rel"]["rel"=51490])->.region; don't works.
Can anyone give me an example, how to make area from OSM element by it's id?


texts: Counting number of plants in area using images from UAV?

texts:
How do I count the number of pineapple plants in a 1-hectare block?

I have heard of these algorithms but I don't know how to do it.

I'm still a newbie of QGIS, ArcGIS and ENVI. I have the NIR, RGB, Red Edge images of the said area. Can anyone teach me how to do it?

I have access to eCognition developer but I'm still new to this software.


texts: Getting from IFeatureClass on annotation feature class to IAnnoClass using ArcObjects

texts:
I am working in Python 2.6 and ArcGIS 10.0 through comtypes, so I have full ArcObjects access.  For example, see the answer I tacked onto Accessing ArcObjects from Python?.  That makes this strictly an object model traversal question.

I have an IFeatureClass on an annotation feature class.  I need to get to IAnnoClass to get the ReferenceScale property.

How do I do that?

I am guessing that it starts with the Extension property inherited from IClass, but I may be way off base.


texts: if-then-else condition in Google Earth Engine

texts:
I have Sample RGB combination from 3 different Bands

// Input RGB imagery
var image = ee.Image("users/k2yulianto/Exercise/Ras/Imagery");
// Use Tree bands from RGB Imagery
var bands = ['b1', 'b2', 'b3'];
var b1 = ['b1'];
var b2 = ['b2'];
var b3 = ['b3'];
// Display the inputs imagery
Map.addLayer (image, {bands: ['b1', 'b2', 'b3'], max: 1, gamma:1, opacity:0.80} , 'RGB imagery',true);
Map.setCenter (101.4476, 0.4925, 12);

And I have also sampling point distribution on top of RGB Imagery.

// Sampling point on top of RGB imagery to identify 4 diferent class
var Sampling = ee.FeatureCollection("users/k2yulianto/Exercise/Sampling/Sampling");
Map.addLayer (Sampling, {opacity:0.80},'Sampling point', true);

Based on the sampling point, I am to produce decision tree using WEKA to distinguish the type of class on RGB imagery.

//var exercise = ee.Algorithms.if(condition,truecase,falsecase)
var exercise = ee.Algorithms.If(b2 <= (0),(ee.Algorithms.If(b3 <= (0),"Sample 1","Sample 2")),(ee.Algorithms.If(b1 <= (0),"Sample 3","Sample 4")));
print('Imagery output: ',exercise);

I want to display if-then-else condition into an image classification.

// Display the 4 Class exercise output imagery
var RGBPalette =
['#5f932c', // Class 01
'#bdff09', // Class 02
'#fdf613', // Class 03
'#abfffb', // Class 04
];
Map.addLayer(exercise,{min: 0, max: 4, palette: RGBPalette}, 'Exercise imagery');

After I run the scripts, I got an error warning:

  Cannot add an object of type  to the map.

Here is my full code:

https://code.earthengine.google.com/77b2238d1a628c59dfa440a27adfa378

How can I solve this?