texts: Manually edit shapefile coordinates (polygon) QGIS 3

texts:
In the QGIS2.18 version I was used to see the point editor table when I would be in editing mode and have clicked a point. A list would appear, with all the points in the selected polygon and their x and y coordinates. Even better, and most importantly for my job, I was able to type in the coordinates in the table and the point would move to that exact location. This was a very useful interactive feature.

In the latest 3.0 and higher versions, the point editor table is no longer visible. Does anyone know if it is still featured somewhere in the program?

If not, does anyone have an idea for me how to do the following:
I have a list of coordinates that should make up the polygon. They have to be exactly those coordinates, so clicking the points on the map with the mouse is not accurate enough. I just want to type them and preferably be able to manually edit them later too.

Looking for online help points me to the Numerical Vertex Edit plug-in for QGIS, which would install a new button in the shape menu bar and allows you to numerically edit the vertex of a polygon typing the coordinates. But this plugin does not seem to be available for version 3.4.


texts: Accessing attribute table within shapefile and replace values?

texts: I would like to find a Python script that will access an attribute table within a shapefile and either update or replace records inside a field column. More specifically, I would like to use the replace expression to convert a table record like "123-456-789" to "123456789". I have many records to process and would like the best approach to do this task using Python.


texts: Using grayed out distance/direction tool in ArcGIS for Desktop?

texts:
I have a set of five points (all with specific latitudes and longitudes) and I want to create a set of points which will be about 80 kms from the first set. I understand that I can use the distance-direction tool in the edit toolbar of ArcGIS, but although I tried several times, I am unable to activate the distance-direction tool. It is grayed out even when I am selecting the original point from which I want to measure the direction and distance and locate another point.

Can someone advise on where I am going wrong?

I'm using ArcGIS 10.2.2 for Desktop.


texts: QGIS hyperlink to specific row/column excel

texts: Is there a way to set a hyperlink to a specific sheet/row/column in an excel document? I am able to link to an excel file on the server and open it from QGIS, but would like to open the document to a specific sheet or row within the excel document. Is there a specific modifier on the file path needed?


texts:
I'm trying to populate a field from another. The catch is that it isn't 1 or 2 values. If my field A = 4 then I want it to return field B to 0. But if Field A is 1, 2 or 3 then I want the user to be able to use the domain drop down list to choose what field B is--because it could be one of 3 other values. Is that possible? Below is the code I think is a start. But I need it to return 1, 2 or 3 if WaterType = anything but 4.

def Reclass (IrrigType):
  if WaterType = 4:
     return 0
  else:
     return 1


texts: KML imports as a polygon without keeping any attributes in QGIS

texts:
I am extremely new to GIS mapping.

I'm doing RF mapping with QGIS 2.18.9, and whenever I try to import a KML (which views just fine in Google Earth as an RF signal and tower), importing the polygon just shows up as a block. I've tried looking on this site and others, but nothing suggested seems to help.

Here's what I start with (a KMZ in Google Earth):

The steps I do are:

