texts: Dividing points into Polygons based on an Attribute of that point and its relative distance to other points and a facility

texts:
I have random point data (school locations) that I need to divide into 'n' polygons (regions) based on: their distance to a facility, their relative distance to each other, and an attribute of that demand point (school impedance value).

This is really a simple Location-Allocation problem with the 'schools' as demand points, and 'Instructors' as facility points.  However, the built in ArcGIS location-Allocation impedance value is only based on either TIME or Distance.  I want to base the impedance value on distance (or time) AS WELL AS on a "school impedance value" (a simple integer field of the Schools attribute table).

The purpose:
There are 'n' number of Instructors (facilities) that have to visit 'n' number of schools (demand points). I need to divide up these schools into regions (polygons). However, some schools take multiple days for the instructor to complete their task at that school (larger schools require more time to complete the task).

For example:
there are 300 schools and 3 instructors. Location Allocation tool would divide these up according to either: equal number of schools (100 schools/instructor) or lowest impedance value based on time/distance.  But of those 300 schools 50 of them could take multiple days to teach and theoretically all 50 multi-day schools could be assigned to 1 instructor; thus making the teaching time uneven.

I do have some programming/scripting background if it is required.


texts: KML/KMZ file <10MB for global shipping lanes?

texts:
I am trying to render an overlay of the world's shipping lanes using google KmlLayer class from the Javascript V3 API but most of the files I find are not rendered because of google's 10MB size limit.

Does anyone know of any helpful KMZ/KML files of a suitable size for world shipping lanes?


texts: postgis finding lines that touch and update

texts:
I am not much of a postgis guy so please bear with me. I have a subset of lines representing exits and roads around those exists. Main roads and exists can be distinguished by their class attribute (exits = class 5).

My final goal is to update each single line that together make up an exit with the class attribute of the main road with the smallest value (say you leave the highway (class 1) and continue to a local road (class 4), the exit lines need to be updated to class 1).

Sounds pretty straight forward but I'm struggling.

I tried to break the problem apart so I came up with the following queries to close in on the problem.

# finds all lines of exits that together form the exit road
SELECT ST_LineMerge (ST_Collect(a.geom))
FROM sr_buffer as a, sr_buffer as b
WHERE ST_Touches(a.geom, b.geom)
AND a.class = 5
AND b.class = 5

This isn't really helping in the end I'm afraid as I'm losing the ids along the way.

# finds all transition points from main road to slip road
CREATE TABLE touch AS
SELECT s2.geom, s2.linkid, s2.class, s1.linkid AS linkid_sr FROM
sr_buffer as s1, sr_buffer as s2
WHERE s1.class = 5 AND s2.class != 5 AND ST_TOUCHES(s2.geom, s1.geom)

This is better already. Now we know that the exists starting with ID 14 and 9 for example should become class 1.

  Class    ID

  1     14

  1 9

Now it should be possible to find all lines that make up the exit starting at ID 14.

Any ideas on how to proceed?


texts: Geometry Generator - Line pattern fill to a buffer of x metres

