texts: What is the difference between pgr_createtopology and pgr_createverticestable?

texts:
I'm using pgRouting and I'm going through the manual (Release 2.1.0) and in chapter five, it says:

• pgr_createTopology - to create a topology based on the geometry.

• pgr_createVerticesTable - to reconstruct the vertices table based on the source and target information.

When I run them, I get the exact same outputs. Both pgr_createtopology and pgr_createverticestable give me a new layer public.network_vertices_pgr with:

              id|cnt|chk|ein|eou|the_geom

As far as I can tell the two functions are doing the exact same thing, and it would be redundant to run both of them, I just wanted to see if there is a difference that I'm overlooking.


texts: Successfully added a Database Server but nothing inside (ArcCatalog)

texts:
I'm just a newbie in using ArcGIS. I'm trying to create a Web editing application in ArcGIS Server Manager like this post does here, but after added a Database Server i realize that there is nothing inside it, so i can't add User and add the SOC account :sad: Is there any way to fix it ?

I'm using Microsoft SQL Server 2008 R2 and ArcGIS Server Enterprise 10.0


texts:
I have imported arcgis rasters into R.  The rasters are binary maps and the legend in arcgis shows this to be binary as the colors are only 0 and 1.  Now, when I import the raster into R it is no longer binary, as there are 5 numbers in the legend.  Does anyone know why this occurred and have a solution to fix this?

My other question is, when I load a raster in arcgis, how do I effectively change the legend and then save the raster as a different name?

This is m R code for importing the raster and then making a 4 paneled plot.  I have no idea how to change the legend in R.

r5<-raster("C:/users/k/desktop/g/A1b00-59/dican00")
r6<-raster("C:/users/k/desktop/g/A1b50-59/dican50")
r7<-raster("C:/users/k/desktop/g/A1b50-59/dican59")
r8<-raster("C:/users/k/desktop/g/A1b00-59/dicancalc")

par(mfrow=c(2,2))
plot(r5,sub="2000 Distribution",col.sub="blue")
points(Dican,col='blue')
title(main="A1b scenario Panic Grass Distributions",col.main="red")
plot(r6,sub="2050 Distribution",col.sub="blue")
points(Dican,col='blue')
plot(r7,sub="2059 Distribution",col.sub="blue")
points(Dican,col='blue')
plot(r8,sub="Net Habitat Change Panic Grass",col.sub="blue")
points(Dican,col='blue')

Thanks in advance


texts: Mapnik + OpenStreetMap + Shapefile with Raster Images + Mod Tile + PostGIS

texts:
I am new to GIS.

I have a Shapefile (.shp, .shx and .pbf). The .pbf file contains path to PNG files located on my machine only. I was able to save the Shp data to PostGis via Shape2Gis tool.

Now that a tables are created with PostGIS. How can I use this table to show data with Mapnik? I am using Mapnik + OpenStreetMap and Mod_tile.I have been previously using osm2pgsql and the maps are being rendered. I have a shapefile, shx file and pbf file. Not sure how to use it to render maps..


texts: Fuzzy memberships with trapezoid shape in ArcMap

texts:
I would like to combine two raster data (temperature and precipitation) using fuzzy modeling in ArcMap. Each raster is defined by three membership functions with trapezoid shape (e.g. cold, warm, and hot for temperature). I am wondering how I can create the membership functions (Fuzzification) for each raster in ArcMap. The problem is that looking to the ArcGIS documentation, I do not see any option for the case of trapezoid shape. The possible membership functions in ArcGIS are: Gaussian, Large, Linear, MSLarge, MSSmall, Near and Small (see the link below).

How can I solve this? Should I use another package, for instance Python NumPy?

Fuzzy membership:
http://resources.arcgis.com/en/help/main/10.1/index.html#/Fuzzy_Membership/009z000000rn000000/

Python NumPy:
http://resources.arcgis.com/en/help/main/10.1/index.html#//002z00000028000000

EDIT - some more information about my question

I edit my post in order to give better idea about my question.
Here I show an example of how is defined the membership functions for precipitation.

Each trapezoid function is define by its four vertices:
Dry (300, 350, 470, 515),
Medium (470, 515, 540, 568),
Wet (540, 568, 580, 665).
And also for the temperature data:
Cold (6, 7, 8.1, 8.3),
Warm (8.1, 8.3, 8.5, 8.6),
Hot (8.5, 8.6, 8.8, 10).

Then I assign the output values for the total combination of all membership functions of the inputs. Because there are 2 inputs (temperature and precipitation) with 3 membership functions each, there are in total 9 combinations.

IF temperature = ‘cold’ AND precipitation= ‘dry’ THEN output = 0.35
IF temperature = ‘cold’ AND precipitation= ‘medium’ THEN output = 0.15
...
IF temperature = ‘hot’ AND precipitation =’wet’ THEN output = 0.85

After defining the membership functions and the output value for each combination, it should be a way to get the ‘truth value’.
For instance,
when temperature is 7.8 and precipitation is 501, there is an overlapping of the membership functions in precipitation. Therefore there are two ‘truth values’:
IF temperature = ‘cold’ AND precipitation = ‘medium’ THEN output = 0.15
IF temperature = ‘cold’ AND precipitation = ‘dry’ THEN output = 0.35
The ‘truth value’ for the case of ‘dry’ precipitation is 0.317 and for ‘medium’ is 0.683. Thus, in order to get my ‘entire output’, the ‘output’ values are multiplied by their ‘truth values’ and then the two new values are added up.
Case for ‘dry’ precipitation:
0.35*0.317 = 0.111
Case for ‘medium’ precipitation:
0.15*0.683 = 0.102
Sum both cases:
0.111 + 0.102 = 0.213

0.213 is my ‘entire output’, this means the output value that I should have when the temperature is 7.8 and precipitation is 501.


texts: How to adjust placement of piechart symbols to prevent overlapping for a hardcopy?

texts: My map is crowded with pie chart symbols that overlap and obscure each other. Is there a way to draw these symbols so they don't overlap, perhaps drawing them not in their actual location but in a nearby empty space, with a line linking them to their actual location on the map? I need to print this map, so solutions that rely on a varying zoom level won't work. I'm using QGIS 2.18.16.


texts: Qgis2Web - Why are exported labels style not displaying correctly as what has just been set?

texts:
I am using qgis2web with QGIS3.0 to make web map.
The problems I encountered were:

Label styles not exported correctly. (e.g. color changed, background disappeared after export etc.)
Zoom level - I set the Label Rendering based on scale-dependent-visibility. It displayed well in QGIS edit panel. However, it just did not want to work correctly after I exported the map into web map.

Just wondering if there is any bug in QGIS/ qgis2web.js or it was because I did not adjust some settings correctly?
This is an example:
picture(1) what I desire to see (in QGIS edit panel)

picture(2) faulty exported webmap: (result)

Scales, color, text of label background style are all wrong!


texts: How to define the symbology used in QGIS print composer without using separate styles for screen and print

texts:
I use scale-dependent symbols for mapping (rendering Openstreetmap vectors as it happens). For instance paths show up when zoomed closer than 1:50000, but only roads show when zoomed out.

The rules for rendering are defined to make mapping usable on a monitor screen.

When printing my needs are different - the paper is bigger than the screen so it doesn't follow that the same rules should be used to define the rendering.

In another application I use it is possible to instruct the print layout to render the map using the rules normal for a particular scale even though actually physically printing at a different scale. For example I may want to render the map using the rules normal for a 1:50000 map, but because I want to squeeze the map onto a smaller piece of paper I need to actually print it to be actually at 1:100,000.

At the moment every time I change map extent or scale everything remains connected (changing map extent changes scale changes rendering style etc etc)... so the only option I see is to re-write the complex qml files I use to define rendering styles, so that I have one set of styles for printing and one for on-screen use. Am I missing something? This will be an extremely inefficient way to work.


texts: Get distances between boundaries from an SpatialPolygon and raster pixels

texts:
I am trying to calculate the distance from several pixels of a raster to a boundary of some polygons in R as it is showed in the image. I have two problems:

How to correctly define a boundary from a SpatialPolygon object. I currently have two options, but I can not decide for which is the best one. First, I can coerce the SpatialPolygon as a SpatialLines object and then to SpatialPoints. Second, I can raster::rasterize the SpatialPolygon. I suppose each option has its pros and cons.
I only want make the distance calculation from one of the sides of the border. In the example, I only want the distance to the border from the pixels inside Colombia to the border between both countries. I tried to "filter" the border using a latitude parameter but the real border is quite irregular, so using a strict cut-off value can lead to error.

The following code can be used to get the a example shapefile:

#Get the world map and select two countries [Colombia and Venezuela]
library(maptools) #To get the shp data
data(wrld_simpl)
colven <- c("Colombia", "Venezuela")
colven_map <- wrld_simpl[wrld_simpl@data$NAME %in% colven, ]

Raster data can be created with this code:

raster <- raster(nrow=100, ncol=100, crs=proj4string(wrld_simpl))
raster_extent <- extent(colven_map)
raster <- setExtent(raster, raster_extent)
raster[] <- 1:length(raster)
raster_colven <- mask(raster, coven)
raster_colven <- crop(raster, colven_map)
raster_colven <- mask(raster_colven, colven_map)


texts: Seeking comprehensive descriptions of tools in QGIS, SAGA, Orfeo, etc. toolboxes?

texts:
There exists a large number of tools for QGIS (e.g. SAGA, GRASS, ORFEO etc.) - however none of these have sound description of the aim of the tool, the algorithm itself, differences to other similar tools and a hints for the interpretion of the result. How can I judge which tool to use if I do not know this essential background information?

This leads to a great problem that many people use sophisticated analyses but apply wrong (or only standard) parameters or thresholds as they do not know which parameter to adjust. ESRI provides simple explanations for their tools (also they seem not to be written by experts but programmers), taudem has useful description - It would be great to have similar descriptions for the other available opensource tools too!

I wish that the programmers, after working so hard could add half a page to describe the approach, add references, and explain (shortly) their parameters and the results and codes.

Where can I get background information on the different morphology/hydrology tools?


texts:
I have been implementing a workflow previously devised ArcGIS 10.2 in R. The process involves a discretization of an integer raster (the cells are counts) through quantile reclassification.

I've discovered that the ArcGIS quantile classify tool produces completely different results to different R packages and GRASS.

For quantile classification where n=6 produces the following breaks(showing break value start points):
ArcGIS:
0,1,3,6,11,97

R (classInt package):
0,0,0,0,0,0,97 - which effectively yields a 2 class map rather than the desired 6 class map.

My raster data does include a large proportion of a single value (0 in my case) which must be skewing the allocation of equal number of values into quantile bins.

So I guess ArcGIS is doing something more to determine the breaks. Does anyone know what this is and how I may implement it myself?

FYI, the statistics of the raster are:
Min = 0
Max = 97
Mean = 0.601
SD = 3.306

Note I've tried various different R packages and functions including classInt, binr and cut2 (table(cut2(rValues, g=nClassBreaks,levels.mean = TRUE)), and also GRASS (r.reclass). All give same or very similar results to each other. But ArcGIS is doing something wildly different.


texts: Use input parameter as variable in expression

texts:
In my project I have a vector layer and a raster layer. I want to remove all features of the vector layer where the feature's location has a height value of 0 in the same location in the raster layer.
This works fine by using "Select by expression" and using the below expression:
raster_value('layerid' ,1,make_point(x($geometry),y($geometry))) = 0

and then replacing 'layerid' with the ID of the raster layer.
I want to automate a lengthy workflow with the graphical modeler, where this is one of the tasks.
It currently looks like this:

I'm trying to use "Extract by expression", with the same expression as above, and it works as long as the layerid parameter is changed prior to running the model, however I want to use the raster layer input heightmap to get the ID. Like so:

This isn't working. The expression returns null every time. I've tried using eval(@heightmap), as well, but no luck.
As an attempt of a workaround I have tried to use "Load layer into project" algorithm, and then using it's output variable in the expression, but this did the same.
How can I put the layerid of the heightmap input into the expression properly?


texts: ArcGIS 10.1: Do I need to deauthorize licences when reinstalling Licence manager on the same PC?

texts:
I have installed ArcGIS Desktop 10.1 Advanced concurrent use and Licence server administrator on the same PC (Windows 7 Ultimate, 32 bit). I would like to reinstall the whole installation of ArcGIS: the Desktop and the Licence server administrator.

Do I need to deauthorize my Licences before the reinstallation of Licence server administrator, if I am using the same PC?

The deauthorization of Licence server administrator from one PC to another is allowed 4 times per year; what is the rule for the same PC/reinstalled Windows?

I don't want to just try to reinstall Licence server because I am afraid to lose my licences.

I have no longer paid customer support from ESRI to ask them directly.


texts:
I am working with MODIS data (MOD14) and I only need the fire mask, which is the first dataset.

The problem is that MODIS data in the sinusoidal projection and I need to reproject into something else.

My code reads/write out the original (unprojected) data and now I am trying to write out the projected data but am having difficulties.

I am pretty new to python and GDAL so I am not sure what I've been doing wrong.  I got the reprojection parameters from this answer on StackExchange but I am having trouble converting it to the python version.

I am not sure if I am not writing it correctly or perhaps the srcSRS option is incorrect.  I got the srcSRS information from SpatialReference.org here SR-ORG:6974

file_name = 'MOD14.A2018304.0525.006.2018304093408.hdf'

hdf_file = gdal.Open(file_name)
#print(hdf_file.GetMetadata())
#hdf_file.GetSubDatasets()  #fire mask is the first subdataset

band = gdal.Open(hdf_file.GetSubDatasets()[0][0], gdal.GA_ReadOnly)
arr = band_ds.ReadAsArray()  #now numpy array

driver = gdal.GetDriverByName("GTiff")
rows,cols = band_array.shape
dtype=gdal.GDT_Byte  # fire mask is an 8-bit unsigned integer
outname = 'Mod14_org.tif'

#write out unprojected geotiff
dataset_out = driver.Create(outname, cols, rows, 1, dtype)   # fire mask is an 8-bit unsigned integer
dataset_out.GetRasterBand(1).WriteArray(arr)
dataset_out.FlushCache()
dataset_out = None

#MODIS is uses the Sinusoidal projection, now reproject it

ds = gdal.Warp('MOD14_warp.tif', 'Mod14_org.tif', dstSRS='EPSG:4326',
               srcSRS='+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs ', outputType=gdal.GDT_Byte)
ds=None


texts: Using input name of non-native raster format to name output files of iterative model?

texts:
I am using ArcGIS Pro 2.0.1. I'm attempting to create a model to save me some time in the long run, just executing some repetitive raster functions. The model runs properly if I put in a static output name (i.e. one output raster that gets overwritten by the subsequent runs).

The problem, as far as I can tell, is that the rasters I am importing are exported from Oasis Montaj as .gxf files, which are able to be read by ArcGIS but are not listed as a raster format in the raster iterator. Because of this, the name from the raster is being read as "**.gxf", so if I try to set the output as %Name%_ex, the model tries to output the name as "**.gxf_ex" and throws an error. What I'm trying to do is have the final name / raster as "**_ex.gxf" and I'm not entirely certain how to do that (not sure how to add strings together properly in model builder). I also tried naming the output %i% as a test, but it only named the output "0" and wrote over it continuously.


texts: Most cells are "null" during join in QGIS?

texts:
I'm really new to QGIS and GIS software in general. I'm trying make a map travel time across different wards in a city, I have the average travel time and range between peak traffic and lowest traffic. Anyway, I joined the shapefile containing the polygons of each ward to a CSV file containing the above parameters (I joined via name ward 1, ward 2.....so on). I added 2 fields, one fore range and the other for average time (both integer values), i cannot find a pattern as to why only three of these rows merged properly.

The tables have joined, however all but 3 of the cells have a value of null and I cannot figure out why. The cells that do have a clue correlate with the CSV file. I'm properly stuck.

Lastly, as I said I'm really new to QGIS and I don't know any code, so I'm sorry I cannot provide that.


texts: Symbol Selector has Save As button grayed out?

texts: I upgraded to 10.3 from 10.2.  I tried to SAVE a symbol but the SAVE AS button is grayed out, rendering saving as not possible.  I've tried rebooting, re-opening the mxd... What solutions are available?


texts:
I have tried both Google and Yahoo's Geocoding services to get the Montreal neighborhood/district from GPS coordinates but while they are both able retrieve the correct (address, city, state, country) they won't retrieve the neighborhood.

Which Geocoding service would offer this service? Otherwise, where can I find the neighborhood data to do it myself?


texts: Zonal computations: computing raster pixel counts within polygons using GDAL versus rasterstats

texts:
I'm computing some zonal statistics, and was comparing the results obtained using the GDAL cookbook example code here with the rasterstats zonal_stats function.

I want to compute the mean value per polygon, as well as the number of raster pixels that contributed to that calculation. I think rasterstats gives you the number of contributing pixels with their 'count' statistic.

