texts:
I have raster data that represents a probability distribution, i.e. each cell has a probability value (in my case the probability that an animal may be found in the cell), and all the cells add up to 100% (I know for certain the animal is within the extents of my raster).  I want to be able to generate vector data for confidence values.  For example, the 95% line/polygon denotes the boundary wherein I am 95% confident that I will find the animal.

Similarly, if I have a kernel density estimate, How do I generate the XX% line/polygon that borders the densest part of the raster containing XX% of the total population?

I am willing to use ArcGIS, or open-source software.  If there isn't a tool to perform this for me, what is an algorithm I can implement?


texts:
I have a Python script to download data in netCFD format. I have managed to extract the information I need in a masked array. This looks like:

print(my_array.shape())
>> (700, 299, 217)

As a working example I have selected only the first array of the first dimension, such as:

data = my_array[0][0]
print (data.shape())
>> (299,217)

This is stored as a masked array by default, this means that some of the values are empty. I am trying to export this to a raster file, so it looks like the image below 1. The array looks likes 2:

However, with the following code (see below) the raster is saved as the second picture 3. The empty values in my array are stored as 2.147483648000000000e+09. See attached filed.

np.savetxt('raster.txt', data)

or

data_1 = np.ma.masked_array(data, fill_value=None)
np.savetxt('raster.txt', data_1)

What I do next is with arcpy

import os
import sys
import numpy as np
import arcpy

arcpy.env.overwriteOutput = True
arcpy.env.outputCoordinateSystem = arcpy.SpatialReference('WGS 1984')
nraster = np.genfromtxt('raster.txt', missing_values='Null')
cell_size = 0.041666667
lat = 18.020833333333343
lon = -104.97916666666666
corner = arcpy.Point(lat, lon)

# Convert array to a geodatabase raster
my_raster = arcpy.NumPyArrayToRaster(nraster, corner, x_cell_size=cell_size, y_cell_size=cell_size)
my_raster.save('working_path/GDB.gdb/my_raster')

So my question is, how can I save this array to a .txt file (in Python) so I can read it back with arcpy.NumPyArrayToRaster and continue my process in ArcMap?


texts: Making directional arrow line symbol in ArcGIS Pro