texts:
I have an image below to help illustrate the style I'm trying to apply to my polygon. I have an idea of what I need to type in the geometry generator but can't get it to work.
Something along the line of:
difference(buffer($geometry,-2)

How do I make the hatching always appear at 45 degrees?

I've tried using a simple line marker but end up with hatching appearing outside the polygon or crossing over each other at times; see below.


texts: Get GDAL up and running with ArcSDE enabled

texts:
I want to migrate rasters from ArcSDE 9.2 to Postgresql.

I have access to a ArcSDE server with several large raster datasets (mosaics) in ArcSDE 9.2. I can use ArcCatalog, ArcMap and ArcToolbox to view and manage the datasets. If I use the Export Raster tool, it fails, since the generated tif is larger that 4Gb.

Unfortunately, I'm not able to compile GDAL with SDE support. I have access to the SDE dlls installed on the server, but I don't have the ArcSDE Development Kit to compile GDAL with SDE support. I was able to compile GDAL, but not with SDE support.

I would like to split the raster mosaics. How can I split the mosaic in several smaller rasters?


texts: Draw streets with bridges using "feature rendering order" instead of "symbol levels" (QGIS3)

texts:
To control rendering orders of symbols, you could use the "symbol levels"-feature to give every geometry a own Z-dimension. If you want it dynamically because of complex situation like bridges over bridges you need more control. The "feature rendering order" enables you to give every feature a own Z-dimension using a field or expression. So far so good, but what if I use more then one symbol level per symbol like streets with borders? The screenshot below shows the situation:

The order of the "main-symbol" (color of the road) is correct, but the borders are all on the same level. That is not the expected result. How could I solve that without using the static "symbol level"-feature?

Added:
Because I sometimes have to work on several scales in one project, my symbols and symbol-levels become confusing:


texts: Creation of conditional statement with Python fails

texts:
I want to calculate new values with the QGIS field calculator, while running the code in the calculator GUI from QGIS it works quite well, but trying to apply this via PyQGIS it just give me in expression_2cl the value after the else statement "Hardwood".

What is wrong with my syntax that QGIS doesn't recognize the full statement?

#Create new fields
pv = layer.dataProvider()
pv.addAttributes([QgsField('Share_SW',QVariant.Double),\
QgsField('2_classes', QVariant.String)])

layer.updateFields()

expression_SW = QgsExpression('("basalareap"+"basalareas")/("basalareap"+"basalareas"+"basalaread")')
expression_2cl = QgsExpression('CASE WHEN \"Share_SW\" >= 0.5 THEN \'Softwood\' ELSE \'Hardwood\' END')

context = QgsExpressionContext()
context.appendScopes(QgsExpressionContextUtils.globalProjectLayerScopes(layer))

System Windows 10 QGIS 3.10


texts: I want to write a script in Python for ArcGIS to search for the size of an activeDataFrame (DataFrame) in meters (scale)?


texts: What are some Domain Specific Languages (DSLs) for GIS?

texts: I'm curious what sort of domain specific languages (DSLs) are used in GIS.  The only one I can think of is the DOCELL language that was provided in ArcInfo Workstation GRID, but is no longer supported.


texts: How to reduce wild swing in Cesium camera.flyTo?

texts:
Cesium's camera flyTo produces a pretty wild swing.  Try it out by going here:

http://cesiumjs.org/Cesium/Apps/Sandcastle/index.html?src=Cesium%20Widget.html&label=Showcases

and paste in this code:

var viewer = new Cesium.Viewer('cesiumContainer');
//Cartesian3.fromDegrees(longitude, latitude, height, ellipsoid, result)

//new york 40.7127° N, 74.0059° W
viewer.camera.flyTo({
    destination: Cesium.Cartesian3.fromDegrees(-74, 41, 500000),
    duration: 0.5
});

setTimeout(function() {
    //Tokyo 35 N, 139 E
    viewer.camera.flyTo({
        destination: Cesium.Cartesian3.fromDegrees(139, 35, 500000)
    });
}, 4000);

We start out looking at New York.  We are trying to fly to Tokyo.  In the middle though we swing way out to look into space.  I feel like the camera should stay pointing to the center of the earth throughout.  Instead it seems like it swings really quickly to look to the direction it will end up looking, then flies to the right place in space to look at the point it was set to fly to.  This seems to happen any time you fly from places that are far enough apart to require a significant camera angle change.

Is there a way to keep it looking at the center of earth throughout, and to not fly quite so far into space?


texts: plpgsql | before update trigger: insert the new pseudo-record in another record and leave the old one as it was

texts:
I've got the following table:

CREATE TABLE alumbrado_puntos_luz
(
  id numeric(8,0),
  "GEOMETRY" geometry,
  revision_actual numeric(10,0),
  revision_expirada numeric(10,0),
  CONSTRAINT alumbrado_puntos_luz_id_key UNIQUE (id)
)
WITH (
  OIDS=TRUE
);

When updating the table i would like to update just the revision_expirada value with the next serial and insert the NEW pseudo-record's values in a new record on the same table. I mean, every time i have a change in the drawing (GEOMETRY) i will create a new record with the changes and leave the old one as it is.

I don't know if a trigger is the best way of doing this or creating a rule would be a better idea.

Regards,


texts: How to install the QGIS Raster Lang plugin?

texts:
I´m trying to install the rasterlang plugin for QGIS. I tried several things on this website but it seems to be quite complicated because I need to install binary files and other files.

The idea would be to make it the way that when I type in my console:
from rasterlang.layers import layerAsArray then It would work nicely. Thanks


texts: ArcGIS Online map freezing at certain locations

texts:
When using the map view of ArcGIS Online I have encountered a number of locations which, when zoomed in to, will cause the web page to become sluggish, freeze, and typically be unresponsive for 30+ seconds. When the view is panned the web page will lock up again.

Limited testing has shown the issue occurs with at least the default 'Topographic' basemap with and without additional data loaded, starting at zoom levels around 14 or greater. Also it appears to affect Google Chrome but not Internet Explorer.

Example coordinates where this behavior has been noticed (long,lat to be pasted into the map view search box):

3.121, 53.869 - link to location at a problem zoom level (17)

3.070, 53.246 - link to location at a problem zoom level (17)

Does anybody have insight into what may be causing these lockups or any solutions?


texts: Understanding unknown spatial reference warning in ArcMap?

texts:
I am updating some layers to arcmap and receing a warning saying "unknown spatial reference". My understanding is that in this case arcmap will try to adopt the data frame coordinate system (in the present case defined by the first layer  uploaded).

The layer seems to be where it is supposed to be, but can someone confirm my thoughts?


texts: Why do different methods give different numbers of pixels in an image in Google Earth Engine?

texts:
I'm trying to figure out why I get different answers when calculating the total number of pixels in an image in different ways.

Calculate a histogram on all possible values, then add together all of the categories to get the total number of pixels. This won't count NA's.
Use the ee.Reducer.count() function. This also doesn't count NA's, so I would expect this to give the same answer as number 1.
Use unmask().ee.Reducer.count() function. This DOES count NA's, so it will be larger than number 2 if there are NA's.
print the image, look at the dimensions, and multiply the 2 dimensions to get the number of pixels. This should count NA's and I would expect it to give the same answer as number 3.

Why do these methods give different counts of the total number of pixels in an image?

Sample code: https://code.earthengine.google.com/2c4767f13001604254fa7d431af14b86


texts: Elevation profile extraction to excel using ArcGIS from DEM in a tiff format

texts: I have an SRTM image in .tiff format downloaded from USGS website and I want to extract elevation along a path to a Excel file in the form of distance vs elevation or lat,long. vs elevation. Can anybody tell me how can I do this in ArcGIS?


texts: Feature rendering issue in Geoserver WMS layer at larger zoom scale

texts:
I am using Geoserver 2.15.1 version to display PostGIS data as WMS layer at client side. There will be a rendering problem in my WMS layer, some dense areas are not visible in larger scale as marked in the below snapshot. Those area are getting generated when zooming closer to the particular location.

Is there any way in GeoServer to make the features visible at all scale levels? I have tried MinScaleDenominator and MaxScaleDenominator, but no luck.


texts: 1-center and k-center problem / maxmin facility problem in Post

texts:
is there anybody who has already solved the 1- and k-center problem within PostGIS (optionally with distance metric obtained by pgrouting)?

References:

1-center problem in Wikipedia and graph center in Wikipedia
Facility location problem in Wikipedia

What I did so far for the 1-center problem (I found a similar solution elsewhere) :

// Generate min. cost table for all pairs
INSERT into catchments
SELECT
    seq, id1 as source, id2 as target, cost
FROM
  pgr_apspWarshall('
    SELECT id, source, target, cost
    FROM network WHERE city = ''New York''', false, false);

// Select max(min(distance))
SELECT source, max(cost) as maxcost
FROM catchments
WHERE cost <= 10 **km** GROUP BY source ORDER BY maxcost ASC;

However, the second step misses the initial choice.

EDIT: Someone has already asked a question about OSS and location allocation which is a nearby topic. However, in the tools listed (Grass v.net.alloc)the functionality divides an area into districts based on predefined hubs/facilities/stations. To emphasize, I current question is about ab-initio optimization without previous definitions of location.

EDIT2: the Floyd-Warshall algorithm returns only the path information which has not previously been calculated- that is called dynamic programming in this sense.
Well, now the question remains how to reestablish the full information to run the 1-center problem?


texts: Delete features from FeatureCollection in GeoJSON?

texts:
I have a FeatureCollection in GeoJSON. I want to delete all features that have an ID that doesn't start with "B".

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "id": 0,
      "properties": {
        "map": "eu",
        "name": "Belgium" // don't delete
      },
      "geometry": {
        "type": "Polygon",
        "coordinates": ...
      }
    },
    {
      "type": "Feature",
      "id": 2,
      "properties": {
        "map": "eu",
        "name": "CzechRepublic" // do delete this Feature
      },
      "geometry": {
        "type": "Polygon",
        "coordinates": [ ...

I can do this manually in Python, but is there a nicer way to do it, using something on the command line?


texts: Problem merging multiple paths of KML files (saved as SHP) in QGIS

texts:
I have several KML files (trail01.kml, trail02.kml etc.) each containing multiple paths (about 30) created in Google Earth. I have opened them in QGIS and saved them as SHP (trail01.shp, trail02.shp etc.), merged the paths so that each SHP contains one path (one feature, as in one line in attribute table instead of 30), and then saved them as KML files and opened them in Google Earth again.

The problem is that when I try to show the elevation profile of a (now) single path, the topology isn't right, as if the previous 30 paths weren't merged in the right order (the profile travels from the beginning of one path to the end of it, and then jumps to another start point, which isn't the next path but seemingly a random one).

I've sorted the paths in Google Earth so that they are in right order in the KML file (if that makes a difference?), and they show up in the right order in the attribute table in the SHP file. I've also tried multistage merging, meaning merging only two paths at a time, but the same problem occurs.

I have very basic skills in Python which I'm willing to try, I just don't have a good idea where to start?


texts: Using update cursor to populate a field based on a Date field in ArcPy

texts:
I'm trying to populate a field (Season, "TEXT") in my attribute table with a season code based on a sampling date (FIXDATE, "DATE"). For example sample collected between 01/01/2016 and 01/03/2016 would get season code "SPRING". I'm new to arcpy and search cursors. Here is my code so far.

import arcpy
from datetime import datetime

caribou = "C:/Python/Test_data.gdb/test_data"
fields = ['FIXDATE', 'Season']
with arcpy.da.UpdateCursor(caribou, fields) as cursor:
    for row in cursor:
        if (row[0] >= datetime.strptime("01/01/2016", "%d/%m/%y") and row[0] <= datetime.strptime("01/03/2016", "%d/%m/%y")):
            row[1] = "SPRING"
        elif (row[0] >= datetime.strptime("02/03/2016", "%d/%m/%y")and row[0] <= datetime.strptime("01/06/2016", "%d/%m/%y")):
            row[1] = "SUMMER"
        elif (row[0] >= datetime.strptime("02/06/2016", "%d/%m/%y") and row[0] <= datetime.strptime("01/10/2016", "%d/%m/%y")):
            row[1] = "FALL"
        elif (row[0] >= datetime.strptime("02/10/2016", "%d/%m/%y")):
            row[1] = "WINTER"

    # Update the cursor with the updated list
        cursor.updateRow(row)


texts: Keeping track of edited location in ArcGIS for Desktop?

texts:
With large datasets, how do you keep track of where you have edited across the map? For instance, I am editing a city's waterline and water valve data. I estimate using a 1:400 average zoom level. This is a large city, and I can quickly lose track of where I have edited, especially after days of working on it.

I'm using ArcMap 10.2 Basic with a file geodatabase. Is there a function to track the location and extent of where my window has been? Or is there a better way to track changes/edited location?

If not, I might write script to make an initial polygon, or draw a feature, based on the current window extents, and then merge new extents to it as I pan across the map. Perhaps stop it from drawing if zoom is beyond 1:1000.

What do you think?


texts: Flipped (Y-axis) Raster from point cloud

texts:
I have been working on a project where I need to create rasters from irregularly spaced point clouds.  I am resampling my point clouds to a regular 0.25 cm grid using a nearest neighbor method pyresample.  I am using the output array to make a raster using GDAL and my raster is inverted on the y-axis.  For a particular dataset I get the extents and apply a geotransfrom using the following parameters:

top_left_x(i.e., xMin)
top_left_y(i.e., yMax)
x_resolution
y_resolution
x_rotation
y_rotation

For a hypothetical dataset, let's say my geo transform parameters are:

top_left_x(i.e., xMin)=223280.0
top_left_y(i.e., yMax)=577877.75
x_resolution=0.25
y_resolution=-0.25
x_rotation=0
y_rotation=0

Using GDAL's GTiff driver I am successfully creating rasters by writing the output array Result to rasterband(1).  The problem is the output raster us upside down w.r.t the y axis.  If I flip the array result using numpy like: np.flipud(result), the raster is plotted in the correct position.  I am at a loss why I need to flip result and want to see if someone can help me understand why.


texts: QGIS Forms - How Do I Display a Read-Only Value?

texts: Older versions of QGIS had an "Immutable" option for form fields to make it obvious that data was read-only.  How do I display values in the current version and make it obvious that they are read-only?  (Useful for showing auto IDs etc.)


texts: Is there an open-source equivalent of ArcGIS "Zonal Statistics as Table" tool?

texts:
With ArcGIS zonal statistics as a table I can easily calculate the mean value from a raster contained within a polygon from a shapefile. I now have a Mac, and therefore I cannot access ArcGIS anymore. Actually, I am trying to use open source software only (e.g. R and QGIS).

I thought that "zonal statistics" from QGIS would do the trick, however, an essential feature for my needs is missing: in ArcGIS, you can specify the column by which you want to aggregate the values. For example, if in your shapefile you have 100 polygons, you may not want to just calculate the mean value from the raster in each polygon. Rather, you may want to know the mean value of a group of polygons. In my case, one column ("veg_class") classifies all polygons as "tropical", "temperate" or "boreal" vegetation.

How can I, in QGIS or any other software, select the column to aggregate the statistics so, in the end, I obtain one mean per veg_class (i.e. 3) instead of per polygon?


texts: Implementing radio buttons in QGIS processing scripts

texts:
Is it possible to implement radio buttons in QGIS 3.x processing scripts?

See also: http://docs.qgis.org/2.8/de/docs/user_manual/processing/scripts.html (supported parameter types for QGIS 2.8)


texts: Arcpy ListFields field types different from Layer Properties

texts:
I'm trying to add some functionality to an ArcPy toolbox to check the column names and types of an input feature class.
It seems that the types that Python uses are different from the types in ArcMap itself?
For example, a column of type float in the layer properties menu in python is type Single
for the other column types:
Text: String,
Long:Integer,
Short:SmallInteger

I haven't found a reference for the full universe of options only this documentation and there's no mention of this distinction in the docs for ListFields, it just uses the Python types, with no mention of the difference between Python and ArcMap.
is there a better way to be working with field types than by name like this? Is there documentation on this that I haven't found?


texts: For street network data lacking z-elevation information

texts:
Can anyone help me in this? Many street network datasets lack z elevation fields needed to model overpass and tunnels (the F_ZLEV and T_ZLEV fields). Although my road network spatial database contains bridges and tunnels as features, I need to create the two fields from it. Is there a quick fix to automate this. I am thinking of something like:

Select all road links and allocate F_ZLEV=0 and T_ZLEV=0
Select all bridges which at their link start connect to road link and allocate F_ZLEV=0 and T_ZLEV=1
Select all bridges which at their link end connect to road link and allocate F_ZLEV=1 and T_ZLEV=0
Allocate all other bridge links F_ZLEV=1 and T_ZLEV=1.

Any ideas will help.


texts: Calculate distance using 2 sets of lat/long values in the same record

texts:
I have a point layer with over 90K records in it. For each record I have lat/long field for incident, and lat/long for engine location. Is there a way to find the distance between Incident point and Engine Location point?

So far, I created points using "Make XY Event layer"  using lat/long for incidents. And again I run the same tool for Engine locations. When I run the Near tool in ArcGIS it will find the nearest engine point around the incident and will give me the distance to the nearest engine but not necessarily the distance from the engine location on the same record. How can I calculate distances using the lat/long values for incident and engine location in the same record?


texts:
I am developing a QGIS plugin. I need to load some polygons from a table in PostGIS. This table is defined with this SQL query :

CREATE TABLE cad.polygon(
    id_polygon serial NOT NULL,
    geometry geometry NOT NULL,
    id_plan integer NOT NULL,
    manually boolean NOT NULL,
    creation_date date NOT NULL,

    CONSTRAINT polygon_pkey PRIMARY KEY(id_polygon),
    CONSTRAINT fk           FOREIGN KEY(id_plan) REFERENCES cad.plan(id_plan)
) WITH (OIDS = FALSE) TABLESPACE pg_default;

and contains 163 entities with two geometry types (curve polygon and polygon).
I wrote a Python script to get these polygons in QGIS :

uri = QgsDataSourceUri()

# Connect to the database on the local host
uri.setConnection(
    "localhost",
    "5432",
    "dbname",
    "username",
    "*****"
)

# SELECT * FROM cad.polygon WHERE id_plan = <plan ID>
uri.setDataSource(
    "cad",
    "polygon",
    "geometry",    # Column name for the geometry
    "id_plan = 2"  # whatever the plan ID
)

uri.setKeyColumn("id_polygon")
uri.setWkbType(QgsWkbTypes.Polygon)

# Create a new QGIS vector layer
layer = QgsVectorLayer(uri.uri(False), "polygons", "postgres")
QgsProject.instance().addMapLayer(layer)

When I run this script, near the layer generated in the legend, there is this message :

  Unusable layer : The data source of this layer has not been found. Click to choose a new data source.

I don't understand where this error comes from and how to fix it.


texts: How can I merge postcode sector polygons into larger territories in QGIS?

texts:
I'm new to GIS so patience might be required :-)

I'm building a basic web-based mapping tool for business use, to be used as a demo to show that we need one but managers don't know it yet!

So far I've built a HTML page displaying an OpenLayers OSM base layer. I've been supplied with some postcode sector polygons in shapefile format that I'm wanting to merge into our business territories (wards, districts etc.).

Aims were to plot these as layers on top of OSM that the user can switch between etc.

Using QGIS, i don't seem to be able to import or link to a table of our geography hierarchy and dissolve the polygons by these fields. Is this best achieved by importing CSV or linking to db etc. I'm a bit lost :-)

our geo hierarchy:

Disrict, Ward, Neighbourhood, PC Sector
11, 11A, 11A15, AA1 1 etc.

Any advice would be gratefully received :-)


texts:
I have been trying to open a map which I downloaded from the ESA website.

There are several files in the folder including .cpd, .csv, .dsr, .qml, .tif and so on. I was trying to open the .tif image in ArcGIS Desktop but it would not display anything. Then I was just trying to identify but I realized something was not right here. Apparently .qml files could be opened in QGIS, but I have no idea how to work around QGIS.

Is it possible to open the QML file in ArcGIS Desktop?


texts: See the result of a query on the map with Arcgis Api JS

texts:
Let's say you enter these parameters. Which method of the "ArcGISImageServiceLayer" class can serve to display the image of the mosaic dataset?

In the attached image the result is an image, but I do not know how to visualize it on the map.

Part of my code is this, I used the "queryVisibleRasters" method to try to visualize it, but it did not work.


texts: Custom Grid axis labels, to only display decimal minutes ignoring degrees

texts:
I have created a chart with a very small grid, I would like the axis labels to only display the decimal minute value, so I have more space to increase font size, by removing the degree value.
I presume I need to use the custom format, but I can't work out the coding.
Picture of my current axis labels at the moment.
other example of how I would like to display the asix labels


texts:
In Arcmap 10.8, I am attempting to draw an orthophoto and several polygons at a high scale (not lower than 1:6000). However, Arcmap simply refuses to draw any layers unless I have a spatial resolution as low as app. 1:8000 or lower.
What could the reason be? I have tried the reverse many time (that you need a quite high scale to display certain layers), but never this.


texts: It seems the Processing options in QGIS 3.0 no longer include R Scripts (they are not listed as available along with SAGA, GRASS, and GDAL). Is this a change in this version of QGIS or have I had a problem with installation?


texts:
Thats my ways.shp file. I loaded this to PostgreSQL DB using PostGIS shape file loader. I ran pgr_anlayzeGraph() on it and found it had a lot of isolated segments.

Next I ran pgr_nodeNetwork() to create a new table 'ways_ways'.
I made the topology  for it and ran the following query in Qgis sql window.

 SELECT * FROM ways_ways JOIN (SELECT seq, id1 AS node, id2 AS edge, cost FROM
pgr_astar('SELECT id, source, target, cost, x1, y1, x2, y2 FROM ways_ways',
1, 108, false, false) )AS route ON ways_ways.id = route.node;

Given below are the screenshots:
Original floorplan and the prescribed route

query ouput

Clearly , I don't have a continuous output. I ran pgr_analyseGraph on ways_ways to find that it has no more isolated segments . It has 51 deadends, which is acceptable.

Did I miss any step in pgrouting?
Is there anyway to see which node has id 1 and id 108 in qgis?
i ran the same query in another map(which i downloaded) and got output for all possible source and target ids. So how can i create a way(a route) that is completely noded?


texts: Cannot install a version of GDAL via Anaconda that permits reading BigTIFFs

texts:
I have some BigTiff files that I need to read into Python using gdal to do some analyses. According to GDAL's documentation:

  When built with internal libtiff or with libtiff >= 4.0, GDAL also supports reading and writing BigTIFF files (evolution of the TIFF format to support files larger than 4 GB).

My version of libtiff is 4.0.9, so there should be no problem. However, if I try

bt = gdal.Open('path_to_bigtiff.tif'),

bt is a NoneType. I tested to ensure that non-BigTIFF files are readable with gdal.Open(), which they are.

Furthermore, if I check the list of gdal's supported drivers' metadata like this:

md = gdal.GetDriverByName('GTiff').GetMetadata()
if md['DMD_CREATIONOPTIONLIST'].find('BigTIFF') == -1:
    return 'BigTIFF is not supported'

BigTIFFs are not supported. My version of gdal is 2.2.2.

Using Anaconda, I created a new environment with

conda create --name test1 gdal

but I get the same result--gdal 2.2.2 is installed along with libtiff 4.0.9 and BigTiff files still cannot be read.

Using Anaconda again, I created a new environment with

conda create --name test2 gdal=2.3.1 -c conda-forge

but when I run import gdal in the test2 environment, I get the error message

  ImportError: DLL load failed: The specified module could not be found.

How can I install a version of gdal that will permit me to load BigTiffs? I would also like it to be compatible with rasterio and fiona if possible.

Here is the output from conda info for the environment containing 2.2.2:
(geoenv) C:\Users\Jon>conda info

     active environment : geoenv
    active env location : C:\Users\Jon\Anaconda3\envs\geoenv
            shell level : 2
       user config file : C:\Users\Jon\.condarc
 populated config files : C:\Users\Jon\.condarc
          conda version : 4.5.9
    conda-build version : 3.0.27
         python version : 3.6.3.final.0
       base environment : C:\Users\Jon\Anaconda3  (writable)
           channel URLs : https://repo.anaconda.com/pkgs/main/win-64
                          https://repo.anaconda.com/pkgs/main/noarch
                          https://repo.anaconda.com/pkgs/free/win-64
                          https://repo.anaconda.com/pkgs/free/noarch
                          https://repo.anaconda.com/pkgs/r/win-64
                          https://repo.anaconda.com/pkgs/r/noarch
                          https://repo.anaconda.com/pkgs/pro/win-64
                          https://repo.anaconda.com/pkgs/pro/noarch
                          https://repo.anaconda.com/pkgs/msys2/win-64
                          https://repo.anaconda.com/pkgs/msys2/noarch
          package cache : C:\Users\Jon\Anaconda3\pkgs
                          C:\Users\Jon\AppData\Local\conda\conda\pkgs
       envs directories : C:\Users\Jon\Anaconda3\envs
                          C:\Users\Jon\AppData\Local\conda\conda\envs
                          C:\Users\Jon\.conda\envs
               platform : win-64
             user-agent : conda/4.5.9 requests/2.18.4 CPython/3.6.3 Windows/10 Windows/10.0.15063
          administrator : False
             netrc file : None
           offline mode : False


texts: How to change attribute value in drop down list with constraint and/or Function Editor in QGIS?

texts:
I'm usig QGIS 2.18.6.

UPDATE2:
I've mangaged to get (assign/find) the actual widget (value Map) for my attribute (field) "ibruk" with calls from the python console in QGIS:

>>>layer = qgis.utils.iface.activeLayer()
>>>index = layer.fieldNameIndex("ibruk")
>>>layer.editFormConfig().widgetType(index)
u'ValueMap'
>>>layer.editFormConfig().widgetConfig(index)
{u'': u'', u'I bruk': u'1', u'Ikke tatt i bruk (under bygging)': u'5', u'(opplagt) ute av bruk': u'3', u'Usikker': u'4'}

But, how do I now set a specific value (like setValue, changeValue, setText or similar) so this value automatically appear in the drop down list and the user may just click OK? (See further down for more explanation.)

--

UPDATE1:
I'm trying to change the value in a Value map widget ("ibruk") in QGIS automatically using python-code in the Function Editor.

I have tried something similar described in Layer with custom UI not working properly, but I'm not able to change the value/text in the combo-box. I want "Usikker" (= value 4) to appear as chosen value in the drop box.

from qgis.core import *
from qgis.gui import *
from qgis.utils import iface
from PyQt4.QtCore import *
from PyQt4.QtGui import *

@qgsfunction(args='auto', group='Custom')
def ingen_adkomst(feature, parent):
    layer = qgis.utils.iface.activeLayer()
    index = layer.fieldNameIndex("ibruk")
    QgsEditorWidgetWrapper.fromWidget( ibruk ).setValue( 4 )
    """return 4"""

I'm having trouble to assign the Value Map (called "ibruk"), then change the value. I have found examples of how to create a QComboBox with PyQt and use setText etc. But in my case, I already have the Value Map (set in Layer Properties in QGIS). I just want a specific option to be chosen based on a constraint, see below. Is it possible?

--

I have one attribute "o1mrk" (Value Map) where the user may choose one of 5 options. If user choose 5 ("o1mrk" = 5), then I want another attribute's ("ibruk") value automatically to be set to 'Usikker' (value 4) saving the user some clicking and thereby some time.

So if "o1mrk"= 5 then "ibruk"='Usikker'.

First I tried this constraint for "ibruk":

CASE
  WHEN "o1mrk"=5 THEN "ibruk"=4
  ELSE "ibruk" IN (1,3,5)
END

But this requires the user to manually choose "ibruk" like 'Usikker' when "o1mrk" is set to 5 - and no less clicking. I want 'Usikker' to appear automatically when "o1mrk" = 5 (i.e. Ingen adkomst).

I have managed to partly achieve what I want by setting a constraint for "ibruk" and using Function Editor (which I'm not mastering, just trying my best by looking at questions and example-code here at gis.stackexchange.com and QGIS Expression):

CASE
  WHEN "o1mrk"=5 THEN ingen_adkomst()
  ELSE "ibruk" IN (1,3,5)
END

where ingen_adkomst() is my own Python function:

from qgis.core import *
from qgis.gui import *
from qgis.utils import iface
from PyQt4.QtCore import QVariant

@qgsfunction(args='auto', group='Custom')
def ingen_adkomst(feature, parent):
    layer = qgis.utils.iface.activeLayer()
    index = layer.fieldNameIndex("ibruk")
    f_id=feature.id()
    layer.changeAttributeValue(f_id,index,"Usikker")

The attribute value in the attribute table changes. However, the visible Value Map-option (value) is still not changed - preventing me from clicking OK.

Originally before changing "o1mrk":

After choosing "o1mrk" = 5 (Ingen adkomst). The value of "ibruk" is appearantly not changed (but the value is changed in the attribute table) and turns red:

Even if I manually choose 'Usikker' from the drop down list, the option stays red:

How may I change my constraint and/or code so that when "o1mrk" = 5 (Ingen adkomst), the value of "ibruk" changes to 'Usikker" also in the drop down list so the red colour disappear and I'm allowed to click OK?


texts: Inserting row into table and then updating it in sum of each column using ArcPy?

texts:
My problem is the following:

I have a table with an unknown number of columns (n)
Use arcpy.ListFields create a list of columns and generates Field name list
Then, using arcpy.da.InsertCursor create a new row at the end of the table ...

for example

  fieldNameList = []
  fields = arcpy.ListFields(fc, "*")

  for field in fields:
     if field.type in ("Double", "Integer", "Single"):
        fieldNameList.append(field.name)

 cur = arcpy.da.InsertCursor(fc, fieldNameList)
 for x in xrange(-1, 0):
     cur.insertRow((SUM_FIELD1, SUM_FIELD2...SUM_FIELD n ))       ?????????!!!!
 del cur

I do not know how to calculate the sum for each column and then the result to update in the created row. The sum should be separately calculated for each column ...


texts: Choosing coordinate system for mapping Mars in ArcGIS Desktop?

texts:
I wish to use ArcGIS Desktop to map planets, in particular Mars.

Where would I start in terms of basemapping and projection?


texts: Is it possible to draw a point interactively and execute an intersection operation (drill down) to identify spatial relationship with other layers?

texts: Is it possible to draw a point interactively and execute an intersection operation (drill down) to identify spatial relationship with other layers using GeoServer? I am a newbie with GeoServer…


texts: Use expression to only keep 1 layer in the legend

texts: When I print my project, that includes different layers, I always end up with all the layers in the legend. What I want when printing is to only show 1 layer, that always has the same name "Werken". The rest of the layers are not necessary to be shown in the legend. What kind of expression would I have to write to become this result?


texts: Applying 'Integrate' tool (data management) using ArcPy?

texts:
I try to integrate two features (shapes) with their ranks using ArcPy.

arcpy.Integrate_management()

I don't get a clue from the explanation/example in ArcGIS help. It's straight forward, however, to apply this tool from the graphical dialog.

I seek a working example having 2 shapes: a.shp (rank 1), b.shp (rank 2).

I'm using ArcGIS 10.2


texts: Creating a cached image server on ArcGIS with EPSG:4326 projection

texts:
I need to create a cached image service for development reasons which doesn't use the ArcGIS Online or Google Maps projection (EPSG:3857). It should use the EPSG:4326 projection.
I couldn't find anywhere on ArcMap to configure it. The Esri documentation suggests to use their tool to export the required XML from an existing cached map service but I don't have one.
It's important to note that I did succeed in publishing an EPSG:4326 image service to my enterprise server but when I attempt to create a cache, it wants me to add tiling schemes manually. Even when I do that, it doesn't cache.
What am I doing wrong?
Bottom Line: How do I publish a cached Image Service on ArcGIS with the EPSG:4326 projection?
P.S. On-the-fly transformation to EPSG:4326 from EPSG:3857 could also work for me. Anyone knows if that's possible? Some query string perhaps?


texts: ArcGIS Pro Label Placement grayed out in Scene

texts: I'm using Scene in  ArcGIS Pro I don't understand why my Label Placement is grayed out even though I've selected my Feature Class in the Contents panel and the Maplex Label Engine is checked? I'm new to Pro and could be missing something obvious.


texts:
I would like to calculate the number of cells from gridded elevation data that drain into a particular cell. First,  I calculate the flow direction by selecting the cell with the steepest slope from one out of eight adjacent cells.

library(raster)
library(gdistance) # package for adjacent function

# make the blank raster, the extents and resolution of which are up to you.
# the cells are numbered 1 to ncell (starting in the top left in R)
r <- raster(xmn=0,xmx=6,ymn=0,ymx=6,res=c(2,2))
r[] <- 1:ncell(r)

# flow direction of digital elevation model
flow <- terrain(subelevation, opt="flowdir", neighbors=8)
# transform to matrix
dim(flow@data@values) <- c(sqrt(length(flow@data@values)), sqrt(length(flow@data@values)))
# yields values from 1 to 128, see Figure 1
flow@data@values

Figure 1

Figure 2 plots flow@data@values transformed to an igraph object to illustrate the problem. For any chosen cell, I would like to calculate the number of incoming cells connected at any distance. For example if the cell labelled A is chosen, it should have the value 3 as B, C and D drain into it. The more efficient the solution, the better as the gridded dataset is large.

Figure 2


texts:
Does QGIS data work in Esri CityEngine?

I have a copy of CityEngine Basic and imported a shapefile which have the street network and lots. Without ArcGIS I can't open data which have building footprints and DEM from USGS/CLICK. Would QGIS open these files so I can edit them to be imported to ESRI CityEngine?

Importing Terrain and Imagery CityEngine Essential Skills #5 is the tutorial that I'm working on which instructed me to use USGS/CLICK to get the data. In this tutorial downtown Redlands CA was the area of interest an orthoimagery and elevation was emailed. The data in the email was extracted using a program named Bitser to a folder then preview in ArcMap. Once in ArcMap it was imported to CityEngine. Is there another way to get this done without ArcMap to use the data in CityEngine?


texts: Extra lines in Sketchup DTM when using stamp tool

texts:
I have a DTM in Sketchup. I am trying to use the stamp tool to make a level area on which to build a building.

However, when I use the stamp tool lots of extra lines appear on the DTM seemingly at random (I realise these are the lines that form the DTM but only some of them appear).

The lines are not  hidden geometry. It may be that this is the reason they are showing - they aren't hidden and the rest of the lines are. If this is the case I'm not sure how to hide them again without exploding the DTM, which I don't want to do.

Can anybody tell me why this is happening? Can anyone tell me how to prevent it?

Example:

Left: Original DTM

Right: DTM after using stamp tool.

Note: The images are of different parts of the DTM! My problem isn't that the DTM is changed (the stamp tool works exactly as it should), it is just that the lines become visible.


texts: Selective masking in QGIS for point markers

texts: According to the changelog of QGIS 3.12 the selective masking will not only work with labels, it should also be available for point features. However, I can't find the function for point features, only labels. Any ideas?


texts: Why is ArcMap leaving connections " in transaction"?

texts:
The TL;DR version: When you use ArcMap Query Layers to access a PostgreSQL table, ArcMap's database connections are almost always sitting idle in unclosed transactions.  Why would this be?

I have PostgreSQL 9.1.x (with PostGIS 2.0.x) with some existing tables and data.  I successfully use ArcMap (10.1 SP1) to connect to one of those tables (using a Query Layer) and render dots on my map.  The ArcGIS enterprise geodatabase tables are not installed.

One day I brought up pgAdmin's server status window while ArcMap was running (and nothing else was accessing the database).  As expected, there are two connections whose application name is "ArcSDE(sv=nnnn;sn={xxx})".  But I was surprised to see that both of those connections are listed as "<IDLE> in transaction" in pgAdmin's "Query" column.

That seemed fishy.  As a rule of thumb, it's generally bad to leave lingering open transactions.  And as a programmer, I've witnessed first-hand the dangers of improper failure handling that leaves unclosed transactions.

So I closed everything down, restarted PostgreSQL, then started ArcMap again.  Same result: two connections, both of which are "<IDLE> in transaction".

I cranked up the PostgreSQL statement logging and watched as ArcMap did its magic.  Sure enough, the very last statement in the log (once the ArcMap window finishes rendering) is BEGIN, with no corresponding COMMIT or ROLLBACK.  And then every time you pan or zoom, the very first thing it does is COMMIT.  Then it gets all the data and ends its rendering with an un-closed BEGIN again.

At first glance, it appears benign.  Those two connections hold no locks other than the locks on their own virtual transaction IDs.  But it sure seems odd.

Any ideas what legitimate reason ArcMap might have for behaving this way?


texts: Possibility to use polygon as global area parameter for overpass query

texts:
I have been trying to provide a polygon with an array of coordinates as a global parameter for my overpass query. But the query is through syntax error.

[bbox:poly:"50.7 7.1 50.7 7.2 50.75 7.15"];
node[amenity=cafe];
out;

Error: line 1: parse error: ',' or ']' expected - ':' found. Error: line 1: static error: A bounding box needs four comma-separated values. Error: line 1: static error: A bounding box needs four comma-separated values. Error: line 1: static error: A bounding box needs four comma-separated values. Error: line 1: parse error: ';' expected - ']' found.

Providing a bounding box instead works as expected.

[bbox:50.7,7.1,50.75,7.25];
node[amenity=cafe];
out;

It also works when the polygon is given in-line,

node[amenity=cafe](poly:"50.7 7.1 50.7 7.2 50.75 7.15");
out;

(as suggested in this tutorial)

As final attemp when below query was tried no error was given but also the result was blank.

area["50.7 7.1 50.7 7.2 50.75 7.15"];
node[amenity=cafe](area);
out;

Is it possible to give a polygon instead of bounding box as a globle parameter in an Overpass query?

I could not find any documentation over the Overpass wiki pages that supports the thing I am trying here.


texts: Plot the spatial join shapefile boundary in Python

texts:
I have two spatial object

Geo_TWS is a shapefile boundary, geo_cell is points and already convert to GeoPandas.

base = Geo_TWS.plot(color='white', edgecolor='black')
 geo_cell.plot(ax = base,color='yellow' )

My problem is: how can I plot the yellow points with the intersects black boundary?

So the final figure shows the yellow shape, and the black boundary line within the shape.

(I have tried gpd.sjoin but seems can't complete it, anyone have better ideas? or proper way to use gpd.sjoin)

Attempt:

import geopandas as gpd

TWS = gpd.read_file('TWP bounadry/TownBoundary2.shp')
TWS = TWS.to_crs({'init': 'epsg:4326'})
Geo_TWS = TWS

#cell file
import pandas as pd

grid_cell = pd.read_csv('ARB_2km_Info.txt', delim_whitespace=True)

from geopandas import GeoDataFrame
from shapely.geometry import Point

#map to point
geometry = [Point(xy) for xy in zip(grid_cell.Lon, grid_cell.Lat)]
new = grid_cell.drop(['Lon', 'Lat'], axis=1)
crs = {'init': 'epsg:4326'}
geo_cell = GeoDataFrame(new, crs=crs, geometry=geometry)

Twship_in_cell = gpd.sjoin(Geo_TWS, geo_cell,  how="left", op='intersects')

Twship_in_cell.plot(color = 'white', edgecolor = 'black')

I have tried many how and op combinations, but still not working, and thats the problem..


texts: How to create .jgw file out of jpg picture with exif header?

texts:
If i have picture with this metadata in exif header:

Size is 4000, 3000
Coordinate System is `'
Metadata:
  EXIF_ApertureValue=(3.625)
  EXIF_ColorSpace=1
  EXIF_ComponentsConfiguration=0x1 0x2 0x3 00
  EXIF_CompressedBitsPerPixel=(5)
  EXIF_CustomRendered=0
  EXIF_DateTime=2015:09:09 15:27:28
  EXIF_DateTimeDigitized=2015:09:09 15:27:28
  EXIF_DateTimeOriginal=2015:09:09 15:27:28
  EXIF_DigitalZoomRatio=(1)
  EXIF_ExifVersion=0230
  EXIF_ExposureBiasValue=(0)
  EXIF_ExposureMode=0
  EXIF_ExposureTime=(0.0005)
  EXIF_FileSource=0x3
  EXIF_Flash=16
  EXIF_FlashpixVersion=0100
  EXIF_FNumber=(3.5)
  EXIF_FocalLength=(4.5)
  EXIF_FocalPlaneResolutionUnit=2
  EXIF_FocalPlaneXResolution=(16393.4)
  EXIF_FocalPlaneYResolution=(16393.4)
  EXIF_GPSAltitude=(91.8)
  EXIF_GPSAltitudeRef=00
  EXIF_GPSDateStamp=2015:09:09
  EXIF_GPSLatitude=(45) (18) (3.564)
  EXIF_GPSLatitudeRef=N
  EXIF_GPSLongitude=(19) (48) (19.842)
  EXIF_GPSLongitudeRef=E
  EXIF_GPSMapDatum=WGS-84
  EXIF_GPSStatus=A
  EXIF_GPSTimeStamp=(14) (27) (29)
  EXIF_GPSVersionID=0x2 0x3 00 00
  EXIF_ImageDescription=
  EXIF_Interoperability_Index=R98
  EXIF_Interoperability_Version=0x30 0x31 0x30 0x30
  EXIF_ISOSpeedRatings=800
  EXIF_Make=Canon
  EXIF_MakerNote=
  EXIF_MaxApertureValue=(3.625)
  EXIF_MeteringMode=5
  EXIF_Model=Canon PowerShot SX280 HS
  EXIF_Orientation=1
  EXIF_PixelXDimension=4000
  EXIF_PixelYDimension=3000
  EXIF_Related_Image_Length=3000
  EXIF_Related_Image_Width=4000
  EXIF_ResolutionUnit=2
  EXIF_SceneCaptureType=0
  EXIF_SensingMethod=2
  EXIF_ShutterSpeedValue=(10.9688)
  EXIF_UserComment=
  EXIF_WhiteBalance=1
  EXIF_XResolution=(180)
  EXIF_YCbCrPositioning=2
  EXIF_YResolution=(180)
Image Structure Metadata:
  COMPRESSION=JPEG
  INTERLEAVE=PIXEL
  SOURCE_COLOR_SPACE=YCbCr
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 3000.0)
Upper Right ( 4000.0,    0.0)
Lower Right ( 4000.0, 3000.0)
Center      ( 2000.0, 1500.0)
Band 1 Block=4000x1 Type=Byte, ColorInterp=Red
  Overviews: 2000x1500, 1000x750, 500x375
  Image Structure Metadata:
    COMPRESSION=JPEG
Band 2 Block=4000x1 Type=Byte, ColorInterp=Green
  Overviews: 2000x1500, 1000x750, 500x375
  Image Structure Metadata:
    COMPRESSION=JPEG
Band 3 Block=4000x1 Type=Byte, ColorInterp=Blue
  Overviews: 2000x1500, 1000x750, 500x375
  Image Structure Metadata:
    COMPRESSION=JPEG

How should i format appropriate .jgw file for this picture?


texts: Changing title automatically according to selected layer using QGIS composer

texts:
I'm working on a project where I have to create different maps for the same area. I have lots of maps and each of them is related to a different theme, that I've grouped into different groups of layers as you can see in the image.
Each group of layers has a name that will be used as the map title. Because I have a lot of them, each time that a select a different group I want to change the map title automatically to have the same name as the layer group that is selected.

How I can do that?


texts: Can I use .aux, .sid, .sid.xml, .sid.aux.xml, or .sdw files as raster in QGIS?

texts:
I am working with data from my boss who uses ArcGIS. I am aware that .aux files accompany raster files, but to my knowledge you can still use .aux files to create rasters in ArcGIS.

I'm working with limited data. Is there a way to use any of these files to create a visible raster in QGIS? I'm getting a visual when I upload the .aux file to QGIS, but it's only a black image. If I fiddle with the style it will disappear and turn into a white image.

Also, my companion was able to make the .sid file work in her QGIS, but she uses a PC and I use a MAC. Would there be any explanation for this?


texts: How to configure custom scale in Leaflet?

texts:
I have a custom WMS tileset generated using ESRI ArcGIS Server.

I have successfully configured an OpenLayers page to display this basemap.  The code is as follows:

var map = new OpenLayers.Map('map', {
    projection : "EPSG:3857",
    resolutions : [ 132.29193125052919, 92.604351875370426, 67.733468800270941, 33.866734400135471, 16.933367200067735, 8.4666836000338677, 4.2333418000169338, 2.1166709000084669, 1.0583354500042335, 0.52916772500211673, 0.26458386250105836 ],
    maxExtent : new OpenLayers.Bounds(-8352920.48698, 4887220.912586, -8082558.309632, 5067712.988343)
});

var basemap = new OpenLayers.Layer.WMS("basic_wm", "http://myhostname/geowebcache/service/wms/", {
    layers : 'basic_wm',
    format : 'image/jpeg'
    }, {
    tileOrigin : new OpenLayers.LonLat(-20037700, 30241100)
});
map.addLayer(basemap);

I am trying to reproduce the equivalent in Leaflet, but have so far been unsuccessful.  I see that I can override the scale function and set a max zoom, but I am not certain how to override the scale function.  I'm also wondering if I shouldn't be able to provide my DPI, tile origin, max extent, etc?

Here is what I have so far in Leaflet:

var map = new L.Map('map', {
     center: new L.LatLng(40.693726386087945, -73.9856945509795),
     zoom: 1
});

var basemap = new L.TileLayer.WMS("http://myhostname/geowebcache/service/wms/",
    {
        layers: 'basic_wm',
        attribution: "Copyright 2012 The City of New York",
        maxZoom: 10
});
map.addLayer(basemap);


texts: Using MAX in SQL Expression of ArcMap ModelBuilder

texts:
I am trying to use a SQL expression in order to select values that are >= 25% of its max value in that column. For instance, VALUE >= MAX(VALUE)*0.25, and those will ultimately be the selected points it pulls out.

Is there syntax to make that happen, or am I going to need to create a new column and use field calculator on it?

Biggest issue with that is I will have to do this for about 250 files. I am also trying to create this in ModelBuilder.

When I try it as above it just returns that there is an error with the expression.


texts: Viewing Esri ASCII Raster file in ArcGIS Desktop?

texts:
i have some problem viewing an ascii file in a structured way in arcGIS.

i have used IDL ENVI too to solve this issue but all data is not displayed.

This link has the files:
https://drive.google.com/?authuser=0#folders/0B8iAKVklcaMUVU5qY2NNVFBRbDQ


texts:
I'm looking for GIS files (shapefiles, geoJSON, KML...) for US special districts, but I can't find them anywhere. Special districts are single-purpose local government distinct from counties and municipalities. For instance, a water management, or a mosquito eradication district.

Some counties have files, but only for districts inside their boundaries. I'm looking for something for the US as a whole.


texts: QGIS shuts down when executing PyQGIS script

texts:
I'm new to PyQGIS. I have some simple PyQGIS scripts.

The only one that I use a bool input doesn't work. Whenever I execute it, QGIS shuts down.

I see that the reason is 'qgis killed by SIGABRT' but I have no idea what it is. Could you please take a look if there're some mistakes in my code or it's about the QGIS/Python package? I'm on QGIS 3.10 sur fedora

class algTest(QgsProcessingAlgorithm):
    OUTPUT_TABLES = 'OUTPUT_TABLES'
    CITY = 'CITY'
    INPUT_LINE_CITY= 'INPUT_LINE_CITY'

    def __init__(self):
        super().__init__()

    def name(self):
        return "download_line_bool"

    def displayName(self):
        return "download lines with city bool"

    def createInstance(self):
        return type(self)()

    def initAlgorithm(self, config=None):
        self.addParameter(QgsProcessingParameterString(
            self.CITY, "Input City"))

        self.addParameter(
            QgsProcessingParameterBoolean (
                self.INPUT_LINE_CITY,
                'Only City'
            )
        )

        self.addParameter(
            QgsProcessingParameterVectorDestination (
                self.OUTPUT_TABLES,
                'Output layer'
            )
        )

    def processAlgorithm(self, parameters, context, feedback):
        my_city = self.parameterAsString(parameters, self.CITY, context)
        my_bool = self.parameterAsBool(parameters, self.INPUT_LINE_CITY, context)
        uri = QgsDataSourceUri()
        uri.setConnection("localhost", "5432", "***", "***", "***")

        if my_bool == True:
            uri.setDataSource("public", "transport_public", "the_geom", "city='"+my_city+"'","id")
        else:
            uri.setDataSource("public", "transport_public", "the_geom", "","id")
        table = QgsVectorLayer(uri.uri(), "table", "postgres")
        return {self.OUTPUT_TABLES:table}


texts: Tilecache + Mapnik + Apache generates incomplete tiles (missing layers)

texts: I have TileCache + Mapnik working with a custom osm.xml stylesheet. When I run tilecache_http_server.py the tiles are served OK. But when I run TileCache under Apache (with mod_wsgi) I see incomplete tiles (with missing layers). What can go wrong? Furthermore, generating tiles is very slow. Can this have something to do with the missing layers? Thanks.


texts: Optimized interpolation method with fixed-location sample data points?

texts:
I'm trying to find an optimized spatial interpolation method that deals with fixed-location sample data points. Say if I have fixed-location sampling sites and collect data daily, and then use these data to create interpolation surfaces.

Is there a method that pre-process the spatial relationships among those sampling sites and then integrate this spatial relationship into the interpolation method to speed up the whole interpolation process?

I searched some information but couldn't find what I want. Can anyone put me on the right direction?


texts: Get name/number of the table currently being edited

texts:
In MapInfo you can edit one table at a time.

How can I programmatically get the name or ID number of the table that is currently being edited?

I tried the following function:

Function GetEditableTable as String
  Dim i, j, WinID as integer
  Dim TableList() as String
  i = NumTables()
  ReDim TableList(i)

  WinID = FrontWindow()
  For j = 1 to i
    Note LayerInfo(WinID, j, LAYER_INFO_EDITABLE)
    If LayerInfo(WinID, j, LAYER_INFO_EDITABLE) = TRUE Then
      GetEditableTable = TableInfo(j, TAB_INFO_NAME)
      Exit Sub
    End If
  Next
End Function

...but it returns the same table name no matter which table is editable.


texts: Debugging ERROR 999999 in Raster to Polygon Scripting?

texts:
I am trying to make a script which selects every .png file in a folder beginning with the letters "LG". I then want the scipt create a shapefile, replacing the "LG" with "SH", and then i want the script to buffer that shapefile and rename the buffer with the first 2 letters being "SB".

I keep getting an error 99999 error message at line 37!
( gp.RasterToPolygon_conversion(INPUT_RASTER, Output_polygon_features, "SIMPLIFY", "VALUE") )

I have posted this on stackoverflow.com but probably here is a better place as scripting for arcmap is probably quite specialised!

Can anyone offer me any advice to why this script doesnt work? I am running ArcMap 9.3.1. and (attempting) to run the script from pyscripter!

# ---------------------------------------------------------------------------
# RASTER2POLYGON.py
# Created on: Wed Dec 01 2010 12:20:48 PM
#   (generated by ArcGIS/ModelBuilder)
# Usage: RASTER2POLYGON <INPUT_RASTER> <Output_polygon_features> <Raster_Buffer_shp>
# ---------------------------------------------------------------------------

# Import system modules
import sys, string, os, arcgisscripting

# Create the Geoprocessor object
gp = arcgisscripting.create()

# Load required toolboxes...
gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Conversion Tools.tbx")
gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Analysis Tools.tbx")

# Script arguments...

folder = "D:\\J04-0083\\IMAGEFILES"
for root, dirs, filenames in os.walk(folder): # returms root, dirs, and files
    for filename in filenames:
        filename_split = os.path.splitext(filename) # filename and extensionname (extension in [1])
        filename_zero = filename_split[0]
    try:
        first_2_letters = filename_zero[0] + filename_zero[1]
    except:
        first_2_letters = "XX"

    if first_2_letters == "LG":
        Output_polygon_features = "D:\\J04-0083\\ShapeFiles.gdb\\" + "SH_" + filename + ".shp"

        # Process: Raster to Polygon...

        INPUT_RASTER = os.path.join(root + "\\" + filename_zero + ".png")
        gp.RasterToPolygon_conversion(INPUT_RASTER, Output_polygon_features, "SIMPLIFY", "VALUE")

        Distance__value_or_field_ = "5 Meters"
        Raster_Buffer_shp = "SB_" + filename + ".shp"

        # Process: Buffer...
        gp.Buffer_analysis(Output_polygon_features, Raster_Buffer_shp, Distance__value_or_field_, "FULL", "ROUND", "NONE", "")


texts: Creating a Loop in R for averaging Raster values to create time series

texts:
I have calculated several different indices for lakes stored as GEOTiff. I now want to write a loop with R to average the values for each lake, so only one value is left for each lake. And then create a table with the global means/medians and the date which is given in the filenames.

Input is raster data with following naming convention:

ndci__LC08__192023__20130515.TIFF

ndci__LC08__192023__20130803.TIFF
.
.
.
I uploaded five of my rasters here:

https://www.dropbox.com/s/fax1a16tvg0iohd/ndci_LC08_192023_201305151.tif?dl=0
https://www.dropbox.com/s/ghre1290pj0mwnl/ndci_LC08_192023_201404161.tif?dl=0
https://www.dropbox.com/s/cu1x3kxnd8koz2d/ndci_LC08_192023_201503181.tif?dl=0
https://www.dropbox.com/s/s2uqsj6f14b5u06/ndci_LC08_192023_201508091.tif?dl=0
https://www.dropbox.com/s/044b5vll1rsrd8a/ndci_LC08_192023_2013080321.tif?dl=0

So the output format I am looking for is:

Date       Mean/Median of Index Satellite
20130515   5.4                  LC08
20130803   8.7                  LC08
.          .                    .
.          .                    .
.          .                    .


texts: Changing mouse wheel zoom settings in QGIS?

texts:
I'm having a problem with the mouse wheel zoom setting on my macbook trackpad. The zoom occurs far too quickly for any level of precision. I noted this comment in the dev log:

2004-08-14 [gsherman] 0.4.0devel9 Added mouse wheel zoom. Moving wheel forward zooms in by a factor of 2.

How can I change the settings so that zooming occurs by a factor of 1 or .5?

I'm using QGIS 2.6 Brighton.


texts: gdalbuildvrt.exe doesn't recognize existing path

texts:
I used this function before but it suddenly stopped working on some updates. The simple command throws an error.

PS C:\> c:\ProgramData\Anaconda3\Library\bin\gdalbuildvrt.exe 'c:/tmp/index.vrt' 'c:/tmp/*.tif'
0...10...20...30...40...50...60...70...80...90...100 - done.
ERROR 4: `c:/tmp/*.tif' does not exist in the file system, and is not recognized as a supported dataset name.
Warning 1: Can't open c:/tmp/*.tif. Skipping it
PS C:\>

I don't know what's wrong. The paths exist 100% and tif files are in it.


texts: How to add CQL filter into GetFeatureInfo request in OpenLayers

texts:
I have made a WMS layer using Openlayers2 and Geoserver and I display the layeer using a specific CQL filter.

Now I want to use a WMSGetFeatureInfo request in order to get specific data on click.
Although I use the CQL parameters in my WMSGetFeatureInfo request nothing seems to change and I get back all the results.

This is my code:

First I create the layer and added to the map:

        cql = "ftth_pop='"+ area+"' AND simul_desc='"+simulation+"'";
        layerName2 = 'ngn_marco:out_ducts';
        wms =  new OpenLayers.Layer.WMS("wms","http://ngn-marco.enomix.gr/geo_wms/",
                {
                    layers: layerName2,
                    format:  "image/png",
                    transparent: true,
                    version: "1.1.1",
                    CQL_FILTER: cql,
                    tiled: true
                },
                {
                    format: "image/png",
                    displayOutsideMaxExtent: false,
                    projection: new OpenLayers.Projection("EPSG:900913"),
                    visibility: true
                }
        );

       map.addLayer(wms);

Then I make the WMSGetFeatureInfo request:

      infoControls = {
        click: new OpenLayers.Control.WMSGetFeatureInfo({
            url: 'http://ngn-marco.enomix.gr/geo_wms/',
            title: 'Identify features by clicking',
            layer: wms,
            CQL_FILTER: cql,
            queryVisible: true
            })
        };

I add the control and activate it as:

    for (var i in infoControls) {
        infoControls[i].events.register("getfeatureinfo", this, showInfo);
        map.addControl(infoControls[i]);
    }

    infoControls.click.activate();

The showInfo() function is:

function showInfo(evt) {
        if (evt.features && evt.features.length) {
            highlightLayer.destroyFeatures();
            highlightLayer.addFeatures(evt.features);
            highlightLayer.redraw();
        } else {
            //console.log(evt.txt);
            document.getElementById('responseText').innerHTML = evt.text;
    }
    }

The problem is that the CQL filter for the WMSGetFeatureInfo request doesnt work.

EDIT

These are the parameters sent to Geoserver:

BBOX    2637105.178295,4616485.137413,2638971.913845,4617560.033123
FEATURE_COUNT   10
FORMAT  image/png
HEIGHT  900
INFO_FORMAT text/html
LAYERS  ngn_marco:out_ducts
QUERY_LAYERS    ngn_marco:out_ducts
REQUEST GetFeatureInfo
SERVICE WMS
SRS EPSG:900913
STYLES
VERSION 1.1.1
WIDTH   1563
X   715
Y   325

And this is part of the response I get back:

     <title>Geoserver GetFeatureInfo output</title>


texts:
I'm having two main problems:

When hovering over a map with multiple layers, even though I only have one layer editable and visible, I am seeing all the nodes of all layers show up as red dots when hovering over the map canvas, making it impossible to select nodes of the layer I am actually working on. Can this behaviour be turned off? It may be related to the option for editing only the selected layer, but I can't actually find this under Settings | Options and maybe it is no longer there? It seems that I can only edit the selected layer anyway, so why do the dots need to show up at all? I have checked snapping settings and have the new topological editing option turned off.
I am finding editing of polygons extremely difficult and clunky now. It is now not possible to move a node with a simple mouse drag. It seems that it is now necessary to drag a new node from a vertex to the new position, then delete the old node. Selection can be done by keyboard or mouse, but deletion can only be done by hitting delete on the keyboard, so the whole procedure takes several seconds. However, what usually happens is that when attempting to create a new node, several connected vertices do move at once even if no nodes are selected. Other strange things are also happening.


texts: Showing all polygons with same ID with Atlas generation in QGIS

texts:
I want to create PDF-maps showing polygons/areas automated by using the Atlas generation.

I have a SHP-file that contains polygons and every polygon/feature (in the attribute table) has an ID that I can use to filter these polygons and only show polygons with a certain ID. So sometimes there is only one polygon assigned to let's say ID "X", but then there are cases where several polygons/features are assigned to ID "Y".

When I generate the atlas using the rule based symbology ($id = @atlas_featureid) it creates an individual map for each polygon with the ID "Y" instead of one map with all "Y"-polygons.

How can I create (with the Atlas tool) maps that show all polygons with a certain ID (be it one or more polygons)? So the goal is to create for each ID a map that shows all the features/polygons with that specific ID.

I can achieve this manually of course by just applying a filter on the SHP-file layer and select each ID manually, but I would prefer an automated way with the Atlas.

I'm using QGIS 3.6.1.


texts: Add python modules so scriptrunner will find them for QGIS on Mac

texts:
I am using QGIS 3.10, and Mac OS 10.(14 & 15)

This issue relates to recent versions of QGIS for the Mac that are distributed as .dmg.

These versions of QGIS bundle their chosen version of Python inside the application so they have control of the version being used.  One consequence of this is that the path that python uses to find modules (sys.path) is relative to the QGIS application directory.  So when one installs a module on your machine with pip it is installed into the system python tree, not the tree rooted in the QGIS application.

I ran into this because I use scriptrunner to run some scripts that include gpxpy.  I initially fixed this by adding the directory where gpxpy was install to sys.path in the script but this made my script machine dependent.

My first thought was to see if there was a version of pip in the QGIS python tree that I could use to install modules but if it is there I did not find it.

Scriptrunnner appears to be unsupported -- there are 3 open issues with no comments or actions over more than a year.

One thing that puzzles me is how do plugins that require other module work around this issue so there must be a known solution.

Can someone point me in the right direction?


texts: Fire frequency Map in ArcGIS for Desktop?

texts: I would like to create a fire frequency map in ArcGIS showing the areas which have burnt 1, 2, 3 times etc (overlapping areas). Currently I have years from 2000 to 2014 of burnt area MODIS data which I have converted into polygons but I'm really struggling to figure out how to do this


texts: PostGIS : Find direction between custom GPS tracks (linestrings) intersected with pgr_dijkstra route

texts: I have many custom GPS tracks(linestrings) and a route generated by pgr_dijkstra by giving source & destination as input. Using ST_Intersects, I able to filter required GPS tracks belongs to the route generated by pgr_dijkstra. Now I wanted to classify intersected GPS tracks direction like how many GPS tracks are from source to destination and how many belongs from destination to source. What are some possible methods to find direction (A->B / B->A) between 2 linestrings?


texts: Iterating rasters using parse path in complex model in ModelBuilder?

texts:
I have a complex model (with many tools/steps chained - named from "a" to "p") and I need iterate two input rasters (p1.img and p2.img), so that all the model outputs returned for the two input files. As there are many chained steps on the model, I applied the Parse Path output as a in-line variable substituition. My attempts and the problem were:

Multivalue iterator: return the outputs correctly but just in relation the file p1.img, doesn't run for the file p2.img. (image below)
Raster iterator: generate a kind of loop, running repeatedly and cumulatively outputs in relation of the file p1.img. (image below)
Turn the raster iterator and parse path tools as a submodel and add it into a main model that contain the steps:  return the outputs correctly but just in relation the file p1.img, doesn't run for the file p2.img. (image below)

How can I run all the model steps for both input files p1 and p2?

I suspect that the problem is the Parse Path tool, that in some way limit the model to run upon the first file (p1) only.


texts: I think this is a fairly easy question so apologies.  Someone sent me a .gdb in Google Drive.  When I download it to my computer, it is a file folder with all the .gdb contents in it.   I can't figure out how to add it as a .gdb in ArcCatalog - I know it must be straightforward but I'd be very grateful if someone could explain it to me.  Thanks in advance!


texts: Validating Values in Multiple Fields in an Attribute Table

texts:
I've been playing around with 'Maintaining Attribute Integrity' functions found here:

http://resources.arcgis.com/en/help/...000000n000000/

I see how I can set domains and subtypes to control the values in the attribute table to maintain quality and how to validate features and relationships using the tools ArcGIS provides.

I don't know if I am missing something, but I have a situation where values in one field in an attribute table, depend on the value in another field in the same attribute table.

Let's have 3 fields (CITY, PROVINCE, COUNTRY)

Let's say I set a domain on the CITY field to populate a list of all the cities in the world. I want to be able to validate that the PROVINCE and COUNTRY fields are correct. So, for example, if I enter 'Toronto' in the CITY field, I want to validate that the PROVINCE field is 'Ontario' and the COUNTRY field is 'Canada'

I'm not sure how to set up the domains or subtypes to do this. I want to be able to run the validate functionality found here:

http://resources.arcgis.com/en/help/...000000v000000/

To throw an invalid error if these conditions aren't met. Can anyone suggest how I can accomplish this?

I can do this through python, but the validation in ArcGIS does a lot of things I was going to script, so I'm hoping it will do this for me too.

Thanks,
Mike


texts: Polygons created with R are corrupted when imported into ArcMAP

texts:
I created some polygons using R and to get a better visualization and some other stuff I need to import those into ArcMap. However, when I imported the polygons I realized there are some holes in some of the polygons,not all. Then I thought there could be a problem with the calculation but when I checked R, there's no problem with the polygons. I'm completely lost about what's wrong. This is my R polygons:
and this is what I got from ArcMap:

Can you help me understand what's wrong with them?


texts: Black vertical lines on exported image from Google Earth Engine

texts:
I tried to export the NDVI obtained from the Landsat 8 Surface Reflectance Collection in Google Earth Engine and I keep on getting black lines in my final export. Everything looks fine when displaying it on the map online, it's just the export to google drive that has the problem.

Here is the code I'm using

var region = ee.FeatureCollection(
  'ft:1A9vRY9l-eBNBFBwAV_mLAEeCIPYS2T14vsKhiEQc');
var crop = region.geometry();

function mask_image(image){
  return image.updateMask(image.select('pixel_qa').bitwiseAnd(ee.Image.constant(2)).eq(2));
}

function ndvi(image){
  return image.addBands(image.normalizedDifference(['B5','B4']).rename('NDVI'));
}

var start = ee.Date('2013-01-01');
var end = ee.Date('2017-12-31');
var image = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
              .filterDate(start,end)
              .filterBounds(crop)
              .map(mask_image)
              .map(ndvi);

var to_export = image.mean().select(['NDVI']);

Export.image.toDrive({
  image: to_export,
  description: 'NDVI-test',
  folder:'Landsat8-NDVI',
  scale: 30,
  region: crop,
  maxPixels:1e10,
  crs: 'EPSG: 32718'
  });

Map.centerObject(image,8);
Map.addLayer(to_export.clipToCollection(region));

This is how it looks like when opening the export image in a GIS software.


texts: Make a WFS-T Insert request to geoserver and add non spatial values to the layer

texts:
My question is somehow related with my previous one Failed to insert data to geoserver layer with WFS-T.
I use python and requests in order to create a new feature in a geoserver layer. I have managed to make a successful request and add the geometry of the feature to the layer.

Now I am trying to make the same post request but also to include other values of attributes that they exist in the layer.

What is the format in order to achieve something like this? As far as I see, the Demo page of geoserver doesn't provide an example like this.

@iant helped me figure this out. This is how my xml looks like:

<wfs:Transaction service="WFS" version="1.1.0"
  xmlns:wfs="http://www.opengis.net/wfs"
  xmlns:gml="http://www.opengis.net/gml">
  <wfs:Insert>
   <bgd_uti_transmissionline_nrel>
     <the_geom>
      <gml:MultiLineString srsName="http://www.opengis.net/gml/srs/epsg.xml#4326">
      <gml:lineStringMember>
       <gml:LineString>
        <gml:coordinates decimal="." cs="," ts=" ">
        39.247354197,21.438097723 89.247354197,18.438097723
        </gml:coordinates>
      </gml:LineString>
     </gml:lineStringMember>
    </gml:MultiLineString>
   </the_geom>
   5
   <POWERLIN_I>1</POWERLIN_I>
   <OBJECTID>2</OBJECTID>
   <GLOBALID>fd</GLOBALID>
  </bgd_uti_transmissionline_nrel>
 </wfs:Insert>
</wfs:Transaction>


texts:
Is there any good, complete API for the GDAL Python bindings?

I am interpolating and rasterizing points with elevation data using the gdal.Grid() method;

output = gdal.Grid('outcome.tif','newpoints.vrt')

In its most basic sense, it's working perfectly. But now, I'd like to try the linear approach;

output = gdal.Grid('outcome.tif','newpoints.vrt', algorithm = 'linear:radius=0')

That does generate a new .tif file (512kb), but is just one big black image. Also, not all options seem to be accepted.. for the invdist:power=4 for example it also results in a useless image. What is going wrong, am I missing something? Below is the .vrt I'm using;

<OGRVRTDataSource>
<OGRVRTLayer name="newpoints">
    <SrcDataSource>newpoints.csv</SrcDataSource>
    <GeometryType>wkbPoint</GeometryType>
    <GeometryField encoding="PointFromColumns" x="long" y="lat" z="depth"/>
</OGRVRTLayer>

EDIT
As a next example I used the gdal.GridOptions(), putting the parameters in that one. It is resulting in a correct file, but still ignores the options. How to actually pass the GridOptions onto the gdal.Grid() ?

edit
I found this question also; GDAL grid produces NoData
Is the linear interpolation just a bug in GDAL?

See my example workfiles;
https://gist.github.com/willemvanopstal/1eae4999fcfa62e89cca39a350b504d8


texts: How is user range error (URE) calculated/estimated in a basic receiver?

texts:
For a basic receiver how can we estimate/calculate URE when the pseudorange observable is used? Assume that it is a single frequency (L1) receiver.

Is there any metric that can be used to estimate this error to some extent? i.e. does the receiver have any estimate of the quality of the signal other than from the elevation of the satellite and DOP. Is the signal-to-noise ratio (SNR) which appears in the RINEX observation file for each pseudorange of any use here?


texts: Error while updating feature information to geoserver using OL3

texts:
I am trying to update feature information of wfs feature to geoserver using WFST in ol3. I am getting the below error in the response

<ows:ExceptionReport xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:ows="http://www.opengis.net/ows" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.0.0" xsi:schemaLocation="http://www.opengis.net/ows http://localhost:8080/geoserver/schemas/ows/1.0.0/owsExceptionReport.xsd">
<ows:Exception exceptionCode="InvalidParameterValue">
<ows:ExceptionText>
org.geotools.referencing.operation.projection.PointOutsideEnvelopeException: 1178724.6835985 outside of (-180.0,180.0)
Parsing failed for LineString: org.geoserver.wfs.WFSException:
org.geotools.referencing.operation.projection.PointOutsideEnvelopeException: 1178724.6835985 outside of (-180.0,180.0)
org.geotools.referencing.operation.projection.PointOutsideEnvelopeException:
1178724.6835985 outside of (-180.0,180.0)
1178724.6835985 outside of (-180.0,180.0)</ows:ExceptionText>
</ows:Exception>
</ows:ExceptionReport>

And my request is

<Transaction xmlns="http://www.opengis.net/wfs" service="WFS" version="1.1.0" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd"><Update typeName="undefined" xmlns:feature="http://localhost:8080/geoserver/web/geowebmap_ws"><Property><Name>the_geom</Name><Value><MultiLineString xmlns="http://www.opengis.net/gml" srsName="EPSG:4326"><lineStringMember><LineString srsName="EPSG:4326"><posList>1178724.6835985 8865818.20524869 1179426.69118129 8864805.19788247 1180462.9028373 8863781.32573395 1185077.23512065 8860641.8489268 1186685.37344036 8858025.82976121 1188893.85395699 8856311.50960299 1192110.55978526 8852515.52609889 1192563.64100789 8851714.02576518 1193798.31753863 8845947.67614209 1194366.28358949 8845081.20975355</posList></LineString></lineStringMember></MultiLineString></Value></Property><Property><Name>TYPE</Name><Value>primary</Value></Property><Property><Name>NAME</Name><Value>Thanjavur - Mannargudy - Thiruthuraipoondi - Vedaranyam - Kodiakkarai Road</Value></Property><Property><Name>ONEWAY</Name><Value>test</Value></Property><Filter xmlns="http://www.opengis.net/ogc"><FeatureId fid="tn_highways.15747"/></Filter></Update></Transaction>


texts:
UPDATE: 7/26/2013

The issue I'm encountering is a known bug with the netCDF-Java API described in ticket TDS-391:

A user would like to convert radar data to netCDF format. Unfortunately, the data files have structure data, which fails badly for writing netcdf-3 (illegal data type, as netcdf-3 does not support structures) or netcdf-4 (unimplemented) files.

The issue is targeted for a fix in API ver. 4.4, which will hopefully will be out this summer, as per an email from Unidata netCDF Java Support received on July 8:

We are going to implement a solution for this in version 4.4, but it won't be out before our training workshop later this month.  We may limit the implementation to NetCDF-4, and offer to skip trying to write structures in NetCDF-3.

I have many level-III NEXRAD precipitation netCDF files which I created from the raw files using the netCDF-Java API (ver. 4.3.17). I am now trying to compute monthly average rainfall rates and returning a new netCDF file which I will then need to convert to ESRI GRID format. I am new to working with netCDF files, but have read up a bit on creating dimensions and variables using the scipy.io.netcdf module.

I am using the first netCDF file in the directory as a sort of template file from which I'm grabbing the x and y locations and projection information. I have 3 dimensions:

In [4]: ncfile.dimensions
Out[4]: {'textStringSize1': 1, 'x': 131, 'y': 131}

I also have 4 variables:

In [3]: ncfile.variables
Out[3]:
{'FlatEarth': <scipy.io.netcdf.netcdf_variable at 0x67527d0>,
 'PrecipArray_0': <scipy.io.netcdf.netcdf_variable at 0x6752690>,
 'x': <scipy.io.netcdf.netcdf_variable at 0x67526d0>,
 'y': <scipy.io.netcdf.netcdf_variable at 0x6752430>}

I have successfully created my x and y dimensions and variables and populated the precip_array_0 variable with my monthly average values. The last piece of the puzzle, I think, is to copy the projection info which I believe is what the FlatEarth variable holds. I'm not sure how to look at what's in the variable:

In [5]: ncfile.variables['FlatEarth'][:]
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
C:\Users\jbellino\Desktop\<ipython-input-5-aec6027d7a77> in <module>()
----> 1 ncfile.variables['FlatEarth'][:]

C:\Python27_epd32\lib\site-packages\scipy\io\netcdf.pyc in __getitem__(self, index)
    837
    838     def __getitem__(self, index):
--> 839         return self.data[index]
    840
    841     def __setitem__(self, index, data):

IndexError: 0-d arrays can't be indexed.

In [6]: ncfile.variables['FlatEarth']
Out[6]: <scipy.io.netcdf.netcdf_variable at 0x67527d0>

In [10]: ncfile.variables['FlatEarth'][()]
Out[10]: ' '

Calling help on the object returns help for the entire variables object, not just the FlatEarth variable:

In [11]: ncfile.variables['FlatEarth']??
Type:       dict
Base Class: <type 'dict'>
String Form:{'y': <scipy.io.netcdf.netcdf_variable object at 0x06752430>, 'x': <scipy.io.netcdf.netcdf_variab <...> e object at 0x067527D0>, 'PrecipArray_0': <scipy.io.netcdf.netcdf_variable object at 0x06752690>}
Namespace:  Interactive
Length:     4
Docstring:
dict() -> new empty dictionary
dict(mapping) -> new dictionary initialized from a mapping object's
    (key, value) pairs
dict(iterable) -> new dictionary initialized as if via:
    d = {}
    for k, v in iterable:
        d[k] = v
dict(**kwargs) -> new dictionary initialized with the name=value pairs
    in the keyword argument list.  For example:  dict(one=1, two=2)

I'm assuming that the FlatEarth variable uses the textStringSize1 dimension, but I'm not sure how to create this variable in the new netCDF file or populate it once it's been created. Any help is much appreciated!


texts: sf: SpatialGridDataFrame vs SpatialPixelsDataFrame: Processing, Visualization and Storage

texts:
I'm trying to understand what are the differences between an SpatialGridDataFrame and SpatialPixelsDataFrame. My final goal is to process some raster (readed with raster package and tif format) and save it using sf writing tools.

I'm trying just this and, in large datasets, the conversion takes to long in both cases:

r <- raster('./raster/raster.tif')

# Change raster format to sf
r.sf1 <- st_as_sf(as(r, "SpatialGridDataFrame"))
r.sf2 <- st_as_sf(as(r, "SpatialPixelsDataFrame"))

Which one performs better in processing (reprojections, raster calcs)?
Which one is less ram consuming?
Which one has more capabilities for visualization?


texts:
is there a way to use/transform from the custom grid in this WMTS Capabilities to the osm grid to use the wmts with leaflet ?

I'm using the leaflet-tilelayer-wmts plugin and i cannot display any tile from the wmts linked above due to other tilenames and grid. Both my map and my requests are in EPSG:3857.

var ortho = L.tileLayer.wmts('https://geodienste.sachsen.de/wmts_geosn_dop-rgb/guest?', {
    layer: 'sn_dop_020',
    style: 'default',
    format: "image/png",
    tilematrixSet: 'grid_3857',

});


texts:
Okay so I'm using QGIS 2.10 (a couple of versions out of date, but shouldn't matter much, eh?) on Windows 7, and I can't figure out this supposedly a very basic geoprocessing tool DIFFERENCE.

I'm trying to use a small polygon to find the DIFFERENCE in a larger polygon that completely covers the smaller polygon, so the result should be the larger polygon minus the smaller polygon, and the shape would look like a donut or bagel, but every time the result comes back as empty.

I'm not sure what I'm doing wrong here, I tried switching the layers and nothing worked either.

I mean, there's only 2 inputs and 1 output for this tool, how difficult could it be?


texts: image size has been increased after i.atcorr

texts:
I have compiled i.atcorr and run for different data sets.  The problem is, it is writing the output 4 times of input image. I checked for data types, they are fine. I am not able to catch where exactly problem is. Can anyone suggest me?

I also want to mention that image first processed for toar using grass gui 6.4.4 and then the image passed for i.atcorr (customized for our purpose) actually in toar itself it is getting doubled. Is there any problem in gui version of i.landsat.toar? I have used for landsat5 tm data path141row43. Please suggest me.


texts: Slope analysis from NED data gives crazy histogram?

texts:
I'm attempting a simple slope analysis in QGIS2 using the plugin at Raster -> Terrain Analysis -> Slope. The slope output layer has a very strange histogram, basically spikes at 0 and at 89.9098 with one big valley in the middle:

Here's the initial DEM:

And the resulting slope layer (note that Contrast Enhancement is already at 'stretch to minmax'):

The place I got the DEM is the NED. I tried the gridfloat format, and the arcgrid format, and the IMG format. I downloaded from nationalmap and also from atlas.ca.gov - all seem to have the same resulting histogram.

The procedure I used to make the slope layer was Raster -> Terrain Analysis -> Slope, GeoTIFF, Z factor = 1.0 (what does that do anyway?)

Seems like simple pilot error; however: I found this tutorial online (for QGIS 1.something) that comes with a sample DEM, and it worked just fine using the same procedure. So, it's not a fundamental QGIS2 bug.

So, it must be slightly more complex pilot error: there must be a conversion step or something I'm missing?


texts: Add metadata to dataset (feature class, coverage or shapefile) from a .csv using arcpy?

texts:
I am looking for a script that will take metadata information provided in an excel table, convert it to a .csv and import the values to populate some of the dataset's ArcGIS metadata elements.

I came across this post today and it is essentially the reverse of what I need (it grabs metadata elements and exports them to a .csv) and thought it might help:

Creating a table containing all filenames (and possibly metadata) in a File Geodatabase

I am a novice programmer (with some python under my belt) and am new to this site.


texts:
Are there any good resources on WFS support in ArcGIS?

The only documents I found so far are these: How_to_work_with_WFS_using_ArcGIS_Desktop and Consommer des services/données WFS dans ArcGIS 9.2 But I want to know how it works!

I'm interested in how ArcGIS can be a client to a WFS Server.

Seems like the Interop conecction is downloading all data on the WFS server into a local Cache, is there any way to get it to fetch the features on-demand?

Update:
I voted for this idea Make ArcMap an effective WFS Client to get better support for WFS.


texts:
I have two shapefiles - Parcels2016 and Base_Zoning_B

Parcels2016 has about 660,000 records with over one hundred attributes (about ~2GB) , so performing spatial joins in QGIS is very inefficient, which is why I am trying to use PostGIS. I would like to find out what is wrong with the SQL command that I used to join attributes by location.  If we look at the images below, the Base_Zoning_B layer is in green:

My goal is to join the attributes of Base_Zoning_B to the individual land parcels in Parcel2016 (in blue) and keep  both matching and non-matching records:

Both commands I have used in PostGIS failed to produce the correct results as not all of the polygons were joined and it seems only matching records were kept, and not all of them either:

As you can see, in the testlayer, which I exported from PostGIS as a shapefile, there are many parcels excluded.  I have executed both of these commands in PostGIS, one with a left join, one without, but both produce the same incomplete results:

CREATE TABLE test AS
SELECT t.gid, t.geom, m.zone_dist, m.pd_num, m.council_da, m.long_zone_, m.globalid, m.ord_num,
m.notes, m.res_num, m.effectived, m.case_numbe, m.common_nam, m.districtus, m.cd_num
from parcel2016 AS t, base_zoning_b AS m
where st_within(t.geom, m.geom) AND ST_IsValid(t.geom) AND ST_IsValid(m.geom)

CREATE TABLE test AS
SELECT t.gid, t.geom, m.zone_dist, m.pd_num, m.council_da, m.long_zone_, m.globalid, m.ord_num,
m.notes, m.res_num, m.effectived, m.case_numbe, m.common_nam, m.districtus, m.cd_num
FROM parcel2016 AS t
LEFT JOIN base_zoning_b AS m
ON ST_Within(t.geom, m.geom)
WHERE ST_IsValid(t.geom) AND ST_IsValid(m.geom)

I also received the following warning messages:

NOTICE:  Ring Self-intersection at or near point 2489538.1553773778 7045488.56387801
NOTICE:  Ring Self-intersection at or near point 2458823.8434288907 7006314.9953470686
NOTICE:  Ring Self-intersection at or near point 2497326.4063920071 7047098.5640874337
NOTICE:  Ring Self-intersection at or near point 2509123.1707222643 6916240.9188464647
NOTICE:  Ring Self-intersection at or near point 2485148.4026713888 6957535.547359148
NOTICE:  Ring Self-intersection at or near point 2468946.5905914963 6947132.4664672073
NOTICE:  Ring Self-intersection at or near point 2483211.8534999937 6963347.7294999957 0
NOTICE:  Ring Self-intersection at or near point 2500677.5199999958 6921843.0100000054 0
NOTICE:  Ring Self-intersection at or near point 2505740.7239122842 7002511.0560633009
NOTICE:  Ring Self-intersection at or near point 2476637.4922218379 6988535.7446923936
NOTICE:  Ring Self-intersection at or near point 2539170.4127625637 6940842.0433131075
NOTICE:  Ring Self-intersection at or near point 2497617.4059891542 7004481.5565984966
NOTICE:  Ring Self-intersection at or near point 2483918.2574999928 6975295.9455000013 0
NOTICE:  Ring Self-intersection at or near point 2479238.1964862673 6956274.0468138419
NOTICE:  Self-intersection at or near point 2430764.1748952293 7027762.7618717942
NOTICE:  Ring Self-intersection at or near point 2508302.6419177013 7006131.0566753978
SELECT 183826

Query returned successfully in 1 min.

For some reason, only 183862 records were returned, when there are actually a total of 665518.


texts: The screenshot shows a calculate field tool in modelbuilder.  The error I'm getting is invalid syntax line 1.  I can't see any issues w/the syntax so I'm hoping someone out there can help. It calculates field A to equal field B if field A is not equal to field B.


texts:
I have an old map over a part of Laos and wondering how do I georeference the map with the help of the Grid lines and other Longitude and Latitude information on the scanned map? I have searched for information on the internet and tried many times. My scanned map wont overly correctly on the base map. Any base map would be fine. As long as I can check the result.

Do I add X/East and Y/North values from the grid line crossings as marked and written in some parts of the map? How do I begin georeferencing this map? Do I have to add zeros after the two digits? For example if the Easting is 24 (its in KM right? Do I add like three/four zeros to convert it to meters?

I georeferenced the map with the option "From Map Canvas" in the georeferencer tool and I succeeded but I want to know the coordinate of a point on my old scanned map and write them as X/East and Y/North to georeference the map.

I hope I am not babbling and someone understands what I mean!

We can set the CRS to WGS84, EPSG:4326.


texts: Splitting polyline into points at certain intervals

texts:
I am trying to split a polyline into points with 2m intervals in ArcGIS but I don't have the necessary license.

Is there a way that I can do this in QGIS?

I have drawn a hypothetical lake surface based upon a contour line and need to divide it.


texts: Splitting polygon in equal parts based on polygon area?

texts:
So, I came across this relevant question, but I think my problem is a bit different.

I have two tables i.e. polygon poly and points pts in my Postgres db. I am on Windows 10 (x64) based machine and using PostgreSQL 9.5.12/PostGIS 2.3. Some of the polygons contain points inside them while, some do not. The sample scenario is given as follows:

For each poly, the following query computes the area and counts number of pts inside polygon.

Select
    a.gid poly_id, count(b.geom) pt_cnt,
    st_area(a.geom)::int poly_area
FROM
    poly a
LEFT JOIN
    pts b
ON st_contains(a.geom, b.geom)
GROUP BY a.gid
Order by a.gid;

The area of the polygons (from left to right) is 1079, 744, 340 meter-square respectively. Now, I would like to split these polygons' geometry based on the area i.e. if area is above 100 meter-square, then each polygon needs to be divided/split into approximate equal parts. Can somebody suggest me a workaround to fix the issue?

UPDATE

Here, is another example of the polygons in my data. I am aware that there can be much more complex examples (e.g. irregular polygons). In that case, I think exceptions can be added.


texts: QGIS: relation to Excel exported files (updates)

texts:
Layers exported from QGIS as Excel format are the topic here.
I see when I edit this files in Excel(exported from QGIS) are simultaneously edited in equivalent file (non-spatial table) which was created during export (and lies in project self).
What I really want is to merge,  (better to say overwrite) updated columns in "original" spatial layer from equvalent non-spatial table ...
I tried with -Layer Properties -Joins functionality in QGIS but it adds a column which I selected and it doesn't update/overwrite selected Column/Field.
Is this possible with standard functionalities in QGIS at all? Or Python plugin is only option (or even maybe PostGIS Database)?
P.S. I have also tried with Excel Sync plugin but seems this updates from  QGIS to Excel and I want reversely


texts: PostGIS, create centroid column from polygon column

texts:
I have a table like this:

id|      shape [geometry]     |  centroid [geomety]
---------------------------------------------------
 1|MULTIPOLYGON(((-70, etc...)|
 2|MULTIPOLYGON(((-71, etc...)|

I need to create the centroid of each shape, how can I do this? The final result should be this table:

id|      shape [geometry]     |  centroid [geomety]
---------------------------------------------------
 1|MULTIPOLYGON(((-70, etc...)|POINT(something)
 2|MULTIPOLYGON(((-71, etc...)|POINT(something_else)

Edit: should be some update with st_centroid(shape) I guess


texts: Exporting MapImageLayer, however getting a blank png

texts:
I am pulling data from Wayback Imagery and attempting to export an image after making it into a arcgis.mapping.MapImageLayer. When I run my code the image it produces is blank and I'm unsure why. The bounding box should be plenty large enough to see something.
import arcgis
import re
from arcgis.gis import GIS
from arcgis.mapping import WebMap
from IPython.display import display

gis = GIS("https://www.arcgis.com", "arcgis_python", "P@ssword123", verify_cert=False)
result = gis.content.search(query="title:Wayback",sort_field="numViews",
                            outside_org=True,max_items=1,item_type = "Map Service")
ids = '903f0abe9c3b452dafe1ca5b8dd858b9'

bbox=bbox=-141.19530416221985,-62.217823180545146,139.27427961579508,84.15317625109763
size = 600,550
format = "png"

print(len(result))
count=0
for item in result:
    print(item)
    map = gis.content.get(item.id)

    if(re.match(r'(World Imagery \(Wayback ).*',map.title)):
        print(map.title)
        print(type(map))
        print(map)
        display(map)
        wm = WebMap()
        wm.add_layer(map)
        print(wm.layers[0]['title'])
        print(wm.layers[0]['url'])
        mapImage = arcgis.mapping.MapImageLayer('https://metadata.maptiles.arcgis.com/arcgis/rest/services/World_Imagery_Metadata_2014_r01/MapServer')
        returnJson = mapImage.export_map(bbox,image_format='html')
        print(returnJson)
        count+=1
    else:
        #will do more stuff later
        pass

Below is the output I am receiving:
<Item title:"World Imagery (Wayback 2014-02-20) Metadata" type:Map Image Layer owner:esri_imagery>
World Imagery (Wayback 2014-02-20) Metadata
<class 'arcgis.gis.Item'>
<Item title:"World Imagery (Wayback 2014-02-20) Metadata" type:Map Image Layer owner:esri_imagery>

World Imagery (Wayback 2014-02-20) Metadata
This World Imagery Metadata layer provides detailed information for each image source in World Imagery (Wayback 2014-02-20).Map Image Layer by esri_imagery
Last Modified: December 11, 2019
0 comments, 59,777 views
1.9cm Resolution Metadata
https://metadata.maptiles.arcgis.com/arcgis/rest/services/World_Imagery_Metadata_2014_r01/MapServer/0
{'href': 'https://metadata.maptiles.arcgis.com/arcgis/rest/directories/arcgisoutput/World_Imagery_Metadata_2014_r01_MapServer/_ags_mapc034e8d33e09406d9405e44efe3de69c.png', 'width': 600, 'height': 550, 'extent': {'xmin': -36067515.01650002, 'ymin': -17088899.391186584, 'xmax': -5180428.986109771, 'ymax': 11224262.803337809, 'spatialReference': {'wkid': 102100, 'latestWkid': 3857}}, 'scale': 405341525.67215466}

I take the https//* and paste it into Google and get a blank image.


texts:
I've installed QGIS 3.4 on my Mac. I don't have any vector or processing menu. In the "plugins" menu under settings, I have the following:

"The settings on this tab are only applicable for Python Plugins. No Python support detected, thus no settings available".

I've tried to troubleshoot this problem on my own in the following ways:

re-installing QGIS.
Under QGIS > Preferences > System, adding the following custom variable under the environment section:

This is as directed by several other threads I have referenced but still no results. Can anyone help?