texts: Pgrouting: problem with pgr_withpointsDD

texts:
I'd like to calculate a propagation from sources using my road graph. In extenso, I'd like to answer the following question: where are the limits of my isochrones for every sources I got. My sources (point) are not on the graph but outside.
No problem to create my graph.
Here's how I create a table with the sources and elements attaching them to my graph:
create table exo.psource as
(with poi as (
    select p.gid as gid_pt,
        p."q (m3/h)" as debit,
        p.volume as volume,
        p.geom
    from exo.epi  p
    )
SELECT
    poi.gid_pt::int4,
    poi.volume::int8,
    poi.debit::int4,
    closest_tr.gid::int4,
    closest_tr.dist::double precision,
    st_linelocatepoint(closest_tr.geom, poi.geom)::double precision as fraction
FROM poi
LEFT JOIN LATERAL
    (SELECT
        t.gid,
        ST_Distance(t.geom, poi.geom) as dist,
        t.geom
    FROM exo.troncon t
    ORDER BY poi.geom <-> t.geom
    LIMIT 1
    ) AS closest_tr
ON TRUE
)

After that I run the following query to get the vertices at 200m or less from my sources:
select
    pgr.*,
t   vp.the_geom
from
pgr_withPointsDD(
    'select gid AS id, nsource::int4 AS source, ncible::int4 AS target, st_length(geom)::float8 AS cost, st_length(geom)::float8 as reverse_cost
    from exo.troncon'::text, 'select gid_pt as pid, gid as edge_id, fraction from exo.source'::text,
    (
        select
            array_agg(exo.source.gid_pt order by gid_pt asc) from exo.source
    ) , -- un array de tous les id de mes sources
    200
,
    details := true
    ) pgr
left join exo.troncon_vertices_pgr tvp on tvp.id = pgr.node

Everything seems correct in my (favorite) GIS:

Except: for some sources, the closest edge hasn't been browsed (see source 2 to vertex 142 by example)
What did I miss?
No vertex is exactly at 200m from a source, how do I calculate the remaining part of edges to get 200m exactly?
My sources are distants of the edge, how do I integrate this distance into the calculation?
I'm quite sure the solution to this probleme may benefit a lot of people.
I orignally post this problem on a french forum.
Here are the datas
I've read those articles and used them partially or fully to solve (not totally) my problem:
présentation foss4g
une question sur gis stackeexchange
un post de bostongis
article sur geoafrica


texts: Can Mapnik render spatial data from SQL Server?

texts:
This is a very simple question - I just want to know if anybody can confirm that they have successfully used Mapnik to render tiles directly from data held in SQL Server. (And, if so, what did you do to make it work?) Please do not reply unless you have personally confirmed that this is possible... not simply that it should be possible

I've seen people recommend Mapnik for rendering data from SQL Server in threads such as this one: Any Good Map Rendering Engines for SqlGeometry/SqlGeography

Also, the mapnik website states that "all OGR formats are supported": http://mapnik.org/faq/#data

However, I have been completely unable to use a connection to my local SQL Server database as a datasource for a Mapnik layer. As a test, I'm using a straightforward test OGR VRT, which executes a SQL statement on a SQL Server to select a single inline geography Polygon as Well-Known Text, as follows:

<OGRVRTDataSource>
  <OGRVRTLayer name="test">
    <SrcDataSource>MSSQL:server=zangief\DENALICTP3;database=TempDB;</SrcDataSource>
    <SrcSQL>SELECT geography::STPolyFromText('POLYGON((0 50, 2 50, 2 53, 0 53, 0 50))', 4326).STAsText() AS geomWKT</SrcSQL>
    <GeometryField encoding="WKT" field="geomWKT"/>
    <LayerSRS>EPSG:4326</LayerSRS>
    <GeometryType>wkbPolygon</GeometryType>
  </OGRVRTLayer>
</OGRVRTDataSource>

(Note I've also tried this with more conventional SELECT * FROM Table syntax). OGRINFO succesfully opens this VRT and reports the following:

ogrinfo mssql.ovf -so -al
INFO: Open of `mssql.ovf'
      using driver `VRT' successful.

Layer name: test
Geometry: Polygon
Feature Count: 1
Extent: (0.000000, 50.000000) - (2.000000, 53.000000)
Layer SRS WKT:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]
geomWKT: String (0.0)

However, setting this same VRT file as a datasource for a Mapnik layer reports the following error:

ERROR 1: Failed to open datasource `MSSQL:server=zangief\DENALICTP3;database=Tem
pDB;'.
Traceback (most recent call last):
  File "rundemo.py", line 48, in <module>
    provlines_lyr.datasource = mapnik.Ogr(file="mssql.ovf", layer="test")
  File "C:\OSGeo4W\apps\Python25\lib\site-packages\mapnik\__init__.py", line 418
, in Ogr
    return CreateDatasource(keywords)
RuntimeError: Failed to open datasource `MSSQL:server=zangief\DENALICTP3;databas
e=TempDB;'.

It seems that there's some problem between the way that the OGR driver creates the SQL Server connection and Mapnik, but I don't know exactly where it is. I've not yet found any evidence or examples to suggest that anybody else has successfully been able to achieve this either...


texts: Using server data store with Python toolbox?

texts:
ArcGIS Server lets you set server side data stores:

This is nice. I can use one connection on my local machine and another on my server. (Right? Please tell me I'm right.) But for the life of me, I can't find a way to actually set the data store on my geoprocessing service. My service is defined using a Python toolbox. It looks something like this:

import mygptool

class MyGPTool(object):
    def __init__(self):
        """Define the tool (tool name is the name of the class)."""
        self.label = "MyGPTool"
        self.description = "Does my stuff!"

    def getParameterInfo(self):
        """Define parameter definitions"""
        return [arcpy.Parameter(displayName="State (Abbreviation)"
                               ,name="state_abbreviation"
                               ,datatype="GPString"
                               ,parameterType="Required"
                               ,direction="Input"
                               )]

    # Boilerplate isLicensed, updateParameters, updateMessages methods

    def execute(self, parameters, messages):
        """The source code of the tool."""
        state_abbreviation = parameters[0].valueAsText
        mygptool.run_my_tool(state_abbreviation)
        return

I have a feeling I need to add something to this definition, but I'm having an awful lot of trouble figuring out what I need to add to give my service a connection to work with. (At the moment, I'm generating a temp connection file based on connection parameters in a config file, but that is something I'm looking to change as it doesn't make moving across environments very easy.)

A workspace parameter doesn't seem to help:

arcpy.Parameter(displayName="My Database"
               ,name="smart_db"
               ,datatype="DEWorkspace"
               ,parameterType="Required"
               ,direction="Input"
               )

This results in a constant parameter on my service, defeating the purpose of being able to specify the connection.

Also, I really do need a connection, not a table or feature class or what have you. The service creates multiple query layers using this connection; we're trying to avoid an MXD as this makes deployment much more difficult and tedious.

How can I make use of server configured data stores with a Python toolbox?


texts: How to define an output parameter in Python script exported from ModelBuilder

texts:
I have a modelbuilder parent model with two nested models in it that ouputs a dwg file.  I exported the parent model to a python script which is shown below.  In modelbuilder, the output dwg file is a parameter so I can define its name and location when I run it.  In the python script, it does not let me set the name and location.  It shows the message "this model has no parameters" when its run in Python.  I would like to be able to define the output name and location but am not sure how to do this in Python.  The output file is "OutputFile.dwg" in the code below.

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# GIStoCAD.py
# Created on: 2014-02-27 13:15:51.00000
#   (generated by ArcGIS/ModelBuilder)
# Usage: GIStoCAD <Output_File_Name>
# Description:
# ---------------------------------------------------------------------------

# Import arcpy module
import arcpy

# Load required toolboxes
arcpy.ImportToolbox("/cadd1/gisusers/CADGIS/Version2014/GIStoCADv10_1.tbx")

# Script arguments
Output_File_Name = arcpy.GetParameterAsText(0)
if Output_File_Name == '#' or not Output_File_Name:
Output_File_Name = "E:\\Workspace\OutputFile.dwg" # provide a default value if unspecified

# Local variables:
Delete_succeeded = Output_File_Name
GIStoCAD_gdb__4_ = "C:\\GIStoCAD.gdb"
V_PROP_LINE_PGCO = "C:\\GIStoCAD.gdb\\V_PROP_LINE_PGCO"
V_PROP_LINE_MOCO = "C:\\GIStoCAD.gdb\\V_PROP_LINE_MOCO"
V_BLDG_OUTL_PGCO = "C:\\GIStoCAD.gdb\\V_BLDG_OUTL_PGCO"
V_BLDG_OUTL_MOCO = "C:\\GIStoCAD.gdb\\V_BLDG_OUTL_MOCO"
V_ROAD_ASPH_PGCO = "C:\\GIStoCAD.gdb\\V_ROAD_ASPH_PGCO"
V_ROAD_ASPH_MOCO = "C:\\GIStoCAD.gdb\\V_ROAD_ASPH_MOCO"
V_SSWR_PIPE = "C:\\GIStoCAD.gdb\\V_SSWR_PIPE"
V_WATR_PIPE = "C:\\GIStoCAD.gdb\\V_WATR_PIPE"
V_WATR_STRC = "C:\\GIStoCAD.gdb\\V_WATR_STRC"
V_SSWR_MHOL__3_ = "C:\\GIStoCAD.gdb\\V_SSWR_MHOL"
V_WATR_INST = "C:\\GIStoCAD.gdb\\V_WATR_INST"

# Process: Setup
arcpy.gp.toolbox = "//cadd1/gisusers/CADGIS/Version2014/GIStoCADv10_1.tbx";
# Warning: the toolbox //cadd1/gisusers/CADGIS/Version2014/GIStoCADv10_1.tbx DOES NOT have an alias.
# Please assign this toolbox an alias to avoid tool name collisions
# And replace arcpy.gp.Setup3(...) with arcpy.Setup3_ALIAS(...)
arcpy.gp.Setup3(GIStoCAD_gdb__4_, V_PROP_LINE_PGCO, V_PROP_LINE_MOCO, V_BLDG_OUTL_PGCO, V_BLDG_OUTL_MOCO, V_ROAD_ASPH_PGCO, V_ROAD_ASPH_MOCO, V_SSWR_PIPE, V_WATR_PIPE, V_WATR_STRC, V_SSWR_MHOL__3_, V_WATR_INST)

# Process: Export
arcpy.gp.toolbox = "//cadd1/gisusers/CADGIS/Version2014/GIStoCADv10_1.tbx";
# Warning: the toolbox //cadd1/gisusers/CADGIS/Version2014/GIStoCADv10_1.tbx DOES NOT have an alias.
# Please assign this toolbox an alias to avoid tool name collisions
# And replace arcpy.gp.Export3(...) with arcpy.Export3_ALIAS(...)
arcpy.gp.Export3(Output_File_Name, "C:\\GIStoCAD.gdb\\V_PROP_LINE_PGCO;C:\\GIStoCAD.gdb\\V_PROP_LINE_MOCO;C:\\GIStoCAD.gdb\\V_BLDG_OUTL_PGCO;C:\\GIStoCAD.gdb\\V_BLDG_OUTL_MOCO;C:\\GIStoCAD.gdb\\V_ROAD_ASPH_PGCO;C:\\GIStoCAD.gdb\\V_ROAD_ASPH_MOCO;C:\\GIStoCAD.gdb\\V_SSWR_PIPE;C:\\GIStoCAD.gdb\\V_WATR_PIPE;C:\\GIStoCAD.gdb\\V_WATR_STRC;C:\\GIStoCAD.gdb\\V_SSWR_MHOL;C:\\GIStoCAD.gdb\\V_WATR_INST")