Create KMZ with CloudRF
Download as KML
Import to QGIS by Layer>Add Layer>Add Vector Layer (Type: File, Encoding: System, Source: Sample_GoogleEarth_KML_Name
Select vector layers to add (Geometry type: point and polygon)

And this is the output I get (in QGIS):

How do I properly import the KML so it keeps the RF shape/etc without becoming a solid block with no attributes?

I was actually able to get it configured and it showed up! CloudRF was finally able to export as an accurate SHP file, so I had the attributes I needed and it worked perfectly.


texts: Calculating percentage of overlap of two layers in QGIS 3

texts:
I am trying to compare two different vector layers on QGIS 3. They are species distribution maps.
I would like to calculate the percentage of overlap of two vector layers. My goal is to determine how similar the two ranges are.
Anyone know how I can do this in QGIS 3.8?


texts: Shortcuts for QGIS python console commands

texts:
I would like to define a shortcut for "Save" in the python console, is there a way to achieve that ?

Is it possible to define shortcuts only for the given context of the python console and make them available (for other tasks) in other contexts when console is closed (or ... inactive ?) ...


texts: I have a points dataset, which contains what are essentially vertices of polygons. They represent a number of separate, discrete polygons, and there is a field delineating these by name. Using this field, upon executing a 'Points to Line' tool, I am able to specify which points belong to which polygon. However, I can't specify which points are to be drawn in which order. The only other field is FID, and the points are not in the correct order such that it creates the proper polygon shape. For example, a polygon that is supposed to be a circle comes out with lines drawn between vertices on opposite sides, rather than the circle it is supposed to be. Is there a way to get the polygons to generate properly?


texts: Error loading heatmap layer using cartodb.js

texts:
I'm trying to load a heatmap layer using cartodb.js, but it is failing whether I use createVis or the manual way.

Using createVis:

cartodb.createVis('mapDIV', 'https://PATH/viz.json');

Building a Leaflet map and then adding it:

var layerUrl = 'https://PATH/viz.json';
var cartoDB_lyr = cartodb.createLayer(map, layerUrl)
cartoDB_lyr.addTo(map);

Either way I am getting the following error message:

I am not sure how to debug this.

If I disable torque in CartoDB (i.e. visualizing the dataset as simple), then the dataset displays without a problem, so it must be specifically related to the heatmap:


texts: Authenticating ArcGIS Desktop to web tier authenticated ArcGIS Server?

texts:
I'm using a webadaptor with ArcGIS Server 10.4 to which I added an authentication filter (CAS client) via the web.xml.  Everything works perfectly for the web application.  However, when I try to connect via the ArcCatalog it will not authenticate my user.  If I use the built-it siteadmin account, I can connect, but as that user does not exist in my enterprise identity store, it is clearly not authenticating at the web tier.  I captured the traffic from ArcCatalog in my filter and can see that it initially requests:

GET /cas-web-adaptor/admin/publicKey?f=json

This leads me to believe that it is going to use the public key to encrypt the username/password on and send that as a second request.  If it is doing that, then only ArcGIS Server would be able to decode the username/password so it would not be using web tier for authentication...   Am I missing something here (I hope so...)? Or is ArcCatalog incapable of authenticating at the web tier?


texts: Chart JS chart not rendering in ArcGIS popup when using navigation arrows

texts:
I'm putting a chart in the popup of my points. On the initial click, the popups display the chart and if I use the forward navigation arrow on the popup, the next chart renders. However, if I use the back arrow the chart will not render again. The method that generates the charts doesn't seem to be firing the second time.

My code is below. Popup:

            popupTemplate: {
            title: "{Location}",
            content: [{
                type: "text",
                text: "Sample Location: {Type}
 Survey: {Survey}
 {Location:getSurveyInfo}
<canvas id='chartArea'>{Location:createChart}</canvas>
"
            },
            {
                type: "attachments"
            }]
        }

Chart:

    createChart = function (location) {
    var date = new Date();
    var chartArea = "chartArea";
    var sub = location.substring(0, 2);

    getData(date.getFullYear(), [[location]], function (data) {
        var maxScale = Math.max(...data[0].tData);
        var colors = [["#34eb58", "#0000ff"]];
        var chartData = buildChartData(name, maxScale, data, colors, null);
        var ctx = document.getElementById(chartArea);
        var chart = new Chart(ctx, chartData);

    }, getChartDataError)

}


texts:
Is there a way of calculating the distance to the nearest points in QGIS? In ArcGIS there is a function called Euclidean Distance, which can take an input vector or raster, and generate a raster of distances.

Is there a similar function in QGIS? I could not find one in the documentation. I did explore the Raster functions in QGIS, and also the Hub Distance calculation fromMMGIS (for a non-raster solution) but I was not able to figure out how do this using either. If a vector approach is the only solution, that could be fine also.


texts: Calculating Voronoi Diagrams for polygons

texts:
There are a variety of good algorithms for generating the Voronoi polygons or their complement, the Delaunay Triangulation for a set of points.

My question is simply, is there an algorithm for generating the Voronoi diagram for a set of input polygons, rather than points?

One technique I've explored is breaking my polygons into sets of vertices, and creating the Voronoi diagrams for those, then combining the resulting shapes for each set of vertices belonging to a particular input polygon. The results are not totally accurate, however. Does anyone have an alternate technique?

EDIT:

Here is a super rough hand-drawn example of what I'm after. I have a set of polygons with gaps. I'm trying to create output polygons with no gaps between them. Ultimately I want to use this to say whether any two nearby polygons can be considered "adjacent" to one another, even if they are not touching.


texts: Mapfile errors, unable to access my shapefiles

texts:
I am getting the following errors in my mapfile

[Thu Oct  3 11:40:23 2019].299949 msWMSLoadGetMapParams(): enabling non-square pixels.
[Thu Oct  3 11:40:23 2019].303483 msDrawMap(): kicking into non-square pixel preserving mode.
[Thu Oct  3 11:40:23 2019].303495 msDrawMap(): rendering using outputformat named png (AGG/PNG).
[Thu Oct  3 11:40:23 2019].303499 msDrawMap(): WMS/WFS set-up and query, 0.000s
[Thu Oct  3 11:40:23 2019].303563 msShapefileOpen(): Unable to access file. (/usr/local/var/www/data/ne_10m_admin_0_countries/ne_10m_admin_0_countries)
[Thu Oct  3 11:40:23 2019].303615 msShapefileOpen(): Unable to access file. (/usr/local/var/www/data/ne_10m_admin_0_countries)
[Thu Oct  3 11:40:23 2019].303624 msDrawMap(): Image handling error. Failed to draw layer named 'Countries'.
[Thu Oct  3 11:40:23 2019].303682 mapserv request processing time (msLoadMap not incl.): 0.004s
[Thu Oct  3 11:40:23 2019].303690 msFreeMap(): freeing map at 0x7ffd85006000.

The mapfile is really simple

MAP
  NAME "MAPSERVER QUICKSTART"
  EXTENT -137 29 -53 88
  UNITS DD
  SHAPEPATH "/usr/local/var/www/data/ne_10m_admin_0_countries"
  SIZE 800 600
  CONFIG "MS_ERRORFILE" "/tmp/ms_error.txt"
  DEBUG 5
  IMAGETYPE PNG24

  PROJECTION
    "init=epsg:4326"
  END

  WEB
    METADATA
      ows_enable_request "*"
    END
  END

  LAYER
    NAME "Countries"
    GROUP "default"
    STATUS ON
    TYPE POLYGON
    DATA "ne_10m_admin_0_countries"
    CLASS
      STYLE
        COLOR 246 241 223
        OUTLINECOLOR 0 0 0
      END
    END
  END

END

Then when I list the contents of the directory I get this

ls -latr /usr/local/var/www/data/ne_10m_admin_0_countries
total 18448
-rw-r--r--@ 1 tommie  staff     2140 May 21  2018 ne_10m_admin_0_countries.shx
-rw-r--r--@ 1 tommie  staff  8800376 May 21  2018 ne_10m_admin_0_countries.shp
-rw-r--r--@ 1 tommie  staff      147 May 21  2018 ne_10m_admin_0_countries.prj
-rw-r--r--@ 1 tommie  staff   599487 May 21  2018 ne_10m_admin_0_countries.dbf
-rw-r--r--@ 1 tommie  staff        5 May 21  2018 ne_10m_admin_0_countries.cpg
-rw-r--r--@ 1 tommie  staff        7 May 21  2018 ne_10m_admin_0_countries.VERSION.txt
-rw-r--r--@ 1 tommie  staff    24083 May 21  2018 ne_10m_admin_0_countries.README.html
drwxr--r--@ 9 tommie  staff      288 Oct  3 11:03 .
drwxr-xr-x  4 tommie  admin      128 Oct  3 11:44 ..

Why would it be unable to access my files.


texts:
I have been having issues deleting large numbers of features (> 1,000,000) from a feature class in ArcMap. I can't seem to get higher than a few thousand before it freezes indefinitely, sometimes overnight, before crashing.

Selecting the features I need and exporting is not an option here because I have to work within a specific gdb with a predetermined structure, and I'd prefer to just delete the features and not mess with deleting/replacing entire layers.

The layer contains approx. 1.2mil records. It's a file geodatabase.

Any ideas?


texts: In a land cover raster (categorical map) that is a product of classified satellite image, what is the best way to remove "salt and pepper" effect from the image after classification. I am not interested in performing classification again, I am rather interested in a procedure that would use some kind of local functions. I use GRASS 7.


texts: How to perform spatial computation in a Mercator projection in turf?

texts:
TL;DR: Spatial computation with turfjs does not work well in a EPSG:3857 map.

In a OpenLayers map, I use turf for transforming a line into a polygon by computing the offset line (aka side-buffer for instance in PostGIS).

Here's the process:

Draw a line using OpenLayers
Compute an offset line by a given distance offsetDistance using turf.lineOffset()
Gather all points and convert to a polygon.

And here's the code

const offsetLine = turf.lineOffset(line, offsetDistance / 1000);
const line2polygon = turf.lineString([
    ...turf.getCoords(line),
    ...turf.getCoords(offsetLine).reverse()
]);
const polygon = turf.lineToPolygon(line2polygon);

The process works well except that there is a projection issue. In the image below I drew a line along the soccer pitch and I enter an offset distance of 53m, which is the width of the pitch. The problem is that the offset line is computed in the WGS84 coordinate system so the distance is actually too small (it is 53 m * cos(latitude)) and it does not result in a square polygon. The OpenLayers map is in EPSG:3857 and I would not like to change that.

I manage to do something OK by converting first the line to Mercator and by multiplying the offset distance by a scale factor (1.783 * 100000) that I got by trial-error:

const offsetLine = turf.lineOffset(turf.toMercator(line), offsetDistance / 1000 * 100000 * 1.783);
const line2polygon = turf.lineString([
    ...turf.getCoords(line),
    ...turf.getCoords(turf.toWgs84(offsetLine)).reverse()
]);

const polygon = turf.lineToPolygon(line2polygon);

And I got this:

So the result is OK but is there a better way to do it? And how can I compute my scale factor properly (I guess it is dependent on the latitude but I could not find the formula)? Surprisingly, my scale factor is not equal to cos(latitude).


texts: React Native - react-native-mapbox-gl vs react-native-maps

texts:
I have to develop apps with React Native for mobile that concern with mapping and found out there are two libraries for Native on this - react-native-mapbox-gl and react-native-maps. From github pages and tutorials I see no big difference on what they can achieve. I can't decide which one to choose and can't afford to change later because one does not meet what I need. What I want the apps to do is to:

show markers and routes
work great with GeoJson and PostGIS
better performance on both Android and iOS

Which one would have advantage over another?


texts: Is it possible to perform spatial join in PostGIS and keep all features?

texts:
I am attempting to create an employment change map using a hex grid. I have individual establishment data as points with attached employment figures. When I ran the spatial join using the postGIS extension in QGIS 2.6 it returns all of the hexes where the join is true, as expected. But for cartographic sake I would like the whole net return where I could symbolize the null cells myself and not just have them drop out.

Is there a way to do this through postGIS similar to the return all features option in ArcMap's spatial join? My sample code right now.

--view joining employment data to hexagon grid
   CREATE VIEW food_employment_2yr_hex AS
SELECT a.gid AS gid, sum(b.avgemp) AS Emp_2010, sum(c.avgemp) AS Emp_2012,   a.geom AS geom
    FROM hex_grid_pdx_nad83 a, qcew_2010_nad83 b,
    qcew_2012_nad83 c
    WHERE ST_Within(b.geom, a.geom) AND ST_Within(c.geom,a.geom)
    GROUP BY a.gid;


texts: Rendering new stratigraphic classification using QGIS Midvatten plugin on plot

texts:
I've a borehole with a new stratigraphy classification: lignite.
First I've inserted (add) the new strata item into the zz_statigraphy_plot table and into the zz_strat table.

zz_strat table:

geoshort | strata
-----------------
Lignite  | lignite

zz_stratigraphy_plots table:

strata  | color_mplot | hatch_mplot | color_qt  | brush_qt
------------------------------------------------------------
lignite | DarkGreen   | -           | darkGreen | HorPattern

obsid   stratid depthtop    depthbot    geology geoshort    capacity    development comment
---------------------------------------------------------------------------------------------
TEST01  1       0           8           peat    peat        0
TEST01  2       8           12          clay    clay        2
TEST01  3       12          16          lignite lignite     1

Then I have inserted (add) the stratigraphy information into the stratigraphy table. I am able to view the stratigraphy plot but the fill color in the logplot is just white in the new stratigraphy item.

How can I fix it?


texts: Estimate spatial kernel density from raster and output kernel values in a grid

texts:
I would like to use a raster map (find below) to estimate or create a kernel density distribution. The objective is to extrapolate the kernel density function and re-use it in order to distribute a traits in space according to such kernel density function.

Do yo have any suggestion?

I tried using the sp.kde function in SpatialEco R package. However, it estimates a kernel density in continuos terms, while I would like to distribute the kernel density values in a 23 x 23 grid (space). In other words, I would write in each cell of the picture below the corresponding kernel value.


texts: Slow rendering with layers when making WMS call on GeoServer

texts:
I am making WMS calls (L.tileLayer.wms) to GeoServer and rendering the layers on a Leaflet map for a production environment, so I have maybe about 100+ concurrent users that might be constantly rendering layers
The issue that I am running into is that the WMS call to GeoServer can be super slow that the layers do not render at all. I can see in the Network calls that the WMS call is still pending. I was wondering how I can optimize GeoServer (or Tomcat) to speed up the layer rendering?
For background, I am using GeoServer version 2.17.2, hosting it on a AWS EC2 instance, and using Tomcat to run GeoServer. As for the CPU and Memory for the EC2 instance, the CPU is at 8 and 16GB for memory. For the store, I am using a PostGIS Database that is AWS RDS and all of my data is spatially indexed. Here is some info about the connection parameters:
Max Connection: 150
Min Connection: 1
Fetch size: 1000
Batch insert size: 100
Connection timeout: 20
Max Connection: 150
Max Connection: 150

I have already done GeoWebCaching for all of my layers on GeoServer already and all of my layers are vector data. I am not using raster data at all for this.
I am not sure what else I can do to optimize GeoServer to render my layers faster.. Should I upgrade the EC2 instance to have a higher CPU and memory to ensure faster speed? Or maybe do improvements on Tomcat? Some of my layers consist of 10,000+ rows so they can be fairly large datasets.
Another thing I would like to add is that I cannot access the Geoserver UI at all when the rendering of the layers are super slow. It takes forever for it to load the Geoserver webpage / it doesn't load whatsoever unless I restart the tomcat


texts: QGIS v/s MapInfo - Differences in Distance Measurements

texts:
Distance calculated in different software (MapInfo and QGIS) has given me differences in distance values. To be more precise, it's mostly happening when the distance is more than 25 miles. I doubt this is because QGIS is always including the earth curve (3D) while calculating distances, while in MapInfo, it's directly measuring the distances in 2D form.
Why I am sure about 3D and 2D impact is because the distance I am getting is always huge in QGIS when compared.
If this is the case, how can we calculate the distances in MapInfo with the earth curve included in measurements and should come similar to QGIS?


texts:
I am attempting to Zoom from an Access record to a specific polygon in ArcMAP but I can't seem to get the function to work.
The field in Access is identical to the field in arcmap. I.E. they both have a unique polygonID (custom field) identifier.

I have been able to open the specific record in Access from ArcMAP but would like the functionality to work both directions.

Does anyone know if there is any examples of this online and if it is even possible/worth doing with VBA being phased out.


texts: I have a road network in QGIS composed of the nodes layers and the links layer. I need to extract a node link table from this which would give the start node, end node and the distance for every link in the network.


texts: PostGIS: Assign ID of point in layer A to closest point in Layer B

texts:
This should have been an obvious precursor (that I didn't ask) to my other question:
How to create spider diagrams (hub lines) in PostGIS?

If I do not know the relationship between a point in layer A (stores) and a point in layer B (customers), I would like to generally say "Customer 1 is serviced by the nearest store."  While I realize this fact may not be true, it can be a decent surrogate.

Using PostGIS, what is the most efficient way to assign the ID of the nearest point in layer A (stores) to each point in layer B (customers).  The output I am looking for is something like below.

Customer | Store
    1    |   A
    2    |   A
    3    |   B
    4    |   C


texts: Buffering continuously throws error "Background Server Threw an Exception" (Error 001143). Why?

texts:
I have written a code to create 1 meter buffers for features of converted feature classes (from rasters). I want to dissolve by the field "Class" so that I will have significantly smaller datasets to work with because the converted feature classes are far to large to do much with. I was able to do the first one using the GUI so it is possible but when I try to run my code I get error 001143. The error does not pop up right away, instead I see it about an hour after execution (the first file is pretty big). I do get the message "[file] is ready for buffer" (Line 13). The error message says the problem is in line 17 which is when I execute the buffer. Are my buffer parameters incorrect? I disabled background processing but nothing changed. I am happy to hear any suggestions and answer questions that you may have.

import os
arcpy.env.workspace = r"C:\Data\ConvertedPolygons.gdb"
arcpy.env.overwriteOutput = True
arcpy.env.addOutputsToMap = False
polygdb = arcpy.ListFeatureClasses()
for fc in polygdb:
    name = fc.split(".")[0]
    if "Test" in name:
        print fc + " was skipped"
    elif "merge" in name:
        print fc + " was skipped"
    else:
        print fc + " ready for buffer"
        outws = r"C:\Data\BufferedPoly.gdb"
        outfc = name + "_buff"
        outfile = os.path.join(outws, outfc)
        arcpy.Buffer_analysis(fc, outfile, "1 Meter", "FULL", "ROUND", "LIST", "Class" )


texts: Referencing CSV file in Select By Attributes SQL query?

texts:
I have 2 databases:

Database 1: Point Shapefile with adress and APN (parcel numbers) as attributes. (Huge number of points, addresses for an entire metropolitan area).

Database 2: CSV file with a smaller database of addresses and APN numbers, and a few other fields. This one is not georeferenced.

I would like to construct a query in Select by Attributes where I reference database 2 and search for APN's within database 1 corresponding to all the ones in database 2. Ultimately I would like to isolate all corresponding points and display them on the map.

This is like a batch form of the following query:

SELECT*FROM ADDRESS_APN WHERE "APN" = '###' AND "APN"='####' AND ......

Only thing is there are about 565 values so it is unfeasible to have an AND statement for each value. I would much rather reference the CSV file.


texts:
A dummy question. I need to change all values in my rasters into ones. I've been able to do this fine in the Rater Calculator just by multiplying the layer by zero and adding one (my_layer*0+1), but when I'm doing the same in arcpy, for some reason it works for a part of my raster but for some it changes all the nodata cells to ones as well. The nodata values in different rasters are different because my rasters have been generated from different source data. Does it not recognize some nodata values? If so, this is odd as it displays them correctly when I just import them on the mxd. Any ideas on why this works sometimes, but sometimes not?

Here's the exact script:

bin_test = arcpy.Raster("F:/datasets/map1.img")*0+1


texts: Using variable in SelectByAttribute_management

texts:
I am new in Python. I was trying to use a variable for selecting a feature instead of giving a direct input. My code sample will explain it better:

import arcpy
from arcpy import env
env.overwriteOutput = True
env.workspace = "H:/shp"

#variables
parcel = "parcel.shp"
parcel_Layer = "parcel_Layer"

Make Feature Layer
arcpy.MakeFeatureLayer_management(parcel, parcel_Layer, "", "", "FID FID VISIBLE NONE;Shape Shape VISIBLE NONE;OBJECTID OBJECTID VISIBLE NONE;ACCOUNTID ACCOUNTID VISIBLE NONE;PropID PropID VISIBLE NONE")

parcel_Layer_selattr = arcpy.SelectLayerByAttribute_management(parcel_Layer, "NEW_SELECTION", "\"PropID\" = 142855")

This code works fine as long as I enter the property ID directly ("\"PropID\" = 142855"). But when I declare the PropID as a variable, it doesnot work. Here is the code sample:

import arcpy
from arcpy import env
env.overwriteOutput = True
env.workspace = "H:/shp"

#variables
parcel = "parcel.shp"
parcel_Layer = "parcel_Layer"
x = 142855
Make Feature Layer
arcpy.MakeFeatureLayer_management(parcel, parcel_Layer, "", "", "FID FID VISIBLE NONE;Shape Shape VISIBLE NONE;OBJECTID OBJECTID VISIBLE NONE;ACCOUNTID ACCOUNTID VISIBLE NONE;PropID PropID VISIBLE NONE")

parcel_Layer_selattr = arcpy.SelectLayerByAttribute_management(parcel_Layer, "NEW_SELECTION", "\"PropID\" = x")

This is not working. it gives me this error msg:

 Runtime error  Traceback (most recent call last):   File "<string>", line 1, in <module>   File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\management.py", line 6435, in SelectLayerByAttribute     raise e ExecuteError: ERROR 000358: Invalid expression Failed to execute (SelectLayerByAttribute).

What I am doing wrong?


texts: How to make Stored Queries for WFS 2.0.0 in Geoserver

texts:
I am a new user of Geoserver. I have created a WFS and also could connect it with the QGIS WFS  2.0.0 Client plugin available now. But as like the default WFS link in that plugin I am not able to make a query in my service. Also the metadata is not coming for my service. Can anyone help me in this regards?

What do I actually have to do to make the stored queries and also to get metadata?
Where to start doing that?


texts: Image is South up while converting NetCDF layer to GTiff with gdal_translate

texts:
Issue: I'm trying to extract a NetCDF layer to GeoTIFF. The output GeoTIFF image has the South up.
What I need: The output GeoTIFF image with North up.
GDAL version: 3.2.2, released 2021/03/05; OS: Windows 10
The test NetCDF file (~42 mb) is available from DropBox:
Code:
gdal_translate -of GTiff -ot float32 -a_nodata none -strict -q -unscale -co TFW=yes -a_srs EPSG:4326 NETCDF:test.nc:chlor_a_diff test.tif

The output GeoTIFF image has the South up:

So far, I have tried:

creation option: -co GDAL_NETCDF_BOTTOMUP=NO. This creation option is not recognized
Used NumPy (np.rot90, with k=2): it did not help.

I might be missing something more fundamental.
How to get the image with North up?


texts: Changing CRS for basemap from OpenLayers plugin

texts:
I'm working on plotting data from NYC (Queens, specifically) and I would like to overlay my heatmap on a basemap from the OpenLayers plugin to show the location of streets, parks, etc. Because I'm plotting city level data I want to use an appropriate CRS for NYC: NAD83 / New York Long Island (ftUS) EPSG:2263. However, the base map from OpenLayers comes WGS 84 / Pseudo Mercator EPSG: 3857. Also, I have on the fly reprojection enabled.

When I try to change the OpenLayers layer to EPSG: 2263 it plots my data from Queens off of the West coast of Africa. I've tried adjusting the layer crs settings and it seems that I can only get my data to plot correctly on the basemap if my openlayers basemap layer is in EPSG: 3857.

Is it possible to project an OpenLayers layer to a CRS other than 3857 without this issue and if so, how?


texts: Handling mountain shadows when classifying Landsat 8

texts:
I'm trying to classify some area in northern Mongolia using Landsat 8 satellite images from 2013. Since the recordings were made in winter, the sun is very low at the acquisition time. Therefore there are very long and dark shadows from the mountains.

I'm able to identify this shadow using a DEM, as described in this Question: How to Remove Topographic Effects Like Shadows from Landsat Image

How can I perform a supervised classification for the clipped shadowed area? Is it possibile to enhance this areas? I tried some band ratios, but I don't know which is the best for my task.

In this image you can see that there are some vegetated areas in the dark shadows, but it is not possible to classify them.


texts: Arcpy add and calculate fields for features in a list

texts:
I have a script I am working on and am hitting a runtime error. I am attempting to create a list of all of the shapefiles within a root directory, then add a field to each shapefile and calculate it with the shapefile name.   I am able to generate the list (converted from a dictionary) but when I attempt to add the field I encounter the error.

def main():
    try:
        import arcpy, sys, traceback, os, glob
        arcpy.env.overwriteOutput = True
        masterFolder = r"Q:\\GIS\\Field_Data\\MT"
        outputFolder = r"C:\tmp\Shp_merged"

        #collect a list of subfolders in master folder
        subfolderLst = os.listdir(masterFolder)

        #declare a dictionary where a key will be shapefile name
        #... and value a list of pathes to shapefile with this name in all subfolders
        shpDict = {}

        #loop through all subfolders
        for subfolder in subfolderLst:
            #check current subfolder and make a list of pathes to each .shp file
            shpLst = glob.glob(os.path.join(masterFolder,subfolder,'*.shp'))

            #add each shapefile path to dictionary
            for shpPath in shpLst:
                shpName = os.path.basename(shpPath)

                if not shpName in shpDict:
                    shpDict[shpName] = []
                    shpDict[shpName].append(shpPath)

                else:
                    shpDict[shpName].append(shpPath)

        shpDict = shpDict.values()

        print shpDict
        for fc in shpDict:
          arcpy.AddField_management(fc, 'shpname','text')
          arcpy.CalculateField_management(fc, 'shpname', '"'+fc+'"'+ time.strftime('%m_%d_%y'))
    #arcpy.Merge_management(fcs, 'out.shp')

    except:
        print arcpy.GetMessages()
        # Get the traceback object  '"' + wildcard + '"'
        tb = sys.exc_info()[2]
        tbinfo = traceback.format_tb(tb)[0]

        # Concatenate information together concerning the error into a
        #   message string
        pymsg = tbinfo + "\n" + str(sys.exc_type)+ ": " + str(sys.exc_value)

        # Return python error messages for use with a script tool
        arcpy.AddError(pymsg)

        # Print Python error messages for use in Python/PythonWin
        print pymsg

if __name__ == '__main__':
    main()

The list generated looks like this:

[['Q:\\\\GIS\\\\Field_Data\\\\MT\\023N052E\\Points.shp', 'Q:\\\\GIS\\\\Field_Data\\\\MT\\023N053E\\Points.shp', 'Q:\\\\GIS\\\\Field_Data\\\\MT\\024N052E\\Points.shp']]]

and the error:

File "P:\Scripts\scratch\scratch.py", line 45, in main
    arcpy.AddField_management(fc, 'shpname','text')

<type 'exceptions.RuntimeError'>: Object: Error in executing tool

I think it has to do with cursors but being a novice I have no idea how to resolve the problem

Thanks


texts:
I need to calculate a metric for similarity of shapes (Hammond or Hausdorff) between polygons created from a vectorized map and official district polygons.

Now, I've done some searching and found this:

Polygon Shape Similarity using QGIS but it apparently applies to QGIS2, so I don't know how to get the script mentioned there into the toolbox of QGIS3.

It also seems like there is some built-in functionality for Hammond or Hausdorff similarity in QGIS itself, judging from this: https://github.com/qgis/QGIS/tree/master/resources/function_help/json .
But I have no idea how to access this or make it run.

My capabilities to perform programming or developing tasks is very limited, I'd either have to run something within the software, by installing a plugin or copying files to a certain location, or maybe through running a python script from outside, but nothing more unless I have a clear guide, so advice like "You can run this [link] in PostGIS" might be good but not help me (sorry for that).


texts: Nesting iterators using ModelBuilder in ArcMap

texts:
I need to iterate all files in one folder/geodatabase, and perform those operations for all files in a different folder/geodatabase i.e. a matrix, of all entries in one, to all entries in another.

ModelBuilder won't permit using two iterators in one model, and I can't seem to get nesting one iterator in a sub-model to work, or two iterators in two sub-models.

Can someone provide a working example?


texts: Downloading point coordinate data from a website

texts:
I would like to copy the red dots of the following web page:
https://www.inspectorulpadurii.ro/#/
to points in QGIS. The accuracy can be bad. I wanted to just place the points in QGIS with a background map for help.
But is there any way to download all these points?
(The project is about saving virgin forests.)


texts:
After installed QGIS 2.12 I don't have SAGA in Processing Toolbox anymore and I need it bad. Am I blind, can't find it or something went wrong? Any suggestions please?
My OS is Win 7 64 bit.


texts: ArcEngine in WPF -> publish -> creates .msi install file -> installer fails

texts:
I am trying to install my custom ArcEngine application on another machine.  I just uninstalled and reinstalled ArcGis Engine Runtime 10 and authorized it.

"Unable to install or run the application.  The application requires that assembly ESRI.ArcGIS.Framework Version 10.0.0.0 be installed in the GAC first"

The weird thing is that a previous version of my application still runs.

Even the most long-shot type advice is welcome.  I'm hitting a wall.

I removed the reference to ESRI.ArcGIS.Framework, I guess I wasn't using it.  So I thought that would fix the problem, but actually the problem persists.  So I assume that one of the assemblies I use is dependant on Framework.  On the target machine, I went to the assembly directory and there is no Framework assembly on that machine and there is on my machine.  So short term quick fix is to get the assembly on that machine, I am installing devkit on the target to see if that gives me the assembly I need.


texts:
I need to create global reference raster (100km grid, Eckert IV). Later, I will project other various rasters (with different resolutions and projections) into the resolution/projection of the reference grid.

I am using the following code:

First to create a global grid in lat-long 1 degree:

ref_grid<-extent(-180, 180, -90, 90)
ref_grid<-raster(ref_grid)
res(ref_grid)<-1
values(ref_grid)<-1#dummy values
projection(ref_grid)<-CRS("+proj=longlat +datum=WGS84 +no_defs    +ellps=WGS84 +towgs84=0,0,0")

Project the reference layer into an equal area grid of 100km

ref_gridProj<-projectRaster(ref_grid, crs="+proj=eck4 +datum=WGS84", res=100000)
Warning message:
In rgdal::rawTransform(projto_int, projfrom, nrow(xy), xy[, 1],  :
3480 projected point(s) not finite

Later project one of my data rasters into this reference grid:

datalayerProj <- projectRaster(datalayer, ref_gridProj)
Warning message:
In rgdal::rawTransform(projto_int, projfrom, nrow(xy), xy[, 1],  :
3480 projected point(s) not finite

My question is whether I am doing this the most efficient way and also whether I need to worry about these warning messages.


texts:
Does anyone know of good tutorials on how to build raster models (using PostGIS) in GeoDjango? I found Geodjango read raster files from 3 years ago, as well as some sparse documentation for django-raster, but no real examples or tutorials of how to use the module.

Could someone point me to example code for how to implement django-raster?

I've successfully gotten raster data into a GeoDjango model/postgis database, but now am having trouble rendering it. The code I have:

In my models.py:

class testRaster(models.Model):
    raster_id = models.TextField(primary_key=True, default=1)
    name = models.TextField(default='test_layer')
    raster = models.RasterField()

In my html:

new L.tileLayer('/raster/tiles/1/{z}/{x}/{y}.png');

I think my understanding of the layer_id as referenced in the django-raster docs is slightly off.

Also, I am unsure if:

I should even be using tileLayer, since the test raster I am using is quite small and may not even have been tiled.
If I should be passing the tiles link through a Django URL like so:

new L.tileLayer('{% url '/raster/tiles/1/{z}/{x}/{y}.png' %}');


texts: Looking for a live feed of tides in Australia

texts:
I'm looking for a feed of predicted tide heights and times, specifically for NSW, Australia.

Tide sites such as the Bureau of Meteorology only provide static images of predicted tide height. NSW Maritime publish a calendar of tide heights in PDF format.

Ideally I'd like to display a dynamic chart for a series of locations along the coast, in a similar manner to the tide chart on this site.

Does anyone know where I can get a feed of the raw numbers for the coming week?

Thanks


texts: Displaying proportional circles with different colors in QGIS?

texts:
I am new to GIS and just loaded QGIS this morning.

I would like to display data on a map of Arizona showing power generation type (depicted by color; e.g. coal-black, oil/gas - brown...) and capacity (depicted by a proportionally sized circle).  I have the latitude and longitudinal coordinates of the stations.  I'd also like to have a label for each circle showing the name of the station.

Agua Fria   33.5542 -112.2161   oil/gas 614 MW

Coolidge    32.9175 -111.5041   gas 575 MW

Coronado    34.5778 -109.2717   coal    785 MW

Craig   40.4627 -107.5912   coal    250 MW

Desert Basin    32.9042 -111.7889   gas 581 MW

Four Corners    36.69   -108.4814   coal    154 MW

Gila River  32.9761 -112.694    gas 1650    MW

Hayden  40.4856 -107.185    coal    130 MW

Kyrene  33.3544 -111.9364   oil/gas 410 MW

Mesquite    33.344  -112.8535   gas 625 MW

Navajo  36.9047 -111.3886   coal    969 MW

Palo Verde  33.3869 -112.8626   nuclear 688 MW

Santan  33.3333 -111.751    gas 1219    MW

Springerville   34.3186 -109.1636   coal    417 MW

For the map of Arizona, I would like an outline Arizona.  An added feature would be rivers, lakes, and reservoirs.  In QGIS, I looked at the maps under the web tab and really like Stamen Terrain-USA/OSM map.

Is there a way to control the size of the labels and the visibility of certain labels on the baseline map?

I know nothing about QGIS and it's been ages since I've done any programming.


texts: Georeferencing a whole heap of raster imagery

texts:
I have a whole heap of satelite imagery coming in that is not georeferenced. They have a grid overlaid so I do know the bounding box/coordinates of the images (and projection). I just need them georeferenced for use in ArcGIS.

I have manually done a few by referencing to a point dataset I traced over the top - and it works great. I just need to automate this process.

Each image is exactly the same, so the same points overlaid will always match up (even if it doesnt, the application of this does not need to be accurate).

Basically I need to apply the 'link table' information from the example I've already done, and apply it to all the images...

Possible in python/model builder?


texts:
I am trying to import some shapefiles into postgis. Now I have an example database within my download of postgis. In this database I can import the shapefile by changing the DBF file character encoding to LATIN1, from the default UTF-8.

When I create my own database and schema the same pattern does not work. The error message received when I use LATIN1 is:

Importing with configuration: gem_2014, public, geom, C:\Users\Sven\Folderstructuur_GIS_Viewer\Brondata_CBS\gem_2014.shp, mode=c, dump=1, simple=0, geography=0, index=1, shape=1, srid=28992
Shapefile type: Polygon
PostGIS type: MULTIPOLYGON[2]
Failed SQL begins: "SET CLIENT_ENCODING TO UTF8;
SET STANDARD_CONFORMING_STRINGS TO ON;
BEGIN;
CREATE TABLE "public"."gem_2014" (gid serial,
"gm_code" varchar(6),
"gm_naam" varchar(60),
"water" varchar(4),
"aant_inw" numeric(10,0),
"aant_man" numeric(10,0),
"aant_vrouw" nume"
Failed in pgui_exec(): ERROR:  function addgeometrycolumn(unknown, unknown, unknown, unknown, unknown, integer) does not exist
LINE 46: SELECT AddGeometryColumn('public','gem_2014','geom','28992',...
                ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.

Shapefile import failed.

When I use UTF-8 it gives:

Importing with configuration: gem_2014, public, geom, C:\Users\Sven\Folderstructuur_GIS_Viewer\Brondata_CBS\gem_2014.shp, mode=c, dump=1, simple=0, geography=0, index=1, shape=1, srid=28992
Shapefile type: Polygon
PostGIS type: MULTIPOLYGON[2]
Failed SQL begins: "SET CLIENT_ENCODING TO UTF8;
SET STANDARD_CONFORMING_STRINGS TO ON;
BEGIN;
CREATE TABLE "public"."gem_2014" (gid serial,
"gm_code" varchar(6),
"gm_naam" varchar(60),
"water" varchar(4),
"aant_inw" numeric(10,0),
"aant_man" numeric(10,0),
"aant_vrouw" nume"
Failed in pgui_exec(): ERROR:  function addgeometrycolumn(unknown, unknown, unknown, unknown, unknown, integer) does not exist
LINE 46: SELECT AddGeometryColumn('public','gem_2014','geom','28992',...
                ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.

Shapefile import failed.

For anyone trying to reproduce the situation I am using the following dataset: gem_2014 which can be downloaded from the Dutch national bureau for statistics: http://www.cbs.nl/nl-NL/menu/themas/dossiers/nederland-regionaal/publicaties/geografische-data/archief/2015/wijk-en-buurtkaart-2014-art.htm

I have postgres/gis 9.2

I am especially curious to why it does not work. Because I am planning to work more with Postgis I'd like to understand the problem I have now.


texts: Fix broken SDE layers with Arcpy - unexpected error

texts:
We are in the middle of our ArcGIS 9.3 to 10.1 migration project and we have a lot of custom map documents (from different users) that among others, point to our current 9.3 SDE database. After our migration, we'll disable this database to the users so they can't connect anymore. That means they'll end up with broken links in their MXDs (that's basically what we want) and need to repoint to the new SDE.

I prepared a little script that should automatically repoint all broken layers to our new 10.1 SDE (using a connection file), but unfortunately I always run into the same error. I already checked out some other posts on this site (e.g. here and here ), but also on the ESRI forums e.g. here and the ESRI resource center of course. Also tried to repoint to e.g. a file geodatabase, but I'm getting the same error.

However, none of the suggestions apply to my case or worked for me, and after playing around almost the whole day now, I thought maybe someone here knows a solution to that problem. The code itself works fine, it finds the MXDs, spits out the broken layers if I ask it to, it's just the actual repointing command(s) that give me a headache..

Below my (simplified) code and the error I get. Thankful for any tips!

# Import system modules
import sys, string, os, arcpy, datetime, shutil

# source folder(s) to look for MXDs
sourcePath = r"D:\scripts\fixBrokenLinks\test"
# new 10.1 SDE connection
sde = r"\\path to *.sde connection file

# loop through folders and find all MXDs (put into a list)
for root, dirs, files in os.walk(sourcePath):
    for f in files:
        if f.endswith(".mxd"):
            mxdPath = root + "\\" + f
            mxd = arcpy.mapping.MapDocument(mxdPath)
            print "converting " + str(mxdPath)
            brknList = arcpy.mapping.ListBrokenDataSources(mxd)
            for brknItem in brknList:
##                print brknItem
                brknItem.replaceDataSource(sde, "SDE_WORKSPACE")#, brknItem, True)

##                # ALTERNATIVE:
##                find_workspace = brknItem.dataSource
##                replace_workspace_path = r"\\path to *.sde connection file
##                print "find_workspace " + str(find_workspace)
##                brknItem.findAndReplaceWorkspacePath(find_workspace, replace_workspace_path)

            newMXD = mxdPath[:-4] + "_NEW.mxd"
            mxd.saveACopy(newMXD)
            del mxd

And here the error I get:

Traceback (most recent call last):
  File "D:\scripts\fixBrokenLinks\fixBrokenLinks_simple.py", line 19, in <module>
    brknItem.replaceDataSource(sde, "SDE_WORKSPACE")#, brknItem, True)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\utils.py", line 181, in fn_
    return fn(*args, **kw)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\_mapping.py", line 680, in replaceDataSource
    return convertArcObjectToPythonObject(self._arc_object.replaceDataSource(*gp_fixargs((workspace_path, workspace_type, dataset_name, validate), True)))
ValueError: Layer: Unexpected error

I'm using ArcGIS 10.1 SP1 (Build 3143) and Python 2.7.2


texts: Creating "average" line from multiple adjacent polygons

texts:
I have a set of stream channels that are currently being represented by multiple adjacent polygons rather than by actual lines. See the green polygons in the image below. The polygons are essentially a pixelated version of the stream channels.
So I'm looking for a way to basically create a line to replace the polygons. The way I picture it, the lines would go through the vertex of two polygons that are kitty-corner from one another. In theory (I think)  any straight section of channel is just represented by a single elongated pixel, so as long as the lines are taking a sort of diagonal slope through each single pixel and crossing through the vertices, it should end up following the "true" stream channel.
I've tried converting features to lines and it just ends up giving me all the borders of each pixel, rather than a single line going through them.


texts: How to test if a layer exist in a MapBox GL JS map?

texts:
I load a GeoJSON layer in a MapBox GL JS map in this way ....

    map.addSource("route", {
        "type": "geojson",
        "data": {
            "type": "Feature",
            "properties": {},
            "geometry": theCoords
        }
    });

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

... but I need to check if the "route" source and layer exist in my map to remove them before to load the new layer with the same source and name

Does exist in MapBox GL JS any way to test if "route" source and layer exist in my map?

Suggestions / examples / useful links?


texts: QGIS WMS mapserver not showing all layers

texts:
I'm publishing a map via qgis_mapserv.fcgi and the WMS behaves strangely, without giving me any errors or other clues as to what might be wrong.

I am wondering:

is the problem with using GeoJSON as a storage format?
is the problem with my .qgs file?
is the problem with the mapserver?

First, it doesn't advertise all the layers. I have defined a bunch of layers:

but when I access the WMS in QGIS to test it, it only shows me a few of them:

interestingly, one of the layers it does advertise isn't even published in properties:

I am doing nothing inbetween. I have QGIs on my desktop and qgis_mapserv.fcgi on my server. I simply upload the entire directory, .qgs file and all .geojson files and then load it in QGIS as a WMS layer to test.

I do get proper responses to Capabilities, etc. so the server itself appears to be working:

My second problem is that some of the layers, especially polygon layers, don't show up at all, even if I add them. It might be related.


texts: Customizing 'Feature Edit Popup' in Community OpenGeo's GeoExplorer?

texts:
I am running with Boundless Opengeo Suite on Linux based server to create a WebGIS application in forestry management. Everything is working fine (at least to this step that I am being stuck here).

In GeoExplorer, the default popup window for Editing (modify or create) is working well, but I want to customize it to attach another PostgreSQL referenced table within this popup windows, to attach files to the feature (something like the functionality in ArcGIS) using PHP (or Javascript). I have googled for several weeks, but still no luck.

Can anyone show me which file(s) (in GeoExplorer and/or in Geoserver) to be customized, or let me know if other open source alternatives are available on the earth?


texts: XML Parsing Error: error in WFS request with Openlayers?

texts:
I have written a script using OpenLayers library to get some layers from Geoserver and display them on the map. I use both WMS and WFS layers.
Everything worked fine in my localhost (MAMP) in windows OS.
Now I have to migrate in an UBUNTU 12.10 system. Although I get the WMS correctly, I can not make it to get the WFS.

I used firebug and this is the POST request:

<wfs:GetFeature xmlns:wfs="http://www.opengis.net/wfs" service="WFS" version="1.1.0" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
 <wfs:Query typeName="rural:foc_network" srsName="EPSG:4326" xmlns:rural="http://www.opengeospatial.net/rural">
 <ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
 <ogc:BBOX>
 <gml:Envelope xmlns:gml="http://www.opengis.net/gml" srsName="EPSG:4326">
 <gml:lowerCorner>21.831815896484 37.948519811767</gml:lowerCorner>
 <gml:upperCorner>22.346800027344 38.120181188719</gml:upperCorner>
 </gml:Envelope>
 </ogc:BBOX>
 </ogc:Filter>
 </wfs:Query>
</wfs:GetFeature>

And this the POST response (which is actually the content of my cgi file:

 #!/usr/bin/env python

"""This is a blind proxy that we use to get around browser
restrictions that prevent the Javascript from loading pages not on the
same server as the Javascript.  This has several problems: it's less
efficient, it might break some sites, and it's a security risk because
people can use this proxy to browse the web and p

In the XML tag of firegu I get this error:

 XML Parsing Error: not well-formed Location: moz-nullprincipal:{79ed7c81-0daf-4525-a315-808f0894befd} Line Number 1, Column 2:

 #!/usr/bin/env python

Do you have any idea what is going wrong here? Is the path: #!/usr/bin/env python wrongly set? Or what else could be the problem?

And this is my WFS request:

  wfs_layer_komvoi = new OpenLayers.Layer.Vector("Κόμβοι", {
        strategies: [new OpenLayers.Strategy.BBOX()],
        protocol: new OpenLayers.Protocol.WFS({
                version: "1.1.0",
                url: "http://localhost:8080/geoserver/wfs",
                featurePrefix: "rural", //workspace from geoserver
                featureType: "komvoi_real", //layer name from geoserver
                        featureNS : "http://www.opengeospatial.net/rural", //namespace from geoserver
                        styles: "point"
                    })
                })


texts: Converting to Linestring using Dataframe Column

texts:
I'm attempting to convert tuples from a dataframe into a linestring. This is part of my dataframe imported from a csv file.

   Unnamed: 0       name        route                                             decode
0          0        Funshine!  ofosF|mqaShJ@?rLh@d@veCIVd@LbEJfJ^f@lE?Rp@^L~g...  [(-105.279, 40.060), (-105.279, 40.058...

If I manually copy and paste the contents of the decode column into the LineString() condition aka

line = LineString([(-105.279, 40.060), (-105.279, 40.058), (-105.282, 40.058), (-105.282, 40.058), (-105.282, 40.03), (-105.28, 40.036), (-105.283, 40.036), (-105.285, 40.036), (-105.285, 40.036)]

It converts it. But if I try to do the same thing by specifying the column and row, it fails.

line = LineString(df.decode[0])
print(line)

This is the error I receive. It is posted below

Traceback (most recent call last):
  File "shapely\speedups\_speedups.pyx", line 86, in shapely.speedups._speedups.geos_linestring_from_py
AttributeError: 'str' object has no attribute '__array_interface__'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "C:/Users/taylo/PycharmProjects/PermitProj/Polyline Decode.py", line 20, in <module>
    line = LineString(df.decode[1])
  File "C:\Users\taylo\Miniconda3\envs\ExcelPermitExe\lib\site-packages\shapely\geometry\linestring.py", line 48, in __init__
    self._set_coords(coordinates)
  File "C:\Users\taylo\Miniconda3\envs\ExcelPermitExe\lib\site-packages\shapely\geometry\linestring.py", line 97, in _set_coords
    ret = geos_linestring_from_py(coordinates)
  File "shapely\speedups\_speedups.pyx", line 166, in shapely.speedups._speedups.geos_linestring_from_py
AssertionError

How do I write this so that I avoid this error and can convert the tuples to a column in my dataframe?


texts:
I have a three column dataset with LAT, LON and Temperature data I would like to produce a raster image that predicts the temperature of the landscape based on 24 data logger data points. The dataset can be accessed here: DATA
This is what I have attempted this far:
Lets try to interpolate the data onto a raster
library (raster)
library (gstat)
library (sp)
#Temp and XY data
temp<-read.csv ('test_temp.csv')
#create a blank raster to the extent of the system
r<- raster (nrows=300, ncols=100, xmn=-84.95, xmx=-84.936, ymn=45.7, ymx=45.74)

#build a prediction model temp ~ LAT*LON
loc<-temp [,c(3,5,6)]
loc<-na.omit (loc)
#TPS model
tps<-Tps(loc, loc$TemperatureC)
#gstat model
mod1<-gstat (data=temp, formula=TemperatureC ~ 1, locations=loc )
summary (mod1)

r2<-interpolate(r, model=tps)

r2<-interpolate(r, model=mod1)

Ultimatley I want to create a series of rasters of interpolated data to show temperature variation at different times of day. Any ideas on how to do this?


texts: Is there a way to display SVG Symbols in an Attribute Table in a Print Layout?

texts:
I'm working on a trail map that contains both hiking and mountain biking trails. In the Attribute Table for the Trails, I'd like to have a column that displays a bicycle SVG if mountain biking is allowed on the trail.
Is there a way to do this?
I'm also open to the idea that this is not the best way to show this information and am open to other ideas (different trail styling, etc.)


texts: Extracting bounding box from the polygon geometry in OpenLayers


texts: How is raw OSM data processed for openstreetmap.org

texts:
Can anyone provide insight into how OSM data is processed or rendered for www.openstreetmap.org?

A specific example... I extracted data from a recent planet.osm PostGIS dataset for an area in Missouri. The OSM data needs a lot of cleaning before it can be rendered using the correct styles. Many water bodies are stored as line strings that don't close properly, so I have to use FME for snapping and then polygon building so that I can have blue filled rivers / lakes.

If I look at the same data here the water bodies are rendered as expected.

I'm having trouble identifying all the cases where snapping is required (e.g. which 'Natural' types require it and what the tolerance should be). Also I suspect there are many other data issues that I will never see as I am dealing with all of North America.

Does everyone who downloads and uses OSM data go through their own cleanup process? Does anyone know how this cleanup is handled by www.openstreetmap.org? It seems like their process would be the best informed and most tested.

Any insight much appreciated.

EDIT : Here is more information on my workflow

A planet.osm file is downloaded and loaded into PostGIS, using Osmosis, into the pgsql schema. I then extract OSM xml from PostGIS for lots of small areas, again using Osmosis. Each of these small xml files is then converted into Shapefiles using FME and its broad feature categories. It is this stage (OSM xml -> Shp via FME) that I am expecting to convert lines into polygons and perform other cleanup on the data.

These Shapefiles are served up through GeoServer (and cached using GWC).


texts: How to implement 'SELECT column_default' in a QGIS processing script?

texts:
When creating a PostgreSQL serial value using the QGIS DB-Manager, the sequence is named /*schema*/./*table*/_/*column*/_seq.

I've created a processing script that drops the default value/sequence if there is no need for it:

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

#--------- define Interface
##PostgreSQL/PostGIS Tools=group
##Drop default value/sequence=name
##PostgreSQL_Table=vector
##Column=field PostgreSQL_Table

from PyQt4.QtCore import *
from PyQt4.QtGui import *
from PyQt4 import *

from qgis.core import *
from qgis.gui import *
from qgis.utils import *

import psycopg2

#get the parameters for the tpg table into a directory
#get the table
pg_table = processing.getObject(PostgreSQL_Table)
#create empty dictionary for key/value pairs of the tables connection parameters
db_params = {}
db_params['clmn'] = Column
#iterate over connection string
progress.setInfo(20*'-' + '  Connection parameters')
for param in pg_table.dataProvider().dataSourceUri().split(' '):
    key_val = param.split('=')
    progress.setInfo(str(key_val))
    try:
        #set key/value pair
        db_params[key_val[0]] = key_val[1]
    except:
        pass

#generate the sql statement string
#the text in round brackets are the keys from the db_params dictionary created above
#the values belonging to the keys are inserted into the string
progress.setInfo(20*'-' + '  SQL statement')
sql = """ALTER TABLE %(table)s
ALTER COLUMN %(clmn)s DROP DEFAULT;
DROP SEQUENCE IF EXISTS %(table)s_%(clmn)s_seq;
COMMIT;""" % db_params
#remove double quotes
sql = sql.replace('"','')
progress.setInfo(sql)

#make connection string
constr = """dbname=%(dbname)s host=%(host)s port=%(port)s user=%(user)s password=%(password)s""" % db_params
progress.setInfo(20*'-' + '  DB Connection string')
progress.setInfo(constr)
#make db connection
con = psycopg2.connect(constr)
cur = con.cursor()
#execute the above created sql statement
progress.setInfo(20*'-' + '  Executing SQL statement ...')
cur.execute(sql)
progress.setInfo(20*'-' + '  ... done.')

I also want to drop the default value if the sequence is named different than /*schema*/./*table*/_/*column*/_seq? How can I implement the following SQL query in the processing script?

SELECT column_default
FROM information_schema.columns
WHERE (table_schema, table_name, column_name) = ('/*schema*/', '/*table*/', '/*column*/');


texts: Login form for loding layer from MS SQL Server

texts:
I created database and featureclass and table in SQL server. I loaded layer in QGIS, my UI editing and veiw from, designed by qt designer and saved to QGIS project.

Now I need to create login form: after opening QGIS, click 1 button and open this form and type user and password for connect to SQL server and load layer.


texts: Extending polylines via basemap reference in ArcGIS Desktop?

texts:
I have a polyline feature of a pre-existing bus route (let's call it "Route A"). I want to extend the polyline of Route A, using a basemap as reference. I don't want to totally freehand it, my first preference is for an at least semi-programmatic method (the existing polyline features are sourced from GIS shapefiles). Maybe a picture would make more sense:

Now, if that's not possible in ArcGIS Desktop 10.5, as an alternative, consider Route A and the polyline of another certain route (let's call it "Route B"). How would I be able to snip Route A at a certain intersection, snip Route B at this same intersection, then merge the northern part of Route A with the southern part of Route B?

If you're confused, here's a picture:

I tried rasterizing the features, but that splits them into hundreds of constituent polylines, whereas I want each feature to stay as one contiguous polyline.


texts: ArcObjects joins and non-persistent fields after joining to multiple tables

texts:
I have inherited a truly awful and very large ArcObjects project that was originally written in VB and then ported to C#. Apparently the source code I have been given is not the latest version and I am chasing down errors everywhere. The most challenging problem goes like this...

A number of tables are being joined (one after the other) to a Feature Layer via IDisplayRelationshipClass.DisplayRelationshipClass. Let's call these tables A and B (there are more than 2). Each of these tables has fields 1 and 2. The ITable objects are created by using IFeatureWorkspace.OpenTable with a DBF file.

After the first join if I check the fields in my Feature Layer and at the end I can see A.1 and A.2 (field name prefixed with the name of the DBF file the table came from) - Great! However, after the second join A.1 and A.2 are gone, and in their place I have B.1, B.2 - not so great...

Do I need to do anything to persist the additional fields across these joins, so that I end up with A.1, A.2, B.1, B.2 etc. ? If I do, either this behaviour has changed since the code was first written, or the code never worked. Given that the fields are properly prefixed they shouldn't be overwriting each other.

The joins are happening like this:

IMemoryRelationshipClassFactory pMemRelFact = new MemoryRelationshipClassFactoryClass();
IRelationshipClass pRelClass = pMemRelFact.Open("TabletoLayer", (IObjectClass)pTTable, dbfJnField, (IObjectClass)pFeatLayer.FeatureClass, layerJnField, "forward", "backward", esriRelCardinality.esriRelCardinalityOneToMany);

IDisplayRelationshipClass pDispRC = (IDisplayRelationshipClass) pFeatLayer;
pDispRC.DisplayRelationshipClass(pRelClass, esriJoinType.esriLeftInnerJoin);

Any assistance with this issue would be hugely appreciated.

EDIT: If this helps, I have observed something quite strange. After the join, pFeatLayer.FeatureClass.Fields.FieldsCount gives the original number of fields before the join, but ((ITable)pFeatLayer.Fields.FieldsCount gives the expected number (with joined fields). Unfortunately the pFeatLayer variable is then being passed back in to the join function, so it appears that tables are repeatedly being joined against the feature layer's initial state, hence only the last join making a difference.

Any suggestions for why an ITable cast would give a different field count to the FeatureClass? Is this a bug or just silly design??


texts: Can you upload maps from QGIS to ArcGIS Online?

texts:
I am making a storymap in ArcGIS Online where I want to show some maps I made in QGIS.

It seems that you can only show maps that are in the ArcGIS Online gallery.

I tried making a map in ArcGIS Online very briefly, it is very limited and
I can't make the same map unfortunately.

Is there a way to upload your maps from QGIS in my storymap from ArcGIS Online and make them interactive ?


texts: php for "selected feature" - Combination of Openlayers and php-Form

texts:
I use a postgres db, geoserver and openlayers to display a map with a wfs-t Layer that displays the geometry column of my table.
Now I use control select to call a popup which contains a form to edit the other attribute-columns, e.g.:

    //add Popup
var select = new OpenLayers.Control.SelectFeature(garden);
        map.addControl(select);
        select.activate();

garden.events.on({
   featureselected: function(event) {
    var feature = event.feature;
    feature.popup = new OpenLayers.Popup.FramedCloud
    ("pop",
    feature.geometry.getBounds().getCenterLonLat(),
    null,
    <form action="myphp.php" method="post">

 Crop:
<input name="crop" type="text" value=+feature.attributes.crop+>
    <input type="submit" name="write">
    </form>);
 map.addPopup(feature.popup);
},

In my php-script I wrote:

$crop= $_POST['crop'];

and:

$res = pg_query ("UPDATE $tabelle SET crop = '".pg_escape_string ($crop)."' WHERE id = 6;");

This works fine, editing a specific line, in this case the line where the "id=6".

My Question: What is the php-Code for sth. like "where the id= selected feature id".
How can I replace the "WHERE id = 6"-part with a variable that says: always update the attributes of the selected polygon"?

Thanks for your help!


texts: Surveys showing GPS accuracy in different models of smartphones?

texts: Are there any surveys that give accuracy in meters for different smartphones (e.g,  iPhone 6: 8 meters, Samsung Galaxy S8: 7 meters) regardless of what the statistical method is being used?


texts: Use SVG-file as point symbol in Geoserver

texts:
I want to make a map of traffic signs by using .svg image files of the traffic signs as point symbols. I put all the svg-files in the folder 'styles' in the data_dir of Geoserver (on my localhost). But for some reason, only grey squares appear, instead of the images. I am certain the files are placed in the correct folder. This is an excerpt of my SLD:

<?xml version="1.0" encoding="UTF-8"?>
<StyledLayerDescriptor version="1.0.0" xmlns="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc"
  xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
  xsi:schemaLocation="http://www.opengis.net/sld http://schemas.opengis.net/sld/1.0.0/StyledLayerDescriptor.xsd">
    <NamedLayer>
        <Name>Traffic</Name>
        <UserStyle>
            <Name>trafficsign</Name>
            <FeatureTypeStyle>
                <Rule>
                    <Title>08970150-8c69-4c99-98d3-c2bd86035ac6</Title>
                    <ogc:Filter>
                        <ogc:PropertyIsEqualTo>
                            <ogc:PropertyName>img</ogc:PropertyName>
                            <ogc:Literal>955544_aanzicht1.svg</ogc:Literal>
                        </ogc:PropertyIsEqualTo>
                    </ogc:Filter>
                    <PointSymbolizer>
                        <Graphic>
                            <ExternalGraphic>
                                <OnlineResource xlink:type="simple" xlink:href="955544_aanzicht1.svg" />
                                <Format>image/svg+xml</Format>
                            </ExternalGraphic>
                            <Size>
                                <ogc:Literal>60</ogc:Literal>
                            </Size>
                        </Graphic>
                    </PointSymbolizer>
                </Rule>

I must also say, that it did work when I uploaded a png-file to imgbb.com to create a URL, and copy/pasted the URL in the SLD. But it really seems like a lot of work to individually upload images and copy/pasting the URL in the SLD file (especially when dealing with close to a thousand different traffic signs).

Is there any reason why geoserver won't show my svg-files when using the relative path?


texts: Exporting the attribute tables of rasters in ArcPy

texts:
Like the title suggests I have this script that worked once before to export the attribute tables of some rasters as CSV files. However when I run this script again, I get an error:

Runtime error
Traceback (most recent call last):
   File "<string>", line 28, in <module>
NameError: name 'row' is not defined

I pieced this script together from some thread a while back performing a similar task. Can someone please help?

import arcpy, os
ws = r'INPUT_RASTERS_LOCATION'
outPath = r'OUTPUT_RASTERS_LOCATION'
outExt = "OUTPUT_NAME.csv"
arcpy.env.workspace = ws
rasters = arcpy.ListRasters("*")
for raster in rasters:
    rasloc = os.path.join(ws, raster)
    fields = '*'
    try:
        flds = arcpy.ListFields(rasloc)
        header = ','.join([fld.name for fld in flds])
        if len(flds) != 0:
            outCSV = os.path.join(outPath, '{0}{1}'.format(raster, outExt))
            with open(outCSV,'w') as f:
                header += ',RasterName\n'
                f.write(header)
                curs = arcpy.SearchCursor(rasloc)
                for row in curs:
                    lst = [row.getValue(fld.name) for fld in flds]
                    lst.append(raster)
                    line = ','.join(str(a) for a in lst)
                    f.write(line + '\n')
    except Exception as e:
        print "Error processing", raster
        print "Error", e
        print "Is raster not integer or is there no attribute table?"
del row, curs

EDIT: I just realized the data is 32-bit floating point. Is this an issue?


texts: Deleting first two letters from string field usimng Python Parser of ArcGIS Field Calculator?

texts:
I'm looking for a field calculator Python statement to delete only the first two letters of string entries.
For Ex:

!Name!.strip([0:2])

I know this syntax is incorrect.


texts: How do I make line symbols intersect cleanly without showing edges in QGIS?

texts:
I digitised a map showing electric grid lines. I merged a number of lines into a single feature. Now, I'd like these lines to be displayed so that the intersections appear continuous. At the moment, the line edges on intersections are visible and it's not pretty, as visible on the screenshot below:

I could swear I've seen something about this somewhere on the Internet but I can't find it anymore.


texts: Printing a 1" scale from 2.8.2 Print Composer

texts:
I'm working on a small-scale site survey in qGIS, and for engineering purposes I need to print physical copies of the map that can be used with an engineering ruler.

Basically, I need a way to print a scale bar that has intervals that will actually be 1". The map is at 1":500", and I'm using NAD83 CA State Plane V (ft). For the life of me I can't find an option that lets me dictate the size of the scale bar on a paper map. I am printing as a raster, rather than scaling to fit the page, but still no luck.

I'm having the issue that's detailed here but that solution is for ArcMap Layout, and I need to use QGIS Print Composer.


texts:
I am working with a simple map application (basemap and some points/markers) that is using ArcGIS services. The map was developed about 2 years ago but this year (due to security measurements, firewalls etc) it is only accesible by a VM without internet connection (the web application is on an internal server).

Is there any way to run an ArcGIS application (using JS API) without any connection to Esri servers or should I look for the alternatives?

Leaflet looks promising but I would not want to re-engineer the whole application.


texts: Getting name of group that layer is in using PyQGIS?

texts:
I've been using the "Layer List" plugin (https://plugins.qgis.org/plugins/layerList/) that does a great job of creating a text file with all the layer in a qgs project file.

I want to extend it to get the group and subgroup that the layer is in.

The desired result will be

  Provider,Group,Subgroup,Name, Path, Filename.ext

  ogr,Administration, Land,CL_TENURE_POLYGON_DC,D:/GIS/Data/DSE/Vicmap/VMCLTENURE/,CL_TENURE_POLYGON_DC.tab

I assume the change would be

def loadActiveLayers(self):
self.dlgCreate.listWidget.clear() # Clear the list widget before it is filled again
canvas = iface.mapCanvas()
activeLayers = canvas.layers() # Create list with all active layers
layers = u"# Layer List - QGIS Plugin by Klas Karlsson\n" # Create list header and "test" line in a text string
for layer in reversed(activeLayers): # Repeat for all layers in the list
    layerType = layer.type() # Is it Vector or Raster?
    layerSource = layer.publicSource() # path or command for the layer source
    provider = layer.providerType() # Example org, gdal, wms, wfs, postgres, etc
    if layerType == QgsMapLayer.VectorLayer:
    layers = layers + (u"%s,%s,%s\n" % (provider, layer.name(), layerSource))

after layerSource add something for group and subgroup (maybe .layerTreeRoot() from https://qgis.org/api/classQgsLayerTreeGroup.html)

    self.dlgCreate.listWidget.addItem(layer.name()) # Add the layer name to the list

self.dlgCreate.listWidget.addItem(layer.?

 if layerType == QgsMapLayer.RasterLayer:
    layers = layers + (u"%s,%s,%s\n" % (provider, layer.name(), layerSource))
    self.dlgCreate.listWidget.addItem(layer.name()) # Add the layer name to the list
return layers   # Send back all the layers in a text string ready to be saved to a file

Would other changes be required?


texts: Finding the distance between two coordinates (longitudes/latitudes) which are not found as nodes in OSM

texts:
I intend to make my own OSM Router, and I do not want to use third party OSM routers.

WHAT I HAVE DONE: I have downloaded OSM data in XML format and applied different algorithms to find the shortest path between two nodes in map using their coordinates (longitude/latitude). Everything is working fine.

THE PROBLEM: The only problem is OSM XML data has only limited number of coordinates (longitude/latitude). I want to find the distance between two coordinates (longitudes/latitudes) which are not found as nodes in OSM. I could not find any thing related to this on internet.

FOR INSTANCE, I have a random longitude/latitude for point A, and I have a random longitude/latitude for point B. These longitudes and latitudes are not available in OSM XML data. How can I find the distance between these two points?


texts: How to calculate the number of cells with a given value?

texts:
I have a raster file where all cells have values between -3 and 7, no data value is -9999. How can I calculate the number of cells with a specific value, e.g. 6?
Can I use the calculator?


texts: Difference between server and geospatial server

texts:
I have been working with Geoserver and ArcGIS Server but I am still not sure what make these geospatial server handle the spatial data transfer compared to the normal server, like apache & jetty, on which these are installed?

What is that factor that makes Geoserver and ArcServers to be spatially enabled?


texts: Create shapefile connecting points "as the crow flies" in R

texts:
I have a set of locations (with latitude and longitude). I know how to get the distance between these points, but I want to create a line that connects them with the shortest path ("as the crow flies"). I would like the line to be a SpatialLines object.

How do I do this? I tried googling a bit and couldn't find it.


texts: Command works in Console, but not in Python Toolbox

texts:
I want to create a memory layer with points in ArcGIS with Python Toolbox. I found an in_memory function. Here is the question: Why this command works in Python Console in ArcGIS, but not in Python Toolbox?

arcpy.CreateFeatureclass_management('in_memory','memory_layer', "POINT")

Not works, I mean no errors and no effects. Function:

x = arcpy.CreateFeatureclass_management('in_memory','memory_layer', "POINT")
arcpy.AddMessage(str(x))

Returns in_memory\memory_layer

Any ideas?

# -*- coding: utf-8 -*-
import arcpy
import os

class Toolbox(object):
    def __init__(self):
        self.label = "YY"
        self.alias = ""
        self.tools = [XX]

class XX(object):
    def __init__(self):
        self.label = "XX"
        self.description = "YY"
        self.canRunInBackground = False

    def getParameterInfo(self):
        parameters = []

        x = arcpy.Parameter(
            displayName="X",
            name="x",
            datatype="String",
            parameterType="Required",
            direction="Input")
        parameters.append(x)

        y = arcpy.Parameter(
            displayName="Y",
            name="y",
            datatype="String",
            parameterType="Required",
            direction="Input")
        parameters.append(y)

        return parameters

    def isLicensed(self):
        return True

    def updateParameters(self, parameters):
        return

    def updateMessages(self, parameters):
        return

    def execute(self, parameters, messages):
        x = float(parameters[0].valueAsText)
        y = float(parameters[1].valueAsText)
        newpoint = arcpy.PointGeometry(arcpy.Point(x,y), arcpy.SpatialReference(4326))
        a = newpoint.firstPoint.X
        b = newpoint.firstPoint.Y
        arcpy.env.overwriteOutput=True
        arcpy.env.workspace = "in_memory"
        x = arcpy.CreateFeatureclass_management('in_memory','memory_layer', "POINT")
        arcpy.AddMessage(str(x))

Console log:

  Executing: XX 10 10 Start Time: Fri Apr 22 11:37:26 2016 Running
  script XX... in_memory\memory_layer Completed script XX... Succeeded
  at Fri Apr 22 11:37:26 2016 (Elapsed Time: 0,09 seconds)


texts: no viewable feature when creating line with ogr in python

texts:
the python code below produces the expected line shapefile viewable in QGIS

import os
from osgeo import ogr

drv=ogr.GetDriverByName('ESRI Shapefile')
if os.path.exists('test.shp'):
    os.remove('test.shp')
ds=drv.CreateDataSource('test.shp')
lyr=ds.CreateLayer('test',None,ogr.wkbUnknown)
flddef=ogr.FieldDefn('Name',ogr.OFTString)
flddef.SetWidth(32)
lyr.CreateField(flddef)

feat=ogr.Feature(lyr.GetLayerDefn())
feat.SetField('Name','george')
line=ogr.Geometry(ogr.wkbLineString)
line.AddPoint(1116651.439379124, 637392.6969887456)
line.AddPoint(1188804.0108498496, 652655.7409537067)
line.AddPoint(1226730.3625203592, 634155.0816022386)
line.AddPoint(1281307.30760719, 636467.6640211721)
feat.SetGeometry(line)
lyr.CreateFeature(feat)

feat.Destroy()
ds.Destroy()

however, when attempting to use the sqlite driver, the database and layer is created with no errors. QGIS can connect to the database and load the layer - though nothing is displayed.

drv=ogr.GetDriverByName('SQLite')
if os.path.exists('test.sqlite'):
    os.remove('test.sqlite')
ds=drv.CreateDataSource('test.sqlite',options=['SPATIALITE=yes'])
lyr=ds.CreateLayer('test',None,options=['FORMAT=SPATIALITE'])
flddef=ogr.FieldDefn('Name',ogr.OFTString)
flddef.SetWidth(32)
lyr.CreateField(flddef)

feat=ogr.Feature(lyr.GetLayerDefn())
feat.SetField('Name','george')
line=ogr.Geometry(ogr.wkbLineString)
line.AddPoint(1116651.439379124, 637392.6969887456)
line.AddPoint(1188804.0108498496, 652655.7409537067)
line.AddPoint(1226730.3625203592, 634155.0816022386)
line.AddPoint(1281307.30760719, 636467.6640211721)
feat.SetGeometry(line)
lyr.CreateFeature(feat)

feat.Destroy()
ds.Destroy()

on win7, python 2.7.5 (python (x,y) distribution)
BTW - creating a point feature rather than a line works for both cases.
also tried removing the spatialite option from the datasource and layer creation - no change.
am i missing some obvious step?


texts: Which FOSS label engine that supports shrink to fit label in polygon?

texts: I've not found a way to label features using shrink to fit label in polygon, i've tried Qgis, and searched on Geoserver, Mapserver and Mapnik, none of them supports this feature, i'm looking for any FOSS that could label features with that option.


texts: ignore Null Values when iterating through fields with python

texts:
So far I have a code which iterates through Fefature Classes and ignores all fields of type String, Geometry and OID.

I now try to add some lines which tell python to ignore all NULL values when looping through an integer/float/double field. I work with ArcGIS 10.1
Any ideas how to do so? I believe it must be with SearchCursor but I do not know how.

for FCLASS in FCs:

#Skip all Text. OID and Geometry Fields - unfortunately WID, Shape Area etv is still included
AllFields = [f.name for f in arcpy.ListFields(FCLASS)]
LeaveOutField = [f.name for f in arcpy.ListFields(FCLASS,field_type='String')]
LeaveOutField.extend([f.name for f in arcpy.ListFields(FCLASS,field_type='Geometry')])
LeaveOutField.extend([f.name for f in arcpy.ListFields(FCLASS,field_type='OID')])
attrlist = [x for x in AllFields if x not in LeaveOutField]

# text element that has to change (e.g. title...)
title = arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT", "map_title")[0]
for attr in attrlist:
    rows = arcpy.SearchCursor(attr, "row IS NOT NULL")
    for row in rows:
         if attr["*"] is not None:
             #Do Stuff here


texts: High-resolution, hyperspectral satellite imaging?

texts: I am hoping to expand my (biology-focused) research to include satellite data. For this, I'd need satellite imaging that includes multiple values in the visible spectrum (specifically, above 490nm) and extends slightly into the infrared (710nm). I have tried MERIS and Hyperion, but the resolution was insufficient (my area of interest is a lake ~300m x 55m). Does anything like this exist, and is it freely available for research use?


texts:
I am currently labelling a roads map and need help removing the labels for a specific "type" of road.
I apologies for my terminology but I am new to the program.
In my attributes table I have columns "ROAD_NAME" & "ROAD_TYPE" with the latter differentiating between whether the road is a Lane, Highway, Freeway, Street, etc. Some roads are labelled "UNNAMED" in both of these columns. On our other system these roads are present on the map but have no naming.
I currently have applied the rule - "ROAD_NAME || ' ' || ROAD_TYPE" and merged the lines, which has worked well for the other roads but has the UNNAMED roads labelled with "UNNAMED".
I would like the UNNAMED roads still visible on the map but I want the label to disappear.


texts: How to create an offset polyline in arcpy?

texts: I am trying to create an offset polyline using arcpy, identical to the function in ArcObjects (constructoffset). After trying to access ArcObjects in arcpy, which I found difficult for those not well-versed in COM languages, and after much searching for functionality in the regular geoprocessing toolboxes, I discovered that this functionality seems to only exist in ArcObjects. There may be ways to program this functionality in arcpy using trig, but that seems really hard. Essentially, the product I'd like is an offset polyline that preserves all original geometry and curvature. Simply using a SearchCursor on the input and building a new feature class based on its geometry requires some algorithm I don't know. This is because without an algorithm based on slope, etc., a one-size-fits-all addition or subtraction will not preserve an offset shape for all slopes. Any help?


texts:
When I modify the code (which works) from:

vertices = fProject.geometry().asMultiPolygon()
for v in vertices:
    points.append(v)
fet.setGeometry(QgsGeometry.fromMultiPolygonXY(points))

to:

pt = QgsPointXY(92.5,22.5)
points.append(pt)
pt = QgsPointXY(93.5,22.5)
points.append(pt)
pt = QgsPointXY(93.5,23.5)
points.append(pt)
fet.setGeometry(QgsGeometry.fromMultiPolygonXY(points))

throws an error: TypeError: QgsGeometry.fromMultiPolygonXY(): argument 1 has unexpected type 'list'


texts: Renaming print task URL of ArcGIS API for JavaScript?

texts:
I have simple print task using ESRI javscript API PrintParameters, PrintTemplate and printTask. Its working fine. Its return URL like this:

http://sampleserver6.arcgisonline.com/arcgis/rest/directories/arcgisoutput/Utilities/PrintingTools_GPServer/x_____x2Aoqp2xDc-wF9_eV5CU5-w..x_____x_ags_8e12cae972fe45ee848cc74d5e96178b.pdf.

But I need URL like this:

http://sampleserver6.arcgisonline.com/arcgis/rest/directories/arcgisoutput/Utilities/PrintingTools_GPServer/{CurrentDateTime}.pdf

format : CurrentDateTime = 20150806110602 (yyymmddhhmmss)

It is possible?


texts: Is it possible to edit domain coded values in ArcGIS Online Feature Layers?

texts:
I'm getting started on the recently released Citizen Problem Reporter AGOL App. One of the steps is to set up an AGOL hosted featurelayer. I've done that, the app is up and running in a test state, but I'd like to edit some of the coded values to be more in tuned to what our organization handles. Is there a way to do this in AGOL?

The alternative is to download the content to a file geodatabase and host it on my SDE and then serve it using ArcGIS Server, making the necessary edits in Catalog, but that is a lot of work if I can just make the corrections on AGOL...


texts: How to open .mdf file from HAZUS data to check its contents?

texts:
I want to work with HAZUS-MH (US tool for natural hazard modeling) data outside of HAZUS software. I do not have access to ESRI tools and cannot install HAZUS on my machine. So, I decided to use HAZUS data with my own python scripts.
I downloaded state data from here. When I extracted the .exe file I found **.mdf file and a folder with **_HUR_50.bin. Do you know how I can use the file and extract geodata? What is the .mdf file and which open-source tool can open that extension?


texts: When I ran "fit2img" in TIMESAT3.3,I got problematic result image.Contrasted with original image, this result image seemed to be cut into two parts and was badly deformed.I did not know what happened.Have you ever been encountered this kind of problem?I am very puzzled.


texts: Cannot access coords property of LineString in GeoDataFrame

texts:
How can I access the coords property of a LineString in GeoDataFrames? I can access other properties, such as area, length, etc, as documented in the  Shapely documentation. See reproducible example below:

# create example GeoDataFrame with LineString`
line_new = LineString([(0.3, 0.6), (1.5, 1.6), (0.9,2.5)])
df = gpd.GeoDataFrame()
df.loc[0, 'geometry'] = line_new

This works:df.geometry.length

0    2.643715
dtype: float64

But this does not: df.geometry.coords. It returns the following error:

    return object.__getattribute__(self, name)

AttributeError: 'GeoSeries' object has no attribute 'coords

My understanding from the geopandas documentation was that the 'geometry' GeoSeries in a GeoDataFrame was a Shapely object, and thus the respective properties/methods would be available:

  The GeoSeries class implements nearly all of the attributes and methods of Shapely objects.

Is the coords attribute simply one of the few that are not implemented?

Is this the same as what is being asked here?


texts: Creating my own slippy map from GeoTIFF images

texts:
I have a Linux server, about 16 Gb of aerial photography files in GeoTIFF format and orders to set up a slippy map web application. This doesn't need to display anything from other data sources at the moment, but later on we might want to add other layers.

So, the questions: what software packages should I look at, and where could I find a decent tutorial? I know I'll need a client side and a server side. On the client side I am experienced with Leaflet, and using it with Google Maps and OpenStreetMap. On the server side, I have some little experience with PostGIS.


texts: Georeferencing a map image with no precise projection info

texts:
I am trying to georeference a map that has meridians and parallels etc. but no projected grid displayed.  I can get the coordinates of the meridians/parallels fine, so georeferencing should not be an issue. HOWEVER the only projection information is "Universal Transverse Mercator" and no zone information is given.  The map in question spans western Canada so zones 10N-15N.

It is probably projected in zone 11 or 12 but how to be sure for georeferencing?


texts: Is inter-sensor calibration for Landsat 7 and Landsat 8 required?

texts: Can you help me to understand whether I need to calibrate between the ETM and OLI sensors for the surface reflectance tier 1 satellite image obtained from the Google Earth Engine?


texts: Subsequent Python Script Processing Doesn't Adjust Zoom to Selected Features

texts:
I have a script that I load into ArcMap 10.5's Python window and I run it and it works fine. If I go back to the script and change a data value for the a selection to look for, and then load it again, it will run, but won't zoom to my desired location. I want to be able to load the code to search for 'NY' and then load the code again after changing 'NY' to 'OH'. When it tries to zoom to Ohio, it zooms to the layer extent, not the state extent like it does the first time.

import arcpy

#DEFINE VARIABLES
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]
lyr = arcpy.mapping.Layer("NP_Allowances")
lyr2 = arcpy.mapping.Layer("National_NPs")

#CREATE SELECTION LAYER AND WHERE CLAUSE
arcpy.SelectLayerByAttribute_management(lyr,"NEW_SELECTION"," STATE_ABBR = 'OH' ")

#ZOOM TO STATE SELECTION THEN CLEAR SELECTION
df.zoomToSelectedFeatures()
arcpy.SelectLayerByAttribute_management(lyr, "CLEAR_SELECTION")

lyr.definitionQuery = "STATE_ABBR = 'OH'"
lyr2.definitionQuery = "VAR32 = 'OH'"

#UPDATE MAP ELEMENTS
for elm in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT"):
    if elm.name == "Title":
        elm.text = elm.text.replace("Map of NPs", "Practice Locations of NY Nurse Practitioners")

#CLEAR DEFINITION QUERIES AND UPDATE MAP VIEW
arcpy.RefreshActiveView()

#DELETE LOCAL VARIABLES
del mxd, df, lyr, lyr2


texts:
I'm trying to create a OD Cost Matrix layer programmatically.  The problem that I have, however, is that the layer doesn't get added to my Table of Contents.

If I manually run the Network Analyst->Analysis->Make OD Cost Matrix Layer tool a new layer is created in my Table of Contents; however, trying to accomplish the same from arcpy yields no layer.

Likewise, I'm trying to programmatically download, extract, and symbolize OSM (OpenStreetMaps) data.  When utilize this OSM tool programmatically a new layer is not added to my Table of Contents.

Here is the relevant code from my toolbox.

# Load the OpenStreetMap toolbox.
arcpy.ImportToolbox(r"C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcToolbox\Toolboxes\OpenStreetMap Toolbox.tbx")

##
# Get the window extents.
##

# This is the current map, which should be an OSM base map.
curMapDoc = arcpy.mapping.MapDocument("CURRENT")

# Get the datafrom from the map (see the DataFrame object of arcpy).
# The DataFrame object has an "extent" object that has the XMin, XMax, YMin, and YMax.
dataFrame = arcpy.mapping.ListDataFrames(curMapDoc, "Layers")[0]
extent    = dataFrame.extent

messages.addMessage("XMin: {0}, XMax: {1}, YMin: {2}, YMax: {3}".format(extent.XMin, extent.XMax, extent.YMin, extent.YMax))

# Download the data from OSM.
# The "OSMLayer" layer does not get added!
arcpy.DownloadExtractSymbolizeOSMData2_osmtools(extent, True, "OSM_SAC_TEST", "OSMLayer")

##
# Convert the OSM data to a network dataset.
##
arcpy.OSMGPCreateNetworkDataset_osmtools("OSM_SAC_TEST", r"DriveGeneric.xml", r"ND")

###
# Create the OD Cost Matrix layer.
###
#network = "OSM_SAC_TEST_ND"
network = r"C:\Users\avejidah\Documents\ArcGIS\DB\ImportTest.gdb\OSM_SAC_TEST\OSM_SAC_TEST_ND"
# The "OD Cost Matrix" layer does not get added!
arcpy.na.MakeODCostMatrixLayer(network, "OD Cost Matrix", "DriveTime")

I don't get any error or anything, the tool runs and says "Completed."  What am I missing here?


texts: Is export of QGIS layer styles to mapnik possible?

texts:
Is it possible to export layer styles to mapnik xml from QGIS?

As I understand it, it used to be possible once by means of quantumnik plug-in, which is discontinued now.

I see this issue has already been discussed two years ago but to no avail.

Is that it? No other ways to do it? I refuse to believe it is not possible.

It's rather frustrating to have a map ready in QGIS and to be unable to export the styles to mapnik sitting on the same data!

You may call it a "duplicate question", but 4 years have passed since it was asked last time, so maybe something has changed?


texts: Link Russian region codes and iso3166-2 code

texts:
Looking at the Federal subjects of Russia page on wikipedia, i see that Russian regions are listed with a region "Code".
Looking at the ISO 3166-2:RU i see that they also have iso codes.
The names are translated differently, and neither page contains both "code" and iso code.  Looking at regions at random, i see that both the code and iso code are listed (Adygea for example, see the "License plates 01" and "ISO 3166 code  RU-AD").
Do the regions correspond?  and is there a resource to lookup there difference between the two?


texts:
The trigger below don't works and I need that when a new pipeline (adutoras) is inserted or updated, the procedure check what are the municipalities (municipios) that intersects the pipes and put the names of the municipalities inside the field "municipios_abastecidos".

CREATE OR REPLACE FUNCTION topology_fields()
RETURNS "trigger" AS
$$
BEGIN
NEW.municipios_abastecidos:=(
                      SELECT string_agg(municipios.nome,', ')
                      FROM municipios, adutoras
                      WHERE ST_Intersects(adutoras.geom, municipios.geom)
                      AND adutoras.gid = NEW.gid
                    );
RETURN NEW;
END;
$$
LANGUAGE 'plpgsql';

CREATE TRIGGER fill_topology_fields
BEFORE INSERT OR UPDATE ON adutoras
FOR EACH ROW EXECUTE PROCEDURE topology_fields();


texts: Download custom tiles from OSM with specified zoom levels?

texts:
I would like to download openstreetmap tiles which would be used in leaflet.js to implement in a hobby application.
I have been looking for the internet for 2 hours, but cant seem to find a solution for my usecase.

The main thing i would need is to select the areas that i wish to download (a route from A to B and about 20 kms from the route). I does not need to be perfectly 20 kms but i would like to minimise the needed data.
I do not need all the zoom levels like from 0-16. I would only need 3-4 so it also helps to minimise the data. And these zoom levels are not necessarily following each other, so not 1-2-3-4 but something like 1-2-4-8.

What service or an app can do these?


texts:
From the Play Store I have installed QGIS Experimental on my Samsung Galaxy Tab 2. All seems to work fine but I cannot set the coordinate reference system.

From the SETTINGS menu I then select OPTIONS then CRS.

Under section for new layers I select USE A DEFAULT CRS but after clicking search icon NO list of world CRS are displayed for selection.
I want to use EPSG:27700 - OSGB1936.

Have done this without any problem on my desktop QGIS 2.18


texts: ESRI Polygon cutouts (neg. space) become multiple layers in PostGIS

texts:
I have a single MULTIPOLYGON from ESRI ArcGIS with 10,000+ sub-polygons in a Shapefile, some having > 1.4 million vertices. I'm importing the data into PostGIS (ogr2ogr or shp2pgsql, same results) so I can query it for intersection/difference with a LINESTRING. A giant polygon can't be indexed and takes hours to query directly as a single row/column, so dividing it (using default max 255 vertices per polygon) -- lateral join is for unique row IDs:

CREATE TABLE region_subdiv AS
SELECT row_number() over() AS gid, label, segment.geom AS geom
FROM region,
LATERAL (SELECT ST_Multi(ST_SubDivide(geom), 255) AS geom) AS segment;
CREATE INDEX region_subdiv_geom_idx ON region_subdiv USING GIST(geom);

When I export the polygons to a Shapefile to view, the cutouts (small negative spaces inside larger polygons) end up as filled polygons stacked on the larger polygons instead. Also, some small shapes become large rectangles, probably because of errors.

I want to calculate the ST_Difference of a bounding box to find parts of a LineString through the region which are not covered, but the subdivided region multi-polygons have spurious overlapping polygons and cause topology exceptions (GEOSDifference: TopologyException: Input geom 1 is invalid: Self-intersection at or near point ...). I've tried a number of methods to simplify the incoming data, including ST_GeomFromTWKB(ST_AsTWKB(geom, <xy_precision>)) with a few different precision levels. It did not appear to change the behavior.

I've tried the main data cleanup approaches:

UPDATE region_subdiv SET geom=ST_Multi(ST_MakeValid(geom));
UPDATE region_subdiv SET geom=ST_Buffer(geom, 0.0);
UPDATE region_subdiv SET geom=ST_SimplifyPreserveTopology(geom, 0.0001) WHERE NOT ST_IsValid(geom);

These cleanup attempts result in geometries for which ST_IsValid is true, but when I visualize my output data, there are still all kinds of overlapping polygons, up to 3 layers deep, and calculating intersections results in topology exceptions. I use the following materialized view to make an inverse of coverage to find the distance and times traveled outside the polygon regions:

CREATE MATERIALIZED VIEW IF NOT EXISTS public.inverse_coverage AS
    SELECT 1 AS gid, region.name AS region_name,
    ST_Difference(ST_MakeEnvelope(124.41,-16.6717, 125.05, -17.1609, 4326), coverage.geom) as geom
    FROM path, LATERAL (SELECT cover.name, ST_Union(geom) AS geom FROM region_subdiv cover) AS coverage;

I'm using the latest PostgreSQL 10.2 and PostGIS 2.4 on macOS High Sierra.

For more context, here is a question. about the path intersection query I'm using.

Here is a related question about the combination part I'm attempting.