texts:
How can I create a directional arrow line symbol in ArcGIS Pro?
On this page, they explain how to create a directional arrow line symbol in ArcMap. However, I'm using ArcGIS Pro and its documentation doesn't seem to mention this kind of symbology at all.
Here's an image of what I'm looking for (copied from the ArcMap's documentation).


texts: I want to create a Voronoi diagram while retaining the data in the data frame

texts:
We will use st_voronoi() to create a Voronoi diagram.

I have three questions.

Point one.

I understand that the argument to st_voronoi() should be the following data.

■SFC format (an error occurs in SF format)

■MULTIPOINT format(Not in POINT form.)
Is this understanding correct?

Point two.

If you run the following code, it will be converted to GEOMETRYCOLLECTION.
Is there a way to convert this to MULTIPOINT?

Point three.

The original "nz_height" contains the word "elevation".
However, the conversion to SFC format has erased that column.
I would like to figure out the elevation of each area of the Voronoi diagram.
What kind of processing should I do?

If there is anyone who knows about it, please tell us about it.

library(sf)
library(sp)
library(spData)

print(class(nz))

print(class(nz_height))
print(head(nz_height))
plot(st_geometry(nz_height))

#POINT to MULTIPOINT
nz_height_mp <- st_cast(nz_height,"MULTIPOINT")
print(head(nz_height_mp))

#sf to sfc
nz_height_mp_sfc <- st_union(nz_height_mp)
print(class(nz_height_mp_sfc))
print(head(nz_height_mp_sfc))

nz_height_v <-st_voronoi(nz_height_mp_sfc)

print(class(nz_height_v))
print(head(nz_height_v))
plot(st_geometry(nz_height_v),col=0)


texts: qgis utils.py appends to a dictionary / hash

texts:
I am working with the QGIS code in the QGIS/python/utils.py code (that file path on GitHub).

In line 280 the method append() is used on a dictionary (initialized in line 231), which I believe is not possible.

The line is a part of a method which searches for and updates the list of available plugins.

I am getting an error when I run the method updateAvailablePlugins() (line 270):

Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/qgis/utils.py", line 283, in updateAvailablePlugins
    plugins.append(pluginName)  # original
AttributeError: 'dict' object has no attribute 'append'

Is this an error in the code?  If 'pluginName' should be added to 'plugins' as a dictionary entry, what should the value of the key?

Something to note here is that I am able to use plugins just fine in the QGIS GUI, however I am getting this error when running the code directly in the Python console.

The end goal is to run a plugin, Time Manager, from the Python console, and this appears to be one step in the process.

Here is the github page with the full code:
https://github.com/qgis/QGIS/blob/3205c96af629194fa485e072050b677dacfc7f23/python/utils.py


texts: Using Local Routing Tool in ArcGIS ModelBuilder?

texts:
When I create a model in ArcGIS ModelBuilder using the Network Analyst Tools, I am able to create Locations and make the Route Layer and Solve.  However, I need to be able to add barriers as well as locations, before solving for the Route.

Is there a way to do this in ArcGIS ModelBuilder?


texts:
I have recently been programming in QGIS plugin with python using an extension from pip package management. I installed this extension "browserhistory" manually via the OSGeo4W shell.
However, when I now want to pass this plugin to someone else, he gets the error message: "No module named browser-history". How can I modify the code of the plugin in such a way that this pip extension is automatically installed when the plugin is installed on another computer? So that the other person does not have to install pip manually.


texts: Error in converting SpatialPolygonsDataFrame to Data Frame, R Script in QGIS

texts:
I am trying to convert a SpatialPolygonsDataFrame layer to a DataFrame to create plots using ggplot2 in an R Script within QGIS. I've used spatial join to append averages of BioClim mean temperature tiles within each of many polygon land class features.

When I pull up the attribute table for the layer, all values are filled for each of the features; however, R in QGIS keeps throwing the error:

Error in data.frame(x=Layer[[X]], y=Layer[[Y]]):
arguments imply differing number of rows: 9382, 0
Execution halted

For some reason it doesn't seem to be picking up the second field of the SpatialPolygonsDataFrame object.

The R script I've composed thus far is here:

##Vector Processing=group
##showplots
##Layer=vector
##X=Field Layer
##Y=Field Layer

library("ggplot2")

ggplot(data=data.frame(x=Layer[[X]], y=Layer[[Y]])) + geom_histogram(aes(x=x)) + facet_wrap(~y)

This question is a follow-up to one of my earlier ones, which is here: Use facet wrap in ggplot2 with QGIS?

Does anyone have an ideas of why this error is being thrown?


texts:
Another labelling question. I have shape-points which are shown as circles in QGIS. The bigger the data, the bigger the circle. Now i want to label the points with the data defined properties, so that the label of each point is directly shown on top of the circle. How to do that? I have to move the y-coordinate...

Basically like this?

CASE

WHEN "ID"= 10 THEN y-coordinate= +5

WHEN "ID"= 20 THEN y-coordinate= +10

and so on...

ELSE y-coordinate should stay normal

END


texts: QGIS behaves different when i have the "python console" active compare to not active

texts:
I'm developing my QGIS3 plugin and I've realized that if I have the Python console active (so I can see it) the plugin behaves as expected, while if I don't see it doesn't. So this is the code:

for field in fields_db:
    field = field[0]
    layer = self.db.addPostGISLayer('fields', 'polygon', 'public',
                                    extra_name=field + '_',
                                    filter_text="field_name='{f}'".format(f=field))
    set_label(layer, 'field_name')
    QgsProject.instance().addMapLayer(layer)
    print('Added')
set_zoom(self.parent.iface, 1.1)
self.defined_field = ''

When I run the code with the python console present the layers adds one by one, and zoom to the fields.
However when I don't have the Python console (or print('added') is comment away with a #). The code wait until all layers are ready to show and adds them at the same time in the end and it doesn't manage to zoom to the fields either(?). Any suggestions what I may edit?

This is the set_zoom code :

def set_zoom(iface, extra_extent):
    """
    Sets the zoom level to include all layers (excluding tiles layer) with some extra extent
    Parameters
    ----------
    iface: QGIS interface
        The QGIS iface module.
    extra_extent: float
        How much extra space around the layers eg. 1.1 is 10% extra
    Returns
    -------

    """
    zoom_extent = QgsRectangle()
    for layer in QgsProject.instance().mapLayers().values():
        if 'xyz&url' not in layer.source():
            zoom_extent.combineExtentWith(layer.extent())
    if zoom_extent.center().x() != 0.0:
        wgsCRS = QgsCoordinateReferenceSystem(4326)
        QgsProject.instance().setCrs(wgsCRS)
        zoom_extent.scale(extra_extent)
        iface.mapCanvas().setExtent(zoom_extent)
        iface.mapCanvas().refresh()
        wgsCRS = QgsCoordinateReferenceSystem(3857)
        QgsProject.instance().setCrs(wgsCRS)


texts:
I'm trying to convert a RasterBrick to polygon, across four different resolutions. It works for the coarsest three resolutions, but for the finest resolution, I get an error that I can't find much information on.

Here's the RasterBrick Info:

    > rass_agg
class       : RasterBrick
dimensions  : 1883, 2926, 5509658, 14  (nrow, ncol, ncell, nlayers)
resolution  : 13710, 12210  (x, y)
extent      : -20051271, 20064189, -2722745, 20268685  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +units=m +no_defs
names       : Wood_Thru//is_Bare_06, Wood_Thru//rubland_06, Wood_Thru//ropland_06, Wood_Thru//oadleaf_06, Wood_Thru//dleleaf_06, Wood_Thru//oadleaf_06, Wood_Thru//dleleaf_06, Wood_Thru//assland_06, Wood_Thru//_Forest_06, Wood_Thru//rubland_06, Wood_Thru//n_Built_06, Wood_Thru//s_Water_06, Wood_Thru//assland_06, Wood_Thru//oodland_06
min values  :           -0.06590122,           -0.03429975,           -0.13205945,           -0.03906334,           -0.03001432,           -0.28689602,           -0.03063185,           -0.28541537,           -0.11917765,           -0.01400324,           -0.06292345,           -0.15828412,           -0.02550580,           -0.09214195
max values  :          5.590892e-03,          9.034161e-02,          5.403338e-02,          9.276621e-01,          2.288492e-03,         -1.000000e-04,          4.427548e-03,          1.393522e-03,          1.852588e-01,          2.967395e-03,          5.623372e-02,          9.160649e-03,         -1.000000e-04,         -9.698019e-05

Then when I call rasterToPolygons(), I get an error.

rass_poly <- rasterToPolygons(rass_agg, n=4, na.rm=T)

The error:

  Error in rbind(xyv, xyvr) :
    number of columns of matrices must match (see arg 2)

Is this a bug?


texts: ArcGIS Javascript: upgrade ComboBox to FIlteringSelect using ItemFileReadStore

texts:
I would like to upgrade a web app I am working on from v2.1 to v3.1 of the Javascript API (I have tried 3.2 but I am getting some very odd rending problems, so I'm settling with 3.1 for now...)
One issue I have encountered is that a ComboBox populated with a store does not seem to fully function with 3.1.  I was told I can fix this by replacing the ComboBox with a FIlteringSelect widget, but so far I haven't succeeded.
Below is my code.  I was hoping I wouldn't need to do much more than replace the ComboBox object with the FilteringSelect object, but that apparently isn't sufficient.  Can somebody help me see what I am missing? Thanks, Jason

JS (old code:)

      function populateSchoolMenu(results)
{
    var schoolName;
    var values = [];
    var testVals = {};

    var features = results.features;
    dojo.forEach(features, function (feature) {
        if ( (schoolLayerFlag == 0) || (schoolLayerFlag == 2) )
        {
            schoolName = feature.attributes.Name;

        }
        else if (schoolLayerFlag == 1)
        {
            schoolName = feature.attributes.Site_name;

        }

        if (!testVals[schoolName])
        {
            testVals[schoolName] = true;
            values.push({ name: schoolName });
        }
    });

    var dataItems = {
        identifier: 'name',
        label: 'name',
        items: values
    };
  var menuStore = new dojo.data.ItemFileReadStore({ data: dataItems });

  dijit.byId("schoolMenu").store = menuStore;
   }

HTML (old:)

            <select id = "schoolMenu"
            dojotype="dijit.form.ComboBox"
            value = "Select School"
            autoComplete="true" class="comboBoxClass"
            forceValidOption="false"

            onchange="executeSelectedQuery(this)">
            </select>


texts:
I want to project the shadows from a DSM raster.

Is it possible to iterate the grass function r.sunmask.datetime
over a list from a csv file composed with interest dates(month/day) and time (hours from 9:00 to 18:00) along the year?


texts: QGis Zonal Statistics Plugin NoData included in Calculations

texts:
I have been trying to use the QGIS zonal statistics plugin. My raster has a nodata value of -99999 and all other values are positive. After computing the mean zonal statistics there are a number of features that have a negative mean indicating that the Zonal Stats plugin must be included nodata values in the calculations. Is there a workaround for this or am I doing something wrong:

from qgis.analysis import QgsZonalStatistics
vectorLayer = qgis.utils.iface.activeLayer()
rasterLayer = "path/to/my/raster/et_30.tif"
prefix = 'Test'
zoneStat = QgsZonalStatistics(vectorLayer, rasterLayer,prefix)
zoneStat.calculateStatistics(None)


texts: How to fix "Unable to initialize ODBC connection to DSN for ALC_AV_CO"?

texts:
So I am new to the QGIS community, and I have been trying to bring in some rather large access databases with little to no success. I have windows 7, so I've gone in through administrative tools/Data Sources (ODBC) to create a new user DSN.  I made it a Microsoft Access Driver (.mdb,.accdb), and linked it to the database I wish to use.

Then in QGIS, I have used add vector layer, selected Database, left type ODBC, given it a name, specified localhost and for Database, written in an identical name to the User DSN I created... BUT... every time I try to test it, it gives the following error:

Connection failed - Check settings and try again.

Extended error information:
Unable to initialize ODBC connection to DSN for ALC_AV_CO,
[Microsoft][ODBC Driver Manager] The specified DSN contains an architecture mismatch between the Driver and Application

Why is it doing this?

Interestingly enough the only time I had success, I managed to somehow create a User DSN which was a Microsoft Access Driver (*.mdb) without the .accdb (I saved a version of my DB in .mdb to make this work too).  However, for reasons unknown to me, this option is no longer available to me in the create new data source menu, so I can't use it again!

I would really appreciate if someone could give me advice on what I am doing wrong!


texts: How can I define hierarchical levels in eCognition

texts: I'd like to define classes with 2 levels. Firstly I'm gonna classify some broad classes such as forest, grassland in a level 1. Furthermore, I'd like to classify these classes into some sub-classes such as Evergreen forest, deciduous forest etc. Grassland will also be classified into 3 sub-classes as well. How can I define this hierarchical levels in eCognition.


texts: Copy single feature into another feature class

texts:
I do a lot of editing, in this case areas. Sometimes I add an area to one shapefile, only to realise a bit later that it should have been in another.

Is there an easy way of copying only that polygon into the proper shapefile? Something along the lines of rightclick the row in attribute table, copy and then paste in the next attribute table.

I do know that I can select the polygon and use Append, but this seem unnecessary complicated?


texts: Label Expression If Statement remove Zeros

texts:
I am trying to perform a semi-complicated labeling expression that shows one field divided by another field and some additional fields with no division. The issue I ran into was that the denominator sometimes shows up as a 0. So, I need to build an if statement to perform the division if [TOT_HH_1] <> 0, do __ and if [TOT_HH_1] = 0, do __.

Also, there are three denominators that may or may not be 0. TOT_HH_1, 5 and 10.

Here is what I have so far.

Function FindLabel ([ELEM_1], [TOT_HH_1], [ELEM_5], [TOT_HH_5], [ELEM_10], [TOT_HH_10] )
  if ( [TOT_HH_1] <> 0 OR [TOT_HH_5] <> 0 OR [TOT_HH_10] <> 0) then
    FindLabel = Round ([ELEM_1] / [TOT_HH_1] , 2)  & " , " & Round ( [ELEM_5]/ [TOT_HH_5] , 2) & " , " & Round ( [ELEM_10] / [TOT_HH_10] , 2)  & "/" & [TOT_HH_1]  & " , " & [TOT_HH_5]  & " , " & [TOT_HH_10]
  elif ( [TOT_HH_1] = 0 ) then
    FindLabel = 0  & " , " & Round ( [ELEM_5]/ [TOT_HH_5] , 2) & " , " & Round ( [ELEM_10] / [TOT_HH_10] , 2)  & "/" & [TOT_HH_1]  & " , " & [TOT_HH_5]  & " , " & [TOT_HH_10]
End Function

Edit:
Correct, the code in the answer worked great, here is a picture of what the labeling looks like.


texts: Error parsing GeoJson - parse error - invalid geometry


texts: Convert 1-band .nc file to GeoTiff with RGB

texts:
How can I convert 1-band .nc file to GeoTiff with RGB?

I have tried using gdal translate but it won't let me apply color interp of RGB to the .nc as it has 1 band.


texts:
I'm really new to QGIS

1 - How can I get the radius for each Hub so it fills out my entire country?
Is there a plugin for that. Because with Puffers I have to guess myself and it doesn't really work out with the distances which brings me to my second question

2 - I have chosen DHDN / Soldner Berlin as my CRS. Is that right to get distances in metres for Germany?


texts:
I always wondered where cartographers get their bathymetric data from to create maps like the Zanzibar map on wikipedia for instance. I could download polygons for the 1000 meter polygons on naturalearthdata.com but I want to show maximum 100 meter levels or lower. Thus the question: What ist the easiest, user-friendliest way? Like naturalearthdata, but in 100s...

There was a similar question asked in this forum, I tried the named sources, but I couldn't figure out where exactly to download and how to transform the format.  Therefore I need a better description on HOW to download the data and convert it to an appropriate GIS format.

I work with ArcMap 10.4 or QGIS 2.16. A SRTM raster tile (geotiff) would be fine (like for the land areas) or Polygon/Line Featureclasses/Shapefiles.


texts: How to split a single feature in my polyline when it changes direction?

texts: I'm wondering if there is a plugin in QGIS where i could split my lines in my poly line when it changes direction, as show in the picture below:


texts:
I have tried to install shapely on Ubuntu 12.10. Therefore I installed libgeos-c1, pip and shapely with $ pip install Shapely. The shapely files are now in /home/user/build/Shapely. When I try to import it in my python script with

from shapely.geometry import Point

gives the following error:

from shapely.geometry import *
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 309, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
ImportError: No module named shapely.geometry

Any idea why it is not working?


texts: Extracting mean of multiple raster layers using R?

texts:
I don't know how to deal with big data that I had from radars.

I've have 144 images .tif per day and I need to extract a mean daily value, for each pixel, and export it in one image .tif.

I thought to stack all the 144 images, but i don´t know how to continue and how to export the result in one image.

I'd started this:

set("D:/R/Radar/Rasday1")
f <- list.files(getwd())
ras <- lapply(f,raster)
STACK1 <- stack(ras)   ## with dimensions :1006, 804, 808824 (nrow, ncol, ncell)

Then I don't know if its OK to create a raster and do something like:

r1<- raster(ncol=804, nrow=1006)
media<-mean(STACK1, r1)


texts: QGIS 3.4 crashing when loading two different files with different coordinate systems

texts:
I'm having trouble loading some files in QGIS 3.45 "Madiera". I have a CSV loading as WGS84, and a BMP loading in as UTM 12N NAD 83. Both the files load fine individually but as soon as they are loaded in the same map, it crashes. I've done the same thing in QGIS 2.18 and it works perfectly fine.

I've tried setting the project CRS to both of the CRS and it hasn't worked.


texts: Offline usage of OpenStreetMap or alternatives

texts:
I'm planning to implement an offline map solution for windows phone using open street maps. To achieve this, would I need to download complete OSM database first and than use it? Or is possible to fetch it on-demand and store it locally?

If OSM is not the best solution what else could be used to have an offline maps application? Offline data would be required only for about 3-5 cities, not more than that.


texts: Calculate max of all cells of a raster with gdal_calc.py

texts:
I want to calculate max of all cells of a raster, i.e. the result should be a scalar.

Regardless of what I call amax, max, numpy.max(A, axis=0) I always get a raster with the same WxH as the original raster. I suppose I should get a single scalar value. Is this possible with gdal_calc.py?

gdal_calc.py -A 0009.TIF --calc "max(A)" --debug
gdal_calc.py -A 0009.TIF --calc "numpy.max(A, axis=0)" --debug

gdal_calc.py starting calculation numpy.max(A, axis=0)
file A: 0009.TIF, dimensions: 7800, 7900, type: UInt16
Generating output file gdal_calc.tif
output file: gdal_calc.tif, dimensions: 7800, 7900, type: UInt16
using blocksize 7800 x 1
0 .. 10 .. 20 .. 30 .. 40 .. 50 .. 60 .. 70 .. 80 .. 90 .. 100 - Done

gdalinfo gdal_calc.tif
Driver: GTiff/GeoTIFF
Files: gdal_calc.tif
Size is 7800, 7900
Coordinate System is:
PROJCS["WGS 84 / UTM zone 33N",
    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"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",15],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32633"]]
Origin = (271485.000000000000000,4743015.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  271485.000, 4743015.000) ( 12d12'19.09"E, 42d48'20.13"N)
Lower Left  (  271485.000, 4506015.000) ( 12d17'46.77"E, 40d40'24.15"N)
Upper Right (  505485.000, 4743015.000) ( 15d 4' 1.63"E, 42d50'22.88"N)
Lower Right (  505485.000, 4506015.000) ( 15d 3'53.75"E, 40d42'18.09"N)
Center      (  388485.000, 4624515.000) ( 13d39'30.36"E, 41d45'52.71"N)
Band 1 Block=7800x1 Type=UInt16, ColorInterp=Gray
  NoData Value=65535


texts:
I have downloaded osm2po 4.7.7 but when trying to follow the instrutions in either the documentation but I'm not able to access the map in my browser on http://localhost:8888/Osm2poService. What I have done is:

to download and extract osm2po 4.7.7 and the map denmark-latest.osm.pbf.
to run (in cmd.exe) java -jar osm2po-core-4.7.7-signed.jar prefix=dk "C:\users\mgn.mes\documents\denmark-latest.osm.pbf" from the directory where i unpacked osm2po 4.7.7.

Then osm2po starts ad I'm prompted for accepting terms by hitting enter. I hit enter and the program ends. Then I try to enter http://localhost:8888/Osm2poService in the browser, but nothing i returned.

Do I need to change the .config file?

Thanks


texts: Better way to change the color of a graphic without affecting future graphics - ArcGIS JSAPI 3.X

texts:
I am trying to see if there is a better way to highlight a graphic in the map when you click on it. I currently have a layer of railroad tracks in the US and the project is to come up with a way to calculate the distance in the highlighted path. Right now, I have it set up so that you can click on a segment of track, click a button to save it, and it will store the information about the segment while keeping it highlighted. In that same list, I want to be able to click on a row in the list and have it change that track segment to a different color. As far as I am aware, the graphics are being stored on the main GraphicsLayer when using map.graphics.add(). The problem is that when it changes the color on the one graphic, any future graphics made get set to that same color. I'm able to reset the color on the graphic (without it visibly changing on the map) as a workaround, but there has to be a better way to do this.

var g = map.graphics.graphics[idx];
var originalColor = g.symbol.color;
map.graphics.remove(g);
g.symbol.setColor(new dojo.Color([0, 0, 0, 1]));
map.graphics.add(g);
g.symbol.setColor(originalColor);


texts: Is there an agreed standard to associate a time stamp with a Shapefile?

texts:
There are many reasons one might like to store temporal information in association with a shapefile. In my case I want a manually accessible timestamp so I can request temporally relative updates to a dataset.

Is there a proper way to do this in the shapefile standard? (or perhaps another QGIS accessible file standard)

I've had a look and the most appropriate place I could find seemed to me to be the geospatial metadata file (suffix .shp.xml) but I can't see how or if one can actually read or write a time stamp to that file.

The next thing I thought was to store a time stamp directly into the file name but that seems a little hacky for such an important factor.

To be clear I don't want to store a time stamp attribute for features in my shapefile as in GIS file format for space and time data? That would be straight forward but not appropriate. I want to store a timestamp for THE shapefile.

Using a timestamp assigned to the files by the OS (e.g.: file creation date or file access date) is also possible but that doesn't give me enough control and may be changed inappropriately by the OS later.

I use QGIS, but if there's no option available to QGIS, then writing my own XML (or similar) reader in Python is a viable solution too.


texts: Add date and time to attribute table using ArcPy insertcursor

texts:
ArcPy under ArcGIS Pro 2.3.2.

I am adding a Python datetime.datetime object to a shapefile attribute table using an arcpy insert cursor. The problem is that when I add the datetime instance, the date is preserved in the attribute table, but the time is set to 0:00.

Here is how I create the field:

arcpy.AddField_management(outfile_points, field_name='timestamp',
                          field_type='DATE',
                          field_is_nullable='NON_NULLABLE')

Here is how I add the data:

    with arcpy.da.InsertCursor(outfile_points, field_list) as cursor:
        for a_storm in history.all_storms.values():
            for a_rec in a_storm.record_dict.values():
                data_list = [None] * len(field_list)
                data_list[field_dict[designation]] = a_storm.designation
                data_list[field_dict[timestamp]] = a_rec.record_time <===
                #       Eliding irrelevant code
                cursor.insertRow(data_list)

The line with <==== is where I assign the value. Note, the type of .record_time is datetime.datetime.

When I check the value of data_list in this module as it runs, the time of day is preserved in the list. But when I view the data from ArcGIS Pro, the time of day is reset to 0:00.

Can anyone recommend what I should do differently at the <=== line?

The whole module is at
CreateHurricanePolylines.py


texts: Buffer around Point is elliptical, not circular

texts:
I'm trying to create a circle around a point using the code from this answer here. Adding below for better visibility as well
from functools import partial
import pyproj
from shapely.ops import transform
from shapely.geometry import Point

proj_wgs84 = pyproj.Proj('+proj=longlat +datum=WGS84')

def geodesic_point_buffer(lat, lon, km):
    # Azimuthal equidistant projection
    aeqd_proj = '+proj=aeqd +lat_0={lat} +lon_0={lon} +x_0=0 +y_0=0'
    project = partial(
        pyproj.transform,
        pyproj.Proj(aeqd_proj.format(lat=lat, lon=lon)),
        proj_wgs84)
    buf = Point(0, 0).buffer(km * 1000)  # distance in metres
    return transform(project, buf).exterior.coords[:]

# Example
b = geodesic_point_buffer(45.4, -75.7, 100.0)

print(b)
# [(-74.42290765358695, 45.39286001598599),
#  (-74.43102886629593, 45.304749544147974),
#  ...
# (-74.42290765358695, 45.392860015985995),
# (-74.42290765358695, 45.39286001598599)]

However when I try plot it on a map it looks more like an Ellipse and less like a circle. The coordinates I am trying to create the buffer around are (31.1482371,72.6859082). Is this some EPSG issue?


texts: Google Earth image projection transformation question

texts:
I have taken a JPEG image snapshot of an area on google earth. Within it I marked and record the coordinates of the four corners. I then added and georeferenced the image to a map I was already working in. However, the layers in my map are all NAD83. The Google earth image does not line up correctly. I want it to line up correctly with my other layers as I need to draw a new polygon using it.

How can I make the google earth image match my other layers?


texts: Understanding unit of area sizes calculated via rgeos::gArea?

texts:
I'm interested in finding the unit for the obtained area size of the polygons without refererring to the externally available documentation.

Example

A suite of shapefiles is imported as follows:

tmpDir <- tempdir()
tmpFle <- tempfile(fileext = ".zip")
download.file(url = "https://borders.ukdataservice.ac.uk/ukborders/easy_download/prebuilt/shape/infuse_rgn_2011_clipped.zip",
              destfile = tmpFle,
              method = "wget")
unzip(zipfile = tmpFle, exdir = tmpDir)
shpsUK <- rgdal::readOGR(dsn = tmpDir,
                         layer = tools::file_path_sans_ext(dir(tmpDir)[3]))

Area sizes

oldScipen <- getOption("scipen")
options(scipen = 999)
format(summary(rgeos::gArea(shpsUK, byid = TRUE)),
       big.mark = ",")
options(scipen = oldScipen)

which produces:

            Min.          1st Qu.           Median             Mean          3rd Qu.             Max.
" 1,573,530,067" "13,003,803,228" "15,411,419,165" "14,492,421,008" "19,085,426,561" "23,851,569,172"

Can I determine the units of the obtained figures without referring to external documentation that may be available on the sources shapefiles?


texts: How to delete a single Feature Class from a geodatabase using FME

texts:
I am trying to produce an FME script that is a replication of an existing python script in order to compare "apples to apples". The point of this project is to prove that FME can process "big data" faster than traditional methods.

Anyways, in order for me to try to accomplish this, I must
1. Not use Python in my FME script and
2. Try to match the logic in the original script verbatim (including unnecessary steps).

My question is then: How can I delete a single feature class from a geodatabase using FME? I understand that the recommended way is to use Python, but since this is not an option for me, and there is no transformer for deleting a single feature class, and there is no way that I can see to do it within a writer, then how can this be accomplished?


texts: Google Earth Engine Mosaic of Individual images

texts:
I would like to make a mosaic of 3 specific Date and Time Landsat 8 images corresponding to thermal band 10, where overlapping pixels between images in the mosaic will get the mean value.
I could select the images from the collections and display them, but separately not as a mosaic-image with the overlapping pixels as average, my code is below.
I would like to also export to drive the mosaic obtained as GeoTiFF file.
How could all this process be done?
Thanks a lot.
// point of interest correspond to LANDSAT tiles
var point1  = ee.Geometry.Point(35.0062, 32.7828);
var point2  = ee.Geometry.Point(34.8863, 32.0921);
var point3  = ee.Geometry.Point(35.0425, 30.6345);

// display points of interest
Map.addLayer(point1, {color: 'FF0000'}, 'Point1')
Map.addLayer(point2, {color: 'FF0000'}, 'Point2')
Map.addLayer(point3, {color: 'FF0000'}, 'Point3')

// image collection filter for each image
var dataset1 = ee.ImageCollection('LANDSAT/LC08/C01/T1')
                  .filterDate('2017-07-01', '2017-08-01')
                  .filterBounds(point1);

var dataset2 = ee.ImageCollection('LANDSAT/LC08/C01/T1')
                  .filterDate('2017-07-01', '2017-08-01')
                  .filterBounds(point2);

var dataset3 = ee.ImageCollection('LANDSAT/LC08/C01/T1')
                  .filterDate('2017-07-01', '2017-08-01')
                  .filterBounds(point3);

// explore image collection
print(dataset1)
print(dataset2)
print(dataset3)

// find images id
//dataset1 id  = LANDSAT/LC08/C01/T1/LC08_174037_20170731
//dataset2 id  = LANDSAT/LC08/C01/T1/LC08_174038_20170731
//dataset3 id  = LANDSAT/LC08/C01/T1/LC08_174039_20170731

// image of interest as variables
var dataset1 = ee.Image("LANDSAT/LC08/C01/T1/LC08_174037_20170731");
var dataset2 = ee.Image("LANDSAT/LC08/C01/T1/LC08_174038_20170731");
var dataset3 = ee.Image("LANDSAT/LC08/C01/T1/LC08_174039_20170731");

// select thermal bands

// thermal bands detaset 1
var thermalBand10_1 = dataset1.select(['B10']);
var thermalBand11_1 = dataset1.select(['B11']);

// thermal bands detaset 2
var thermalBand10_2 = dataset2.select(['B10']);
var thermalBand11_2 = dataset2.select(['B11']);

// thermal bands detaset 3
var thermalBand10_3 = dataset3.select(['B10']);
var thermalBand11_2 = dataset3.select(['B11']);

// display thermal  bands
Map.addLayer(thermalBand10_1,{min:0,max:65535},'thermalBand10_1');
Map.addLayer(thermalBand10_2,{min:0,max:65535},'thermalBand10_2');
Map.addLayer(thermalBand10_3,{min:0,max:65535},'thermalBand10_3');


texts: Is there a way to show the direction of a route?

texts: I use pgrouting and openstreetmaps. I would like to know if there is a way to show the direction of the route on the map (e.g. something like this ---->---- or any other way, or idea)


texts: ogr2ogr failed to reproject shapefile error exceeds limits

texts:
I cannot get correct transformation from WGS84 to the Hungarian EOV EPSG:23700 (Hungarian Datum 72) using ogr2ogr.

I am using ogr2ogr from the fwtools shell. I have a shapefile in GCS WGS84 (epsg 4326) and I want to convert it to EOV, a projected coordinate system using the Hungarian Datum 1972.
The contents of the prj file are:

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]

When I run this command:

ogr2ogr -t_srs EPSG:23700 EOV_HD72.shp WGS84.shp

I get failure and the following error:
Error 1:latitude or longitude exceeds limits. Failed to reproject feature 0 (geometry probably out of source or destination SRS).

Any suggestions on why this re-projection fails?


texts: ArcGIS ModelBuilder iterator reports only last input shapefile?

texts:
Using Arc 10.3.1, if I manually enter the names of multiple shapefiles (all from the same folder) into the Check Geometry tool, the (single) output table lists the errors that exist in all of the shapefiles.

That is, only one output table is created, even if multiple input shapefiles have errors.  For example, if two shapefiles have errors, both will be listed in the single output table. Here's a screen print of the input:

However, if I iterate those same shapefiles (using the Feature Classes option), the output table only lists the last shapefile that went through the iteration.  Here's a screenshot of the model:

I suspect that the iterator is running Check Geometry not once, but three times, with a single input shapefile going through each iteration, overwriting the previous output.

How do I get the iterator to list all of the shapefiles with errors, either in separate output tables or aggregated in a single output table?


texts: Calculate Polygon Centerline with a LeastCostPath approach

texts:
I'm working on deriving the centerline of a polygon, similar to this question or this one, using a different approach. I'm 90% done with the task, but cannot reach the final step. My approach is to calculate a LeastCostPath through the channel center.

In order to do this, I first calculated a surface within the polygon with the distance to the polygon edge (see image 2). In order to get values from 0 - 1 with low values being in the polygon center, I used the following operation:

rastervalueNormalizedInverted = 1-(eucideanDistance/maxRasterValue)

Now I can just calculate the least cost path between the start and end of my polygon and I get a beautiful line without dangles and without the issue of having to specify threshold values (see image 3).

I have the feeling that I'm 90% done at this point. All I have to do now is consider the side channel in polygon, which were not considered by the initial LeastCostPath calculation.

I am completely stuck at this point! I've tried various approaches to get the diverging and converging line for the two side arms seen in the image, with no success. Can anybody provide a clue in which direction I could start looking?

Here are the images to illustrate:

The original Polygon:

Normalized, inverted euclidean distance to the polygon edge:

Least Cost Path from the start to endpoint of the channel returns a nice line:


texts: How to change/identify Linear units in Near Analysis Attributes Table

texts:
I need to get the closest distance from multiple points of a multipoint feature to a polyline that delineates a land contour. To do this, I used the Near(analysis) feature. It worked, but I can't tell what units it has produced. I right clicked on the layer, but it doesn't have a Linear Units category under the Source Tab of Properties. Also opened up a .prj of the layer to see what the units are -also unspecified. Checked the dataframe properties and it said my Linear Units are in meters (which is what I want them to be in), but that is not correct given the measurements the Near tool produced.

Any suggestions for how to calculate the distance in meters OR change whatever units are in the NEAR_DIST column of my attributes table?


texts: ArcGis Layers Not Loading on Xamarin.iOS App - Invalid JSON Exception

texts:
Am currently using the Esri.ArcGISRuntime.Xamarin.iOS v100.0.0 version on the following version of Xamarin Studio for Mac(6.1.2 build 44) with Xamarin.IOS Version: 10.2.1.5 & XCode Version 8.1.

I get the following Exception when attempting to load a Feature Layer by doing a LoadAsync on the ServiceFeatureTable queried from a given URL. “Esri.ArcGISRuntimeException: Invalid JSON at System RuntimeExceptionServices.Exception” following is the code used.

// create a service feature table of census block groups
var blockGroupsUri = new Uri(Config.MapServerUrl + "/" + layer.ID);
var featureTable = new ServiceFeatureTable(blockGroupsUri);
featureTable.FeatureRequestMode =
FeatureRequestMode.OnInteractionCache;

// load the feature table
await featureTable.LoadAsync();
featureTables.Add(featureTable);
var lyr = new FeatureLayer(featureTable);
myMap.OperationalLayers.Add(lyr);

This only happens for certain layers, have tried validating the JSON returned from querying the layer, there seems to be no issue there. The attributes on the Layers that work and don’t work are pretty much the same.

But the thing to note here is that, the above layer in question loads up on the map fine when its saved in a GeoDatabase (for offline querying), and when queried using the GeodatabaseFeatureTable subsequently.


texts:
I am developing with ArcGIS REST API Vehicle Routing Problem service.

When I sent the request using HTTP GET method, I got the 414 error as the URI was too long. (I requested with 11 routes and 238 orders data.)

So I changed to POST method. Then I could get jobId.
However, when I used the jobId to request the solution with output parameters (out_stops, out_routes, out_directions, out_unassigned_stops and solve_succeeded), I got the errorMessage

  "Unable to complete operation"

Then I POSTed the request with less data that could be solved successfully using GET method before, but I still got the errorMessage "Unable to complete operation".

Does it mean ArcGIS REST API Vehicle Routing Problem service does not accept HTTP POST method?

Here is my HTTP request URL:

  https://logistics.arcgis.com/arcgis/rest/services/World/VehicleRoutingProblem/GPServer/SolveVehicleRoutingProblem/submitJob?token=&orders={...}&order_pairs={...}&depots={...}&routes={...}&breaks=&time_units=Seconds&distance_units=Kilometers&uturn_policy=NO_UTURNS&default_date=1465394400000&time_window_factor=High&f=json

I use java Jersey2.x client to do the POST request.

Client client = ClientBuilder.newClient();
WebTarget target = client.target("https://logistics.arcgis.com/arcgis/rest/services/World/VehicleRoutingProblem/GPServer/SolveVehicleRoutingProblem/")
            .path("submitJob");
    Form form = new Form();
    form.param("token", problem.getToken());
    form.param("orders", ordersRequestParameter);
    form.param("order_pairs", orderPairsRequestParameter);
    form.param("depots", depotsRequestParameter);
    form.param("routes", routesRequestParameter);
    form.param("breaks", breaksRequestParameter);
    form.param("time_units", "Seconds");
    form.param("distance_units", "Kilometers");
    form.param("uturn_policy", "NO_UTURNS");
    form.param("default_date", "1465394400000");
    form.param("time_window_factor", "High");
    form.param("f", "json");

    String jobIdResponse = target.request(MediaType.APPLICATION_JSON_TYPE)
                                .post(Entity.entity(form, MediaType.APPLICATION_FORM_URLENCODED_TYPE), String.class);


texts: Calculating multiple fields using ArcPy?

texts:
I have found a script that can add multiple fields to my existing shapefefile. Now, I need to add to my existing script so that it will calculate the new fields that have been added. In most cases, these calculations will derive from existing values in other fields ex. Crdrt1 = [Field1], Crdrt_min1 = [Field2]. I know that I need to use "fields.arcpy.CalculateField_management" but not sure how or where to add to existing script to get it to loop through and calculate.

import arcpy

#Select Geodatabase
arcpy.env.workspace = "C:.. FOLDER\CASINGS_DATA\casings.gdb"
#Select FeatureClass
shapefile = "Casings"

#List of fields to add, all with same type, precision, and scale
#Fields added to attribute table in reverse order
fields = {
    "Crdrt1": "[DES]",
    "Crdrt_min1": "[GRADE]",
    "Testing": "[NAME]"
}
alias = {
    "Crdrt1": "ok",
    "Crdrt_min1": "test"
}

existing_fields = arcpy.Describe(shapefile).fields

for field in fields:
    if not field in existing_fields:
        arcpy.AddField_management(shapefile, field, "TEXT", 0, 0, 50, alias [aliass])
    arcpy.CalculateField_management(shapefile, field, fields[field])


texts: Import shp to Postgis using Python and ogr

texts:
I just export a Postgis Table to shp using
this tips
but I'm not able to import a shp to Postgis using the same library(ogr).
Any idea?
Thanks a lot
f.


texts: Why do I get 'KeyError: 0' when using Line Intersection tool?

texts:
I'm using Quantum GIS version 1.8.0-Lisboa and trying to create intersection point using Intersection tool. (Vector -> Analysis tools -> Line Intersections )

But when I do that for my shapefile, It pops up the following error in the dialog box:

An error has occured while executing Python code:

Traceback (most recent call last):
  File "/usr/share/qgis/python/plugins/fTools/tools/doIntersectLines.py", line 96, in accept
    self.compute(line1, line2, field1, field2, outPath, self.progressBar)
  File "/usr/share/qgis/python/plugins/fTools/tools/doIntersectLines.py", line 174, in compute
    outFeat.addAttribute(0, atMap1[index1])
KeyError: 0

Python version:
2.7.3 (default, Aug  1 2012, 05:27:35)
[GCC 4.6.3]

QGIS version:
1.8.0-Lisboa Lisboa, exported

Python path: ['/usr/share/qgis/python', '/home/sibi/.qgis//python', '/home/sibi/.qgis//python/plugins', '/usr/share/qgis/python/plugins', '/usr/local/lib/python2.7/dist-packages/identi-0.1.3-py2.7.egg', '/usr/local/lib/python2.7/dist-packages/pyipmsg-0.0.1.0-py2.7.egg', '/usr/local/lib/python2.7/dist-packages', '/usr/local/lib/python2.7/dist-packages/rope-0.9.4-py2.7.egg', '/usr/local/lib/python2.7/dist-packages/ropemode-0.2-py2.7.egg', '/usr/local/lib/python2.7/dist-packages/ropemacs-0.7-py2.7.egg', '/usr/lib/python2.7', '/usr/lib/python2.7/plat-linux2', '/usr/lib/python2.7/lib-tk', '/usr/lib/python2.7/lib-old', '/usr/lib/python2.7/lib-dynload', '/usr/local/lib/python2.7/dist-packages', '/usr/lib/python2.7/dist-packages', '/usr/lib/python2.7/dist-packages/PIL', '/usr/lib/python2.7/dist-packages/gst-0.10', '/usr/lib/python2.7/dist-packages/gtk-2.0', '/usr/lib/pymodules/python2.7', '/usr/lib/python2.7/dist-packages/ubuntu-sso-client', '/usr/lib/python2.7/dist-packages/ubuntuone-client', '/usr/lib/python2.7/dist-packages/ubuntuone-control-panel', '/usr/lib/python2.7/dist-packages/ubuntuone-couch', '/usr/lib/python2.7/dist-packages/ubuntuone-installer', '/usr/lib/python2.7/dist-packages/ubuntuone-storage-protocol', '/usr/lib/python2.7/dist-packages/wx-2.8-gtk2-unicode', '/usr/share/qgis/python/plugins/fTools/tools']

Is it some bug in the Quantum GIS program or I'm missing something here?


texts: How to fix: update of view in private qgis.db failed error

texts:
After a change in permissions, I started to get this error message:

  Update of view in private qgis.db failed.
  table vw_srs already exists

I know that changing the ownership of the db or deleting it and letting QGIS create a new one could solve the problem, but in my case, this is not an option. An option would be not to update qgis.db. Does anybody know, how to "switch off" this checking and updating?

And why is this db checked every time I start QGIS anyway?


texts: GEE large export of pixels using shapefile polygons: Is it possible to automate a division of the polygons asset?

texts:
I need to export pixels from a processed image of Sentinel-2 using GEE Javascript with a 10 pixel size. The export is for pixels within 120,000 land parcel polygons.
I am way above the limit for pixel export (I even increased the MaxPixels). I am trying to divide the original shapefile with the polygons into parts to import to GEE one piece at a time. But, it is very tedious work to do that one by one and run the code and export accordingly.
Is it possible to work with the original asset with all polygons and iterate over subsets of it in order to export data as csv automatically for each subset? I can then merge the subsets once all are exported.


texts: Continuous colour gradient within line segments in a vector layer

texts:
I am new to QGIS (using version 3.13) and my main objective with it is easy to explain: creating a line vector (multiple points, representing a river) in which the colour of the line changes based on a Continuous Attribute (colour gradient within line segments). It is about plotting the changes in the concentration of chemicals in a river (from hereon the "proxy"). I tried with R, SAGA and QGIS. I think what I'm trying is actually not possible or difficult so help would be appreciated:

In QGIS I tried the following:
1. Layer/Add layer/Delimited Text Layer, attributes: Long/Lat/River(name)/SequenceNumber in the river (1-7 points)/pollutant concentration (continuous variable, not integer). Gives points plot.
2. Convert points to a line, with Points to Path and Saga-version Points-to-Line(s), same output, order field SequenceNumber, Group Field RIver. Everything OK, a line vector with lines of different colour.

See image in which I try to explain some things.

The problem is that I don't know how to couple the line vector to the attribute Proxy (continuous variable). In the attributes form of the line vector (symbolism/colour categorized option), I cannot select this proxy (only the river name, a discrete parameter), which is normal. I would need the "graduated colour" of the vector but this brings to me a "value" field that I don't know how to handle. Anyway, the desired Proxy is not in the attributes form of the line vector, so it seems that this is not the way to go.

On page (https://docs.qgis.org/testing/en/docs/gentle_gis_introduction/vector_attribute_data.html) it is explained how to use continuous colour symbols. Here Symbology includes Legend Type: Continuous Color and allows to select the desired proxy (in the example, HEIGHT). I tried using an earlier QGIS version (2.12) but I was not able to reproduce this. Maybe this is because the example uses a colour gradient to indicate contour lines where hence every line has only one colour. It is stated that this approach is useful as "The GIS Application will use a numerical attribute value from a feature (e.g. contour heights or pollution levels in a stream) to decide which colour to use." but the example given (HEIGHTS) has a really different nature that the pollution example, so it doesn't help.

In a probably very relevant discussion on this forum "Styling line layer based on two attributes in QGIS ? (1 discrete and 1 continuous)" the answer of user MarHoff is: "If you only got one discrete attribute (road type [IN MY CASE RIVER NAME]) and one pseudo-continuous attribute (number of lane) I think you don't even need to use the rule-based renderer which is the way to go if you got more complex situations."

So, the way to go for me is probably the "rule-based renderer" but I don't know how to do this or prepare the data so that it can be done.

And with the Geometry Generator (creating polygons from the line) I think it is not possible to link colours to the Proxy.

I would very much appreciate any assistance.

With kind regards, Joeri


texts: Given two polygons, find the the one with more detailed/accurate shoreline?

texts:
I have a polygons for lakes from different sources. Some polygons are more detailed in one source, and some are more detailed in the other.

Given two polygon of the same lake, I want to find which of these two polygons captures the lake shoreline more accurately. In the below example you can see the second one captures the shoreline boundaries better than the first one.


texts:
No GIS Certification Programs in my area...Are Online Certs worth their weight?

Background Information:

Currently a Computer Science student, with 3 semesters remaining. My college doesn't offer any GIS courses, and there aren't any local certificate programs available. Would it be worth trying an Online Certification Program?

Or should I try to teach myself as much information as possible and try to get into a position after graduation without any formal GIS info on my resume?


texts: Obtaining sublayer of layer for calculating shortest path algoritm in PyQGIS?

texts:
I am trying to calculate the shortest path over a series of points via the road layer. For this, I am using the shortest path algorithm in the network analysis tool via the Python console.

However, the algorithm takes a while to be solved because the area of the layer roads that I'm using is quite large. I was wondering if there is a faster way to do this, maybe obtaining a sublayer of the road graph, and then running the algorithm over that sublayer.

To give an example, for the figure in below I'd like to calculate the shortest path between the green and orange point, over the area marked in red and not the entire road layer.

Is there a simple way to obtain the area marked in red, any suggestions on how can I implement it?


texts:
I'm currently writing some code which creates a new layer from the selected features of another layer.
layer_4 = QgsProject.instance().mapLayersByName('layer_1')[0]
iface.setActiveLayer(layer_4)
layer_2 = iface.activeLayer()
selected_2 = layer_4.selectedFeatures()
crs_2 = QgsProject.instance().crs().authid()
selection = QgsVectorLayer("Point?crs=" + str(crs_2), "selection", "memory")
new_features_2 = []

for feat in selected_2:
    new_features_2.append(feat)

selection.dataProvider().addFeatures(new_features_2)
QgsProject.instance().addMapLayer(selection)

This seems to work, however the points are not displayed in the new layer and when I'm checking the attribute table it says that all elements are filtered. Why? And how do I remove the filter?


texts: Failure in printing OSM as a basemap through GeoServer Mapfish Plugin

texts:
I am trying to print OSM as a basemap through GeoServer. The WMS typed layer gets printed but the OSM doesn't. The projection coordinates I am using is ESPG:4236. But I found that the majority of OSM basemap through GeoServer uses ESPG:900913 coordinates. However I somehow managed to pull some data from online research for 4326 cordinates as shown below. I am not sure whether I can do like this or not.

spec.json

{
  "layout":size,
  "srs":"EPSG:4326",
  "units":"degrees",
  "dpi":dpi,
  "layers":[
    {
      "baseURL":"http://c.tile.openstreetmap.org",
      "extension":"png",
      "opacity":1,
      "singleTile":false,
      "type":"OSM",
      "maxExtent":[80.05844110511946,26.34796713166333,88.20152186778043,30.447429596886728],
      "tileSize":[256,256],
      "resolutions":[5.625,2.813,1.406,0.703,0.352,0.176,0.088,0.044,0.022,0.011,0.005,0.003,0.001,0.0005]
    },
    {
      "baseURL":"http://localhost:9090/geoserver/BernHardt_Practical/wms",
      "opacity":1.0,
      "singleTile":false,
      "type":"WMS",
      "layers":["Nepal"],
      "format":"image/png"
    }
  ],
  "pages":[
    {
      "bbox":[80.05844110511946, 26.34796713166333, 88.20152186778043, 30.447429596886728],
      "scale":scale,
      "strictEpsg4326": true
    }
  ]
};

I am not receiving any errors, but no OSM basemap layer is being drawn except the WMS layer. This is the output I am getting.

Imcomplete Output:

The output should have been like this,

Do we have any other proper resolutions?


texts:
I have created an OpenStreeMap layer of the indoor of a bulding by following this tutorial: https://wiki.openstreetmap.org/wiki/IndoorOSM, since that is only a proposed tagging schema for indoor mapping when I have tried to load it into QGIS the only thing I got was the building and not the indoor nodes.

How can I turn those rooms into something that can be interpreted by more advanced tools like QGIS or a postgis database?

PS: My project consists on providing routes between two different rooms of a building, any help or suggestion on how to do it or how to render osm's indoor layers is welcome.


texts: Fetching the nearest location of points, while accounting for bodies of water

texts:
I've got a database of points (in this case, schools) and in my app, users search for the nearest ones to them. Under the hood, we're currently using ElasticSearch to filter by latlng (using Geo Distance, which gets the distance as the crow flies. For the majority of places, this works fine, but in some coastal areas, this will pick up places that are impossible to get to, in the example below, a radius of 20 miles will pick up schools in Weston-Super-Mare, in reality 55 miles:

I initially decided to use the Google Maps Distance Matrix API to filter my inital as the crow flies search, but there's a limit of 25 destinations per query, and as the requests will be dynamic and user-facing, it's not practical to parcel these requests up into small pieces and pop in a background job.

Is there any way to carry out these calculations while accounting for bodies of water on a database level? The schools are stored in a Postgres database, so I thoughts about using PostGIS and some kind of point in polygon query, but I have no idea where to start looking.


texts: QGIS Standalone C++ Application with Qt: "The program has unexpectedly finished"

texts:
I'm trying to develop an Standalone Application with Qt and QGIS API.

There are several code examples that should work, but they don't.

The Problem is that everytime I use "QgsApplication" the application Output says: "The program has unexpectedly finished."

Even this simple code does not work:

#include "mainwindow.h"
#include <QApplication>

#include "qgsapplication.h"

int main(int argc, char *argv[])
{
    QgsApplication app(argc, argv, true);
}

My Project.pro file:

QT       += core gui

greaterThan(QT_MAJOR_VERSION, 4): QT += widgets

TARGET = Project
TEMPLATE = app

SOURCES += main.cpp\
    mainwindow.cpp

HEADERS  += mainwindow.h

FORMS    += mainwindow.ui

LIBS += -L"C:/QGIS Building/lib" -lqgis_core -lqgis_gui

INCLUDEPATH += "C:/QGIS Building/include"

DEFINES +=      CORE_EXPORT=""
DEFINES +=      GUI_EXPORT=""


texts: Leaflet GeoJSON with Jquery UI autocomplete

texts:
How do a autocomplete when I search for something?

I have this:

<form>

        <input type="text" name="busca" id="txtSearch" class="form-control" placeholder="Buscar..." />

        <button class="btn btn-default" id="btnSearch"></button>
</form>

Nothing more, just a form with a input text and a buttom. When write anything I need that autocomplete, when find the word click in buttom.

And my Js have:

stComerciaisLayer = L.geoJSON(setoresComerciais, {
                style: function (feature) {
                    ...
                },
                onEachFeature: onEachFeature
            }).addTo(map);

The variable stComercialLayer is a object that generate the mapping layers. setoresComerciais is the same variable in my json (geoSetoresComerciais)

    setoresComerciais={
      "type": "FeatureCollection",
      "totalFeatures": 116,
      "features": [
        {
         (...)

and

$("#txtSearch").autocomplete({
                source: setoresComerciais
            });

is the autocomplete, but the source, I need a function or a array for work, I think.

The problem is that stComerciaisLayer is a Object and the source just accept type Array, String or Function, conform jquery UI's documentation :(
What i need do?

Resuming: I just need that the autocomplete works! My code is on the git:

https://github.com/eltonsantos/leaflet-tests/tree/master/teste2


texts: ArcGIS 10.1 Aggregate Points cannot create output

texts:
I'm trying to use the tool Aggregate Points with output created using in_memory but the tool keeps failing with: "ERROR 000210 Cannot create output \Output_Tbl"

The tool works fine if I use a GDB feature class as output. Does anyone know how to use output for in_memory with this tool?


texts:
I need to use qgis:union in pyqgis algorithm and to merge 2 layers, layer1 and layer2 and "see" the output layer
But i don´t understand how to work with the "output layer"

this "output layer" have the layer1 and layer 2 features?
how to save a "output layer"?


texts:
I have a problem with ArcGIS 10.1:
I have two TINs and I want to execute Cut Fill (3D Analyst) tool :

Data:

Adding data to ArcScene:

Tool:

Result: Problem

I did the same thing in ArcGIS 9.3 and it's works, I don't know where is the problem, can somebody help?


texts:
I am trying to generate an atlas with python in qgis 3.0. I have :

projectInstance= QgsProject.instance()
projectLayoutManager = projectInstance.layoutManager()
for comp in projectLayoutManager.printLayouts():
    if comp.atlas().enabled():
        comp.atlas().beginRender()
        nbobj = comp.atlas().updateFeatures()
        comp.atlas().first()
        for i in range(0, comp.atlas().count ()):
            comp.atlas().refreshCurrentFeature ()
            exporter =  QgsLayoutExporter(comp.atlas().layout () )
            pdf_settings = exporter.ImageExportSettings () #dpi?
            exporter.exportToImage (<filename>, pdf_settings)
        comp.atlas().endRender()

But I would like the filename to have some attributes of the atlas feature.

How do I get the value of the attribute?


texts:
When I try to clip a raster layer with a polygon vector layer. QGIS throws and error:

  ERROR 1: The transformation is already "north up" or a transformation
  between pixel/line and georeferenced coordinates cannot be computed
  for C:\Users\bhoke\Desktop\gisFiles\uydu1.tif. There is no affine
  transformation and no GCPs. Specify transformation option
  SRC_METHOD=NO_GEOTRANSFORM to bypass this check.

I am new to GIS and I don't know what it means. CRS of both vector and raster layers are same. I tried -to src_method=no_geotransform but it results with a blank image also warns about SRS of vector layer.

What can be the problem here?


texts: GeoServer GeoCSS generates too many symbols

texts:
Can anyone help me with the GeoCSS below, I don't understand why this generates 5 legend categories when I was expecting only 4 categories?

/* @title Monitor */
* {
    mark: symbol(circle);
    mark-size: 12px;
    :mark {
        stroke: black;
        fill: #00CC22;
    };
}

/* @title Warning */
[repair = 0][reading >= site_normal*1.4][reading < site_normal*1.8]:mark {
    fill: #cc9c00;
}

/* @title Alert */
[repair = 0][reading >= site_normal*1.8]:mark {
    fill: #cc0000;
}

/* @title Under Repair */
[repair = 1]:mark {
    fill: #ffffff;
}

That CSS generates this legend, I was expecting only 4 symbols and I don't want the "Monitor, Warning, Alert" symbol.

Its almost a copy of this example, so I'm not sure why there is the extra symbol: http://docs.geoserver.org/latest/en/user/styling/css/tutorial.html#refining-the-style


texts: pyqgis - post-processing of layer output of QgsProcessingAlgorithm

texts:
While extending a plugin with a Processing Provider and a set of algorithms I mostly followed the defaults of Plugin Builder 3. So, the basic code is smth along the lines of

class ORSisochronesAlgo(QgsProcessingAlgorithm):

    def initAlgorithm(self, configuration, p_str=None, Any=None, *args, **kwargs):

        # add all input parameters

        self.addParameter(
            QgsProcessingParameterFeatureSink(
                'OUTPUT',
                'Output layer'
            )
        )

    def processAlgorithm(self, parameters, context, feedback):

        source = self.parameterAsSource(parameters, self.INPUT, context)
        (sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT,
                context, source.fields(), source.wkbType(), source.sourceCrs())

        features = source.getFeatures()

        for feature in features:
            new_feature = QgsFeature()

            # Some API calls and processing to populate new feature

            sink.addFeature(new_feature, QgsFeatureSink.FastInsert)

        return {self.OUTPUT: dest_id}

So, the output layer is automatically handled and added to the map. Very nice!

But now: How do I post-process that output layer within the Processing Algorithm? Sounds weird maybe, but here are 2 scenarios:

Depending on the value of a Boolean input parameter, I'd like to be able to calculate geometric differences of the output features. For this to work, the layer has to be built first. And then handled in some sort of post-processing.
Also, I'd like to be able to customize the style of the output layer.

There is the postProcessAlgorithm() method for QgsProcessingAlgorithm, but I'm not so sure how to work with this.. It doesn't seem to take parameters as input, no idea how to access the output layer.

Being able to do this would make the plugin so much more user-friendly.


texts:
I am using leaflet for creating maps in web and geoserver for serving maps as WMS layer.Leaflet has box zoom-in tool when we click shift and drag on the map.

Similarly I need box zoom-out in leaflet like in QGIS

Is there any way to achieve this using leaflet?


texts:
I'm trying to take a field calculator equation and do it in arcpy directly through IDLE, but am running into a problem where the process runs correctly, but the results in the field show up as 0. I've tried a variety of ways of trouble shooting this, including a massive amount of time and energy trying it in field calculator itself, and playing around with the field types , but not luck so far. Any constructive advice would be much appreciated.

import arcpy, sys

arcpy.env.workspace = "C:\Temp"

shapefile = "USCancer2000.shp"

field_name = "Check8"

print "start your engines"

#Add the field
arcpy.AddField_management(shapefile,field_name,"FLOAT", "5","5","","","NULLABLE","NON_REQUIRED","")

cursor = arcpy.UpdateCursor(shapefile)

for row in cursor:
    #Get the value of each of the fields, and put them into variables
    c1 = row.getValue("Cnt1")
    c2 = row.getValue("Cnt2")
    c3 = row.getValue("Cnt3")
    p1 = row.getValue("Pop1")
    p2 = row.getValue("Pop2")
    p3 = row.getValue("Pop3")
    #check for missing values, assuming that non-missing values are all greater than or equal to zero
    try:

note I've tried changing this field_value=0 part to 7 to see if it populates the field with 7 but it doesn't

        if((min(c1,c2,c3,p1,p2,p3) < 0)):
            field_value = 0
        else:

have tried putting in float in both cases, and the .0 attached to the 100,000 to try to #force float

             field_value= 100000.0* (float(c1)/float(p1))
    except:

note I've tried changing this field_value=0 part to 8 to see if it populates the field with 8 but it doesn't

###If something goes wrong, set the value to 0 still
        field_value = 0

#Set the value of field for just this row; note have tried inserting the cursor. #name for both ".setValue", and for just .updateRow, and also just for one. Still #doesn't work
    cursor.setValue(field_name,field_value)
    cursor.updateRow(row)
del row
del cursor
print "done"


texts: How do I project OSM tiles to a custom coordinate system?

texts:
Background: I'm planning to use OSM data to generate tiled maps and then serve them up as an internal service (hosting service still TBD). I want to be able to consume the tiled maps in both ArcMAP and QGIS. For my day job, I work with ESRI products. For my night job, I like QGIS.

I've tinkered with Maperitive and TileMill for generating some tiles. Both seem to be decent products. I have no preference at this time but I haven't beat on either product long enough to find weaknesses. One weakness I may have found, TileMill outputs to Web Mercator.

My question: How do I generate tiles from OSM Data and project them into a map projection created for my state? The map projection is 3071 - Wisconsin Transverse Mercator.


texts: How to apply a filter to multiple vector layers in QGIS?

texts:
I am trying to apply an attribute filter to multiple layers in QGIS at once. The filter is to be applied to .osm data brought in through PostGIS to QGIS.

I have tried selecting multiple layers -> right-click --> filter but it only applies it to the last layer I selected of the group. I want to be able to apply a filter through a query like "osm_user" = 'donlaser' to multiple layers at once, rather than having to repeat the process many times.

Is there a standard way to do this within QGIS that I am missing? Is there a plugin that would provide this functionality?

A Python-based solution would also be useful for me.


texts: Calculating and setting length value for each row in a shapefile returns empty column

texts:
I have a Linestrings' shapefile to which I added a column named "length". I want to set the length of each row (geometry) in that column. Here's my code :

import geopandas as gpd
from shapely.geometry import *
from shapely.ops import transform
from functools import partial
import pyproj

def addLengthField(shapefile, output):
shp = gpd.GeoDataFrame.from_file(shapefile)
shp['length'] = ''
shp.to_file(output)

addLengthField('./Shapefiles/curbs.shp', './Outputs/curbsWithLength.shp')

def setLengthValue(shapefile, output):
shp = gpd.GeoDataFrame.from_file(shapefile)
for index, row in shp.iterrows():
    geom = shape(row['geometry'])
   # since the result of geom.length is not in Meters, I wanted to make a transformation
    # Geometry transform partial function based on pyproj.transform
    project = partial(
        pyproj.transform,
        pyproj.Proj(init='EPSG:4152'),
        pyproj.Proj(init='EPSG:4152'))
    projectedCurbs = transform(project, geom)

    row['length'] = int(projectedCurbs.length)
    print row['length'] # length values are well displayed
    # Save the final output in another shapefile
shp.to_file(output)

setLengthValue('./Outputs/curbsWithLength.shp','.Outputs/finaloutput.shp')

When I run the script the printed values are well displayed but when I open the output shapefile in QGIS the length's column is empty!
Any ideas on what could be wrong?


texts: How do I set the projection and coordinate system for a TIFF file?

texts:
I have a topographic map of California (obtained from CAL-Atlas) as a TIFF file. I tried importing it into QGis, but it displays the map using a pixel-based coordinate system rather than latitude/longitude (the .prj file indicates that it uses NAD83 Albers projection).

How can I rectify the coordinate system for this raster file in QGis? I'm using QGis 1.8 for this project, but also have 2.0 if necessary.

What if I am using QGis >=2.7v what should I do for the same..?


texts: ArcGIS Javascript API - Looking For Styled Navigation Plugin/Widget

texts:
I am looking for a fancy navigation plugin to use with the ArcGIS Javascript API.  Something along the lines of the Navigation Widget (see image below) available in the Flex Viewer would be great.

Anyone have any ideas/suggestions?


texts: populating one field value based on another field value

texts:
I have a feature class in a personal geodatabase with 2 fields. They both have domains set up to ease the editing process for my users. If my WaterType field is set as "Not Irrigated" (or code 4) then I need my IrrigType to be set as "Not Irrigated" (code 0). I'm not sure if I want to use the description or the code for the domains.

I would also like to set this up so that it runs automatically...so that I don't need to calculate field and run the Python script every time.

Any suggestions? I tried the following code to no avail:

def calIrrigType(WaterType, IrrigType):
  if WaterType=='Not Irrigated':
    return 0
 else:
  return IrrigType


texts: Downloading shapefile from OpenStreetMap?

texts:
I am trying to download the shapefile (orange line) from
https://www.openstreetmap.org/relation/271888#map=12/48.7810/9.2247

Can you explain how I can extract the shapefile from here?


texts:
I am trying to pass the results of a list as a parameter in a SelectLayerByAttribute expression to output and excel spreadsheet of the attribute table. Here is the code, the results create an empty excel spreadsheet.  I am little lost.

# Make feature layer
   arcpy.MakeFeatureLayer_management(OwnerList, "OwnerTable")
   with arcpy.da.SearchCursor("OwnerTable", "ClubName") as cursor:
        myList = sorted({str(row[0]) for row in cursor})

        # Create Club lists
        for i in myList:
            expression = '"ClubName" = \'str(i)\' '
            arcpy.SelectLayerByAttribute_management("OwnerTable", "NEW_SELECTION", expression)
            arcpy.TableToExcel_conversion("OwnerTable", out_xls + "_" + i + ".xls")

Solution:
 # Make feature layer
   arcpy.MakeFeatureLayer_management(OwnerList, "OwnerTable")
   with arcpy.da.SearchCursor("OwnerTable", "ClubName") as cursor:
        myList = sorted({str(row[0]) for row in cursor})
        print myList

    # Create Club lists
    for i in myList:
        expression = ' "ClubName" = \'%s\' ' % i
        print expression
        print arcpy.AddMessage(expression)
        arcpy.SelectLayerByAttribute_management("OwnerTable", "NEW_SELECTION", expression )
        arcpy.TableToExcel_conversion("OwnerTable", out_xls + "_" + i + ".xls")


texts: Arcpy Append Management Function Copying Blank Records

texts:
Below is the code I'm struggling with.

import arcpy

import datetime

#Create copy of GDB 2.0 Template
date = time.strftime('%Y%m%d%H%M%S')

emptyLLDGDB = "\\\\Bvweb104draco\\e\\08_Scripts\\GDB21_LLD_VendorChanges_05Sep_V5.gdb\\GDB21_LLD_VendorChanges_05Sep_V5.gdb"

newGDB = "\\\\Bvweb104draco\\e\\08_Scripts\\Scratch\\Survey" + date

arcpy.Copy_management*(emptyLLDGDB, newGDB)

#Turn ADDRESS feature class from service into a table to be consumed by the append function
ADDRESS = r'\\Bvweb104draco\e\03_MapServices\SDE_CONNECTIONS\bvtampaprod.sde\bvtampaprod.sde.Telecom\bvtampaprod.sde.ADDRESS'

addressTableName = "Table" + date

Address_Table = "\\\\Bvweb104draco\\e\\08_Scripts\\Scratch"

arcpy.TableToTable_conversion*("ADDRESS", Address_Table, addressTableName)

#Append ADDRESS table to survey attribute table in copied GDB 2.0

Survey_Attribute = "\\\\Bvweb104draco\\e\\08_Scripts\\Scratch\\Survey" + date + ".gdb\\Survey_Attributes"

arcpy.Append_management(addressTableName, Survey_Attributes, "NO_TEST", "LOCATION_TYPE \"LOCATION_TYPE\" true true false 50 Text 0 0 ,First,#;BLDG_TYPE \"BLDG_TYPE\" true true false 50 Text 0 0 ,First,#,Table20180918141600,bv_loc_typ,-1,-1;BLDG_SUB_TYPE \"BLDG_SUB_TYPE\" true true false 50 Text 0 0 ,First,#,Table20180918141600,bv_bldg_su,-1,-1;BLDG_LAT \"BLDG_LAT\" true true false 50 Text 0 0 ,First,#,Table20180918141600,latitude,-1,-1;BLDG_LONG \"BLDG_LONG\" true true false 50 Text 0 0 ,First,#,Table20180918141600,longitude,-1,-1;ALLOCATION_STATUS \"ALLOCATION_STATUS\" true true false 50 Text 0 0 ,First,#,Table20180918141600,bv_alloc_s,-1,-1;NOBUILD_REASON \"NOBUILD_REASON\" true true false 50 Text 0 0 ,First,#,Table20180918141600,bv_nobuild,-1,-1;CURRENT_FAC_ENTRANCE \"CURRENT_FAC_ENTRANCE\" true true false 50 Text 0 0 ,First,#,Table20180918141600,bv_curr_fa,-1,-1;PROPOSED_ENTRANCE \"PROPOSED_ENTRANCE\" true true false 50 Text 0 0 ,First,#,Table20180918141600,bv_prop_en,-1,-1;DROP_TYPE \"DROP_TYPE\" true true false 50 Text 0 0 ,First,#,Table20180918141600,bv_drop_ty,-1,-1;VIRTUAL_TERM \"VIRTUAL_TERM\" true true false 50 Text 0 0 ,First,#,Table20180918141600,bv_virtual,-1,-1;PROPOSED_TERM_HHLAT \"PROPOSED_TERM_HHLAT\" true true false 50 Text 0 0 ,First,#,Table20180918141600,bv_prop_te,-1,-1;PROPOSED_TERM_HHLONG \"PROPOSED_TERM_HHLONG\" true true false 50 Text 0 0 ,First,#,Table20180918141600,bv_prop__1,-1,-1;PROPERTY_ID \"PROPERTY_ID\" true true false 4 Long 0 0 ,First,#,Table20180918141600,bv_propert,-1,-1;FIM_ADDRESS_ID \"FIM_ADDRESS_ID\" true true false 4 Long 0 0 ,First,#;FIM_SUBADDRESS_ID \"DATE_CREATED\" true true false 4 Long 0 0 ,First,#;OPERATION_TYPE \"OPERATION_TYPE\" true true false 50 Text 0 0 ,First,#;ENGG_INTERVAL \"ENGG_INTERVAL\" true true false 50 Text 0 0 ,First,#", "")

Arcpy.Append_management is appending blank rows over into the table, but when I do this in model builder it works correctly. Why would the append function append blank rows instead of actually filling it out with information?


texts:
SO I am looking to do a split on a attribute that looks like 45-4.  Now all the attributes have values.

here is my code so far:

import fme
import fmeobjects

def processFeature(feature):

    f_type = feature.getAttribute('F_TYP')
    if f_type:
        f_type.split('-')
        Height = f_type[0]
        Class = f_type[-1]

        feature.setAttribute("Height", Height)
        feature.setAttribute("Class", Class)

It seems to work but only gives me as an example the height = 4 and the class 4.  I guess I need to take everything in front of the dash when it exists for height and everything after the dash for class.  I do need the if statement to deal with no values.


texts: Canceling a Geoprocessing Tool (Arcobjects)

texts:
I am trying to cancel all the ArcMap geoprocessing tools executions (GPTools) using an ArcObjects Extension in C#, but it doesn't seem to work or do anything for any GPTool that I am trying to cancel.
The code I am using to cancel a GPTool is this:

public void Startup(ref object initializationData)
{
   Geoprocessor GP = new Geoprocessor();
   GP.ToolExecuting += GP_ToolExecuting;
}

private void GP_ToolExecuting(object sender, ToolExecutingEventArgs e)
{
   IGeoprocessorResult2 result = e.GPResult as IGeoProcessorResult2;
   result.Cancel();
}

The application does run the Cancel procedure and the 'result.IsCanceled' property does change to 'true', yet the execution of the GPTool doesn't seem to stop. The ArcMap version is 10.2.2.

Any ideas what am I doing wrong or there is another way to prevent GPTools from running?


texts: Changing color of symbol in UniqueValuesSymbology object using ArcPy?

texts:
I have an template mxd I reference with several pre-defined layers. The only variable I need to modify is the color of certain symbols in the layer. I have tried to modify and dig around the Symbology Object but I cannot access the colors of each Unique Value Symbol.

This is how I am accessing the the Symbology of the individual layers. There doesn't seem to be and option to change the color a value.

I am using ArcGIS 10.2.1

mxd_file = r"some\path\file.mxd"
mxd = arcpy.mapping.MapDocument(mxd_file)

# getting data frame
df = arcpy.mapping.ListDataFrames(mxd)[0]

# getting street layer
street_layer = arcpy.mapping.ListLayers(mxd, 'Street', df)[0]

if street_layer.symbologyType == 'UNIQUE_VALUES':
    # trying to access all possible properties of the symbology object
    street_symbology = street_layer.symbology
    desc = street_symbology.classDescription
    labels = street_symbology.classLabels
    values = street_symbology.classValues
    show_other_values = street_symbology.showOtherValues
    value_field = street_symbology.valueField
del mxd, df

layer_file = r"some\path\file.lyr"
arcpy.SaveToLayerFile_management(path_layer, layer_file)


texts:
I am serving a set of GeoTiffs through GeoServer.
In order to save disk space I have multiplied my original 32 bit floating point data by a factor 1000 and then rounded to the nearest integer, so I could save the data as 16 bit integer GeoTiffs.
This is all fine for a WMS service, as GeoServer applies an appropriate color map.

But later on I want to set up a WCS service as well, allowing the user to download the data itself. Is there a way to let the WCS service rescale the data (so in my case divide by 1000.) when replying to a WCS requests?

I have looked for scaling options online, but only found things related to scaling the spatial resolution rather than the data values.


texts: QGIS and Intel Tablet Processor Compatibility

texts: Will any Windows based tablet with an intel processor be able to run QGIS, provided it meets the requirements, or is it only the Surface Pro that is capable? Looking for a cheaper solution to the pro.


texts: To see the geodata in a PDF is necessary to process something first in QGIS?

texts:
I want to see if a pdf file has geographic data and how to access them.
In the properties of the PDF document I see the pdf file was made in "AutoCAD Map 3D 2013," and have 6 layers.

In QGIS 2.14, when I add a vector layer, appears only 3 layers and the attribute tables only show the ID field. I read somewhere that something must be processed first to see more data but I can´t find it anymore.

It is true that to see more data is required to use some command in QGIS or process something?


texts: Densifying polyline according to specific length using QGIS?

texts:
I am using QGIS 1.8.0
I would like to densify a polyline according to a particular length of segment (ex. 200 m). I know it is possible with other GIS software.

Is it possible to do this in QGIS?


texts:
My  question is about Postgres/PgRouting syntax, when I wrote query like this:

SELECT * FROM shortest_path('SELECT gid  AS id,
      source::int4 AS source,
      target::int4 AS target,
      cost::float8 AS cost,
      reverse_cost::float8 AS reverse_cost
from network where country='USA'  ',
257027,
276521,
true,
true)

I got  error:

ERROR:  syntax error at or near "USA"
LINE 6:  from network where country='USA'  ',
ERROR: syntax error at or near "USA"
SQL state: 42601
Character: 226

So, I create new field in the network table countryInt and I updated  field like this:

Update network set countryInt =1 where country='USA'

So my query is:

SELECT * FROM shortest_path('SELECT gid  AS id,
    source::int4 AS source,
    target::int4 AS target,
    cost::float8 AS cost,
    reverse_cost::float8 AS reverse_cost
from network where countryInt=1 ',
257027,
276521,
true,
true)

And now it works. But what if I have to put some string in condition, quotation marks are problem?

How to resolve this error?

Thank you in advance.


texts: Projecting GPS point to point in closest pgrouting road (map matching)?

texts:
I'm quite new to PostGIS and pgRouting.

I'm interested in doing a simple map matching. Precisely, finding the closest point in a given OSM road network.

I’ve successfully imported the OSM network with osm2pgrouting.
After some googling and lots of read I was able to find how to get the closest road vertex from a given point P (-49.239099 -16.667802) using the following query:

SELECT * FROM ways_vertices_pgr ORDER BY the_geom <->
ST_GeometryFromText('POINT(-49.239099 -16.667802)',4326) LIMIT 1

However, what I'm interested in on the project point of P that should be located on the closest edge in the road network.

How can I do that?


texts: Any way to create polygons with brush-like tool?

texts:
I’m making a hiking map for a small range of mountains. I’d like to show surface coverage, including a few broad categories of vegetation type, rock, sand, etc. I’d like to do this as vector data.

I haven’t managed to get satisfactory results by processing remote imagery, although I’m very inexperienced with that and would welcome any suggestions.

I’ve resigned myself to the fact that I’m going to need to trace it manually from Google Earth or something. But, for covering large areas quickly with regions that are not always obviously contiguous or not, creating polygons by drawing outlines feels cumbersome and not intuitive. I’d much rather be able to use something like a brush tool from a drawing application to “paint” in the areas that I want for each coverage type.

Are there any tools or apps that would help me do this? I’m using QGIS, but am open to incorporating other software into my workflow.

I’d also welcome suggestions for alternative methods of generating vector data for surface coverage.


texts: Having 'can't adapt type 'Point' error when try to load data from local to PostgreSQL database

texts:
I am trying to load this shapefile into the PostgreSQL database.
Here is what data looks like:
first, read the file
data_path = "./RFSNSW_BFPL"
rfsnw_bfpl= gpd.read_file( os.path.join(data_path, "RFSNSW_BFPL.shp") )

then, I am trying to  create schema in the database
conn.execute("DROP TABLE IF EXISTS rfsnw_bfpl")
rfsnw_bfpl_schema = '''CREATE TABLE rfsnw_bfpl (
                     category NUMERIC,
                     shape_leng NUMERIC,
                     shape_area NUMERIC,
                     geom GEOMETRY(POINT,4283))'''

conn.execute(rfsnw_bfpl_schema)

Here, I try to change one column's name.
srid = 4283
rfsnw_bfplWkCpy = rfsnw_bfpl.copy()
rfsnw_bfplWkCpy['geom'] = rfsnw_bfpl['geometry']
rfsnw_bfplWkCpy = rfsnw_bfplWkCpy.drop(columns="geometry")

Lastly, when I try to load my data to database, it shows me 'can't adapt type 'Point''
rfsnw_bfplWkCpy_table_name = "rfsnw_bfplWkCpy"
rfsnw_bfplWkCpy.to_sql(rfsnw_bfplWkCpy_table_name, conn, if_exists='append', index=False,
                         dtype={'geom': Geometry(geometry_type='POINT', srid= 4283)})

Here is the screenshot of the error


texts: Drape multiple lines on a DEM with Postgis

texts:
I am trying to drape lines on a DEM by following this great tutorial:
http://blog.mathieu-leplatre.info/drape-lines-on-a-dem-with-postgis.html
When I do it on a single line by giving points coordinates it works perfectly:

WITH line AS
    -- From an arbitrary line
    (SELECT 'SRID=32632;LINESTRING (348595 4889225,352577 4887465,354784 4883841)'::geometry AS geom),
  linemesure AS
    -- Add a mesure dimension to extract steps
    (SELECT ST_AddMeasure(line.geom, 0, ST_Length(line.geom)) as linem,
            generate_series(0, ST_Length(line.geom)::int, 50) as i
     FROM line),
  points2d AS
    (SELECT ST_GeometryN(ST_LocateAlong(linem, i), 1) AS geom FROM linemesure),
  cells AS
    -- Get DEM elevation for each
    (SELECT p.geom AS geom, ST_Value(mnt.rast, 1, p.geom) AS val
     FROM mnt, points2d p
     WHERE ST_Intersects(mnt.rast, p.geom)),
    -- Instantiate 3D points
  points3d AS
    (SELECT ST_SetSRID(ST_MakePoint(ST_X(geom), ST_Y(geom), val), 32632) AS geom FROM cells)
-- Build 3D line from 3D points
SELECT ST_MakeLine(geom) FROM points3d;

However I am at a loss when it comes to applying the query to a whole table. The author mentioned that 'Your function would accept a parameter line, that you would use in the code! Just remove the WITH line AS part'
So the function would look like

CREATE OR REPLACE FUNCTION drape(line geometry) RETURNS geometry AS $$
.......

WITH linemesure AS (SELECT ST_AddMeasure(line, ...

But since there is no more temporary table (like line was before) I do not know what to select linemesure FROM.


texts: Openlayers squeezes placemarks after changing map projection

texts:
I have a map with OSM as a base layers and several layers on top of that. Two layers (red polygons with gray fill and black placemarks) are instances of TileLayer (source TileWMS -> geoserver). Both layers are stored in postgis database and published on geoserver.

If I set map projection 'EPSG:4326' than placemarks are displayed correctly (as it looks when looking at png representing placemark). Image below:

EPSG:4326 doesn't fit my needs as it distorts northern parts of map. That's why I tried to change to EPSG:3857.

And I understand that polygons are squeezed because of projection change. However, is there possibility to prevent placemarks from being squeezed/streched after change of projection.


texts:
SELECT sum(st_length(st_intersection(POLY.geom, LINE.geom))) AS sum_line_length...
Is this function OK to calculate the sum of lines length? Wonder if it do the clip line where one line go from one polgon to another or it just sum the line length to polygon from witch the line started? Thanks for answer. Best regards.


texts:
I am working with Landsat 8 image acquired in 2014. I am trying to convert the raw image (DN) to Top Of Atmosphere Reflectance/Radiance (i.landsat.toar) in Grass 7.1 (still under development). The reason for using this version is that it is caters for landsat 8. I specify the metadata file and load the information etc. I get the same error every time:

  ERROR: Unable to open header file for raster map

Base Name of input raster band: LC81700832014160LGN00_B1@PERMANENT
Prefix for output raster maps: LC81700832014160LGN00_B1.toar
Name of metadata file: C:\Users\carinswart\Documents\GIS_DataBase\Landsat8\GHT\09June14\LC81700832014160LGN00_MTL.txt

I have tried to be careful not to leave gaps in my file names.


texts: Tagging polygons within larger polygons (spatial join)

texts:
another QGIS newbie here. I have a number of smaller polygons which perfectly aggregate to a larger polygon. I want to "tag" the smaller polygons with the ID of the larger polygon they make up.  Unfortunately when I used spatial join many of the smaller polygons were incorrectly joined to the large polygon they were adjacent to - rather than the large polygon they were contained by (all polygons are contiguous). The larger polygons were created using "dissolve" from the smaller polygons, so there shouldn't be a problem in compatibility of the geometries.

In MapInfo you would use a spatial join term like "contained by" or "within" for this. What is the equivalent in QGIS to restrict the spatial join to large polygons the smaller polygon is entirely, or at least primarily, within?

PS I have managed to partially work around this by creating centroids and then tagging them. However, that is somewhat unsatisfactory, and doesn't work where the polygons are strange shapes, meaning the centroid is actually outside the polygon.


texts: Why is 'Line Dissolve' producing different output than 'Dissolve' from Ftools?

texts:
I'm curious if anyone can explain the strange effects I'm seeing from the SAGA 'Line Dissolve' tool. I have a ton of line features that had previously been 'exploded' (which is a necessary step), and now need to be rejoined together so that 'Create Points Along Lines' works properly.

Using the Vector > Geoprocessing > Dissolve tool takes an ungodly amount of time to process (I'm talking hours for a single dataset, which is much too long), but creates expected results after adding the points along the dissolved lines. When I run the 'Line Dissolve' tool from the processing toolbox, it takes about 5 seconds to process, but the points along the lines are sporadic and seemingly placed at random (see attached screenshots).

Points along lines after Vector > Geoprocessing > Dissolve:

Points along lines after "Line Dissolve":


texts: Drawing a line with the same length as the layer extent in QGIS Modeler

texts: How can I use the height of a polygon  (layer extent) to construct a line slightly longer than the polygon? I’m in the Modeler and my input is s polygon layer. I will be constructing the line from a point With an azimuth derived from another line. I guess I will be using the geometry from expression tool, but I don’t know the expression.


texts:
I use gwt-openlayers lib in my GWT app. I have added a drawFeature control and now I want to show a radius of circle while I draw a feature. Now it looks like this:

vectorLayer.addVectorFeatureAddedListener(new VectorFeatureAddedListener() {
        @Override
        public void onFeatureAdded(FeatureAddedEvent eventObject) {
            final VectorFeature feature = eventObject.getVectorFeature();
            if (feature.getGeometry().getClassName().equals("OpenLayers.Geometry.Polygon") && !hasTemplate) {
                if (currentAreaFeature != null) {
                    vectorLayer.removeFeature(currentAreaFeature);
                }
                if (templateAreaFeature != null) {
                    removeTemplateArea();
                }
                currentAreaFeature = feature;
                Geometry geometry = currentAreaFeature.getGeometry().clone();
                geometry.transform(new Projection(map.getProjection()), DEFAULT_PROJECTION);
                geometry.calculateBounds();
                Bounds bounds = geometry.getBounds();
                currentCenter = bounds.getCenterLonLat();
                currentRadius = calculateDistance(currentCenter, new LonLat(bounds.getUpperRightX(), currentCenter.lat()));

                style.setLabel(currentRadius.intValue() + "m");
                feature.setStyle(style);
                feature.redrawParent();
            }
            hasTemplate = false;
        }
    });

In this case, a label with radius represents after the feature added. I want to represent a radius always if(drawFeature.isActive). The first idea to realize this - using the MousePosition control. Thus I can take a second LonLat parameter to calculate a radius using my custom method calulateRadius().  I hope that I can take the first Lonlat (= future center of circle) using mapMouseClickEvent. But if drawFeature.isActive it is impossible to handle this event. So, I am at a deadlock.
This is my solution that doesn't work:

   map.addMapClickListener(new MapClickListener() {
                    @Override
                    public void onClick(final MapClickEvent mapClickEvent) {
                        MousePositionOutput mpOut = new MousePositionOutput() {
                            @Override
                            public String format(LonLat lonLat, Map map) {
                                String out = calculateDistance(mapClickEvent.getLonLat(),lonLat)+"m";

                                return out;
                            }
                        };

                        MousePositionOptions mpOptions = new MousePositionOptions();
                        mpOptions.setFormatOutput(mpOut); // rename to setFormatOutput

                        map.addControl(new MousePosition(mpOptions));
                    }
                });


texts: How many interior and exterior rings can a polygon have in a standard ArcGIS shapefile?

texts:
working with GDAL/OGR library and studying these two slides about geometries and vector data in the OGR library:

slide 1
slide 2

I found that a polygon can have an exterior ring and lots of interior rings for an arbitrary format that OGR supports.
But I'm developing a program that will support only ESRI's point, line and polygon shapefiles.
For testing my program, I decided to draw a polygon shapefile in the Arcmap with the Editor that will have ten polygons and each polygon will have different numbers of interior rings.
But as much as I play with the editor trying to create a polygon that has one interior ring ( for example I created two polygons - one inside another - and then merged them together ), I couldn't create such polygon.

So could you tell me :

How many exterior and interior rings can a polygon have in an ESRI's
polygon shapefile ?
How can I create a polygon that has 1 exterior and 2 interior rings
with the editor in Arcmap?


texts: GeoJSON layer shows/enables only one row in attribute table in QGIS

texts:
I have a GeoJSON polygon layer which I converted from a shapefile layer. The layer opens up in geojson.io, for example, with no problem, displays the table correctly, etc:

However, in QGIS the Attribute Table only shows one polygon, with "Features Total:146, Filtered:1, Selected:0", see the image below

I've tried to remove this "filter", using SQL to select all the polygons or just a couple of the polygons, then export selected, and still, I only get this one unit, 146.

How can I get QGIS to display the whole attribute table in GeoJSON and other formats?

One caveat: the original shapefile had a "lock" file on it. I don't have this problem when I export to .shp, though if I try to select polygons and export a subset of the original shapefile in any format, I run into the same problem.


texts: How to do loops on raster cells with python console in QGIS?

texts:
I'm new here and beginner with python. I'm looking for a solution to my problem for some hours but now you are my last chance...
I'm using the python (2.7.5) console in QGIS (2.4.0). Is there a solution to do loops on  raster cells? I would like to create a raster in which each value cell will depend on the position (column and row number) in the raster.

I'm looking for something like this python code in Arcgis (I found it here: How to get X Y coordinates and cell value of each pixel in a raster using Python?)

for row in range(0,height):
     for col in range(0,width):
        #now you have access to single array of values for one cell location

Thanks!


texts: Using OSTN15 to convert British National Grid data to KML

texts:
I'm trying to save 2D data on a layer (currently to British National Grid co-ordinates) to an accurate KML version using QGIS 2.14.19 Essen. Saving directly to KML using the WGS84 CRS doesn't give an accurate result.

I have noted that accurate transformation needs to use the OSTN15_NTv2.gsb files provided by Ordnance Survey, together with an updated SRS database file - both of which I have saved down into the appropriate folders (as recommended here, for example, https://dragons8mycat.com/2016/09/26/add-ostn15-to-qgis-2-16/) but when I try to save the transformed data in QGIS (right click on the layer > Save as...) the OSTN15 CRS option isn't listed.

I have restarted QGIS etc. but I can't fathom how to try the file/coordinate system conversion using this transformation CRS.

Can anyone help, please?


texts:
I have an attribute table with over 2500 fields. About 2000 fields start with the name "All_txt_" (eg "All_txt_Watermelon").

Is there a way to remove this "All_txt_" from all the fieldnames automatically? (so the name would be "Watermelon").
Or does anyone know how to do so in Python? I work with the 10.1 version and the Alter Field tool is not an option.


texts: OpenLayers Modify interaction on line - use different styles for "Move Vertex" and "Add Vertex"

texts:
I use OpenLayers 6.5 for my mapping application.
With the Modify Interaction, it is possible to add or move vertices for a line layer. However, the symbol (circle) is the same whether the mouse is over an existing vertex (to be moved) or a non-existing (potential) vertex laying on a line (to be created), see the pictures below.
The is only a single style option for the Modify interaction. Any hint on how to setup a different style for these two types of nodes?
Existing vertex (can be moved)

Non-existing vertex (can be created)

Here is the code for my current Modify interaction, where vectorLayerSource is a LineString vector source:
let modify = new ol.interaction.Modify({
    source: vectorLayerSource,
    deleteCondition: ol.events.condition.shiftKeyOnly,
    insertVertexCondition: ol.events.condition.shiftKeyOnly,
});


texts: PyProj: Transform from LLA to ORTHO with WGS84 gives sphere instead of ellipsoid

texts:
I want to convert my coodrinates from lon, lat to x, y of orthographic projection. I tried to use PyProj to do it, but it uses sphere, not ellipsoid. Here it is:

lonlat_proj = pyproj.Proj(proj="longlat", ellps="WGS84")
lat = 90
lon = 0
alt = 0

# ECEF IS OKAY
ecef_proj = pyproj.Proj(proj="geocent", ellps="WGS84")
ortho_proj = pyproj.Proj(proj="ortho", lat_0=0, lon_0=0, ellps="WGS84")

# ORTHO IS NOT!!!
x, y, z = pyproj.transform(lonlat_proj, ecef_proj, lon, lat, alt)
print(x, y, z)
# Gives: 3.918e-10 (okay), 0.0, 6356752.314 (okay, it's earth semiminor)

x_ortho, y_ortho = pyproj.transform(lonlat_proj, ortho_proj, lon, lat)
print(x_ortho, y_ortho)
# Gives: 0.0 (okay), 6378137,0 (earth SEMIMAJOR instead of semiminor)
# WHY?!!

How can I cope with it?


texts:
The query is:

Select COUNT(*) from highway_ref where (ST_DWithin(ST_Transform(highway_ref.geom,2163),ST_Transform(ST_GeomFromText('POINT(73.919595 18.5622609)', 4326),2163),500))

I am trying to find the total number of roads which are within 500 metres from a particular lat long.

It's taking too long. How can I optimize it?