I'm not sure how to get that with the GDAL code, though. I first thought that's what their 'numpoints' variable was counting, but their numpoints corresponds to the number of points within a ring, and I think that means the number of vertices of the ring (please correct me if I'm wrong).

Can anyone explain how to compute the number of raster pixels contributing to the GDAL cookbook calculation? Counting the number of pixels of their 'zoneraster' array variable doesn't come close to agreeing with the rasterstats count, either.


texts:
I've been experiencing slower than usual performance (noticeably slow time to draw as well as query, compared to similar feature classes) with a specific feature class in my organizations enterprise geodatabase.  The feature class consists of tax parcel polygons in the Local Government schema (vector_data.DBO.ParcelPublishing/vector_data.DBO.TaxParcel).  It is a non-versioned feature class with archiving enabled.

I've tried the following to speed it up (with no success):

We reconcile, post, and compress our GDB nightly
Feature class has an attribute index as well as a spatial index
I have run "rebuild indexes" as well as "analyze datasets" on the feature class with just about every combination of options possible

Any suggestions as to what I should try next?

Our enterprise geodatabase is on version 10.3 on top of SQL Server 2008.


texts: Leaflet-WFST and custom markers and markercluster

texts:
I want to add the ability to edit a point layer in my leaflet-map.

When I displayed objects in this way, using custom markers and the markercluster plug-in did not have any problems:

var treeIconList = L.icon({
        iconUrl: '/site/images/list16-3.png'    });
var treeIconHvoya = L.icon({
        iconUrl: '/site/images/hvoya16.png'     });

var owsrootUrl = 'http://localhost:8080/geoserver/ows';

var defaultParameters = {
service : 'WFS',
version : '2.0.0',
request : 'GetFeature',
typeName : 'postgis01:flora',
outputFormat : 'text/javascript',
format_options : 'callback:getJson',
SrsName : 'EPSG:4326'
};

var parameters = L.Util.extend(defaultParameters);
var URL = owsrootUrl + L.Util.getParamString(parameters);

var WFSLayer = null;

var markers = L.markerClusterGroup();

var ajax = $.ajax({
url : URL,
dataType : 'jsonp',
jsonpCallback : 'getJson',
success : function (response) {
WFSLayer = L.geoJson(response, {
pointToLayer: function (feature, latlng) {

    switch (feature.properties.type) {
        case 'F': return L.marker(latlng, {icon: treeIconList});
        case 'T': return L.marker(latlng, {icon: treeIconHvoya});
                }
        },
onEachFeature: function (feature, layer) {
var input = L.DomUtil.create('input', 'my-input');
    input.value = feature.properties.type;
    L.DomEvent.addListener(input, 'change', function () {
      feature.properties.type = input.value;
    });
    layer.bindPopup(input);
}
}).addTo(markers.addTo(mymap));
}
});

After using the https://github.com/Flexberry/Leaflet-WFST plugin, no work custom markers, markercluster:

var treeIconList = L.icon({
        iconUrl: '/site/images/list16-3.png'   });
var treeIconHvoya = L.icon({
        iconUrl: '/site/images/hvoya16.png'    });

var markers = L.markerClusterGroup();
var points = new L.WFS({
   url: 'http://localhost:8080/geoserver/ows',
    typeNS: 'postgis01',
    typeName: 'flora',
    crs: L.CRS.EPSG3395,
    geometryField: 'shape',

    pointToLayer: function (feature, latlng) {

    switch (feature.properties.type) {
        case 'F': return L.marker(latlng, {icon: treeIconList});
        case 'T': return L.marker(latlng, {icon: treeIconHvoya});
                }
        }
  }).addTo(markers.addTo(mymap));

What do you think is the problem? Or how else can you solve the task of editing a point layer in a leaflet or at least editing the attributes of individual objects of a point layer.


texts: convert a raster (.tif format) into a mesh

texts: I want to use the cloud to mesh registration tools in Cloud Compare to align a raster (.tif format) to a point cloud (.las format) for surface differencing. what is the best raster to mesh conversion to use?


texts: Calculating Detours based on given blocked roads with Leaflet?

texts:
Is it possible to use a given route of "blocked roads" and calculate a route from A to B, considering that the given route is off limits using Leaflet?

Where I'm from, there is an annual bicycle sports event and many roads will be blocked for the racers to drive on. So I want to calculate detours based on blocked roads.

I want to calculate routes with Leaflets Routing and if the calculated route crosses a blocked road I want to calculate a detour if possible.


texts: spectral profile for multilayer/time series

texts:
I am missing a spectral profile in QGIS which enables me to see the values for all bands of a given pixel. (absolutely necessary for multispectral and multitemporal/time series work)

It seems, there has been such a plug in called "earthbeat". http://hub.qgis.org/projects/earthbeat/wiki

This is exactly what I need, but the plugin seems to not work for QGIS 2.2 . It should be in the official repository. Am I getting anything wrong? Thanks for help and clarification.


texts: Guideline/scheme for naming properties in GeoJSON

texts:
Are there any guidelines/requirements for naming properties in GeoJSON? Should I add a vendor-prefix for my properties which only make in use with a specific application?

I already read the specification however I didn't found any further information on this matter.

Examples:

{
    "foobar-id": 21
}

or

{
    "com.foobar": {
        "id": 21
    }
}


texts: What could cause text under 255 characters to import as BLOB format?

texts:
I'm importing a table from an .accdb file in Access 2016 into ArcMap 10.6 via an OLE DB connection and one of my nine text fields is importing in BLOB format. It and all the other fields are formatted as short text with a 255 character limit. The longest record in the problem field is ~ 100 characters. I thought ArcMap only used BLOB format for text if it was over 255 characters.

Below is a summary of the structure of the table in the Access database and how it is formatted in ArcMap after it is added by field type and number of fields:

Any suggestions for what I am missing and how I might get ArcMap to recognize this field as a string rather than a BLOB?


texts: How to calculate te area between contour lines on a DEM with QGIS 1.8?

texts: Does anyone know how to calculate the area in km² for each altitudinal band (between 100 and 200, 200 and 300...) based on a DEM? I added the contour lines with a distance of 100 metres but could not find a way to modify them in order to allow area calculation. I am working with qgis 1.8.0. DEM is taken from ASTER.


texts: AMD processor has compatibility with QGIS and Grass GIS?

texts: I would like to ask if the processors AMD are compatible with Q GIS, GRASS GIS as well as Rstudio. Until now I had Intel i5 and I would like to go to AMD Ryzen™ 7


texts:
I have a series of tile sources equating to different map providers at specific resolutions and I would like to set an extent so that when the relevant layer is selected, the map pans and zooms to the extent of the tiles.

Following this example I've set the extent for each source. I'm also able to get the extent of the selected layer using this code in my layerswitcher control:

curExtent = map.getView().calculateExtent(map.getSize());
            if (extent != null) {
                extent = layers[s].getExtent();
            }
            else {
                extent = map.getView().calculateExtent(map.getSize());
            }
            layerZoom = layers[s].getProperties().minZoom;
            map.getView().fit(extent, map.getSize());
            layers[s].setVisible(true);

I'm able to get the resolution/zoom at this point but I'm not sure how to set the view to the minResolution/zoom property of each source. Can anyone please help?

Further to this I'm now able to get the map to pan and zoom but the layer tiles don't appear. I thought it was a projection issue but I'm in the right location and zoom. Do I need to delay loading of the tiles maybe?

Here is the amended code:

var curExtent = map.getView().calculateExtent(map.getSize());
if (title == selectedValue) {
            if (extent != []) {
                extent = layers[s].getExtent();
            }
            else {
                extent = curExtent;
            }
            console.log(extent);
            layerZoom = layers[s].getProperties().minZoom;
            var sourceCentre = ol.extent.getCenter(extent);
            console.log(sourceCentre);
            console.log(layerZoom);
            console.log(curExtent);
            console.log(extent);
            map.getView().setCenter(sourceCentre);
            map.getView().setZoom(layerZoom);
            layers[s].setVisible(true);
        }
        else {
            layers[s].setVisible(false);
        }

I've set up a simple fiddle which now seems to work although I'm struggling to get the extent of the other standard map layers and the functionality seems to break after the initial successful action.


texts: How to load color attribute from shapefile in MapServer 5.0.x?

texts:
according to the docs of mapserver 6.0.x it is possible to load the color of a polygon shape from the shapefile. This is a great feature cos otherwise we would have to define a class for each different colored shape.

But when I create the following mapfile and read it with mapserver 5.0.4, I get a

[MapServer Error]: getSymbol(): Parsing error near (SHAPE_COLOR):(line 1192)

The class in my mapfile is:

CLASS #default coloring
    OUTLINECOLOR 0 255 0
    COLOR [SHAPE_COLOR]
END

What is the proper syntax?


texts: Programmatically modifying ArcGIS Measurement object's units property

texts:
I am using ArcGIS's JS API for Measurement tool to measure and display Area, Distance, and Location points, but I cannot seem to be able to change the units programmatically. There's a method for GetUnits(), but no corresponding SetUnits(), and modifying the values in the DOM doesn't work because the dropdown that their widget creates isn't a dropdown at all, its a div and some span elements modified to look like a dropdown.

I just started a new job with a mapping company that uses this, and I'm rather new to GIS, so I apologize if there is an obvious solution that I'm not aware of.


texts: REST GeoServer services for anonymous users gives me 404 status

texts:
I am trying to access to a REST service by an anonymous user on GeoServer 2.9.4. I modified my resp.properties file to allow anonymous role free access to GET:

/**;GET=IS_AUTHENTICATED_ANONYMOUSLY
/**;POST,DELETE,PUT=ADMIN;

Once I did this, the service gives me 404 status instead 401, and browser mode gives me an empty result, for example for service /geoserver/rest/layers/layer_name, I get:
No such layer: layer_name
But for service: /geoserver/rest/styles everything works great.
What am I missing?


texts: Administrative level and coordinates from Yahoo GeoPlanet WOEID

texts:
The Yahoo GeoPlanet data can be freely downloaded from:

http://developer.yahoo.com/geo/geoplanet/data/.

If I understand it correctly, each WOEID corresponds to a unique location on earth, with associated hierarchy indicators (continent, country, state, etc.). The complementary GeoPlanet Explorer site  allows you to search by location name, and provides a nice hierarchy of parent, siblings and children. For example, this link provides an overview for the "Boudh" district in the state of Orissa in India, this link lists all its siblings, and this all the towns in the district. For every location (parent, siblings, children), the site queries and returns hierarchy (country, state, district, town) and coordinates (lat/long).

My question is that is possible to get the same information by running bulk queries on the WOEIDs in the GeoPlanet data, rather than using the explorer website and querying by name and then drilling down?


texts: Finding Raster Overlap within a Directory

texts:
Here is the scenario:

I have a directory with 100+ map projected, single band, rasters that weighs in a 40GB+.  I have used ArcMap's built-in mosaicing tools to create a raster mosaic.  I am using a custom built suite of image processing tools to perform some stretching and model out some camera error.  The next step is to work on contrast matching the overlapping rasters.

I have used the color matching within Arc and am not happy with it.  I want to implement my own modified algorithm using either poisson image gradient reconstruction or using a multiresolution spline technique.  This will be done in python using GDAL, Numpy, and Scipy.

Does anyone have any suggestions on techniques for finding raster overlap within a directory?  The rasters are all variable size and orientation.  They are all map projected to the same projection.

Ideas so far that I have discounted:

Manually create a multiband .vrt, perform some array math, and return an array containing just the overlap.  This unfortunately, fails to be automated as rasters need to be manually selected.
The images are all relatively close in contrast, so I considered performing some image matching to find areas of overlap.  I discounted this based upon the size of the directory to be parsed.
Convert the image to a binary representation (1=data, 0=nodata), then polygonize, and check for overlap with the polygons.

Item three looks to have the most promise, but I am hoping that someone else will have a novel, fast, method.


texts: How to select classes from a GeoTIFF world map and compute the area by regions defined in a shapefile

texts:
I have two files:

fileLUse.tif defines land cover classes in GeoTiff format (from https://www.esa-landcover-cci.org/?q=node/164)
fileRegions.shp defines administrative borders and it is a shapefile (from https://gadm.org)

Both files are EPSG 4326, that if I understood it well it is WTS84 datum with coordinates system (no projected) and cover the whole World.
I need to select the land use classes from fileLUse whose value is in a certain range ([40-130]) and compute the area and area ratio of this aggregate for each region in fileRegions.
My "output" would then be a CSV with region name (from fileRegions), total area, share of aggregate land use.
My main problem is that the files are not projected, so I don't really know how to give them a unit, as all "equal area" projections I saw are specific for a certain area of the World.
Which approach should I take ? I could use QGIS, GRASS, R, Python... (I'm on Linux)
[EDIT]
What I have done so far:

Computed 0/1 binary data of Land Covers in the interested range:

gdal_calc.py -A fileLUse.tif --outfile=fileNatVegAreas.tif --calc="logical_and(A>=40,A<=140)"

I used rasterstats to get the mean/count of each "polygon"

from rasterstats import zonal_stats
stats = zonal_stats("fileRegions.shp", "fileNatVegAreas.tif")

What I am left with:

stats is a vector of dictionaries with the various stats, e.g.

>>> stats[0]
{'min': 0.0, 'max': 1.0, 'mean': 0.5842002754571467, 'count': 739861}

How do I link back my stats to the field "RegionName" in my shapefile, and possibly export everything as a CSV ?


texts: Batch convert geotiff to jpg but ensure nodata is white

texts:
I've been using gdal_translate to convert geotiffs of rotated paper maps to jpg+wld, but nodata sections would always turn up black (i.e., where the original data has been rotated).

So based on other answers on the web I've tried using a combination of gdalwarp and gdal_translate as such:

for %i in (*.tif) do gdalwarp -srcnodata "0,0,0" -dstnodata "255,255,255" %i %i_fix.tif
for %i in (*_fix.tif) do gdal_translate -of JPEG -co "worldfile=yes" %i %~ni.jpg

I've tried many parameters for -srcnodata and -dstnodata, with/without the alpha band, "0", 0, but still get large black edges/spaces in the output jpg.

How can I ensure that nodata ends up white in the resultant jpg instead of black? Using only gdal, preferably. I would like to avoid python scripts for now if possible.


texts:
I am working with QGis, and I have a layer that has some double features. I want to delete one of them, but not both of them.
How can I do that easily?


texts: Using %i%".csv" in ArcPy to name output CSV file?

texts:
I am using arcpy to describe datasets and output the results into a csv file. I am doing this via a script in arctoolbox and would like to save the file in the same folder as the file being analysed. I don´t want to overwrite any existing files, and thus would like to simply add a number to the file name.

In the FieldCalculator we can use %i% or %n% to have a variable suffix - for example if my output txt file is called output.txt and a second output.txt is created, this is automatically changed to output01.txt

I want to do the same with ModelBuilder but if I can´t use %i%.

import arcpy, os
from arcpy import env
import os.path

env.overwriteOutput = True

#Enter FeatureClass path in modelbuilder dialog
fc = arcpy.GetParameterAsText(0)
ischecked  = arcpy.GetParameterAsText(1)

#Get FeatureClass name
fcname = os.path.basename(fc)

#define Output Directory
outputTxtFilePath = "PATH"

#Define filename
outputTxtFilename = "\\"+fcname+ %i% +".txt"

#full path of file
outputTxtFile = outputTxtFilePath+outputTxtFilename

# Create a list of fields using the ListFields function
fields = arcpy.ListFields(fc)

#open csv file
csv = open(outputTxtFile,'w')

#write describe results to txt file
csv.write(fc)
for field in fields:

    #create string
    txtStr = "Feldname: {0} \n Type: {1} \n Länge {2}".format(field.name, field.type, field.length) + '\n'

    csv.write(txtStr)

csv.close()

#Open the txt file if is checked
if str(ischecked) == 'true':

    os.system('notepad '+outputTxtFile)

else:
    exit

it works ok if I remove %i%, but then each file gets overwritten.


texts: Use of QGIS/GRASS/SAGA procedures on Geometry Generator

texts:
Im new at QGIS, started 2 weeks ago. I'm trying to automatize the display of a map and get rid of some "duplicated" layers and all the extra work. I need to know if SAGA's dissolve can be used in the Geometry Generator, using a personalized function. I tried something like the following

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

@qgsfunction(args='auto', group='Custom')
def DisolverSAGA(layer, field, feature, parent):
    asd = processing.runalg('saga:polygondissolvebyattribute', layer, field)
    return asd

Honestly, the fact it didn't work doesn't surprise me, I'm too new at this and I'm not a programmer or anything.

Already tried using buffer(collect($geometry,"REGLABEL"),0) in the geometry generator expression, but I'm using that symbology with a transparency and I noticed that the borders of the resulting polygons have different scales of color, as if that generated one different output for each original geometry and they were stacking, thats why im trying to emulate the process used to generate the existing shapefiles


texts: Do you know a free program that transform *.ssf to *.shp files?


texts: PostgreSQL function return list of geometry

texts:
I have a query that return result as expected - list of points.
I would like to create a function with two parameters:

input polygon geometry
N as number of points to create

WITH
polygon AS
    (select geometry from property_gis where id=1247675),
intervals AS
    (SELECT generate_series (0, 9) as steps)
SELECT ST_LineInterpolatePoint(ST_ExteriorRing(((ST_Dump(polygon.geometry)).geom)), steps/(select count(steps)::float-1 from intervals)) as polygon_points
FROM
    polygon,intervals;

I've created next function:

CREATE OR REPLACE FUNCTION get_points(polygon geometry)
RETURNS geometry[] AS
$$
DECLARE
    geom_out geometry[];
BEGIN
    WITH
        intervals AS (SELECT generate_series (0, 9) as steps)
        SELECT
            ST_LineInterpolatePoint(ST_ExteriorRing(((ST_Dump(polygon)).geom)),steps/(select
count(steps)::float-1 from intervals))
        INTO
            geom_out
        FROM
            intervals;

   RETURN geom_out;
END;
$$
LANGUAGE plpgsql ;

but when I try to call it, like this:

with
    poly AS (select geometry from property_gis where id=1247675)
select  get_points(poly.geometry)
from poly;

I get the following error:

  [22P02] ERROR: malformed array literal: "01010000000A17E1637B1054C0815E5CA915553A40" Detail: Array value must start with "{" or dimension information. Where: PL/pgSQL function get_points(geometry) line 5 at SQL statement


texts: GeoServer 2.5.2 + OpenLayers 3 not working

texts:
I have installed GeoServer 2.5 and by default it had OpenLayers 2.8, now I have replaced the OpenLayers folder with OpenLayers 3 (ol3) contents and it is not working.

The examples in the ol3 folder are not working.

So please suggest me how should I use OpenLayers 3?

Which folders have to be changed?


texts: Exporting buildings and their types with Overpass

texts:
I am trying to extract all buildings in one file within a certain city (Radebeul, Sachsen, Germany) from Overpass Turbo.

My query looks as follows:

[out:json];

// fetch area “Radebeul” to search in
{{geocodeArea:Radebeul}}->.searchArea;

// gather results
(
  // query part for: “admin_level=*”
  node["admin_level"](area.searchArea);
  way["admin_level"](area.searchArea);
  relation["admin_level"](area.searchArea);
);
  // query part for: “building=*”
(
  way["building"](area.searchArea);
  relation["building"](area.searchArea);
);

// print results
out body;
>;
out skel qt;

The output:

Everything looks okay, however, I can not understand why do my previously queried administrative districts disappear in the final output? I just want to visually check if buildings are within city borders.

Moreover, I was wondering what is the most acceptable file format for QGIS import? Is it GeoJSON? While I am not sure if Overpass Turbo will provide me with buildings types as well.

References:

Finding residential buildings with overpass-turbo?
Overpass API: Select all building polygons (excluding POIs and lines) within certain area
exporting to shapefiles


texts: How to show size/length of sketch while creating new feature?

texts:
While adding a new feature to a feature layer I like to show the size (area)/length of the sketch so that the user has a better feedback while placing the vertexes.
So I hook up MapControl's OnMouseMove and get the IEngineEditSketch. But that one only gives me the size/length without the current mouse position beeing a vertex. So I tried to add that point to a clone but the size of the clone still is the size without the actual mouse position.

So how can I show the user the size/length the feature would have, if he/she would add a vertex at the current mouse position?

private void mapControl_OnMouseMove(object sender, IMapControlEvents2_OnMouseMoveEvent e)
{
  if (m_engineEditor.EditState == esriEngineEditState.esriEngineStateEditing)
  {
    IEngineEditSketch editSketch = m_engineEditor as IEngineEditSketch;
    if (editSketch != null)
    {
      switch (editSketch.GeometryType)
      {
        case esriGeometryType.esriGeometryPoint:
          // TODO
          break;
        case esriGeometryType.esriGeometryPolygon:
          IPolygon polygon = editSketch.Geometry as IPolygon;
          if (polygon != null)
          {
            IClone clone = polygon as IClone;
            if (clone != null)
            {
              IPolygon copyPolygon = clone.Clone() as IPolygon;
              IPointCollection copyPoints = copyPolygon as IPointCollection;
              if (copyPolygon != null && copyPoints != null)
              {
                IPoint thisPoint = new PointClass() { X = e.mapX, Y = e.mapY };
                copyPoints.AddPoint(thisPoint);
                copyPolygon.Close();
                IArea area = copyPolygon as IArea;
                if (area != null)
                {
                  StatusText = String.Format("Area: {0}", area.Area);
                }
              }
            }
          }
          break;
        case esriGeometryType.esriGeometryPolyline:
          // TODO
          break;
        default:
          break;
      }
    }
  }
}


texts: Setting legend in QGIS to only show values that actually occur?

texts:
I am using QGIS 2.18.3 (I know, it's an old version now) on Windows 7. I have a raster of forest groups for Oregon, saved as the .img format and is rendered as Paletted. It was clipped from a raster of forest groups for the entire U.S. The legend in the Layers Panel for the OR layer has all of the forest group pixel values (a few hundred different values), even though there are only about 25 different pixel values in the OR layer.

How do I get the legend to not show all the values that don't actually occur in the layer?

Related to this, if I save as .tiff, the layer's legend only has min and max values, which leaves out a bunch of the pixel values.

EDIT: I've added a screenshot of QGIS to show the legend. All the values that appear in the screenshot legend aren't used in the raster. Those, plus the other unused values, are the ones I want to remove.


texts: How do I configure the "Processing R Provider"-Plugin in QGIS 3.10 on Linux Mint

texts:
I am using the "Processing R Provider"-Plugin in QGIS 3.10 running on Linux Mint 18. I have R installed, but QGIS is not able to find it:

R is installed in "/user/bin":

I have already tried to adjust the "R Folder" link in the Processing Options, but to no avail:

According to the QGIS manual QGIS should know where R is, if I can run "R" by typing R in the terminal and that works.

But QGIS still prompts "This algorithm cannot be run".

If i run an empty R-Script I get:

Any ideas?


texts:
I want to use mapscript with rubyonrails

how should i start doing this?
can anyone point to a good tutorial on mapscript?


texts: Check two geometries intersects in PostGIS?

texts:
I've created a function called 'ST_LineFishnet' which makes a grid network of line string geometry. I wanted to split a polygon shapefile with my grid. I've just used the first grid as "LINESTRING(0 0,0 10,10 10,10 0,0 0)" to split my shapefile according to this snippet code, but the problem is it always return geometry even two geometries don't have any intersection.

How can I check if my grid intersects with polygon before split function?

SELECT ST_AsEWKT(wkt) As result FROM
(SELECT ST_AsText((ST_Dump(ST_Split(circle, line))).geom) As wkt
FROM (SELECT
     (SELECT ST_AsEWKT(geom) FROM (
     SELECT geom FROM ST_LineFishnet(4,6,10,10) As cells WHERE col=1 AND row=1) As koo) As line,
     (SELECT ST_AsEWKT(geo) FROM (SELECT geom As geo FROM block WHERE gid = 1) As k) As circle) As foo) As m;


texts: Python 2.6 Module for GDAL 1.9.2, ImportError

texts:
I work on an internal network (i.e. I can't get the latest 1.11 stuff) and have GDAL 1.9.2 set up on a Red Hat Linux Enterprise Server 6. My Linux box has python 2.6.6 however the Python binding I received in the swig directory appears to be for either Python 2.7.5 or Python 3.2.

I ran the setup.py but when I try to import OSGEO or GDAL I get this message.

Traceback (most recent call last):
    File "<stdin>:, line 1, in <module>
    File "osgeo/__init__.py", line 21 in <module>
        _gdal = swig_import_helper()
    File "osgeo/__init__.py", line 13, in swig_import_helper
        import _gdal
ImportError: No module named _gdal

Anyone got a link to the proper Python module or know of how I can get this to work. Got a lot of raster data that needs ripping apart.


texts: Openlayers - Cloned features get style from original feature

texts:
I am having a problem with the styling of a WFS vector layer from GeoServer after cloning features from it and adding them to  a new feature.

First I create a new empty vector layer and add it to the map object:

var polyAdHoc = new OpenLayers.Layer.Vector("Search Results");
map.addLayer(polyAdHoc);

Then I use a selectFeature control to get a feature from another WFS layer:

var mLayer = map.getLayersBy('name', 'Test Stands 4thQ 2012' );
var name = mLayer[0].name;
var layerFeatureArray = mLayer[0].features;
ctrlSearch.select(layerFeatureArray[x]);

Then I clone the feature from the WFS vector layer:

var cloneFeature = layerFeatureArray[x].clone();

Then I add the cloned feature to the vector layer :

polyAdHoc.addFeatures(cloneFeature);

Then I set the style of the feature in the new vector layer so that it will appear hilighted against the original layer:

cloneFeature.style.strokeColor = '#FF0000';
cloneFeature.style.strokeWidth = 3;
cloneFeature.style.fillOpacity = 0.2;
cloneFeature.style.fillColor = '#FF0000';

The new layer refreshes and the map looks fine until I turn the original WFS layer off and then on again. When it comes on again the features from the WFS that were cloned now share the same style as the cloned features and I want them to stay the same as they were originally.

I do not want the style to apply to the original WFS shape, only the features in the new layer.

Even more frustrating, if I try to reset the original WFS feature style it changes the new vector layers style so I end up with the same problem of both layers having the same style again.

I am using OpenLayers 2.12 and GeoServer 2.1.4.

Any suggestions would be appreciated. Thanks!

@CaptDragon - EDIT:

Here is the complete code I am using now. Everything here works great except for the fact that this changes the style of the features in the original WFS layer.

I have modified this slightly from my original post. Instead of cloning the features I am just doing a straight-up addFeature to the new vector layer.

//FIRST I GET THE WFS LAYER I WANT TO COPY FEATURES FROM
var mLayer = map.getLayersBy('name', 'Test Stands 4thQ 2012' );
/THEN I CREATE A VARIABLE FOR THE ARRAY OF FEATURES IN THE LAYER
var layerFeatureArray = mLayer[0].features; //the array of features in a layer
//THEN I CREATE THE NEW EMPTY VECTOR LAYER
var polyAdHoc = new OpenLayers.Layer.Vector("Search Results");
//I ADD THE NEW EMPTY VECTOR LAYER TO THE MAP OBJECT
map.addLayer(polyAdHoc);

//THE selectItems IS AN ARRAY OF VALUES I RETRIEVE FROM A DATABASE
//I LOOP THROUGH THESE ITEMS
for (a = 0; a < selectItems.length; a++) {
//THEN I LOOP THROUGH THE FEATURES IN THE WFS LAYER
    for (x = 0; x < layerFeatureArray.length; x++) {
//IF THE UID FROM THE WFS FEATURE MATCHES THE VALUE FROM THE selectItems VAR
        if (layerFeatureArray[x].attributes['UID'] == selectItems[a]) {
            //I ASSIGN A VARIABLE TO THE FEATURE
            var searchFeature = layerFeatureArray[x];
            //AND THEN ADD THE FEATURE TO MY NEW VECTOR LAYER
            polyAdHoc.addFeatures(searchFeature);
            //AS I ADD EACH FEATURE I GIVE IT STYLING
            searchFeature.style.strokeColor = '#FF0000';
            searchFeature.style.strokeWidth = 3;
            searchFeature.style.fillColor = '#FF0000';
            searchFeature.style.fillOpacity = 0.2;
            break;
        }//END if
    } //END feature array FOR LOOP
} //END selectItems FOR

//ZOOM TO THE EXTENT OF THE NEW VECTOR LAYER
map.zoomToExtent(polyAdHoc.getDataExtent());


texts: Properly writing floats to feature properties of a new shapefile in Fiona

texts:
The following piece of code returns a shp following the driver and crs of a model .shp.

import numpy as np #just to test I correctly handle data types
import fiona
from shapely.geometry import Polygon, mapping
with fiona.open('model.shp') as source:
    source_driver = source.driver
    source_crs = source.crs
    #source_schema = source.schema
    source_schema = {'geometry': 'Polygon',
                 'properties': {'field_name': 'float'}} #Beware: field names cannot be longer than 10 characters!
        #writing a new file
            with fiona.open('output.shp',
                            'w',
                             driver=source_driver,
                             crs=source_crs,
                             schema=source_schema) as collection:
                rec = {}
                rec['geometry'] = mapping(Polygon([(0,0),(1,0),(1,1),(0,0)]))
                rec['properties'] = {'field_name': float(np.float(0.3))}
                collection.write(rec)

The written record can be checked in this way:

collection = fiona.open('output.shp')
rec = collection.next()
print rec

Which returns: {'geometry': {'type': 'Polygon', 'coordinates': [[(0.0, 0.0), (1.0, 1.0), (1.0, 0.0), (0.0, 0.0)]]}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'field_name', 0.3)])}

I followed @genes' hints from this post and @sgillies' recommendations below.


texts: Is it recommended to have Esri ArcMap MXD files under version control e.g. stored on SVN or TFS?


texts:
Utilizing arcpy, I need to reclassify some rasters.  Once they're reclassified, I need to multiply all the rasters together and then clip them to a set feature. Can anyone assist?

import arcpy
from arcpy import env
from arcpy.sa import *

env.workspace = "C:\\Data"

outReclass1 = Reclassify("reclassstream", "Value",
                     RemapValue([1,0], [0,1]))
outReclass1.save("C:\\Data\\reclassStream1")

outReclass2 = Reclassify("pondbuffer", "Value",
                     RemapValue([1,0], [0,1]))
outReclass2.save("C:\\Data\\pondBuffer2")

outReclass4 = Reclassify("soilanalysis", "Value",
                     RemapValue([1,3,1], [3,6,2], [6,9,3], [9,12,4], [12,15,5]))
outReclass4.save("C:\\Data\\soilAnalysis4")

Raster1: reclassstream
Description: Streams with 150ft buffer
Reclassify Rules;
Current Value = 1 -----> New Value = 0
Current Value = NoData -----> New Value = 1

Raster2: pondbuffer
Description: Ponds with 200ft buffer
Reclassify Rules;
Current Value = 1 -----> New Value = 0
Current Value = NoData -----> New Value = 1

Raster3: PropertyUse
Description: Property Use Raster
Reclassify Rules;
No reclassification necessary

Raster4: soilanalysis
Description: Soil Analysis Raster
Reclassify Rules;
Current Value = 1-3 -----> New Value = 1
Current Value = 3-6 -----> New Value = 2
Current Value = 6-9 -----> New Value = 3
Current Value = 9-12 -----> New Value = 4
Current Value = 12-15 -----> New Value = 5


texts: problem to use openstreetmap data file with pgrouting2.0 on PostgreSql9.3

texts:
I'm a newcomer on pgrouting, and I want to test its functionalities.

I've installed and configured my computer with Postgresql9.3, pgrouting2, and PostGIS.

I've downloaded from openstreetmap the .osm file of the zone I need (4.0706 North, 9.7321 East, 3.9972 South, 9.6737 West).

My problem is that it comes with 4 layers: points, lines, multipolygons, and multilinestrings.

I don't know which layer I must use to test the pgrouting routage function.
I've taken the lines layer, inserted it on postgres, and added 2 columns (source:integer, target:integer, cost:double precision).

Now when I use

select pgr_createTopology('littoral_lines', 0.0001, 'geom', 'osm_id');

where my table is littoral_lines. I get the following errors:

NOTICE:  PROCESSING:
NOTICE:  pgr_createTopology('littoral_lines',0.0001,'geom','osm_id','source','target','true')
NOTICE:  Performing checks, pelase wait .....
NOTICE:  Creating Topology, Please wait...
NOTICE:  l'opérateur n'existe pas : character varying = integer
NOTICE:  UPDATE public.littoral_lines SET source = 1,target = 2 WHERE osm_id =  4173643
Durée totale d'exécution de la requête :1751 ms.
1 ligne récupérée.

I don't know what's wrong, if I'm using the correct table, or if I miss something. There aren't any "roads" or "ways" tables. Please, help !!!


texts: QGIS and Shapely : Shell is not a LinearRing

texts:
When trying to create a polygon with Shapely with QGIS python console, i'm getting the error :

ValueError: Null geometry supports no operations

And it prints Shell is not a LinearRing in the console.

 layer = QgsVectorLayer(myshapefile, "testlayer_shp", "ogr")
 geom = next(layer.getFeatures()).geometry()
 multiPolygon = geom.get()
 pointList = multiPolygon.coordinateSequence()[0][0]
 poly = shapely.geometry.Polygon([[p.x(), p.y()] for p in pointList])

I first tried doing this with geopandas, but I get the same error and I assume it comes from the Shapely package. I installed Shapely with the Christoph Gohlke's python packages.

I'm running it with QGIS python36 on windows10 and I don't have any error when running it directly from this python, only when trying to launch it in QGIS python console.

My pip info:

alabaster==0.7.10
astroid==1.6.3
Babel==2.5.3
backcall==0.1.0
bleach==2.1.3
certifi==2018.1.18
chardet==3.0.4
click==6.7
click-plugins==1.0.3
cligj==0.4.0
cloudpickle==0.5.2
colorama==0.3.9
cycler==0.10.0
decorator==4.2.1
descartes==1.1.0
Django==2.0.5
docutils==0.14
entrypoints==0.2.3
ExifRead==2.1.2
fiona==1.7.11
Flask==1.0.2
future==0.16.0
GDAL==2.2.4
geographiclib==1.49
geopandas==0.3.0
geopy==1.13.0
geos==0.2.1
html5lib==1.0.1
httplib2==0.10.2
idna==2.6
image==1.5.20
imagesize==1.0.0
ipykernel==4.8.2
ipython==6.3.1
ipython-genutils==0.2.0
isort==4.3.4
itsdangerous==0.24
jedi==0.12.0
Jinja2==2.10
jsonschema==2.6.0
jupyter-client==5.2.3
jupyter-core==4.4.0
kiwisolver==1.0.1
lazy-object-proxy==1.3.1
lxml==4.2.1
MarkupSafe==0.23
matplotlib==2.2.2
mccabe==0.6.1
mistune==0.8.3
mock==2.0.0
munch==2.3.1
nbconvert==5.3.1
nbformat==4.4.0
networkx==2.1
nose2==0.6.5
numpy==1.14.2
numpydoc==0.8.0
osmnx==0.7.4
OWSLib==0.14.0
packaging==17.1
pandas==0.22.0
pandocfilters==1.4.2
parso==0.2.0
pbr==1.10.0
pickleshare==0.7.4
Pillow==5.1.0
plotly==2.5.1
prompt-toolkit==1.0.15
psutil==5.4.5
psycopg2==2.7.3.2
PuLP==1.6.8
pycodestyle==2.4.0
pyflakes==1.6.0
Pygments==2.2.0
pylint==1.8.4
pyodbc==4.0.9
pyparsing==2.1.10
pypiwin32==220
pyproj==1.9.5.1
PyQt5==5.9.2
python-dateutil==2.6.0
pytz==2018.3
PyYAML==3.12
pyzmq==17.0.0
QtAwesome==0.4.4
qtconsole==4.3.1
QtPy==1.4.0
requests==2.18.4
rope==0.10.7
Rtree==0.8.3
scipy==0.19.0rc2
Shapely==1.6.4.post1
simplegeneric==0.8.1
simplejson==3.10.0
sip==4.19.8
six==1.11.0
snowballstemmer==1.2.1
Sphinx==1.7.4
sphinxcontrib-websupport==1.0.1
spyder==3.2.8
testpath==0.3.1
tornado==5.0.2
traitlets==4.3.2
urllib3==1.22
wcwidth==0.1.7
webencodings==0.5.1
Werkzeug==0.14.1
wrapt==1.10.11
xlrd==1.0.0
xlwt==1.2.0


texts: How to loop through an array of point features?

texts:
I have a set of point features and I wish to find the points that have common coordinates. Now I changed the code as shwon bewlo by using the IPoint and the compare method through a function. I have 2 loops with cursors. The first one should read the first point while the second loop the second point etc. How can I force the cursor of the second loop to go to the second point? I addedd another NextFeature but it did not work.The result is that all the points are the same between the 2 loops.
Thanks
Demetris

Set pFeature = pFCursor.NextFeature
Do Until pFeature Is Nothing
Set pPointA = pFeature.Shape

Set pFeature2 = pFCursor2.NextFeature
Do Until pFeature2 Is Nothing
'Set pFeature2 = pFCursor2.NextFeature

Set pPointB = pFeature2.Shape
CommonPoints = ComparePoints(pPointA, pPointB)
If CommonPoints = True Then
MsgBox "Yes"
End If

Set pFeature2 = pFCursor2.NextFeature
Loop

Function ComparePoints(pointA As IPoint, pointB As IPoint) As Boolean

If pointA.Compare(pointB) Then

    ComparePoints = True

Else

    ComparePoints = False

End If

End Function


texts:
I'm trying to build a python plugin that utilizes the pyqgis class "QgsAttributeDialog". I've only found one website (aside from the pyqgis doc page) that shows an example of how to implement the script correctly, and this is the example I based my own script on.

https://qgis.org/api/classQgsAttributeDialog.html

(example script)
https://webgeodatavore.github.io/pyqgis-samples/gui-group/QgsAttributeDialog.html

from PyQt4.QtCore import QSettings, QTranslator, qVersion, QCoreApplication, QObject, SIGNAL, QTimer
from PyQt4.QtGui import QAction, QIcon, QMessageBox, QColor, QDialogButtonBox
from qgis.gui import QgsMapTool, QgsMapToolEmitPoint, QgsMapTip, QgsAttributeDialog
from qgis.core import (
QGis,
QgsCoordinateTransform,
QgsCoordinateReferenceSystem,
QgsRectangle,
QgsMapLayer,
QgsExpression,
QgsFeatureRequest,
QgsFeature,
QgsVectorFileWriter,
QgsGeometry,
QgsVectorLayer,
QgsRasterLayer,
QgsMapLayerRegistry,
QgsPoint)

def run(self):
    global layers
    layers = self.iface.legendInterface().layers()
    self.pointEmitter = QgsMapToolEmitPoint(self.iface.mapCanvas())
    QObject.connect(self.pointEmitter, SIGNAL("canvasClicked(const QgsPoint, Qt::MouseButton)"), self.selectNow)
    self.iface.mapCanvas().setMapTool(self.pointEmitter)

def selectNow(self, point, button):
    selectedLayer.removeSelection()
    self.iface.mapCanvas().refresh()
    radius = .00001
    rect = QgsRectangle(point.x() - radius, point.y() - radius,
    point.x() + radius, point.y() + radius)
    selectedLayer.select(rect, False)
    selected_features = selectedLayer.selectedFeatures()
    features = selected_features[0]
    attribute_dialog = QgsAttributeDialog(selectedLayer, features, True)
    attribute_dialog_from_iface = self.iface.getFeatureForm(selectedLayer, features)
    attribute_dialog.show()
    attribute_dialog_from_iface.show()

The issue is in the two "attribute_dialog.show()" and "attribute_dialog_from_iface.show()" lines in the script above (my plugin script). Currently the script allows me to click on a layer of road segments and selects whatever segment I've clicked on. That works just fine, and I have the selectNow function clear the selection each time I click (so it will only select one feature at a time). I then took the example script and changed the code slightly so that it takes the selected layer and QgsFeature's attributes and inputs it into the Attribute dialog box. However when I integrated the QgsAttributeDialog lines from the example, it opens the dialog box (with attribute data visible) for 1-2 seconds and then immediately causes my session in QGIS to crash.

I cannot find anything wrong with the input (the layer of interest, QgsFeature object, and boolean value). I was hoping that someone on stack exchange has used this class before and can point me in the right direction/assist with troubleshooting.

My QGIS version is 2.14.2.


texts: Function for Labelling Icons per category in Leaflet

texts:
I have geographic coordinates that I display to a Leaflet Map. Each feature is ether in the Category (2018, Never Covered or Other periods).

I would like to create a simple function to have different colors for my categories.

Whats the most simple way to do that?

Ive joined my blockbuilder: http://blockbuilder.org/renauld94/5678103bddc3c05c5410e1579b62944a

Can i call a function directly after adding my markers ?

var mymap = L.map('mapid').setView([10.762622, 106.660172], 13);
var pharmaciesMarkers = new L.FeatureGroup();

L.tileLayer('http://{s}.tile.osm.org/{z}/{x}/{y}.png', {
        attribution: '&copy; OpenStreetMap contributors'
    }).addTo(mymap);

Any examples will be usefull.


texts: Is there a way to disable right-click in QGIS when drawing

texts:
I would like to know if I can disable the rightclick in QGIS. The issue is that I draw polygons that are kind of big, let's say 4000 polygon points. And as you know rightclick brakes the drawing process and you can save or delete what you draw.

But drawing 4k points is terrible for a hand. So it often happens to me that I accidentially click the right mouse button though I was not finished with drawing.

Is there a way to disable it. I mean instead of right click maybe something in the menu to stop drawing?


texts: How to export the Chart as .csv file automatically without clicking the top right button of the chart displayed in the console of Earth Engine?

texts:
The way that I export the time series of a point of an image collection as .csv file is to click the top right button of the chart displayed in the console of earth engine, which generated from the time series. The following code shows how I extract the time series of a point by its coordinate and generate a chart to display in the console. Then I click the top right button of the chart to export as .csv file.

However, I have 30 points of interest, in my way, I have to input the coordinate of each point manually and run the code every time, and then click the top right button of the chart of each point to export .csv file again and again until I finish the export of time series for all points. So is there any convenient way to accomplish this? Maybe I should input the coordinates of points in a loop, and export the .csv file also in a loop, but I don't know how to do it.

// the coordinate of a point, I enter it manually every time.
    var lng = -92.56894444;
    var lat = 37.21702778;
    var r = ee.Geometry.Point(lng, lat); // r is the point of interest
    // Load Landsat 7 imagery and filter it
     var collection = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')
        .filterDate('2011-01-01','2011-12-31') // the period of interest.
        .filterBounds(r);
Map.setCenter(lng, lat, 13);

//plot a chart of the imagecollection B4 band; then I click the top right button of the chart in the console to export the time series as a .csv file.
print(ui.Chart.image.series(collection.select('B4'), r, ee.Reducer.mean(), 30))


texts: Using Mapsui, how do I create a new transformation class for coordinate system read from shapefile prj file?

texts:
How do I create a new transformation class compatible with MapSui.Projections.ITransformation interface using ProjNet4GeoAPI that reads the source coordindate system from a prj file.

From Mapsui source code, there is a MinimalTransformation that implements ITransformation interface to convert between SphericalMercator and WGS84.

From Mapsui documentation :
The out of the box Mapsui support for projections is limited. The MinimalProjection class only projects between SphericalMercator (EPSG:3857) and WGS84 (EPSG:4326). It is however possible to create your own Transformation. You need to implement the ITransformation interface. Within this implementation you need to use some other projection library. A recommended one is ProjNet4GeoAPI.

I can create a working the transformation class with ProjNet4GeoAPI but it implements GeoAPI.CoordinateSystems.Transformations.ICoordinateTransformation not Mapsui.Projection.ITransformation

            // (FROM SOURCE) prj name: NAD_1983_StatePlane_Massachusetts_Mainland_FIPS_2001"
            ICoordinateSystemFactory csFac = new ProjNet.CoordinateSystems.CoordinateSystemFactory();
            string file = @"C:\DRC_Data\Arcview\USA\Townships\NYTOWNS_POLY.prj";
            string wkt= System.IO.File.ReadAllText(file);
            var csFrom = csFac.CreateFromWkt(wkt);

            //(TO) Prj name: "WGS 84 / Pseudo-Mercator"
            file = @"C:\DRC_Data\Arcview\3857.prj";
            wkt = System.IO.File.ReadAllText(file);
            ICoordinateSystem csTo = csFac.CreateFromWkt(wkt);

            //Step 2) Create transformation class.
            CoordinateTransformationFactory ctFac = new CoordinateTransformationFactory();

            //To 3857
            //var is ICoordinateTransformation
            ICoordinateTransformation ct = ctFac.CreateFromCoordinateSystems(csFrom, ProjNet.CoordinateSystems.ProjectedCoordinateSystem.WebMercator);

How do I use the ICoordinateTransformation class with Mapsui?
Do I create a projection class like SphericalMercator in Mapsui.Projection
(see code below)?

From Mapsui.Projection:

public class MinimalTransformation : ITransformation
    {
        private readonly IDictionary<string, Func<double, double, Point>> _toLonLat = new Dictionary<string, Func<double, double, Point>>();
        private readonly IDictionary<string, Func<double, double, Point>> _fromLonLat = new Dictionary<string, Func<double, double, Point>>();

        public MinimalTransformation()
        {
            _toLonLat["EPSG:4326"] = (x, y) => new Point(x, y);
            _fromLonLat["EPSG:4326"] = (x, y) => new Point(x, y);
            _toLonLat["EPSG:3857"] = SphericalMercato.ToLonLat;
            _fromLonLat["EPSG:3857"] = SphericalMercator.FromLonLat;
        }


texts: "..does not contain a valid set of GRASS 7 modules.", how to fix?

texts:
Using 3.10.8-A Coruña, installation with GRASS tools.
I am trying to run Grass' r.water.outlet tool and getting this error:
This algorithm cannot be run :-(
The specified GRASS 7 folder "C:\PROGRA~1\QGIS3~1.10\bin\bin" does not contain
a valid set of GRASS 7 modules. Please, go to the Processing settings dialog,
and check that the GRASS 7 folder is correctly configured

It is not clear how to correctly configure the Grass 7 folder. My installation only has a GRASS 78 folder.  I have tried the fix Micha explains here:
GRASS not working with QGIS 3.7
I go to Settings > Options > Advanced >  Processing > configuration
The instructions say to go to the GRASS7 folder but my installation does not have a folder, it has GRASS7 log commands, log console, use vexternal and help; non of these allow me to edit the file path, apart from the Grass_help_path.
Can someone update this help for this 3.10.8-A Coruña version of QGIS please?


texts: How to create a QGIS Action which loads a raster?

texts:
I am following an article here:Actions HowTo , specifically the piece about an action to load a raster:

Python type
    Add Raster Layer
            qgis.utils.iface.addRasterLayer('PATH_FILE','NAME_FILE_IN_QGIS')

Can anyone advise what i need to do to open a raster based on a value in the vector table (essentially I have a grid with a field containing a field of Grid_Ref for all items where a raster image exists).

I have been fiddling with the example above and got as far as an Action of Type Python, Named OpenRaster, with the following action:

qgis.utils.iface.addRasterLayer('E:\Plot Sheet Devt\1974-1984\1250k\' + "Grid_Ref" + '.TIF',"Grid_Ref")

The Python error that I am getting on running this is

  File "", line 1
      qgis.utils.iface.addRasterLayer('E:\Plot Sheet Devt\1974-1984\1250k\' + "Grid_Ref" + '.TIF',"Grid_Ref")
                                                                                ^

  SyntaxError: EOL while scanning string literal

How should I go about forming the open 'PATH_FILE' and 'NAME_FILE_IN_QGIS' so that it picks the information from the Grid_Ref field in the underlying vector data?

As always, thanks for your time

Mark


texts: How to reproject raster image from WGS84/Pseudo-Mercator to ECEF to ENU in python?

texts:
I am new to GIS and I am writing a program that should work with a cartesian coordinate system, ENU. Therefore my plan is convert all point coordinates from WGS84/Pseudo-Mercator to ECEF and then to ENU.

I can get an OSM map in WGS84/Pseudo-Mercator as geotiff via an API.

How can I reproject the raster map image from WGS84/Pseudo-Mercator to ECEF and then to ENU using python?

Or can I neglect the error on a square map of 50kmˆ2 and avoid reprojecting the raster image.


texts: How raster data can be converted to KML? I have searched but could not find any answer. Are there any tools in Arc toolbox for such process ?


texts: Channel swap when reading created image with Gdal

texts:
I use a function on python (3.5) to create GTIFF with GDAL:

def array2raster2(fname,   matriz, geot, proj):
    drv = gdal.GetDriverByName("GTiff")
    dst_ds = drv.Create(fname, matriz.shape[1], matriz.shape[0], 3,  gdal.GDT_UInt16)
    dst_ds.SetGeoTransform(geot)
    dst_ds.SetProjection(proj)
    dst_ds.GetRasterBand(1).WriteArray(matriz[:, :, 0])  # write r-band to the    raster
    dst_ds.GetRasterBand(2).WriteArray(matriz[:, :, 1])  # write g-band to the raster
    dst_ds.GetRasterBand(3).WriteArray(matriz[:, :, 2])# write g-band to the raster
    dst_ds.FlushCache()
    dst_ds=None

and works fine, but when I open again the image with gdal.Open() and gdal.ReadAsArray() and check the shape, instead of i.e. (100, 101, 3), it shows (3, 100, 101), can anyone tell me why?

Also when I open the GTIFF on ENVI the orders of band is BGR, when I am saving the image on RGB.


texts: Writing Python code block in field calculator?

texts:
I am a novice in Python as well has having minimal experience with ArcMap. Currently I am using a U.S. Census TIGER file, mapping roads in Harris County, TX. I want to add a new field called SPEED, which will return a speed associated with the MTFCC string (MTFCC is the column which classifies road type. e.g. S1400, S1200, etc.).

I opened my attributes table, added a new field called SPEED, clicked on the Field Calculator of the speed column, and inserted my own Python code. I made sure to click Python (not VB), and have my return value be a number. My code is as follows...

def SpeedCalc(MTFCC):
 if MTFCC == "S1400":
      return 25

 elif MTFCC == "S1200":
      return 45

 elif MTFCC == "S1100":
      return 65

 elif MTFCC == "S1630":
      return 25

 elif MTFCC == "S1730":
      return 15

 else:
      return 9999

I made sure below where 'SPEED = ' is to have "SpeedCalc(!MTFCC!) written. No errors were thrown, but all my columns were listed as 0. I do not have administrative rights to my working computer or else I would insert my code in to Python first and see what is wrong. As far as I'm aware, my syntax is correct and the statement executes but nothing is done. MTFCC are listed as strings in the properties tab.


texts: Osm2pgsql: Error Loading in north-america-latest.osm.pbf

texts:
I've had my OSM Tile server working great for several days now, but just with Indiana. I decided to fill the gis database with all of North America.
I run the following command:

sudo -u user osm2pgsql --slim -d gis -C 2000 number-processes 3 /usr/local/share/maps/planet/north-america-latest.osm.pbf

After Running I get this Error:

  Using built-in tag processing pipeline
  Using projection SRS 900913 (Spherical Mercator)
  Setting up table: planet_osm_point
  Setting up table: planet_osm_line
  Setting up table: planet_osm_polygon
  Setting up table: planet_osm_roads
  Allocating memory for dense node cache
  Allocating dense node cache in one big chunk
  Allocating memory for sparse node cache
  Sharing dense sparse
  Node-cache: cache=2000MB, maxblocks=256000*8192, allocation method=11
  Mid: pgsql, scale=100 cache=2000
  Setting up table: planet_osm_nodes
  Setting up table: planet_osm_ways
  Setting up table: planet_osm_rels

  Reading in file: number-processes
  error while opening file number-processes

I've searched the internet for about 2 hours now and can't find anything useful. Maybe I'm missing some key information. Anyone ran into this?


texts: ArcGIS Server: Limit Access to Port 6080 Based on IP?

texts:
We have an ESRI ArcGIS Server 10x in our intranet which listens on both port 80 (through the Web Adaptor) and port 6080. So the web "Services" are accessible from both ports. I want to limit access to who can access the Services based on IPs. In case of port 80, it was fairly easy: Using IIS Manager, select the virtual directory and then using the "IP Address and Domain Restrictions" feature assign the allowed IPs.

But how to restrict based on IPs for port 6080? I think I could simply disable port 6080 and use 80 but that would require some applications to be adjusted.


texts: Why is my infowindow disabled on the data set I uploaded from NYCOpen Data or the sample data on CartoDB?


texts: ArcFM Session Manager is showing error while OPEN EDIT SESSION

texts:
I am configuring session manager from scratch. I followed below link and executed the steps.
http://resources.arcfmsolution.com/10.1/EngineConfig/SMOracle.html

I configured process framework administrator also. Screen shot attached.

I am using SDE user. when I try to Open edit session after creation I am getting error and not able to open the session:

This session's version could not be found.  You may not be connected to the

database, the version referenced by this session may have been deleted, or

the version's permission may have been changed to private.

I checked with below things:

SDE has all the permission on the mm tables
Compressed the database and checked state is 0 for the version

Not sure what is wrong with database . Please suggest


texts: Python, shapely library: is it possible to do an affine operation on shape polygon?

texts:
I am trying to create a rotated rectangle in lat/lon coordinates.

Basically, I have the coordinates of the center and rotation angle around the center.

Is it possible to do using, shapely or gdal, or something else?


texts: GeoJSON generator with OpenLayers support

texts: Is there a service that allows me to draw on an OpenLayers map, and see the resulting GeoJSON?


texts: arcpy with multiprocessing throws "UnpickleableError: Cannot pickle  objects"

texts:
I'm attempting to speed up a process which is currently running synchronously, by using the python multiprocessing module.

I'm having trouble sending a feature layer to a function which is called by multiprocessing, as demonstrated in this simple script:

import multiprocessing, arcpy

def doProcess(lyr):
    print(lyr.name)

if __name__ == '__main__':

    #Create an array of feature layers
    arcpy.env.workspace = "C:\Program Files (x86)\ArcGIS\Desktop10.2\TemplateData\TemplateData.gdb"
    featureLayers = []
    fcs = arcpy.ListFeatureClasses("*","All","World")
    for fc in fcs:
        arcpy.Delete_management(fc + "_lyr")
        lyrName = fc + "_lyr"
        arcpy.MakeFeatureLayer_management(fc, lyrName)
        featureLayers.append(arcpy.mapping.Layer(lyrName))

    #This works when not using multiprocessing:
    for featureLayer in featureLayers:
        doProcess(featureLayer)

    #This fails with "UnpickleableError: Cannot pickle <type 'geoprocessing Layer object'> objects"
    pool = multiprocessing.Pool()
    pool.map(doProcess, featureLayers)
    pool.close()
    pool.join()

When iterating over the array manually, rather than using multiprocessing, the function has access to the feature layer. But when using multiprocessing, this error message is shown:

  UnpickleableError: Cannot pickle type 'geoprocessing Layer object'
  objects

What is the correct syntax/approach to handle a feature layer within the multiprocessing environment? I based the above script on the example on the Esri blog Multiprocessing with ArcGIS


texts:
I'm trying to plot GIS data, a tiff file using gdal and then imshow() function but I am getting the following error -

My code is -

import gdal
import sys
import matplotlib.pyplot as pyplot
import numpy as np
import matplotlib as ml

try:
    tif = gdal.Open('lulc250k_1415_13861.tif')
    tifArray = tif.ReadAsArray()
except:
    print 'The file does not exist.'
    sys.exit(0)

band1 = tif.GetRasterBand(1)

band1Array = band1.ReadAsArray()

pyplot.imshow(band1Array)
pyplot.show()

The gdalinfo of the file -

Driver: GTiff/GeoTIFF
Files: lulc250k_1415_13861.tif
Size is 32716, 22789
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (72.486328125000000,32.233886718750000)
Pixel Size = (0.000509086484824,-0.000509086072228)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  72.4863281,  32.2338867) ( 72d29'10.78"E, 32d14' 1.99"N)
Lower Left  (  72.4863281,  20.6323242) ( 72d29'10.78"E, 20d37'56.37"N)
Upper Right (  89.1416016,  32.2338867) ( 89d 8'29.77"E, 32d14' 1.99"N)
Lower Right (  89.1416016,  20.6323242) ( 89d 8'29.77"E, 20d37'56.37"N)
Center      (  80.8139648,  26.4331055) ( 80d48'50.27"E, 26d25'59.18"N)
Band 1 Block=32716x1 Type=Byte, ColorInterp=Palette
  NoData Value=0
  Color Table (RGB with 256 entries)


texts: How to create landscape resistance map using path level spatial randomisation in QGIS?

texts: I need to create a landscape resistance map based on path level spatial randomisation, as I have read in some research papers. Are there tools for this?  I need to randomly shift and rotate the the utilised path I have to create available paths with identical spatial topology. I then need to compare the distribution of utilised path with avialable paths from each landscape resistant factor to determine resistance values.


texts: Toggle Layers On/Off using keyboard in Arcmap 10.3

texts:
Is there a way to set keyboard shortcuts to toggle layers on/off while editing?

For example F1 will turn on /off my top layer in my table of contents, F2 will toggle the 2nd layer, ect.


texts: SQL Server Geometry draw rectangle around a line

texts:
I am VERY new to spatial data queries, so please bear with me.
I have a LineString (geometry) in SQL Server and need to "draw" a rectangle around it, based on padding for the Width and the Height, using the existing line.
Is there a direct way to do so using the Spatial API in SQL Server?


texts: What does the "service" box in the PostGIS connection dialog stand for?

texts:
Can sоmеbody explain me what to write in the "service" box when connecting to the local DB.

QGIS Lisboa and Postgres are installed.


texts: Seeking way to tile/pyramid 4 million images using gdal for geoserver WMS?

texts:
My scenario consist of 4 million geolocated images each measuring 1536x1536 pixels in size.

What is the most efficient way to pyramid these images?

My end goal is to serve these images using geoserver via WMS requests.

I have tried creating a virtual dataset of all images and putting them into gdal_retile.py. I have also tried merging each tile with a transparent "null" image that covers the globe then pyramid the result. Both methods seem like they will take forever and make me think there is a better way.

I also have the requirement of adding additional images in the future.

Is there an easy way to do this without having to regenerate everything again?

I am new to geoserver and gdal.


texts: Uncaught SyntaxError: Unexpected token < (Geoserver popup in leaflet using ajax)

texts:
I am trying to show attribute information of features in a popup from my postGIS db tables via Geoserver WMS.
I have successfully been able to add these layers to the map. However, when I try to grab the attribute information for the same layers via Ajax, console prints an error message stating Uncaught SyntaxError: Unexpected token <
On the other hand, the status for URL is 200 i.e. ok  as can be seen here:

Still, the preview of this query is

<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<!DOCTYPE ServiceExceptionReport SYSTEM "http://localhost:8888/geoserver/schemas/wms/1.1.1/WMS_exception_1_1_1.dtd">
<ServiceExceptionReport version="1.1.1" >
    <ServiceException code="LayerNotDefined" locator="layers">
        Could not find layer workspace:
    </ServiceException>
</ServiceExceptionReport>

Below are my query string parameters

1.  Query String Parameters (17)
2.  SERVICE:  WMS
3.  version: 1.1.0
4.  request:  GetFeatureInfo
5.  layers: workspace%3AlayerA,workspace%3AlayerB,workspace%3AlayerC,workspace%3AlayerD,workspace%3AlayerE,workspace%3AlayerF,workspace%3AlayerE,
6.  styles:
7.  bbox: 74.22517776489258,31.487235582017444,74.4019889831543,31.559668260234595
8.  feature_count:  1
9.  width:  1030
10. height:  495
11. format:  image/png
12. srs:  EPSG:4326
13. (empty)
14. x:  247
15. y:  30
16. info_format:  text/javascript
17. callback:  jQuery20309529768737157973_1547708903539_:1547708903542

Everything seems to be working fine here except that it is unable to find the specified layers. I have cross checked in my geoserver’s web.xml file that jsonp is enabled as can be seen below

I just can’t understand what possibly could go wrong here when I’m successfully calling the same layers in a different function with the same URL and here it doesn’t work.

For further reefrence I’m pasting the code snippet here:

URL='http://localhost:8888/geoserver/wasa_assert/wms?SERVICE=WMS
&version=1.1.0
&request=GetFeatureInfo
&layers=workspace%3AlayerA,workspace%3AlayerB,workspace%3AlayerC,workspace%3AlayerD,workspace%3AlayerE,workspace%3AlayerF,workspace%3AlayerE,
&styles=
&bbox=' + BBOX + '
&feature_count=1
&width=' + WIDTH + '
&height=' + HEIGHT + '
&format=image%2Fpng
&srs=EPSG%3A4326&
&x=' + X + '
&y=' + Y + '
&info_format=text/javascript';

$.ajax({
    url: URL,
    crossDomain: true,
    dataType: "jsonp",
    success: parseResponse
});


texts: Problem Merging multiple shapefiles with different CRS in QGIS

texts:
I am trying to merge multiple shapefiles into one, but I encounter one error, not all of my shapefiles end up merged. I am able to load all of my shapefiles but when I attempt to merge them, the merged attribute table contains all the information that it should but not all the shapefiles show-up. I made sure to change the the headings in the attribute tables that I want merged.

I have tried to use a number of tools (Merge Shapefiles to One, MMQGIS Combine etc.) but the problems seems to be that some of my shapefiles are projected differently. I have tried to change the coordinate reference system under Properties but after I save the file, the shapefile disappears. I have also found a guide to use GRASS to change the projection but I was not able to follow it.


texts:
I'm running a process that requires me to change the amount of data input depending on the situation.  With my current system I need to have a separate model depending on the number of files I need to process.  I would like something that's more flexible so I only need one model that can handle different amounts of files.

What I'm trying to do is find any point outside of a buffer around the input layer(s).

Here's what my model currently looks like but it will only process 8 inputs into 'images'.   I need it to be able to handle any amount eg 1-8.  I have them set as optional but it still errors saying I need data.


texts: Using geocode function with ArcPy cursors?

texts:
How do I use a cursor with a geocode address function?

I want to be able to:

read from a table
only geocode addresses(from table1) whose street number is not 0 using a cursor
then update a master table with the results from the geocode

# only try to geocode those addresses whose first street number digit is not zero

  # Select the first char of the street number in the aRecord.ADDRESS1 field
    MyStr = "!aRecord.ADDRESS1!.split(" ")[0]"
    if MyStr != '0':
        address_table = "table1"
        address_locator = r"http://myservername/ArcGIS/rest/services/Composite_Locator/GeocodeServer"
        geocode_result = "geocode_result"

        #Add XY fields and calculate XY
        sr = arcpy.SpatialReference(4100)
        with arcpy.da.UpdateCursor(geocode_result, ["SHAPE@XY", "X", "Y"], "", sr) as cursor:
            for row in cursor:
                row[1] = row[0][0]
                row[2] = row[0][1]

                # geocode is successful only if it returns a valid X (and Y)

            if (row[1] > 0.0):
                cursor.updateRow(row)
                Row_tobe_Inserted.X = row[1]
                Row_tobe_Inserted.Y = row[2]
                XY_Method = "Geocode"

            del row, cursor
    else:
        print "There are no records to be geocoded"


texts: Changing spatial extent of shapefile using ArcGIS Desktop?

texts:
I want to change the extent of one of my shapefiles(lets call it A) to the same extent of another shapefile(B) but I can't seem to do it. I tried using the 'recalculate feature class extent' and changing the processing extent under Environments from Default to shapefile B's extent. I, of course made sure the input feature was shapefile A. I clicked OK and checked the extent of shapefile A but there was no change.

What can I do to change the extent?


texts: ArcPy Utility Network Trace Not Returning Any Result

texts:
Using ESRI Tutorial Data for Utility Network Tracing (Montgomery.gdb) which is a File Geodatabase and following Arcpy script  I am not getting any erro on running the script to find the path between two points in Flag (point) layer

# Import arcpy module
import arcpy
# Local variables:
gnVersionFDS_Net = "E:/Montgomery.gdb/Water/Water_Net"
Flags = "E:/Montgomery.gdb/Water/Flags"
gnVersionFDS_1 = "gnVersionFDS_1_Net"
# Process: Trace Geometric Network
arcpy.TraceGeometricNetwork_management(gnVersionFDS_Net, gnVersionFDS_1, Flags, "FIND_PATH", "", "", "", "", "", "NO_TRACE_ENDS", "", "", "", "AS_IS", "", "", "", "AS_IS")
print "The Proccess Done"

I added the Flag layer manually to the Dataset and it has two points as

again, I am not getting any error but I am not seeing any layer and path to be generated or added to somewhere. Can you please let me know what I am doing wrong or how can I make this rum properly?

  Update

I try to copy the output as a shapefile using this ArcPy

arcpy.Copy_management(gnVersionFDS_1, "C:/data/path.shp")

but I am getting  this error now


texts: Assuming that I have a point on a map and a multipolygon on a map. When I query with ST_Intersects or ST_Within I get the same result. I want to know what is different between these two methods?


texts: I need to export the list of IDs of the satellite images available in Google Earth Engine to work in Excel. I just need the ID

texts: I am working with Google Earth Engine and I need to export the collection as a list, to be able to work it in Excel. What I mean is being able to export the list of the ID of the available satellite images of any Landsat, that is located in the console of Google Earth Engine. What I need is just this Ej: LANDSAT/LT05/C01/T1_TOA/LT05_001083_19860204.


texts: Get Unique Field Values Within Python Script Tool

texts:
Using ArcGIS 10.2.2 and Python. I have a script that currently counts the number of points (representing crimes) in each zip code area, and writes that value to a new column in the FC that contains zip code boundaries.

Basically, the user inputs a point FC (crimes) and a polygon FC (zip boundaries), and the tool writes the number of crimes in each zip code boundary to a new field in the zip boundaries polygon FC.

What I would like to do is improve this code by allowing the user to choose the categories (which is contained in a separate field in the crimes point FC) of crimes that get counted for each zip code boundary.  For example, rather than just counting all crimes, I would like the user to be able to select, say, Aggravated Assaults and Grand Theft (two of many values contained within the Category field of the Crimes FC) and have the tool only count those points.

# Import ArcPy Site Package
# Set environments settings: workspace location
import arcpy, sys

# Define location of farms and census sub division feature classes
crimes = arcpy.GetParameterAsText(0)
zip_parcels = arcpy.GetParameterAsText(1)
workspace = arcpy.GetParameterAsText(2)
new_field = arcpy.GetParameterAsText(3)
new_field_alias = arcpy.GetParameterAsText(4)

# Check if the field named new_field exists in the 'zip_parcels'
# feature class. If it exists, delete it.
if arcpy.ListFields(zip_parcels, new_field):
    sys.exit(arcpy.AddError("Error: Field Already Exists"))

# Create a new 'NumOfFarms' field in the 'zip_parcels' feature class
arcpy.AddField_management(zip_parcels, new_field, 'Long', '', '', '', new_field_alias)

# Create temporary feature layers
crimes_lyr = arcpy.MakeFeatureLayer_management(crimes,'crimes_lyr')
zip_parcels_lyr = arcpy.MakeFeatureLayer_management(zip_parcels,'zip_parcels_lyr')

# Create an update cursor to access and update the temporary
# census subdivision feature layer.  For each sub division (or
# row), the NumOfFarms field will be updated
fields = ['ZIPCODE',new_field]
with arcpy.da.UpdateCursor(zip_parcels_lyr,fields) as cur:
    for row in cur:

        # Create a query to select a census sub division that will dictate which farms
        # are selected
        zip_parcels = row[0]
        where = '"ZIPCODE" = \'{}\''.format(zip_parcels)
        arcpy.SelectLayerByAttribute_management(zip_parcels_lyr,'NEW_SELECTION',where)

        # Select the farms that intersect (or are within) the selected census sub
        # division and count how many are selected
        arcpy.SelectLayerByLocation_management(crimes_lyr,'WITHIN',zip_parcels_lyr)
        number_of_crimes = int(arcpy.GetCount_management(crimes_lyr).getOutput(0))

        # Populate the NumOfFarms field with the counted number of selected farms, which
        # also updates the the field in the csd feature class
        row[1] = number_of_crimes
        cur.updateRow(row)

Please excuse some of the commenting - it sometimes says farms and census subdivisions because I used this for a different project and have been slowly modifying it.

Where I need help is how to generate a multivalue checklist within the script tool GUI (when you double click the script in ArcCatalog sidebar) that lists all of the unique values within the Category field of the crime point FC.  The checkboxes (field values) that the user selects will dictate the crimes that get counted.

Additionally, I would like to add that I have tried the solution listed here: http://blogs.esri.com/esri/arcgis/2012/07/16/generating-a-multivalue-choice-list/

This solution crashes Arc every single time and I feel there must be some easier way accomplish this task.


texts: Rasterize vector to presence-absence raster incorrect in R

texts:
I want to create a presence-absence raster out of a shapefile (vector) with the rasterize command of the R-package in R. Until now I only get wrong results.

The raster should just distinguish between the two options ('shapefile on this cell present', 'shapefile on this cell abscent').

I tried it out with several options of the rasterize-command:

forest_polygons_raster_mask <- rasterize(forest_polygons, ras, mask = TRUE)

forest_polygons_raster_update <- rasterize(forest_polygons, ras, update =
TRUE)

forest_polygons_raster_cov <- rasterize(forest_polygons, ras, getCover =
TRUE)

the command with the getCover option does not stop, also after several hours. With both of the other commands I get a wrong result and my raster has presence-cells in areas where the shapefile is abscent.

The original shapefile with two holes:

That is the raster I get after the rasterize command with the update-option. One hole is missing, and the other is split into two smaller holes. No abscence was recognized wrong, just presence was recognized wrong:

Shapefile layer over raster layer, the holes in the shapefile are not left out as abscence in the raster:

Did anyone have a similar problem? What might the problem be? My shapefile is quite ordinary with different fields of the categories: ID, Shape (Polygon), Classification (text), Shape_Area, Shape_Length.
The cellsize of the input raster for the rasterize process is quite small (so I do not expect that the errors are due to two big cells).


texts: Advantage of Digital Elevation Model (DEM) in Agriculture?

texts: I have given the task to create Digital Elevation Model (DEM)'s of an agricultural area. I am programmer with Computer Science background and very new to GIS can anyone explain what is advantage of using Digital Elevation Model (DEM) in agriculture ?


texts: For-Looping through Bookmarks in Python Script tool?

texts:
I'm using this script tool, in ModelBuilder for ArcGIS Pro, to set up a map layout, reference a bookmark, and export the layout as a PDF. This model iterates through Urgent Care centers and its corresponding data one at a time (in alphabetical order by center name). Since my model loops through each Urgent Care center (via the Iterate Field Values), it repeats the same bookmark through every loop/iteration. How do I tell it to go from one bookmark to the other?? (I have about 90 centers, so that means 90 bookmarks also)

Here is my script tool:

import arcpy

#input layer
lyr = arcpy.GetParameterAsText(0)

# input name of layout
p = arcpy.mp.ArcGISProject("CURRENT")
lyt = p.listLayouts("Layout_King")[0]

# Reposition the scale bar
scaleBar = lyt.listElements("MAPSURROUND_ELEMENT", "Scale Bar")[0]
mf = scaleBar.mapFrame
scaleBar.elementPositionX = mf.elementPositionX + 0.0
scaleBar.elementPositionY = mf.elementPositionY - 0.5

# Reposition the north arrow
northArrow = lyt.listElements("MAPSURROUND_ELEMENT", "North Arrow")[0]
mf = northArrow.mapFrame
northArrow.elementPositionX = mf.elementPositionX + 8.8
northArrow.elementPositionY = mf.elementPositionY + 0.7

# Align the title with the center of the map frame
title = lyt.listElements("TEXT_ELEMENT","Name of Map Text")[0]
mf = lyt.listElements('MAPFRAME_ELEMENT',"Map Frame")[0]
title.elementPositionX = mf.elementPositionX + (mf.elementWidth / 3.7)
title.elementPositionY = mf.elementPositionY + (mf.elementHeight / 0.98)

# Reposition the Legend and fix legend title
legend = lyt.listElements("LEGEND_ELEMENT", "Legend")[0]
legend.title = "Legend"
legend.elementPositionX = mf.elementPositionX + 7.7
legend.elementPositionY = mf.elementPositionY + 7.15

# setting layout to bookmark
aprx = arcpy.mp.ArcGISProject("Current")

# add name of layout
lyt = aprx.listLayouts("Layout_King")[0]
mf = lyt.listElements("MAPFRAME_ELEMENT")[0]

# add name of bookmark
bkmks = mf.map.listBookmarks()
bkmks.sort(key=lambda x: x.name, reverse=True)
for bkmk in bkmks:
  mf.zoomToBookmark(bkmk)
lyt.exportToPDF(r"C:\arcGIS_Shared\Exports" + "\\" + bkmk.name + ".pdf")

I have heard that this could be accomplished via a for-loop, but don't know how to write one.


texts: Adding z-value to the triangulation output from the TIN interpolation in QGIS

texts:
For a project I'm working on I need to do a triangulation over points (csv) or lines (dxf). When running the TIN interpolation from the processing toolbox I save the triangulation as a temp file.
This polygon output has no z-value in its geometry.

With drape, I can get the z-value from the raster-output from the TIN interpolation, but not all points are on the raster so there are lines that have the default value.

Is there a way to fix this, or add the z-value to the geometry from the input source?


texts: Is it possible to create a geometry from a proj4 string and Cartesian coordinates in PostGIS?

texts:
Considering I have a proj4 string such as:
+proj=tmerc +lat_0=12 +lon_0=-12 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m

Is it possible to create a geometry given X and Y offset (in meters)? For example a point 100 meters away from the origin (lat:12,lon:-12) in the X axis:
ST_MakePoint(100, 0, '+proj=tmerc +lat_0=12 +lon_0=-12 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m')


texts:
I currently using Barefoot (https://github.com/bmwcarit/barefoot) for several Route Matching tasks.

The installation and usage of Barefoot just worked fine. The matching already works and I got some good results.

Now I would like to get additional information about the matched roads. Barefoot returns data like this (slimjson mode) :

 ... {
    "route": "LINESTRING (106.89085170694243 -6.142643394836986, 106.89085170694243 -6.142643394836986)",
    "road": 1199220,
    "frac": 0.7514822019311651
} ...

I would like to request information like speed level or road type. How can I retrieve this kind of information. Is there a way to query for example Nominatim?

What kind of Road ID is the returned one?