# Process: Delete
arcpy.Delete_management(GIStoCAD_gdb__4_, "")


texts:
It's posible to automaticly generate tiles from QGIS web server ?

I need export a project (europe size) from zoom 9 to 13 in png Format in order to use it with leaflet like.
I tried many extention that are not working properly.


texts: Dynamically assign style to chloropleth map leaflet

texts:
I'm trying to create a chloropleth using Leaflet, inspired by the official Leaflet documentation, but I've been having some trouble in achieving this, I think due to how my data is structured.
I have a javascript object which has been drawn from a PostGIS database and stored in a dictionary object elsewhere in my project. The following code is used to display the chloropleth when a button is toggled on the website:
function getColor(d,range) {
    return  d === null ? '#FFFFFF':
        d > range[6] ? '#800026' :
        d > range[5] ? '#BD0026' :
        d > range[4] ? '#E31A1C' :
        d > range[3] ? '#FC4E2A' :
        d > range[2] ? '#FD8D3C' :
        d > range[1] ? '#FEB24C' :
        d > range[0] ? '#FED976' :
                    '#FFEDA0' ;
    }

function onDataLayerToggled(params) {
    dataLayer = geoDataLayers.find( layer => layer.type == params['dssOverlay'] );
    var range = [500,1000,2000,3000,4000,5000,6000];
    for (i=0;i<dataLayer.layer.options.data.length; i++) {
        dataLayer.layer.setStyle({
            fillColor: getColor(dataLayer.layer.options.data[i]['total'], range),
            weight: 2,
            opacity: 1,
            color: 'white',
            dashArray: '3',
            fillOpacity: 0.7
        });
    }
    dssMap.addLayer( dataLayer.layer );

The structure of the data section of dataLayer (the bit I'd like to assign colours to) is as follows:
Array(22): [{...}, {...}, {...}, ....], where inside each dictionary is the following:
aimag_code: 181
camel: 8.11
cattle: 224.995
goat: 1338.223
horse: 236.749
period: 2020
sheep: 1799.761
total: 3607.838

When I try to use the above code to assign a colour to the "total" field in each dictionary, it assigns the last RGB code to every value in the array. When I try to implement this in the same way as shown in the leaflet documentation (using the code snippet below), the same result happens.
dataLayer.layer.setStyle({
    fillColor: getColor(dataLayer.layer, range),
    ....
}

I've been searching around for a solution all day and haven't managed to find something to fix this yet.


texts: Outputting temporary feature layer from geoprocessing script tool to current map (ArcGIS Pro)

texts:
This is my first time attempting to create a script tool. I have gotten my script to work when running it in the Python window in ArcGIS Pro, however running from the script tool is slightly different.
Explanation of purpose: We have a dataset detailing crimes. One field called "Label" is the general description (theft, assault, arcson, etc.). There is a secondary field in the dataset called "description" that describes the offense in more detail (assault with a deadly weapon, misdemeanor assault, etc.). The goal is to take main dataset, and explode it into separate feature layers, one feature layer for each unique value in the "label" field. So a feature class for assault, a feature class for arson, etc. These feature classes will be symbolized by unique values for each of the more detailed offenses.
This dataset changes, and there are not always the same number of unique values in the monthly report, so it needs to be flexible in how many output feature layers are created. When I run the script with the parameters set correctly in the Python window, the feature layers are automatically added to the map document. However when I run this from the geoprocessing script tool, it runs successfully, but no output is added.
import os
import arcpy

# Set environment
arcpy.env.overwriteOutput = True

proProject = arcpy.mp.ArcGISProject("CURRENT")
m = proProject.listMaps("Map")[0]

# Get parameters
Input_FeatureClass = arcpy.GetParameterAsText(0)
input_Field1 = arcpy.GetParameterAsText(1)
input_Field2 = arcpy.GetParameterAsText(2)

# Building unique values list...
arcpy.AddMessage("Building unique values list...")
valueList = []
with arcpy.da.SearchCursor(Input_FeatureClass, [input_Field1]) as cursor:
  for row in cursor:
    l = list(row)
    valueList.append(l[0])
uniqueValues = list(set(valueList))

# Creating Feature Layers...
arcpy.AddMessage("Creating Feature Layers...")
for value in uniqueValues:
  arcpy.management.MakeFeatureLayer(Input_FeatureClass, str(value),  str(input_Field1) + " = '" + str(value) + "'")


texts:
I would like to move the whole DEM up of fixed offset (z=+0.613).

I am a beginner in QGIS. How can I do this?


texts: Is it possible to change the current map document from a geoprocessing tool

texts:
I'm writing a BaseGeoprocessingTool with the ArcObjects API and I need to open an Mxd then change the current document by this opened document.

I managed to open the Mxd document:

IMapDocument mapDocument = new MapDocument();
mapDocument.open(filename, "");

And inside a geoprocessing tool, it's possible to get the current map or page layout by using gpUtilities.getMap() or gpUtilities.getActiveView().

But is it possible to get the current IApplication in order to set the current document?

Or is there another way to change the current document?

Else, I need to handle the current document, in order to change data frame settings, Maplex settings, print and export settings... not only layers.

I use ArcMap 10.2.


texts: Join the result of an SQL query from PostgreSQL to QGIS

texts:
I want to send an SQL query from a Python script in QGIS to PostgreSQL/PostGIS and get back the result as a table which I then join in QGIS to a shapefile. How do I do that? At the moment I try to achieve the first step: sending an SQL query to PostgreSQL/PostGIS and getting a joinable table back.
It works manually with DB Manager but I don’t know how to translate that into Python. The function seems to be qgis:postgisexecuteandloadsql. But that fails even for a simple SQL:
Parameters = { 'DATABASE' : 'HistoPG', 'GEOMETRY_FIELD' : '', 'ID_FIELD' : 'eidied', 'SQL' : 'Select eidied FROM testunion5' }
processing.run("qgis:postgisexecuteandloadsql",  Parameters)

The error message is:
Traceback (most recent call last):
File "C:\PROGRA~2\QGIS3~1.10\apps\Python37\lib\code.py", line 90, in runcode
exec(code, self.locals)
File "", line 1, in
File "", line 18, in
File "C:/PROGRA~2/QGIS3~1.10/apps/qgis-ltr/./python/plugins\processing\tools\general.py", line 106, in run
return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)
File "C:/PROGRA~2/QGIS3~1.10/apps/qgis-ltr/./python/plugins\processing\core\Processing.py", line 181, in runAlgorithm
raise QgsProcessingException(msg)
_core.QgsProcessingException: There were errors executing the algorithm.
The table testunion5 actually has a geometry field, but I want just to see a table in QGIS, not a layer. By the way with a geometry field it doesn’t work either:
Parameters = { 'DATABASE' : 'HistoPG', 'GEOMETRY_FIELD' : 'geom', 'ID_FIELD' : 'eidied', 'SQL' : 'Select eidied FROM testunion5' }

And the actual SQL text which I eventually want to use and which works in DB Manager doesn’t have a geometry field anyway.
I tried also ‘PostgreSQL execute and load SQL’ from the processing toolbox.

The PostGIS log shows the following:
ERROR:  zero-length delimited identifier at or near """" at character 8
STATEMENT:  SELECT "" FROM (Select eidied FROM testunion5) AS "subQuery_0" LIMIT 0
What can I do to get a table from PostgreSQL/PostGIS which I then can join to a shapefile?
I’m using QGIS 3.10, PostgreSQL 10 and PostGIS 2.5.
The following works in DB Manager:

If I try to repeat that in the console:
Parameters = { 'DATABASE' : 'HistoPG', 'GEOMETRY_FIELD' : '', 'ID_FIELD' : 'nummer', 'SQL' : 'Select * FROM testunion' }
processing.run("qgis:postgisexecuteandloadsql", Parameters)

I get the following error message in View/Panels/Log Messages/PostGIS:
WARNING    Invalid PostgreSQL layer
View/Panels/Log Messages/Processing has the following messages:
CRITICAL    This layer is invalid!
Please check the PostGIS log for error messages.
CRITICAL    This layer is invalid!
Please check the PostGIS log for error messages.
CRITICAL    There were errors executing the algorithm.
The layer testunion has two fields: nummer (integer, primary key) and geom (geometry).
When I try to define the geometry column, even DB Manager has problems:

View/Panels/Log Messages/PostGIS:
WARNING    Invalid PostgreSQL layer


texts: Find nearest point with the same attribute

texts:
I have a points data set of toponyms with several thousand elements. There are a lot of duplicate names that fall in two categories:

legit ones, i.e. different places with the same name at several km of distance
points that have been entered two (rarely more) times at a somewhat slightly different position but actually refer to the same place

I'd like to filter data to only show the nearest homonym (which likely are in category 2) so that I can decide which one keep, without the distant ones getting in the way. I looked at all the various clustering, matrix and nearest neighbor tools I found in QGIS (including SAGA and GRASS in the toolbox) but they all work on the entire data set, while I want to measure the distance only within points that share the same value in the name column.
Any tip? I can switch to other (free) software if necesssary


texts: Eliminate red polygons and join adjacent ones

texts: How I can eliminate the red polygons and join them; the adjacent ones biggest one.


texts: Adding attributes to a new layer in pyqgis

texts:
I am trying to create a temporary layer in a script but I can't get the attributes defined:

# create layer
    vl = QgsVectorLayer("Points", name, "memory")
    vl.addAttribute(QgsField('name', 10))

    print [f.name() for f in vl.pendingFields()]

  prints []

Clearly I am missing something

I am working off code from replacement-of-qvariant-and-setattributemap-in-pyqgis


texts:
I want to use QGIS 3.4.13 to convert las files into xyz-text files.
So far I installed the LAStools Toolbox in QGIS and I think I can use the las2txt tool for the conversion. However I get the error message

The system cannot find the path specified

According to this link I'm supposed to download the LAStools Software from here. However, the software requires to pay a licence fee (pricing), but I am looking for a cost free solution.
In this chat a user said:

I would suggest using an open source tool.   LASTools is a great open source tool for handling and converting LAS files.

You would need to use the LAS2TXT function and you could get an ascii
text file of the lidar data.

Thus, this particular tool is for free after all? But then why do I get the error message?
Perhabs there is an alternative to LAStools to convert las data into xyz?


texts:
Im trying to get familiar with spatial data using sql server 2008.

I use book: http://www.beginningspatial.com/

and there is sample data: http://www.census.gov/geo/cob/bdy/zt/z500shp/zt06_d00_shp.zip

I uploaded it to my sql server database and Spatial results are presented well. Everything seems to be ok but when I tried to calculate distance between 2 geometries result is really small.

For example 0.27

The same is if I try to calculate STLength of given region. Result also is very small like 0.19 for example.

I wonder why the results are so small. I used SRID 4269 importing data.

Does anyone have any idea why is it like that ? I think that the result whould be in meters.

Thanks for any advice


texts: What does the blocking option do in GRASS r.watershed?

texts: What does the blocking option do to the r.watershed program in grass? According to the manual, it should stop water from being routed past that location, but it appears to have no effect. I'm trying to determine if this is something I'm not understanding, or a bug in the program, but have not had luck with the grass mailing list (no replies at the time of this writing, see thread here)


texts: Calculate the value of two different fields (int), and summarize those two into a new field in attribute table

texts:
I am trying to make a pie chart, that means I have to connect some fields in AT.
How do I sum two integer fields into a new created field. The sum should also be connected to the different object id.

The picture show the type of calculation I have used, and the field bef_0_4 and bef_10_14 wont sum into the new f14 field? Any clue of what has happened? The first row should sum into 70, but it calculated 62 after using the field calculator.


texts: Viewing/creating/editing Google My Maps KML files on Android

texts:
I have two large (several hundred features) maps created by Google My Maps that I use for my field research. They contain important information required by the field crew for their work on site.

Currently these maps are accessed via iPhone by using a free app called My Maps Editor. This allows each individual to view all the features as well as create and edit new features as needed, that the rest of the crew will view.

However, I have been wanting to switch from iPhone to Android for a very long time. Currently the ONLY thing holding me back is that I can't get the same functionality on Android that I can on Iphone. Yes, Google's native maps and Google Earth apps allow you to VIEW the My Maps, but there are several limitations which make this essentially useless. What I need is to be able to:

View multiple maps overlaid (Google only shows one at a time)
View ALL features on each map (Google caps off at around 200 and arbitrarily decides which ones to show)
View and edit these features (Google does not match the names with the features, and does not allow existing features to be edited).

Is anyone aware of an app that allows this on Android?

Alternatively, is anyone aware of another strategy that might allow KML files to be accessed via mobile phones over 3G?


texts: Calculating centroid latitude in ArcMap throws error 000539 in Arcpy?

texts:
I can create a centroid latitude field in ArcMap10.1 using "Calculate Geometry..." in the attribute table, but when I try the same thing using:
arcpy.CalculateField_management(bins_WGS84, "cent_lat", "!SHAPE.CENTROID.Y!", "PYTHON_9.3"),

I get the following error:  arcgisscripting.ExecuteError: ERROR 000539: Error running expression: GPVARIANTOBJECT0.CENTROID.Y.  It only happens occasionally, but kills my script and costs me lots of time/data. On the particular dataset I'm testing, it happens about 10 rows in, although the geometry looks fine and Arcmap doesn't complain when finding the centroid.

I know it is possible, because it's working for the same dataset in ArcMap, however, this action does not show up in the Results window.


texts: Interpreting Sentinel 1 Imagery from Alaska Satellite Facility

texts:
The Alaska Satellite Facility's Vertex data portal provides Sentinel 1A/B data, and for some granules they also provide false coloring. Is there some documentation or idea of how they do this coloration?

Here is an example with the following granule (from March 9, 2018). Is it delineating ice from water here?

S1A_IW_OCN__2SDV_20180309T171806_20180309T171831_020941_023EFA_5CE0


texts: Assigning new values to points from its position?

texts:
I want to classify a set of points as it is shown on the figure, from the upper part of the map (First point), to the lower part (last point). The question is that I want to be the FID=0 the first point in the figure, the FID=1, the second one,etc.

I read something about modifying the FID,but in ArcGIS I don't know how to do it. It classifies the points with apparently no order, and the FID= 0 is on the middle part of the figure. The points were calculated decomposing a vector line.


texts: Displaying a color composite using GEE python API

texts:
I want to display an color composite image from my image collection:
listOfImages = dataset_fxm.toList(dataset_fxm.size())

firstImage = ee.Image(listOfImages.get(0))

Map.addLayer(firstImage,{band:['B5','B4','B3']},"patan")

But the following error is showing up: name 'band' is not defined.
And I have also used 'bands' for defining the bands, but still the same error.
I am running the task on jupyter notebook python API GEE.


texts: Hiding NULL text in attributes using QGIS?

texts:
In QGIS 3, how do I hide the 'NULL' text from empty attributes?

QGIS 2.18 has a setting for that, but can't find it in QGIS 3.4.


texts: How to Exclude Processes from Iteration in ModelBuilder?

texts:
Is it possible to break Iterate Raster after certain set of processes ?

I have a list of Image files which I am processing using Iterate Rasters and performing aggregate and Extract Values to Points and then writing the output to a textfile.

I have three processes:

Aggregate
Extract Values to Points
Writing to textfile

I want to execute process 1 and 2 in Iterate, but 3 Process I don't want to perform in the Iteration.  Instead, I want to perform process 3 after all the Iteration got completed for the 1 and 2 process in the same Model.

Anybody has any idea how can I achieve this in the same model?


texts: QGIS 2.99 dev python import processing fails with ImportError

texts:
I've just installed one of the nightly build versions of QGIS 3.0 to have a play with, but on start-up, it's complaining that:

Python error: Couldn't load plugin 'processing' See message log (Python Error) for more details.

The full details from the log messages panel is like so:

  Traceback (most recent call last):
  File "C:\OSGEO4~1\apps\Python36\lib\code.py", line 91, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "C:/OSGEO4~1/apps/qgis-dev/./python\qgis\utils.py", line 647, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:/OSGEO4~1/apps/qgis-dev/./python/plugins\processing\__init__.py", line 29, in <module>
    from processing.tools.general import *              # NOQA
  File "C:/OSGEO4~1/apps/qgis-dev/./python\qgis\utils.py", line 647, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:/OSGEO4~1/apps/qgis-dev/./python/plugins\processing\tools\general.py", line 39, in <module>
    from processing.core.Processing import Processing
  File "C:/OSGEO4~1/apps/qgis-dev/./python\qgis\utils.py", line 647, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:/OSGEO4~1/apps/qgis-dev/./python/plugins\processing\core\Processing.py", line 45, in <module>
    from processing.script.ScriptUtils import ScriptUtils
  File "C:/OSGEO4~1/apps/qgis-dev/./python\qgis\utils.py", line 647, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:/OSGEO4~1/apps/qgis-dev/./python/plugins\processing\script\ScriptUtils.py", line 34, in <module>
    from processing.script.ScriptAlgorithm import ScriptAlgorithm
  File "C:/OSGEO4~1/apps/qgis-dev/./python\qgis\utils.py", line 647, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:/OSGEO4~1/apps/qgis-dev/./python/plugins\processing\script\ScriptAlgorithm.py", line 40, in <module>
    from processing.core.GeoAlgorithm import GeoAlgorithm
  File "C:/OSGEO4~1/apps/qgis-dev/./python\qgis\utils.py", line 647, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:/OSGEO4~1/apps/qgis-dev/./python/plugins\processing\core\GeoAlgorithm.py", line 51, in <module>
    from processing.algs.help import shortHelp
  File "C:/OSGEO4~1/apps/qgis-dev/./python\qgis\utils.py", line 647, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:/OSGEO4~1/apps/qgis-dev/./python/plugins\processing\algs\help\__init__.py", line 28, in <module>
    import yaml
  File "C:/OSGEO4~1/apps/qgis-dev/./python\qgis\utils.py", line 647, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
ImportError: bad magic number in 'yaml': b'\x03\xf3\r\n'

My question is, is this maybe a problem with my installation (I installed through the OSGeo4W setup tool, selected Advanced Install and chose qgis-dev with automatic installation of other dependencies), as in I missed something out? Or possibly there is some incompatibility issues as I still have QGIS 2.18 installed and it's not liking the fact that I'm mixing both versions?


texts:
I'm trying to display USPS Shipping Zone information but the section of the USPS site I used to scrape this information from has been broken for months and I'm wondering if anyone knows another way to get it in a parsable format?

The page that has been broken for months is here: http://postcalc.usps.gov/Zonecharts/ZoneChartPrintable.aspx?zipcode=931

They have the same chart in a not-so-parsable format here: http://postcalc.usps.gov/Zonecharts/

I do cache this information so it doesn't have to be requested constantly, but the cache is cleared occasionally since the zones do change from time to time.


texts: remove values in list using logical operators

texts:
Just would like to ask how can I masked or remove the values in my list based on logical operators.

Here is the case:

I converted an NDVI.img file into a 1 dimensional array. and as part of the preprocessing, I would like to remove NDVI values in my array that are less than 0.3 and greater than 0.8. Values greater than and equal to 0.3 to 0.8 are only my interest.


texts: How to prevent accidential deleting attributes and geometries by del/backspace in form view?

texts:
When editing the attribute table of a shapefile in QGIS (2.14.8 ltr) using a value relation field (one or multiple choices) in formular view, hitting Ctrl-A selects all attribute table entries; subsequently hitting Backspace or Delete deletes all content - geometry and attributes. There is no warning. There is no way to restore the entry besides turning off the editing mode without saving your changes.

I'm not even sure if this is intended behaviour. Anyone knows a way to prevent that?

FTR, I also lost single entries by hitting Del or Backspace when entering data, and this inspired my post. However, I couldn't reproduce this on purpose, but if I find out how, I'll come back to this post for an update.


texts: Adding "noise" to overlapping X,Y coordinates so no longer in EXACT same place?

texts:
I have a series of points representing animal observations. Many of them overlap, having the exact same coordinates, and this is causing me trouble when I want to run some analytical tools which requite no identical points. I was wondering if there's any tool that can add some "noise" to the coordinates, that's to say, adding some random decimal values at the end in order to obtain slightly different coordinates which are still close enough in space so as to not really modify the observation. I've heard about such a tool, but I can't find it, at least in ArcGis 10. Maybe there's another software that can do it?
Thanks!


texts: Changed attribute order of PostGIS layer in QGIS

texts:
I'm working in a PostgreSQL table (PostGIS) that I loaded into QGIS. This layer has been created by a SQL Statement in pgAdmin.

I noticed that the order of attributes shown in the attribute table differs from that I declared with the SQL Statement. But checking the layer properties in QGIS (tab "fields"), the order is correct.

To make it more weird, I tested it with QGIS 2.14, 2.18 and 3.6 and the wrong order in the attribute table doesn't occur in 2.14, but in the others.

Does anyone have an idea why order changes?


texts: Using ArcGIS Desktop v10.1 with an Advanced level license we have georeferenced our engineering drawings in PNG format. Next we have created a polygon layer for every single drawing to see the Area of Interest depicted in the drawing. Is there an automated way to bring individual PNG layers by clicking on their respective polygon? I want a simpler way to add the georeferenced PNG into the map without having to use the classic "add layer" tool, where a user has to find the respective layer in their location on the network.


texts: Using hydrological equations or DEM map to calculate area inundated for 1m sea level rise scenario?

texts:
I want to know how I can calculate the area inundated for a 1m sea level rise (SLR) scenario in ArcGIS Desktop and also the procedure of creating 1m SLR for a particular area.

Do I need to use hydrological equations to get area inundated or by directly from DEM map I can get it?


texts: How to homogenize NA values within a .gpkg from R

texts:
I have created an SpatialPointsDataFrame in R with some data.
Here is the outputs from R:

> class(Fotos)
[1] "SpatialPointsDataFrame"
attr(,"package")
[1] "sp"

> class(Fotos$Azimuth)
[1] "numeric"

> Fotos$Azimuth
  [1]  NA  NA  NA  NA   0  NA  NA  NA  NA  NA
 [11]  45  NA  90  NA  NA  NA 135  NA  NA  NA
 [21]  NA 180  NA  NA  NA  NA  NA  NA 225  NA
 [31]  NA  NA  NA  NA 270 315  NA  NA  NA  NA
 ...

> class(Fotos$labelpos)
[1] "numeric"

> Fotos$labelpos
  [1] NA NA NA NA  0 NA NA NA NA NA  0 NA  0 NA
 [15] NA NA  0 NA NA NA NA  0 NA NA NA NA NA NA
 [29]  0 NA NA NA NA NA  0  1 NA NA NA NA NA NA
 [43] NA NA NA NA NA NA  0 NA NA NA NA NA  1 NA
 ...

It seems that both Attributes have the same data structure and class but when I export this SPDF to .gpkg some attributes replace NA by NULL and others replace NA by nothing. Like here:

How could I homogeneize those null/na/empty values in the .gpkg exported from R?
I tryed something like this without fourtune:

Fotos@data[is.na(Fotos@data)] <- NA

EDIT:
I tried to homogenize those NA with this:

 # Homonegeneizar NAs
    storage.mode(Fotos$OrdenSel) <- "integer"
    storage.mode(Fotos$Plano) <- "integer"
    storage.mode(Fotos$Cluster) <- "integer"
    storage.mode(Fotos$labelpos) <- "integer"
    storage.mode(Fotos$Azimuth) <- "integer"

But QGIS read those values as integers, but some has NULL and others empty...


texts: Unable to Make PostGIS 2.0 - file not found

texts:
I am struggling to get the latest trunk of PostGIS to install on a fresh Mac running Lion.  I did install PostGIS 1.5.3-2 for Postgres 9.1 from [KyngChaos Mac OSX Installers][1].  This setup worked fine for me and I was able to setup spatial databases against a PostGreSQL 9.1 database instance.  But I realized I need the new GeoJSON support in the PostGIS 2.0 trunk.  I went ahead and downloaded the latest source code(revision 8401) and compiled with the JSON optional argument (downloaded JSON-C libraries):

sudo ./configure --with-gui --with-raster --with-jsondir=/usr/local/

The compile works successfully.  My issue is when issuing the make command.  It appears that the make command is looking for a /usr/bin/postgres/ directory on my system which does not exist:

ld: file not found: **/usr/bin/postgres**
collect2: ld returned 1 exit status
ld: file not found: /usr/bin/postgres
collect2: ld returned 1 exit status
lipo: can't open input file: /var/tmp//cc82cOrG.out (No such file or    directory)
make[1]: *** [postgis-2.0.so] Error 1
make: *** [all] Error 1

I am thus unable to install PostGIS 2.0 on my system.  Anyone have any ideas what specific files (libraries, etc) the PostGIS 2.0 install is looking for under /usr/bin/postgres?

Thanks.


texts: What are the attributes of the lines my line connects with?

texts:
I have a series of lines with different attributes and have used Generate Near Table, Summary Statistics and Join Field to add a column with the number of lines each line is near to.

(I have had to use 'Near' as the lines do not always intersect).

    FID | Attribute | Near Count |
    1   | Red       | 2          |
    2   | Red       | 1          |
    3   | Blue      | 3          |

I want to find out how many times each line connects with a line with a specific attribute.

    FID | Attribute | Near Count | Near to Red Count | Near to Blue Count |
    1   | Red       | 2          | 1                 | 1                  |
    2   | Red       | 1          | 1                 | 0                  |
    3   | Blue      | 3          | 0                 | 3                  |

How would I go about adding these columns to my data?


texts: Identify closed polylines in ArcGIS using Python

texts:
Running ArcGIS 10.2 here and a complete beginner with Python. I'm trying to figure out a workflow to test for closed polylines using the field calculator. Essentially what I'd like to do is to populate a specific field with "closed" or "open". This is what I've come across and if I understand it correctly, it loops through each row and checks if the first x coordinate matches the last x coordinate, if there is a match then it prints the OID.

How do I convert this to a field calculator expression?

Do I need to use an update cursor rather than search to populate my new field?

fc = "Boundary"

rows = arcpy.SearchCursor(fc)
for row in rows:
    geom = row.shape
    if geom.firstPoint.X == geom.lastPoint.X:
        print row.OBJECTID

del row, row


texts: Calculating sum of distances between multiple points in ArcMap?

texts:
I am using ArcMap 10.2.

I have two layers- one is around 2,000 geolocated ballot boxes while the other is around a dozen geolocated facilities.  I would like to calculate the total distance between each ballot box and all dozen or so facilities.  For instance, Ballot Box 1 is 10km from Facility X, 15km from Facility Y, and 20km from Facility Z, and so on.  So Ballot Box 1 would be assigned a score of 45km.  I'd like to calculate this score for all 2000 or so ballot boxes.

So far, I've used the point distance tool where the ballot boxes are the input features and the facilities are the near features, left the search radius blank, and ArcGIS produces a table, with a "DISTANCE" column populated, so it seems to have worked.  My question is threefold:

Is this the correct procedure?
In what unit is this distance reported (Projections are UTM, so I assume it defaults to meters?)
If I leave the "search radius" blank, does ArcGIS default to figure the distances in the layer's full extent?


texts:
I am trying to get labels that look like the image below

I have been able to come up with something similar. however the leader arrow takes on the callout box boundary properties. as seen in the image below

How can I go about this?

I am working with ArcGIS Desktop 10.2.2


texts: Where to store lat/lon data with additional informations?

texts:
I have huge database of lat/lon GPS coordinate points. For every point, I have several other variables (height, surface type, temperature, pressure...).

Also, data will be changed every day and old ones drop, no history.

I want to be able for given point find N-closest points in my database. Should I use PostGIS for this, or classic PostreSQL or MySQL is enough? I dont need any other data transformations, only finding closest (and even not by real distance in km). This is supported also by classic database systems and I am not sure if PostGIS is not like going with cannon to pigeon.


texts:
Using ArcGIS 10.1 for Desktop.

I created perpendicular lines along a polyline.  In the process of doing that, the lines were numbered in order (automatically).  I manually edited some of the perpendicular lines which re-assigned them a number and therefore the ID's aren't increasing numerically in order along the polyline.  I then used the interesect tool to create a point at the intersection of each perpendicular line along the polyline.  I want to create a polyline from these data points, but since the id's are not in numeric order the line is not being created correctly.

Question:
How can I assign a number to a data point that increases as the points move along the polyline?

I have also tried using ET Geowizards but I don't have a full license and am limited to 100 points but this data set is approximately 20,000 points.


texts: HTTP WMS GetMap using default SLD instead of provided one

texts:
Im trying to give a sld with my wms getMap request.
My issue is that instead of using the provided sld it uses the default one.
Am I missing something obvious?
Here is the request:

http://url.dev/geoserver/workspace/wms
?SERVICE=WMS
&VERSION=1.3.0
&SLD=http://url.dev/sld/test.xml
&REQUEST=GetMap
&FORMAT=image/png
&TRANSPARENT=true
&LAYERS=workspace:layername
&transparent=true
&format=application/openlayers
&viewparams=trace:6
&t=355
&WIDTH=256
&HEIGHT=256
&SRS=EPSG:3857
&BBOX=8140237.764258131,5009377.085697312,8766409.899970295,5635549.221409475

with sld body beeing a simple red point :

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0"
 xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
 xmlns="http://www.opengis.net/sld"
 xmlns:ogc="http://www.opengis.net/ogc"
 xmlns:xlink="http://www.w3.org/1999/xlink"
 xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <!-- a Named Layer is the basic building block of an SLD document -->
  <NamedLayer>
    <Name>default_point</Name>
    <UserStyle>
    <!-- Styles can have names, titles and abstracts -->
      <Title>Default Point</Title>
      <Abstract>A sample style that draws a point</Abstract>
      <!-- FeatureTypeStyles describe how to render different features -->
      <!-- A FeatureTypeStyle for rendering points -->
      <FeatureTypeStyle>
        <Rule>
          <Name>rule1</Name>
          <Title>Red Square</Title>
          <Abstract>A 6 pixel square with a red fill and no stroke</Abstract>
            <PointSymbolizer>
              <Graphic>
                <Mark>
                  <WellKnownName>square</WellKnownName>
                  <Fill>
                    <CssParameter name="fill">#FF0000</CssParameter>
                  </Fill>
                </Mark>
              <Size>6</Size>
            </Graphic>
          </PointSymbolizer>
        </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>

Using geoserver 2.8.


texts: Raster projects correctly in ArcMap, Cannot align in QGIS

texts:
I have a geoTIFF with orthorectified aerial imagery in .sid format that does not project correctly in QGIS.  The file originally did not have a .prj file, but projection information is stored in the file metadata, listed as "NAD_1927_Washington_North"

The file opens up and projects fine in ArcMap 10.4 (it lines up fine with other aerials, rasters, and vector data), but I don't have regular access to ArcMap.  I tried exporting a .prj file for the projection from ArcMap with it correctly aligned, but that doesn't make a difference when I open it in QGIS.

My project is in EPSG:102748- Washington StatePlane North NAD 83 US Feet I have on the fly projection turned on.  I have several different layer in different projections (e.g. UTM zone 10N, WGS 84) that all project fine, but I can't make this aerial line up.

I tried opening the raster in a clean workspace, then adding other data, but this doesn't seem to solve the problem.

I tried assigning a number of different projections, but couldn't find anything that worked.  In EPSG: 32148, NAD83/Washington North, I can make it line up about 150 km east of the "correct" location, but can't do better without georeferencing the image from scratch.

I have intermittent access to ArcMap, so if there is something I can do to reproject the file there, it would be fine, but need to be able to use the layer in QGIS.


texts: Network analyst scaled cost barriers in ArcGIS JS API or REST API

texts: I am developping a GIS web application using the ArcGIS JS API. I know that scaled cost barriers are available in ArcMap, but I am wondering if they exist in the JS API or in the REST API and if not I am wondering if there are any alternatives


texts: Subset a shapefile that excludes certain field values by using NOT in Select by Attribute where clause

texts:
I am using ArcPy 2.7 to subset a shapefile. I am trying to use arcpy.Select_analysis to exclude certain features from a shapefile (idle fields denoted as 1401 in the field 'Crop_Type'. my 'Crop_Type' field is of type long.

I'm receiving the following error:
ExecuteError: ERROR 000358: Invalid expression  NOT "Crop_Type" = '1401'

Looking at ESRI's Building a SQL Expression, my code is:

import arcpy
from arcpy.sa import *

arcpy.env.overwriteOutput =  True

in_features = 'Copy_SaltonT4_2017.shp'
out_feature_class = 'CroppedFields.shp'
where_clause = 'NOT "Crop_Type" = \'1401\''

arcpy.Select_analysis(in_features, out_feature_class, where_clause)

arcpy.CopyFeatures_management(in_features, out_feature_class)


texts:
I want to apply a GLCM texture filter to some aerial photographs that I have. I am using the skimage package to compute the GLCM and a package called rios to handle the i/o.

Rios - raster input output simplification - reads imagery in as a 400 x 400 x nbands numpy array and handles the creation of the output dataset paramaters. It is designed to handle very large datasets in a memory efficient way and also enables the user to set a moving window. I understand that I need to use caution when reading an image in in blocks and applying a moving window but I am not at the stage where I am implementing a window yet.

I can get my code to run but the result is an output made up of zeros with the occasional grey pixel. Ultimately I don't understand what the parameters I am using do. I consulted the documentation but it doesn't give any more than a brief description or examples of what different parameter values will achieve.

Thus far, I have:

# Set up input and output filenames.
infiles = applier.FilenameAssociations()
infiles.image1 = "infile.tif"

outfiles = applier.FilenameAssociations()
outfiles.outimage = "outfile.tif"

controls = applier.ApplierControls()
controls.progress = cuiprogress.CUIProgressBar()

# Apply a moving window if required
# 3x3 the overlap is 1, 5x5 overlap is 2 etc
# I have yet to try this as i'm trying to get something basic working first
#controls.setOverlap(4)

# Apply Texture
def doFilter(info, infiles, outfiles, controls=controls):
    g = greycomatrix(infiles.image1[3], [1], [0], 256, normed=True, symmetric=True)
    filtered = greycoprops(g, 'contrast')
    # create 3d image from 2d array
    outfiles.outimage = numpy.expand_dims(filtered, axis=0)

applier.apply(doFilter, infiles, outfiles, controls=controls)

Could someone help explain what the parameters within the greycomatrix mean and how I would expect the GLCM to behave under different parameter values? They are referred to in the documentation as:

greycomatrix(input array, distances, angles, levels, normed=True, symmetric=True)

The parameters distances, angles and levels are the ones I understand the least.


texts: Unable to compile using make test(windows)

texts:
While trying to compile my test plugin I encounter this error. At the moment i can't figure out why it isn't working. I m using a windows and try to execute the command in OSGeo Shell.

I typed in OSGeo Shell: make test and got:

pyrcc4 -o resources_rc.py  resources.qrc
make: pyrcc4: Command not found
make: *** [resources_rc.py] Error 127

From the looks of it, I do not have the make command installed?


texts:
I would like to add a band to each image in an image collection ('median_ndvi' in the linked code). The band should be a sequence of ee.Image() constants from 1 to 12.

https://code.earthengine.google.com/e6f22ac8bc6763aa795a14272c89ce6f


texts:
I'm trying to find the polygons a point is in and then the intersecting polygons of the result of the first query.

I'm trying along the lines of:

WITH RECURSIVE all_the_areas AS (
    SELECT a.name,a.geom FROM areas a WHERE ST_Contains(a.geom, ST_SetSRID(ST_MakePoint(lat,lng),4326))
  UNION ALL
    SELECT b.name,b.geom FROM areas b, all_the_areas c
    WHERE st_intersects(b.geom, c.geom)

)
SELECT * FROM all_the_areas;

The problem is that this query hangs plus I'm not experienced at all with this type of query so it's probably a poor way to solve the problem.

Can you help me find an efficient solution to this problem?


texts: ArcMap: export each x/y coordinates to a separate row

texts:
I am trying to sort out how the coordinates of each point in a multipoint or polyline feature could be exported to a .txt with python so the exported file looks like this:

point1_coordinates_x
point1_coordinates_y
point2_coordinates_x
point2_coordinates_y
.
.
.
pointn_coordinates_x
pointn_coordinates_y

I don't need the elevation / z value in further processing of the data, only x & y coordinates. Is there any way to easily generate this output with python?

I know of the export to csv / txt settings in ArcGIS, also of the standard python code matching that function - however I struggle with putting each value in a separate row.


texts: Earth Engine - returning/passing values between UI widgets

texts:
Runnable code snippet here: https://code.earthengine.google.com/24f4e6e5b961e18c7200fd58884fff17

I am building an app that will display a geometry and some data on the map based on user input from multiple ui.Select dropdown menus. My goal is to have the user choose her parameters from each dropdown, then click a button labeled Submit to add the geometry/data to the map all at once. Currently I can add the geometry to the map as part of the onChange function, but I am having issues returning a value from that function to pass to another widget. Here is a code snippet:

var map = ui.Map();
var states = ee.FeatureCollection('TIGER/2018/States').toList(60);

function makeDropdown() {
  var dropdownItems = states.map(function(el){
    return ee.Feature(el).get('NAME');
  }).sort();
  var placeholder = 'Select a state...';

  // Dropdown menu
  var select = ui.Select({
    items: dropdownItems.getInfo(),
    onChange: function(key) {
      var newGeo = states.filter(ee.Filter.eq('NAME', key));
      var newGeoLayer = ee.Feature(newGeo.get(0)).geometry();
      // Uncomment to add state to map - this works.
      // map.addLayer(newGeoLayer, null, 'Selected state');
      // map.centerObject(newGeoLayer, 6);
    }
  });

  // Submit button - this doesn't work.
  var submit = ui.Button({
    label: 'Submit',
    onClick: function(newGeoLayer) {
      map.addLayer(newGeoLayer, null, 'Selected region');
      map.centerObject(newGeoLayer, 6);
    }
  });

  select.setPlaceholder(placeholder);
  return ui.Panel([select, submit], ui.Panel.Layout.flow('horizontal'));
}

// Initialize the split panel UI.
function init() {
  var panel = ui.Panel({
    style: {width: '30%', border: '3px solid #ddd'}
  });
  var dropdownPanel = ui.Panel();
  panel.add(dropdownPanel);
  var splitPanel = ui.SplitPanel({
    firstPanel: panel,
    secondPanel: map,
  });

  ui.root.clear();
  ui.root.add(splitPanel);
  dropdownPanel.add(makeDropdown());
  return dropdownPanel;
}

var panels = init();

What I am hoping to accomplish in this example is to define newGeoLayer inside the onChange function of the ui.Select widget, but also return it somehow so I can pass it to the ui.Button widget. Right now the newGeoLayer variable is only defined inside the function and is out of scope for the button widget.


texts: Workspace Error type- 000301 for reconcile version

texts:
I have a arcpy script to reconcile and post edits from 1 child to 1 parent 9the parent is not the default version).
my code:

import arcpy

print "Preparing to begin Water Version Reconile and POST..."

workspace = r"C:\Users\path\to\my\GDB.sde"
arcpy.env.workspace = workspace
childVer = "sde.child_version"
parentVer = "sde.parent_version"
print 'Creating list of Database versions to work on...'

arcpy.ReconcileVersions_management(workspace, "BLOCKING_VERSIONS",     parentVer, childVer, "LOCK_ACQUIRED", "NO_ABORT", "BY_OBJECT",     "FAVOR_EDIT_VERSION", "POST", "KEEP_VERSION", "c:/temp/reconcilelog.txt")

And the error I receive:

File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\management.py", line 16614, in ReconcileVersions
    raise e
ExecuteError: ERROR 000301: The workspace is of the wrong type
Failed to execute (ReconcileVersions).


texts: Automatically update the result of a QueryTask

texts:
I have developed a application using the ESRI JS API to allow users to create and edit features in a web map.  I have used the QueryTask to create summary statistics to report back to the user on total length and such.  To create this report back, I am using .innerHTML, which works great, but the statistic is not refreshed when a change is made to the feature.  Thus in order for a user to see updated statistics they need to refresh the page.  This is not ideal - can I add a listener or other technique to automatically update the statistics generated with the Query Task?

The Query and Query Task is:

queryTask = new QueryTask("URL");
    query = new Query();
    statisticDefinition = new StatisticDefinition();
    statisticDefinition.statisticType = "sum";
    statisticDefinition.onStatisticField = "SHAPE.STArea()";
    statisticDefinition.outStatisticFieldName = "total_area";
    query.outStatistics = [statisticDefinition];
    queryTask.execute(query, function (result) {
        total_area = result.features[0].attributes;
        postExecute();

    });

and the postExecute Function that gives the statistic to the DIV:

function postExecute() {
    var roadLength = (total['total'])";
    document.getElementById("length").innerHTML = roadLength;
}


texts:
How do you SVG-marker layers in PyQGIS? I'm importing a vector layer and want to display it with SVG-markers.
It seems that QgsSvgMarkerSymbolLayerV2 does not exist in PyQGIS...


texts:
After running the query

UPDATE parcel
SET wkb_geometry = ST_MakeValid(wkb_geometry)
WHERE ST_IsValid(wkb_geometry) = 'f'

I receive the following error:

NOTICE:  Too few points in geometry component at or near point -121.78764087332524 37.700189714113364
ERROR:  Geometry type (GeometryCollection) does not match column type (MultiPolygon)

However, this is even after I imported the GeoJSON file with the ogr2ogr command -nlt PROMOTE_TO_MULTI, which is supposed to covert polygons to multipolygons  and  running

ALTER TABLE parcel ALTER COLUMN wkb_geometry
    SET DATA TYPE geometry(MultiPolygon) USING ST_Multi(wkb_geometry);

in postgres.  What might be going wrong and how can I remove and/or alter the polygon that has too few points?


texts: What exactly is the memory impact on the database when reducing cells size in a Raster GIS?

texts:
I understand that the size of the memory of a raster database increases when reducing the cell resolution.
I read some about it in the ArcGIS Help.

However, I am wondering about the exact relation on the memory footprint and was hoping that someone with more database-experience can tell me whether my reasoning here makes sense:

Given cells of 10x10 meters, and reducing the resolution to 5x5 meters, like so:

The relation for the increase of number of cells would be n*4, so I thought that the memory increase in the database file would follow a similar pattern.

Here's my attempt to generalize this:

  The amount of cells and the memory size of the database both increase by the square of the number that the
  sides of the initial raster cells are divided by.

E.g. 10/2 (2 being the divisor, thus 2^2 = 4 being the factor in this case) results in *4 the amount of cells - so does it also result in *4 the amount of memory size of the database?

Given that my assumption is correct, dividing a raster cell by 4, e.g. 10/4 (with 4 being the divisor thus 4^2 = 16 being the factor), would  result in *16 the amount of cells, and *16 the amount of memory size of the database!
That would be a large increase of file size pretty quickly!

Is this generalized formula applicable or does the size of the database increase following a different pattern - and if so, which one?


texts: Getting flood depth grid values for road segments?

texts: I have a 2 features, the first contains values for flood depth (Raster, GRID). the other is a road network (vector, polyline) that has been clipped to the boundary of the depth grid. I am racking my brain trying to figure out how to get the maximum flood depth for each road segment where it crosses the Flood Depth GRID and I am unable to come up with a method to achieve this. I am using ArcGIS 10.2 Basic, but do have access to an advanced license if necessary.


texts:
I need to develop a couple of qgis-plugins. I made some promising first steps and completed a few easy plugins. Now its time to move to more complex plugins which need to import several tool-classes I developed.

As this is quite a lot of stuff and most of my plugins need to use these libraries I wonder wether it is possible to copy them somewhere where all plugins could import them from. This needs to be a folder which is accessable for the current user. I thought about using something like "/home/username/.qgis/python/plugins/includes" but I couldn't get this working.

Can anybody give me a tip?

Thanks and regards

Ok - hello again,

on writing it occured to me - I can use "/home/username/.qgis/python" to put my classes in - I guess that's why its called like that!

Thanks again


texts:
I have been looking for this for a while and can't find the answer: when you startup geoserver a webpage is opened at localhost:8080/geoserver/web but if you have changed your port config in jetty.xml then this is the wrong address... still, when you start up geoserver opens the page at 8080, which fails now of course...

here is my question: where is this behaviour set in geoserver?

while i am at it, is it possible for geoserver to autostart?

thanks in advance

update: this is for running geoserver on an OSX Machine


texts:
I'm trying to map an area by using ggmap. The problem is that ggmap only plots a square map. What I'd like to do is to expand just a little bit the map, so the map would fit perfectly. Here is my code. You can download the shapefile data here.
AMB_Contour <- readOGR(dsn = '.', layer='AMB_Contour')
myprojection_latlong <- "+proj=longlat +zone=31 +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
AMB_Contour <-  spTransform(AMB_Contour, myprojection_latlong)
tract_plot <- tidy(AMB_Contour)

CenterOfMap <- geocode("41.396108,2.059348")
Ciutat_NoLabels <- get_googlemap(c(lon=CenterOfMap$lon, lat=CenterOfMap$lat),zoom = 11, maptype = "terrain", source = "google", color="bw", style=c(feature="all", element="labels",visibility="off"))
ggmap(Ciutat_NoLabels)

CityMap_NoLabels <- ggmap(Ciutat_NoLabels, extent = "device")
print(CityMap_NoLabels)

CityMapWithRoutesAndStops <- CityMap_NoLabels + geom_polygon(aes(x=long, y=lat, group=group), fill=NA, size=1.5, color='black', data=tract_plot, alpha=1)
print(CityMapWithRoutesAndStops)

In the following image the eastern part of my map doesn't fit completely. I've also tried adjusting the x-y coordinates manually but it doesn't work. I don't want to change the zoom argument because for me it's perfect, the only problem is the eastern part. Finally, I want that maptype (terrain).


texts:
Does anybody know a package who work in python 3.4 as "LatLon" Package ?

https://pypi.python.org/pypi/LatLon/1.0.1#downloads


texts: Force ArcMap to respect the attribute table entries in the relationship tabs

texts:
Is there any method to force ArcMap to respect the attribute table entries in the relationship tabs?

I have a table named “aaa” which is related to a feature class “communities”, the values in the table “aaa” are entered based on the approval year,

Now when I identify a point I was surprised that the tabs doesn’t respect the attribute table sort which is based on the approval year or the objectID


texts:
About 3-5 years ago I created a Story Map using the Map Journal template.
That Story Map is still available at https://www.arcgis.com/apps/MapJournal/index.html?appid=163b6d4177354fbb93ba9ed4aedd8ee3
However, I now want to review/update it to comply with FAQ: What do I need to know about HTTPS Only and the ArcGIS Platform? and I have not been able to locate where it is stored.
I thought that I would be able to go to ArcGIS.com, login using my Esri Global Account, and see it in my Gallery.  However, it does not appear to be there.
Is there a way to use the URL of a Story Map to find out where it is stored?


texts: Spatially join polygons (target) to polylines (join) on longest intersecting polyline

texts:
I am curious if I can, without having to write a python or an ArcObjects script, create a spatial join that will result in a polygon feature class containing the attributes of the LONGEST intersecting join polylines.   The input polylines I have roughly correspond to the polygons to which I want to transfer the attributes.  I've determined that it is almost always the longest intersecting polyline whose attribute I need to transfer into the polygon.

Currently, using the figure below as an example, the spatial join hits 3 polylines (join_count = 3 for this particular polygon) but instead of joining on the longest intersecting polyline it seems to choose the first one it encounters. (which is also contrary to what the tool's description says about the JOIN_ONE_TO_ONE operation which is suppose to aggregate the join features into the output feature. The description speaks of SUM merge rule in its example which I do not see in this particular tool as an option)

I am using INTERSECT because there is no suitable spatial relationship I can think of in this "polyline to polygon" situation.  I understand what is going on but don't know whether a workaround exists. I just wish there would be another set sub-parameters allow me to control what happens when more than one join features exists; in this case selecting the longest polyline by default.

I am using ArcGIS 10.2.1 Standard.


texts: What would be the Overpass query to export minimum bounding rectangles (MBR) of cities/provinces/countries from OpenStreetMap?


texts: Cannot load image by ID Sentinel-2 Google Earth Engine

texts:
I am trying to start using GEE. However, I cannot visualize an image when using the Code Editor. Here is my code:

var s2 = ee.Image('COPERNICUS/S2_SR/20190314T155521_20190314T205515_T17PLK')
Map.addLayer(s2, {bands: ['B8', 'B11', 'B4'], max: 4000}, 'S2')
Map.centerObject(s2, 8);

I followed the recommendation found in this post:

https://gis.stackexchange.com/questions/200043/how-do-i-find-sentinel-2-dataset-id-on-google-earth-engine

However, it does not work. The image I want to visualize is this one:

https://scihub.copernicus.eu/dhus/odata/v1/Products('d3c42954-1eb8-4fb2-b297-061f1f9c0d19')/$value
```[![enter image description here][1]][1]

  [1]: https://i.stack.imgur.com/mijBJ.jpg


texts: Return raster that intersects another raster in gdal

texts:
I'm trying to programmatically orthorectify a bunch of Rasters (.ntf format), however I need a DEM that covers each raster in order to do so.

I have DEM rasters (.tif) as tiles that cover the entire region I'm trying to orthorectify, but I need an automated way to find which tiled DEM covers each individual scene. The tiles are Projected in Albers Equal Area Conic while the scenes that need to be orthorectified are in GCS84.

Something similar to 'Select by Location' in Arc (but for rasters) would be preferable but after searching the web it doesn't seem like this exists in GDAL.

I can get the Corner Coordinates of the DEM .tif from 'gdalinfo', and I can get the x and y max/min coordinates of the .ntf raster from its associated .xml file...

However the coordinates for the DEM from gdalinfo are Projected and the coordinates for the .ntf from the xml file are Geographic.

Previously I had just mosaicked the tile DEM pieces together and then used gdalwarp to clip the mosaicked DEM to the to-be-orthorectified scene's extent, but for this particular process the area is way too large to do this.

Any thoughts as to how I can (simply) find the intersecting DEM raster for each .ntf scene?


texts:
I have the following very simple code:

@map_center = [-32.951106, -60.669952]
@map = new ol.Map({
        target: 'map-canvas',
        layers: [new ol.layer.Tile({source: new ol.source.OSM()})],
        view: new ol.View({
          center: @map_center,
          zoom: 5
        })
    })

It's in coffeescript, but you will get the idea. The problem is, the map does not center at all. It gets stuck in [0, 0]
Am I doing something wrong?


texts: PostGIS: Merge multiple table with same GID

texts:
I just imported some shapefiles with the PostGIS importer. Because several errors and issues I had to split my shapefile in multiple parts. Then I imported them one by one. I thought I can easily merge them later on.

But now there is a problem. The import tool created a primary key called GID. And these identical PKs prohibit the merge.
Doeas anybody know how to delete these GIDS, merge all the tables and create new PKs?

I'm thankful for any answer.


texts: Are there roadmap datasets available for qgis?

texts: On ArcGIS.com, you can download public datasets such as road maps. Do these work with qgis or are there public datasets available for qgis? By public, I'm referring to datasets that others have created and offered for others to use.


texts: Leaflet map in R Markdown - not rendering when placed after other plots

texts:
I have an R markdown document using shiny with 6 or 7 plots. They all display until the 7th plot, which is a leaflet map, which doesn't display - when I 'run document' there is a big empty space at the bottom.

However, I have verified that the leaflet code works by running it in console, in another markdown document, and by moving it to the top of the current document - in the latter instance it does indeed render properly once placed at the top. But then the plot which is now at the bottom gets cut off instead!

I'm fairly new to markdown, but I wondered for example whether there is a maximum number of lines you can have in your document, and whether this could explain what's going on?

The global environment chunk is set simply as follows. I can post further code if helpful.

{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

To reiterate, the plot succeeds if placed at the top of my page but not the bottom after my vanilla plots. If placed at the bottom,  nothing after the leaflet appears either - there's just a massive gap of nothing.


texts:
This is part of a mathematical challenge which I am solving. In the final deciphered message, I am given coordinates which correspond to some location that I need to check.

Coordinates are: 511716656388765455430016138955706839007890052532, 1622805609316535864254436412730925222158623332074

These are very large numbers and the information obtained is correct. I verified it. Also, I have another set of coordinates similar to these which I need to look up.

I read about different Geographical Coordinate Systems here: https://en.wikiversity.org/wiki/Geographic_coordinate_conversion

However, none of them match it.

Could it be a non standard format of representing coordinates?


texts: How to to get time slice in raster dataset using ArcGIS?

texts:
I have a couple of rasters made in different time. Some of the rasters are overlapping and some of them - not.

I need fast and efficient way to get time sliced values for each pixel.

For example:

    pixel # 1
    value: 255         200         150
     date: 01.01.2015  02.01.2015  02.01.2015

    pixel # 2
    value: 240         250         160
     date: 03.02.2015  04.02.2015  05.02.2015

    and so on ...

There are any tools for this task in ArcGIS, ArcPy?

I am sure first step should be put those rasters in raster dataset or mosaic dataset. But I can't proceed further.


texts: Why is a vertical polygon defined by WKT PolygonZ invalid?

texts:
Can you tell me, why this PolygonZ is invalid?

PolygonZ ((409662 5317369 285, 409662 5317377 285, 409662 5317377 302, 409662 5317369 302, 409662 5317369 285))

The polygon is completely vertical, but due to z coordinate points are not overlaying.

From GEOS validity checker I get the following error message:

Too few points in geometry component

The coordinates are in EPSG:32632 (WGS84/UTM Zone 32N)


texts: Drawing lines from two points in CSV using QGIS?

texts:
I have a CSV file where one line looks like this:

Duration, user, lat-start, long-start, lat-end, long-end
298, Casual,38.9101,-77.0444,38.91554,-77.03818

Each row in the CSV has two points and I need to draw a line between those two points for each row.

How can I do this using QGIS 1.8.0?

I have tried using the delimited text file plug-in, but it just draws points and has no line option from what I can see.


texts: Analyzing precipitation data with conflict data in QGIS

texts:
For my Master thesis, I want to address the question, whether rainfall anomalies can increase the risk of conflicts in Sub-Saharan Africa. As a moderator I have picked climate-related aid programs. The Aid data is a .csv file and comes from http://aiddata.org. The conflict data is a .csv file also and comes from the Social Conflict Analysis Database (https://www.strausscenter.org/scad.html). My precipitation data is from the University of East Anglia (http://wps-web1.ceda.ac.uk/submit/form?proc_id=Subsetter), gridded 0.5° x 0.5°  and is a .nc data. I could also download it as a .csv dataset.
After loading the data to QGIS it looks like this

I would like to combine the data in one dataset and the grid cells as unity of analysis. So I could analyze it with STATA or R. In the journal articles I have read so far, that they always work with grid cells, otherwise they could not combine the point data from conflict datasets to rainfall data. I guess I could use the 0.5° x 0.5° raster of the precipitation data or should I use artificial grid cells. Basically, my question is, how can I combine the data to do regression analysis? I did the QGIS-Tutorial, but did not find any answers to my question. Please help a political science student :). I am happy for any help. Unfortunately, we do not have a geographical faculty at my university.


texts: Is there a parellel line function for slds?

texts:
I am attempting to make an sld style which symbolizes a line with an offset. I have tried with filters, functions online graphics but can achieve only ugly results. The line consists out of 2 line symbolizers. one thick underlying line which is displayed at an offset westerly of the achse. The other thinner line lies above the thicker line overlapping the easterly edge of the thicker line.

an example is here: http://maps.zgb.de/rrop/sld_img/style_hatchedline_offset.png

If anyone knows what I can do to produce this sld I´d be most grateful,

yours,

Rob


texts: I just started using SatScan, but I've been getting this error “Reading the coordinates file” repeatedly and can’t run the program. I don’t know what could be the problem. Someone please help.


texts: Population data to Population Density Raster Data

texts: I'm currently using QGIS 3.16 for desktop. I have the population data of a district in excel ( number of inhabitants of a district) as well as openstreetmap raster. I am wondering if there is a way to make a population density data from the current data i have using QGIS 3.16 tools?


texts: Creating TFW and PRJ files for folder of GeoTIFF files?

texts:
I have a folder which contains a large number of files in GeoTIFF format.

I need to generate both PRJ and TFW files for each of these images.

Is there some way to do this?


texts: Can we take the length of a linestring by calculating ST_Distance from its ST_StartPoint and ST_EndPoint?

texts:
I took the values of my StartPoint and EndPoint by doing this procedure:

a)

SELECT COUNT(
CASE WHEN ST_NumGeometries(geom) > 1 THEN 1 END
    ) AS multi, COUNT(geom) AS total FROM public.roads;

because AS_Text(geom) was giving me multi-linestrings.the multi was 0 so I updated the geometry type to linestring:

b)

ALTER TABLE public.roads
    ALTER COLUMN geom TYPE geometry(LineString, 4326)
    USING ST_GeometryN(geom, 1);

and finaly

c)

INSERT INTO myschema.nodes
SELECT ST_StartPoint(ST_GeometryN(ST_Multi(geometry),1)),
ST_EndPoint(ST_GeometryN(ST_Multi(geometry),1)) FROM public.roads

is that procedure right?

and after that I was looking the way to find the length of my roads.
the only command that worked to me was the:

SELECT  ST_Length(geom::geography)/1000 from public.roads;

but I thought that maybe I can find it by:

SELECT ST_DISTANCE(st_p::geography,en_p::geography)/1000,weight FROM myschema.nodes

the results are different and I'd like to understand the reason.

I have not been working on GIS for very long and I'm a little bit confused.


texts: Load non-spatial table in QGIS by querying an Oracle database

texts:
I am trying to create all (non-spatial) layer in QGIS2.16.2 Which uses a "select / join query" on Oracle database as data source.

I can successfully verify my query by executing it in QGIS DB Manager but It seems not to be possible to "load it as New Layer". After pressing the button "load now" nothing happens.

After some experiments I came to the conclusion that QGIS DB Manager only accepts simple queries (no joins ect ..) and needs a unique key column and a geometry column to load ash layer.

Unfortunately I have "read-only" permissions on the Oracle database, so creating a new view in the database is not possible.

Querying the different Oracle tables locally, using virtual layers takes to long or even makes QGIS crash due to the large data sets.

Are there other options?


texts:
I might have went about this the wrong way, though I thought it may be a good idea to ask before starting over again. I have a picture column in cartodb which has multiple picture ids separated by commas. So I created a script that would take those  values and parse the ids out into an a href tag. My question: Is there a way to grab the json value inside  infowindow template {{ photo_column}}  and run the value to my script?

var str = '7BC414B0-B38E-4A70-9464-65CCAEC39C8E,327CF268-6742-4E20-9C92-41E9945FCCBD';
var myarray = str.split(',');
var url = "https://web.fulcrumapp.com/api/v2/photos/";
var token = ".jpg?token=3231c4d7309903b454ac10da556778cdf775b72192a96b3";
    for (var i = 0; i < myarray.length; i++) {
       document.write("

");


texts: ArcGIS Server Geoprocessing Service: Good Async, Bad Sync

texts:
I'm working a geoprocessing service for ArcGIS 10.1+ written in Python / arcpy.

It runs fine in asynchronous mode, but does not write output to disk when switched to synchronous.  I know that synchronous services are supposed to use the arcgisoutput directory, but nothing appears there, or output is deleted before the client can request it.

In forming the output path, I use the arcpy.env.scratchWorkspace variable, which is set to the jobs directory in both async and sync.

What can explain this instability in an ArcGIS Server geoprocessing service when running synchronously?

Thanks in advance


texts: Showing only point markers when labels are shown in QGIS?

texts:
Is there a chance to define a QGIS style rule to just show point markers when their label is shown?

My aim is to exclude all points which are not labeled because of collision between labels.

Example:


texts: Using r and sf, create an concave envelope POLYGON around MULTILINESTRING

texts:
I am interested in creating the minimum possible polygon around a set of nearly parallel cross sections.
For example:

mls <- st_multilinestring(list(rbind(c(2,2),c(1,3)),
                    rbind(c(0,0),c(1,1),c(2,1)),
                    rbind(c(3,4),c(0,4))))
plot(mls)

I know there's the function st_convex_hull but that doesn't accomplish what I want because it doesn't generate the minimum possible polygon.

I came up with almost a solution, like this:

mls %>% st_cast(to = "MULTIPOINT") %>%
st_cast(to = "POLYGON") %>%
st_simplify() %>%
plot()

Unfortunately, this only produces the desired output for this simple case.  But it quickly falls apart with more lines.  Here's example bad output:

mls2 <- st_multilinestring(list(rbind(c(2,2),c(1,3)),
                               rbind(c(0,0),c(1,1),c(2,1)),
                               rbind(c(3,4),c(0,4)),
                               rbind(c(1,10),c(-1,12))))
plot(mls2)
mls2 %>% st_cast(to = "MULTIPOINT") %>%
  st_cast(to = "POLYGON") %>%
  st_simplify() %>%
  plot()

One idea I had was to pair adjacent lines and generate a polygon for each pair and then merge those polygons, but that feels cumbersome.  Any other solutions or clean coding for my first idea would be appreciated.  If sf doesn't have the necessary methods, I'm open to other packages.


texts: How to put this in a for-loop? I'm recreating ArcMap's 'Tabulate Area' tool in Python

texts:
I have two shapefiles, a province with 39 municipalities and a land cover of the province.
Goal: get the sqm of each land cover classification for each municipality. For example, how much open forest, closed forest, eroded areas, etc. are there in Muni 1.
I managed to write a code that does this but only for one of the 39 municipalities.
gdf1 = gpd.read_file('munis.shp')
gdf2 = gpd.read_file('landcover.shp')

muni1 = gdf1.loc[gdf1['NAME_2'] == 'muni1'] # NAME_2 is field containing all munis
intersection = gpd.overlay(muni1, gdf2, how='intersection')
area = intersection['geometry'].map(lambda p: p.area)
print(area)

The result is the sqm of each land cover classification in muni1. How to insert this in a for-loop so I can do the process for each 39 municipalities? I tried this but not working:
for i in range(gdf1.shape[0]):
    munis = gdf1.loc[gdf1['NAME_2'][i]]
    intersection = gpd.overlay(munis, gdf2, how='intersection')
    area = intersection['geometry'].map(lambda p: p.area)
    print(area)

Result
KeyError: 'Abad' # first municipality if sorted alphabetically


texts: Avoiding loss of text (labels, legend, title, page numbers) when exporting many Map Series pages from ArcGIS Pro?

texts:
This is the first page of an ArcGIS Pro layout, with Map Series enabled, that has been exported to PDF when choosing to do just pages 1-10:

This is the first page of the same ArcGIS Pro layout, with Map Series enabled, that has been exported to PDF when choosing to do all 255 pages:

As you can see, no text is coming through in the labels, legend, static title, or dynamic page numbers, and the same loss of text is seen when I export pages 1-100.

I am using ArcGIS Pro 1.2, on Windows 7 SP1, 64bit, with 8Gb RAM.  The data being used is from Natural Earth 1:10M Cultural Vectors: Admin 0 – Countries and Admin 1 – States, Provinces.

It seems like the Map Series export is unable to handle this modest amount of text, for more than a very modest number of pages.  Has anyone encountered, and then been able to avoid losing text in similar circumstances?

Both layers added from these shapefiles have Page Queries in place to label only the province names from the current country, and to label all countries except the current country, but I do not suspect these to be the cause of the symptoms that I am observing.  When I removed these two page queries as a test the problem persisted.

One thing I am suspicious of is that I am currently running with only 1-2Gb disk space spare, so I will try freeing up some more to see if that has any effect.

The same project (which has nothing more in it than what is needed to do the above task) has given three Serious Application Errors that also seem to be related to exporting its layout.  I plan to recreate the project from scratch to see if that too has the "text dropping problem" and whether it gets any Serious Application Errors.

I now have more than 4Gb of free disk space and I have repeated the steps to recreate an identical project from scratch.  No further Serious Application Errors have been experienced.  I exported the layout to PDF a number of times.  The text stayed for the first few of 255 map pages (the only ones I checked) as I added more complexity to the map (labelling countries, dynamic title, legend, dynamic page numbers) but when I added a second layer for states and provinces, and labelled that, the exported PDF was missing the text on the first three pages BUT appears to have it on the remaining 252 pages!  I did not use any page queries in this test.

I've now reproduced this using far less steps and complexity than in my original test case:

Download ne_10m_admin_0_countries.shp and ne_10m_admin_1_states_provinces.shp from Natural Earth
Create a new Blank project
INSERT a New Basemap
In Project pane Add Folder Connection to the Natural Earth data
Drag and drop ne_10m_admin_0_countries.shp into the Basemap
Label the layer using default NAME field
Use FEATURE LAYER | LABELING to access Layer Placement group and Label Class pane where Conflict Resolution tab can be used to Remove duplicate labels > Remove all
Drag and drop ne_10m_admin_1_states_provinces.shp into the Basemap
Label the layer using default name field
Use FEATURE LAYER | LABELING to access Layer Placement group and Label Class pane where Conflict Resolution tab can be used to Remove duplicate labels > Remove all
INSERT a New Layout - I used ISO - Portrait > A3
INSERT a Map Frame for the Basemap
Use LAYOUT > Map Series to enable a Spatial map series for the Basemap Map Frame, using layer ne_10m_admin_0_countries, and field NAME for both the Name and Sort fields
Use SHARE > Export Layout and do all 255 pages

When I examine the output PDF there are no labels on the first three pages, but they appear to be on the remaining 252 pages.


texts:
For my master’s thesis I’m using ArcMap to create maps based on a large historical dataset. I’m mapping individual movements over long periods of time.
I’m struggling to design the database and figure out how I’m going to actually make the maps.
Here, I’ve made a basic, scaled-down example of my project. I want to map which counties the individuals are living in, in 1880 and in 1910.

I have a spreadsheet of individual data:

It gives the county that each individual is living in in each year. There is no spatial data in it, just the names of counties and county ID

Using a file geodatabe, I made a feature classes of county boundaries as polygons:

Because the boundaries of counties change over the decades, I made one for 1880 and one for 1910. The counties each have a unique id (ID_NUM). I can use this unique ID to join or relate the two tables together.

My question is where do I go from here?

I want to create two maps. Where the individuals are living 1880 and where they are living in 1910.

Do I join or relate the individual table to both county tables?

Once the tables are connected somehow, how do I use the symbology of the county maps to show where the individuals are living?

How do I give spatial data to the individuals based on the spatial data from the county tables?

I don’t know exactly how to ask for what I need. I’m just at the start of this project and I’m feeling overwhelmed in how to design it.

I am looking to find out the number of individuals living in each county at each time period. Sort of like a population density map


texts: How can I interact with PyGrass after a script is launched?

texts:
From Grass GUI 7.4.1, I launch a .py script (File > launch a script) that builds maps.

How can I interact with the maps viewer, selecting the coordinates of some features, giving them as a new input to continue the process? The simple command x=input() is not operating.


texts: SLD to CSS. Geoserver CSS module problem

texts:
I styled Geoserver maps with SLD. And now I have to convert this style to CSS. CSS is easier to support. Well, I have this fragment of SLD:

<Geometry>
    <ogc:Function name="offset">
        <ogc:PropertyName>way</ogc:PropertyName>
        <ogc:Literal>-13</ogc:Literal>
        <ogc:Literal>13</ogc:Literal>
    </ogc:Function>
</Geometry>

I don't know how to use 'offset' function in CSS. Can you help me?


texts: Most Effective Format to Manage Aerial Photography for Viewing Only

texts:
What is the most effective format to manage aerial photography for viewing in ArcGIS?

I manage several datasets of Aerial Photography at a few different companies. I have backups of the original photography, so I'm not worried about keeping the original format in any way. What I need to do is create a small dataset, as too not take up much space on the server, that also draws as quickly as possible in ArcGIS.

So far I had been using GRID's in a File GeoDatabase, but that was just because I has assumed it was the standard. I then switched to TIFF's which were smaller, but seemed to draw slower. Before I do another mass conversion, I'm hoping someone can offer guidance on this. I've heard some great things about JP2's and that's where I am likely heading next.

Also I don't want to reduce cell size, but other than that I have no restrictions in any answer.


texts: Transform data stored in a TIFF file to a CSV with latitude/longitude and data?

texts:
I have obtained a population density map from the European commission's Global Human Settlement Layer. The download contained a .tif file and .tif.ovr file sharing the same name. Their specifications say about the file "Values are expressed as decimals (Float)
 representing the absolute number of inhabitants of the cell."

I need to convert this data to a format where it is listed, in plain text, the latitude, longitude, and that float value from earlier for every cell (csv would probably work best).

So what would be the best way to go about this conversion?

Sorry if this is rather basic or vague as I'm new to this field.


texts:
I'm trying to close the window (by code) that appears right after the "add feature" tool is called. This is the window I'm talking about:

This is my plugin code:
    def run(self):
    """Run method that performs all the real work"""
    self.iface.actionAddFeature().trigger()
    # CLOSE WINDOW CODE HERE

Everything I tried closed my QGIS program. Is what I'm asking possible?
If so, how?


texts: GDAL3 gdal2tiles AttributeError: 'SpatialReference' object has no attribute 'this'

texts:
Windows10 64bit.
I did a clean Python & GDAL install using the instructions on this site, but then for the latest 3.1.2 GDAL version and Python 3.7.9.
I'm trying to run a most basic gdal2tiles.py from CMD on a georeferenced (EPSG 31370) pdf:
gdal2tiles.py myfile.pdf

I got the error code below. I also tried running from the python IDLE with os.system, which of course returns the same error. I also tried with a georeferenced tiff without succes. I noticed the same error was posted recently on the OSGeo github issues page.
ERROR 1: PROJ: proj_create_from_database: Cannot find proj.db
Traceback (most recent call last):
  File "C:\Users\myuser\AppData\Local\Programs\Python\Python37\lib\site-packages\osgeo\osr.py", line 933, in __init__
    self.this.append(this)
AttributeError: 'SpatialReference' object has no attribute 'this'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "C:\Users\myuser\AppData\Local\Programs\Python\Python37\Scripts\gdal2tiles.py", line 2928, in <module>
    main()
  File "C:\Users\myuser\AppData\Local\Programs\Python\Python37\Scripts\gdal2tiles.py", line 2922, in main
    single_threaded_tiling(input_file, output_folder, options)
  File "C:\Users\myuser\AppData\Local\Programs\Python\Python37\Scripts\gdal2tiles.py", line 2849, in single_threaded_tiling
    conf, tile_details = worker_tile_details(input_file, output_folder, options)
  File "C:\Users\myuser\AppData\Local\Programs\Python\Python37\Scripts\gdal2tiles.py", line 2765, in worker_tile_details
    gdal2tiles.open_input()
  File "C:\Users\myuser\AppData\Local\Programs\Python\Python37\Scripts\gdal2tiles.py", line 1535, in open_input
    self.in_srs, self.in_srs_wkt = setup_input_srs(input_dataset, self.options)
  File "C:\Users\myuser\AppData\Local\Programs\Python\Python37\Scripts\gdal2tiles.py", line 737, in setup_input_srs
    input_srs = osr.SpatialReference()
  File "C:\Users\myuser\AppData\Local\Programs\Python\Python37\lib\site-packages\osgeo\osr.py", line 934, in __init__
    except __builtin__.Exception:
NameError: name '__builtin__' is not defined

Any ideas on how to fix this?
/EDIT/
I took a look at python37/lib/osgeo/site-packages/osr.py
It says version 2.0.4, whereas on github the latest version is 3.0.12
I tried replacing the source code with this of version 3.0.12
This tiles are created, but throws the error below. The created tiles have the wrong numbering, so the spatialReference is still faulty.
ERROR 1: PROJ: proj_create_from_database: Cannot find proj.db


texts: How to speed up osm2po when using as java lib?

texts:
I'am trying to utilize osm2po as java lib for calculation of distances matrix between some set of points (lat/lon) and for each pair of points I call the following method that code is based on sample provided at osm2po.de web-site:

public int getDistance(GeoPoint source, GeoPoint target) {
        int res = Integer.MAX_VALUE;
        int sourceId = graph.findClosestVertexId(source.getLatitude(), source.getLongitude(), 1);
        int targetId = graph.findClosestVertexId(target.getLatitude(), target.getLongitude(), 1);

        router.traverse(graph, sourceId, targetId, Float.MAX_VALUE, params);

        if (router.isVisited(targetId)) { // Found!
            int[] path = router.makePath(targetId);
            float distKm = 0.0f;

            for (int segmentId : path) {
                RoutingResultSegment rrs = graph.lookupSegment(segmentId);
                distKm = distKm + rrs.getKm();
            }

            res = (int)(distKm * 1000);
        }

        router.reset();
        return res;

But I noticed that it takes about 2-3 seconds per one point to calculate whole distance matrix (I have 10-25 points per matrix in average) that looks not too fast (30-60secs per matrix).
Could someone advice what could be improved here - especially I am not sure in correct usage of reset() call - when actually it should be done? - it's not too much documentation on osm2po usage within java...
Also I would appreciate any performance tips to improve this code. One more guess I have is that findClosestVertexId() is expensive enough call and cashing its results for subsequent calls could improve the situation.


texts: Iterating polygons in folder as input in QGIS Modeler like ArcGIS ModelBuilder

texts:
How do I iterate polygons in a folder as input in QGIS modeller?

Usually, ArcGIS can iterate polygons as an input for functions like union, merge, intersect, etc.

But I cannot find the way to do this in QGIS modeller, not a pyqgis (for convenience).


texts: Conversion from .shp to .dxf in one unified layer

texts:
I am trying to make a 3D .dxf out of a 2D shapefile for several hundred buildings. Doing this in one fell swoop is proving to be a roadblock.

In ArcScene, I extruded building heights for maybe two dozen city blocks in Tel Aviv then later applied the Layer 3D to Feature Class and Multipatch to Collada tools. The immediate goal was to take this swath of hundreds of cubes of various shapes and sizes representing buildings and import them into Sketchup and/or 3DS Max as a .dae for later conversion to a .dxf.

Unfortunately, the buildings can only be imported individually, one at a time, without spatial alignment to one another in both Sketchup and 3DS Max. Precision is key, so needless to say this is unacceptable. I need one cohesive layer.

The eventual goal is to import these buildings into the gaming software Unity for unwrapping and rendering with images as part of both car accident and parking spot simulations.


texts:
I am using QGIS AND struggling with setting projections. I am attempting to reproject all of my files that I am accumulating into NAD83-UTMzone13 utilizing the reprojection tool.  The tool appears to run correctly, I save the file. I set the project projection under the properties. When I add the shape, 1. the project reverts back to WGS84 and 2. the newly projected shape still does not show up, even when I reset the project back to the desired projection.
Thoughts?


texts:
I'm trying to identify this grid system that consists of two letter followed by two numbers.

If it helps I know GU37 is defined by the following vertices:
  Lat (North)   Lng (West)
1 48° 30' 00"   64° 00' 00"
2 48° 40' 00"   64° 00' 00"
3 48° 40' 00"   63° 50' 00"
4 48° 30' 00"   63° 50' 00"
5 48° 30' 00"   64° 00' 00"

Point references are based on the geodesic system North American Datum 1983 (NAD83)
My ultimate goal is to be able to convert a point using PostGIS to a grid name.
Given that I know an origin - GU37 - and I know this is a 10 minute grid (I think) - is this something I could generate?


texts: Point sampling tool returning no value for for extractions in QGIS

texts:
I am attempting to use the point sampling plugin in QGIS to extract values from phenology layers for the Western US provided by USGS. The data comes a .hdr file  with a prj file specifying a Lambert Azimuthal Equal-area projection.

PROJCS["Lambert Azimuthal Equal area",
GEOGCS["GCS_Sphere_ARC_INFO",DATUM["D_Sphere_ARC_INFO",SPHEROID["Sphere_ARC_INFO",6370997.0,0.0]],
PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],
PROJECTION["Lambert_Azimuthal_Equal_Area"],
PARAMETER["False_Easting",0.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",-100.0],
PARAMETER["Latitude_Of_Origin",45.0],
UNIT["Meter",1.0]]

Points for my sites were generated from a csv with lat long coords in decimal degree. I have tried placing the points into the same crs as the raster, I have tried modifying the false northing of the raster to align the points and the raster, I have tried reprojecting both to a third projection and modifying the false northing of the raster to make them align.

When I use a categorical raster with the same crs args as the point file, the extractions work as expected.

The result with the phenology raster extraction is empty cells in datatable.

Any ideas?

(ignore the fact that the phenology data is likely highly correlated with the landcover data - boss wants it)


texts: OpenLayers 3: does not display all Features of Layer

texts:
I'm using openlayers 3 and geoserver. I want to get and display a wfs layer with openlayers.

this is source of my layer:

vectorSource = new ol.source.Vector({
            format: new ol.format.GeoJSON(),
            url: function (extent) {
                return 'http://myUrl/geoserver/wfs?service=WFS&' +
                    'version=1.1.0&request=GetFeature&typename=myLayer&' +
                    'outputFormat=application/json&srsname=EPSG:4326&' +
                    'bbox=' + extent.join(',') + ',EPSG:4326';
            },
            strategy: ol.loadingstrategy.tile(ol.tilegrid.createXYZ({
                maxZoom: 19
            }))
        });

and style for my vector layer:

new ol.style.Style({
                stroke: new ol.style.Stroke({
                    color: 'rgba(255, 0, 0, 1.0)',
                    width: 2
                }),
                geometry: function (feature) {
                    return feature.getGeometry();
                }

all features are type of MultiLineString and loaded correctly in client. But Map does not display some of features.

wms from above source correctly works.

my geojson and geometries are correct.

please help.


texts: How to make dashed lines between adjacent polygons render correctly in QGIS?

texts:
I have adjacent polygons that I want to be marked out with dashed or dotted lines. This sometimes works but there appears to be a problem where the pattern isn't always synchronised on the shared boundaries which leads to odd patterns or solid lines. Is there a simple way round this?

Using GQIS 2.10.1 Pisa


texts:
I'm trying to deal with an ECW file with QuantumGis 1.8 Lisboa, on Ubuntu 11.10. Obviously it tells me that .ecw is not a supported kind of file on QuantumGis. I have tried to install all the libraries but it still doesn't work.

What sould I do?

Thanks,

Jacopo


texts: How to lock map labels for inset map in Print Composer?

texts:
OK, so I am making a print layout which has two maps, one an inset map with just a skeleton of a line layer and the other is the detail map at higher resolution. In the detail map, I want to show the labels for each arm of the skeleton. For the inset map, I want the labels to be off.

Even when I 'lock layers for map item', the labels still update in the inset map when the map is exported to PDF.

The update preview is set to Cache for both, but I don't think this has an impact of what gets exported.