texts:
I am trying to navigate in an outdoor area from one point to the other, e.g. in an area with hills or mountains, and I would like to declare a hill as an obstacle if the change in elevation is too large relative to the distance. What tools would be best to use to work through this problem, and from what sources I can get DEM data for this purpose. (e.g. Scottish Highlands)

I would like to consider an area of 100 metre radius at a time for the purpose of obstacle identification and avoidance.

Can the hills (Obstacles) be avoided from the gradient descent threshold?

I am new to GIS and don't know much about it.


texts:
I want to create a snap control of my layer. My simple example demo is here

var raster = new ol.layer.Tile({
    source: new ol.source.OSM({
    })
});

var source = new ol.source.Vector({
        format: new ol.format.GeoJSON(),
        url: 'https://raw.githubusercontent.com/openlayers/ol3/6838fdd4c94fe80f1a3c98ca92f84cf1454e232a/examples/data/geojson/countries.geojson'
    })

// vector layer
var vector = new ol.layer.Vector({
    source: source,
    style: new ol.style.Style({
        stroke: new ol.style.Stroke({
            color: '#eedd00',
            width: 1.5
        })
    })
});

// render the map
var map = new ol.Map({
    target: 'map',
    renderer: 'canvas',
    layers: [raster, vector],
    view: new ol.View({
        center: [0, 0],
        zoom: 1
    })
});

var snap = new ol.interaction.Snap({
  source: source
});

map.addInteraction(snap);

snap.setActive(true);

I activated the snap interaction but it does not snap vectors


texts: determining GPS points from a grid in meters

texts: I am new to GIS use and I have a plan (scanned from paper, old fashion) with a grid of 5x5 meters cells and a GPS point in the middle. I wish to georeference my plan, but I only have one known point. I was wondering if there was a calculation to find out what coordinates are associated to the points on the grid?


texts: Converting Labels to Annotations moves position of labels

texts:
Why are my contour line labels moving after I convert them to annotations?

I am using Maplex to place them. But if they still move after I converte them to annotations, the whole maplex placing process doesn't make sense.


texts: Multi-temporal Maximum Likelihood classification?

texts: Here is my basic questions. Does it make sense from a theoretical point of view to use the Maximum Likelihood classifier in a multi-temporal dataset of satellite images (Sentinel-2)? I mean, perform a single MLC classification for the complete multitemporal dataset, not MLC for each image.


texts: Interpolate Line Tool giving different results?

texts:
I have been using 'Interpolate Line' Tool in the 3D Analyst tool pack to create a profile graph so I can analyze interpolated results gathered from an IDW created surface. However, when I use the tool in the same spot (has the same X Y starting point and the same width each time) at different times, I will obtain different profile graphs each time.

For example, the first time I used this tool, I created a line with a starting point of xxx,xxx m E, x,xxx,xxx m N with a length of 34,500 m. I obtained this profile graph from the first trial.

The second time I used this tool, with the same starting point and width as the first line, I obtained this profile graph.

Can anyone explain why this is, or how I can fix it?


texts: "Hole" argument ignored in call to Polygon function

texts:
I am working with a dataset consisting of a list of coordinates for a set of polygons. Some of these polygons are holes. However, the dataset does not follow the convention of listing the coordinates of holes in counter-clockwise direction. Therefore, I have to check for each polygon whether they are a hole, and use the result as argument to the Polygon function. Here is a simple dummy dataset:

library(sp)

# Coordinates for two dummy polygons (squares), one of which is fully contained within the other

x1 <- c(5,9,9,5,5)
y1 <- c(5,5,1,1,5)
x2 <- c(6,8,8,6,6)
y2 <- c(4,4,2,2,4)

tables <- list()    # Will contain the coordinates of each object, each in a data.frame
tables[[1]] <- data.frame(x1,y1)
tables[[2]] <- data.frame(x2,y2)

nobjs <- 2  # Number of objects

Next, I implemented an ugly but functional way of determining whether a polygon is a hole (i.e. fully contained within another polygon), by checking whether all of the nodes of a certain polygon are contained in any other polygon of the dataset:

holes <- logical(nobjs)
# Check for polygons that are holes
for(i in 1:nobjs){
    vectTF <- vector("logical",length(tables[[i]]$V1))
    for(j in 1:length(tables[[i]][,1])){
        vectTFj <- vector("logical",length(tables))
        for(k in 1:length(tables)){
            if(k==i){
                vectTFj[k] <- FALSE
            } else {
            vectTFj[k] <- point.in.polygon(tables[[i]][j,1],tables[[i]][j,2],tables[[k]][,1],tables[[k]][,2])
            }
        }
        vectTF[j] <- as.logical(sum(vectTFj))
    }
    holes[i] <- as.logical(sum(vectTF))
}

In this example, the vector holes correctly tells me that my second object is contained within another polygon, while the first one is not:

holes

  [1] FALSE  TRUE

Next, I try to build a Polygons object using the vector holes as the hole argument to the Polygon function:

objs <- list()
id <- integer(nobjs)
for(i in 1:nobjs){
    objs[[i]] <- Polygons(list(Polygon(tables[[i]][,1:2],hole=holes[i])),as.character(i))
    id[i] <- i
}

Now, when I check the structure of the new Polygons object, I see that none of my polygons are marked as being holes:

  str(objs)
      List of 2
      $ :Formal class 'Polygons' [package "sp"] with 5 slots
    .. ..@ Polygons :List of 1
    .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
    .. .. .. .. ..@ labpt  : num [1:2] 7 3
    .. .. .. .. ..@ area   : num 16
    .. .. .. .. ..@ hole   : logi FALSE
    .. .. .. .. ..@ ringDir: int 1
    .. .. .. .. ..@ coords : num [1:5, 1:2] 5 9 9 5 5 5 5 1 1 5
    .. .. .. .. .. ..- attr(, "dimnames")=List of 2
    .. .. .. .. .. .. ..$ : NULL
    .. .. .. .. .. .. ..$ : chr [1:2] "x1" "y1"
    .. ..@ plotOrder: int 1
    .. ..@ labpt    : num [1:2] 7 3
    .. ..@ ID       : chr "1"
    .. ..@ area     : num 16
   $ :Formal class 'Polygons' [package "sp"] with 5 slots
    .. ..@ Polygons :List of 1
    .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
    .. .. .. .. ..@ labpt  : num [1:2] 7 3
    .. .. .. .. ..@ area   : num 4
    .. .. .. .. ..@ hole   : logi FALSE
    .. .. .. .. ..@ ringDir: int 1
    .. .. .. .. ..@ coords : num [1:5, 1:2] 6 8 8 6 6 4 4 2 2 4
    .. .. .. .. .. ..- attr(, "dimnames")=List of 2
    .. .. .. .. .. .. ..$ : NULL
    .. .. .. .. .. .. ..$ : chr [1:2] "x2" "y2"
    .. ..@ plotOrder: int 1
    .. ..@ labpt    : num [1:2] 7 3
    .. ..@ ID       : chr "2"
    .. ..@ area     : num 4

Can anyone tell me what I'm doing wrong?


texts: Iterating through table for IDW says parameters are not valid?

texts:
I have a table with precipitation data from 5 weather stations and I am trying to iterate through the table to perform IDW on each field (month of precipitation).

import arcpy
from arcpy import env
arcpy.CheckOutExtension("GeoStats")
env.workspace = "C:\Users\Lab R220\Desktop\Riley"
arcpy.env.extent = "Thames_watershed.shp"
fieldList = arcpy.ListFields("stationswithprecip.shp")

for field in fieldList:
    try:
        outraster = str(field.name) + ".tiff"
        arcpy.IDW_ga("stationswithprecip.shp", field.name, outraster, 1, 2)
    except:
        print arcpy.GetMessages()
        pass

The error I am receiving are as follows:

Executing: IDW "C:\Users\Lab R220\Desktop\Riley\stationswithprecip.shp" id id.tiff "C:\Users\Lab R220\Desktop\Riley\1" 2 2 "NBRTYPE=Standard S_MAJOR=25461.6476799568 S_MINOR=25461.6476799568 ANGLE=0 NBR_MAX=15 NBR_MIN=10 SECTOR_TYPE=ONE_SECTOR" #
Start Time: Mon Jul 22 09:36:46 2019
Failed to execute. Parameters are not valid.
ERROR 000361: The name starts with an invalid character

However if I run the months through the IDW tool individually they work. I have 120 of these so I would rather not run them one by one.

My field names look like this with the numbers increasing by 1  "Jan2008" and they are all double


texts: Colouring attribute/field on QGIS data form based on value?

texts:
Using QGIS (2.14) is it possible on the data form for an attribute to have a colour (or other formatting) based on its value?

For example if my data has a field 'quality' which has a value 'good' or 'medium' or 'bad'. I'd like the field in the attribute form to show up green/orange/red depending on this value.

I'd be happy with a workaround - for example a symbol or icon that appears beside the field, coloured background to the field, and so on.

I'm happy to use QT Designer to create the data form. I'm happy to edit the form xml by hand. I'll struggle if the answer is that I need to use more than a line or two of Python code so I'm probably looking for answers which don't require this.

I also need the formatting to show up whether or not the layer is in edit mode (because the main data I'm using won't be editable).

EDIT: Note that I'm only looking for solutions that work for a QGIS data form, not on the attribute table, and not on exported data. The aim is to make data as easy to understand as possible, so a user working with the identify tool can quickly see/understand key information about a feature shown on the map. The real data  I'm working with is complex and front-line users need to be assisted to understand it.


texts:
I have two versions, from same base table, with only one feature class in both that I want to combine in a database view for a presentation. Version feature classes differ both with attributes and geometry but are children of a same parent. I wrote a SQL query that worked:
CREATE OR ALTER sde.CominedViews AS
SELECT *, ROUND(Shape.STLength(), 1) AS [CalcLength]
FROM [sde].[ATableForABaseTable]
WHERE SDE_STATE_ID IN
(SELECT state_id
FROM SDE.SDE_versions
WHERE [name] = 'VERSIONNAME1')
UNION ALL
SELECT *, ROUND(Shape.STLength(), 1) AS [CalcLength]
FROM [sde].[ATableForABaseTable]
WHERE SDE_STATE_ID IN
(SELECT state_id
FROM SDE.SDE_versions
WHERE [name] = 'VERSIONNAME2')

However, on some occasions, when the user does updates, one of the versions can not found state_id in a A delta table. I have checked A and D delta tables and SDE_version. In this example, SDE_versions table have state_id 15360 assigned and update that user did is in A delta with state_id 15356 assigned (One that it should be shown). D table does increases this number but highest value in SDE_STATE_ID is 1537, still not the one I am looking for from SDE_versions.
Where this state_id in SDE_versions could came from?
I thought it might be not saved updates, but shouldn't then that number be higher than SDE_versions number for a particular version?
And if I open that version in ArcGIS, it is shown as it should, updated.
Compress is not run in the meantime. All changes are done by the SDE owner.
Data is registered as traditionally versioned without archiving.
Do you have any tips?


texts: arcpy.ListDatasets returns an empty list

texts:
I can't seem to get arcpy.ListDatasets() to work correctly. My code and output is below. What am I doing wrong?

import arcpy
import os

work_path = os.getcwd() + "\\Lesson4_Data"
arcpy.env.workspace = work_path
print(arcpy.ListFeatureClasses())
print(arcpy.ListDatasets())

Here is the output

[u'Adams.shp', u'Alamosa.shp', u'Arapahoe.shp', u'Archuleta.shp', u'Baca.shp', u'Bent.shp', u'Boulder.shp', u'Broomfield.shp',
u'Cheyenne.shp', u'ConeJos.shp', u'Costilla.shp', u'Crowley.shp', u'Custer.shp', u'Delta.shp', u'Denver.shp', u'Dolores.shp', u'Douglas.shp', u'Eagle.shp', u'Elbert.shp', u'ElPaso.shp', u'Garfield.shp', u'Gilpin.shp', u'Grand.shp', u'Gunnison.shp', u'Moffat.shp']
[]


texts: Using decimals in ArcGIS Raster Calculator tool?

texts:
I would like to create a new raster by using the following regression equation in the Raster Calculator tool in ArcGIS 10.3:

Outputras=0,000826345+(0,000203494*"Roads_density.tif")

("Roads_density.tif" being an existing raster among my set of layers)

However I get the following message error:

ERROR 000989 Python syntax error: Parsing error SyntaxError: invalid syntax (line 1).

When I try :

Outputras=0,8+0,2*"Roads_density.tif"

the message error disappears.

Is it thus a problem related to the presence of too many decimals ? What do I need to do to force ArcGIS to take in account all the decimals?


texts: How to check if a polygon contains another polygon in OL3?

texts:
Is there any possibility to check if a polygon is contained in another polygon in OL3?

I noticed there is a couple of methods in ol.extent e.g. ol.extent.containsExtent(extent1, extent2). But this checks only the extents and not the entire polygon. So in some cases I get "true" although the polygon is way outside the other polygon, but still in the extent...
What I need is a function or algorithm that checks if a polygon entirely contains another polygon like I usually do in any GIS.


texts: Select element switching to pan using ArcMap

texts:
Does anyone know if there's a way to keep the "select element" tool from automatically switching to the "pan" tool when changing from layout view to data view in ArcMap 10.7?
I feel like there should be a config setting somewhere, I just can't find it.  I keep moving my map b/c it's automatically switching to pan and then I have to manually move it back (not always exactly how it was) because undo doesn't work for that.


texts: Why is SQLite python command not working?

texts:
I am fairly new to python and am writing a code to read imagery and insert them into an SQLite database but I get the following error

cur.execute(statement, (index,fullpath,filename))
sqlite3.OperationalError: near "index": syntax error.

What is it that I am getting wrong? The following is the code:

for root, dirnames, filenames in os.walk(basedir):
    print 'root is ' + root
    for filename in filenames:
        ext = os.path.splitext(filename)[1].lower() #extension as lowercase
            if ext == extension:
                index = filename[:9] #tile ID
            filebase = os.path.splitext(filename)[0] #first part of filename
                fullpath = os.path.join(root,filebase) + ext #build up path again
                statement = "INSERT INTO %s (%s) VALUES " % (table,columns) + "(%s,%s,%s);"
                cur.execute(statement, (index,fullpath,filename))
            print 'catalog entry: ' + index + ', ' + fullpath


texts:
I want to import cartopy on Python3, but it keeps showing me
File "shapely\speedups\_speedups.pyx",
  line 408, in shapely.speedups._speedups.geos_linearring_from_py
  GEOSGeom_createLinearRing_r returned a NULL pointer

I've tried reinstalling cartopy, but it won't work at all. Is there something wrong with GDAL or Fiona?


texts: How to correct aspect values in a North Pole Stereographic projection?

texts:
I'm working on a DEM in the Arctic region (initially in a NAD83 coordinate system) and I projected it into a North Pole Stereographic Cartesian coordinate system to compute slope and aspect. Now, because of the distortions, slope and aspect have to be corrected.

I corrected the slope without much trouble (Using SRTM Global DEM for Slope calculation?) but I'm stuck with the aspect correction.

I carefully read the suggestion made in How to use SRTM Global DEM for Aspect calculation? to add the longitude of the location + 90 degrees to the computed aspect value but it doesn't work as expected in my case (some values are negative as far as -87 degrees). I am not sure that the person was speaking about the same polar projection I'm using (he spoke about a Northern Stereographic projection while I use a North Pole Stereographic with Greenwich as prime meridian: 0.0).

Here are the characteristics of the projection I use:

Prime Meridian Greenwich   0.0
False Northing             0.0
Central Meridian           0.0
False Easting              0.0
Linear Unit Meter          1.0
Latitude of Origin        90.0

After much trial and error, I am looking now for any general formula for this aspect correction, taking into account the characteristics of the projection I use (or any other projection if it exists).

EDIT:

I have tried the method mentioned in my comment and I still have doubts about the results. I still have an offset of more or less 70 degrees between the aspects produced using Mercator (which are correct) and the aspects produced with a North Pole projection which are corrected by subtracting the longitude and then adding or subtracting 360 if the obtained result is < 0 or > 360 respectively. Any further ideas welcome...


texts:
I have a gps device that does not output very accurate gps positions.

How can I round those coordonates so that they will match the street course and not buildings etc ?

I want to draw a route on google maps but the line strikes buildings etc. I don't want that i just want a route on the street.

Thanks


texts: 3rd party JSON into ArcGIS Operations Dashboard

texts:
I have a 3rd party JSON response of coordinates and attributes coming from SQL Server via WebMethods and want to get it into an ArcGIS Operations Dashboard.
I can't seem to find a straightforward way to do this.

I can break it apart and get the coordinates to put it into an ArcGIS JavaScript web map, but was hoping to use Operations Dashboard for some of the UI and charting functionality.

This is a simplified version of what the response looks like:

{"records":[{"ID":1,"X_COORD":"50","Y_COORD":"50"},{"ID":2,"X_COORD":"60","Y_COORD":"60"}]}


texts: Geotools: unexpected difference between CRS.decode("EPSG:4326") and DefaultGeographicCRS.WGS84

texts:
I wanted to use geotools (15.3) to reproject data from EPSG:4326 to a local projection. To see if everything works as expected I printed the length of two test linestrings. The expected result is around 200 meters.

However I noticed that the result varies greatly depending on how I define the source CoordinateReferenceSystem. I thought these two CRS should be identical, however they are not:

DefaultGeographicCRS.WGS84 -> works as expected
CRS.decode("EPSG:4326") -> completely wrong results

For DefaultGeographicCRS.WGS84 the defined CRS and the (correct) lengths are:

GEOGCS["WGS84(DD)",
    DATUM["WGS84",
        SPHEROID["WGS84", 6378137.0, 298.257223563]],
    PRIMEM["Greenwich", 0.0],
    UNIT["degree", 0.017453292519943295],
    AXIS["Geodetic longitude", EAST],
    AXIS["Geodetic latitude", NORTH]]
northSouth: 199.70672773686056m
eastWest: 200.33865746495295m

For CRS.decode("EPSG:4326") the defined CRS and the (completely wrong) lengths are:

GEOGCS["WGS 84",
    DATUM["World Geodetic System 1984",
        SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
    UNIT["degree", 0.017453292519943295],
    AXIS["Geodetic latitude", NORTH],
    AXIS["Geodetic longitude", EAST],
    AUTHORITY["EPSG","4326"]]
northSouth: 222.76855610197927m
eastWest: 346.2445236833851m

How is this possible or what am I doing wrong? The result is the same with geotools 14.0.

import org.geotools.factory.Hints;
import org.geotools.geometry.jts.JTS;
import org.geotools.geometry.jts.JTSFactoryFinder;
import org.geotools.geometry.jts.WKTReader2;
import org.geotools.referencing.CRS;
import org.geotools.referencing.crs.DefaultGeographicCRS;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.MathTransform;

import com.vividsolutions.jts.geom.Geometry;
import com.vividsolutions.jts.geom.GeometryFactory;

public class TransformFailMain {

    public static void main(String[] args) throws Exception {
        CoordinateReferenceSystem sourceCRS = CRS.decode("EPSG:4326");
        // CoordinateReferenceSystem sourceCRS = DefaultGeographicCRS.WGS84;
        CoordinateReferenceSystem targetCRS = CRS.decode("EPSG:31256");

        Hints hints = new Hints(Hints.CRS, sourceCRS);
        GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory(hints);
        WKTReader2 reader = new WKTReader2(geometryFactory);

        // two lineStrings in Vienna, both around 200m long
        Geometry linestringNorthSouth = reader.read("LINESTRING(16.368831 48.201796,16.368831 48.2)");
        Geometry linestringEastWest = reader.read("LINESTRING(16.368831 48.2,16.371526 48.2)");

        MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS);
        System.out.println(sourceCRS);
        System.out.println("northSouth: " + JTS.transform(linestringNorthSouth, transform).getLength() + "m");
        System.out.println("eastWest: " + JTS.transform(linestringEastWest, transform).getLength() + "m");
    }

}


texts: How to get the server address of the ArcGIS lyr. file

texts: I have a lyr. file and want use the content in QGIS. So my question is how can I check where this lyr. linked to. I don't need the style just want data itself.


texts: Problem adding feature name to field in attribute table

texts:
I have a small script that should take the name of the file and add it to a newly created field in several shapefiles.

import arcpy

arcpy.env.workspace = r"c:\users\athom\desktop\timber\New4"

arcpy.env.overwriteOutput = True

fcs = arcpy.ListFeatureClasses()

for fc in fcs:
    field = 'NAME'
    expression = str(fc)
    print expression
    arcpy.AddField_management(fc, field, "TEXT")
    arcpy.CalculateField_management(fc, field, "expression", "PYTHON")

del expression
del fc
del fcs

However, the file name of the first feature is written into the attribute table for each following features even though it shows in the python interpreter as the correct file name as "expression" - the attribute for NAME ends up being SA14107.shp in each shapefile.

Why is the expression not being overwritten?


texts:
I'm trying to use arcpy SearchCursor to iterate through an attribute table, and the only fields that I'm interested are line_code, w_min and ang_ataque.
This is not running properly because I got the type error :

  "'Row' object is not iterable"

What's wrong in my code?

import arcpy
featureClass="C:\\Users\\andreia.medeiros\\Documents\\ArcGIS\\Trabalho\\Modelo\\linhas_150_5km.shp"

data=arcpy.SearchCursor(featureClass,"","","Line_code;w_min;ang_ataque","Line_code;w_min;ang_ataque")

import collections
d=collections.defaultdict(list)
for Line_code,w_min,ang_ataque in data:
     d[Line_code].append(w_min*ang_ataque)


texts:
I'm new to GIS. It is very intersting.

I have installed Postgis, insert shapefile, get data in SVG without problem (ST_AsSVG).

Now, I'd like to export the same svg but in pixel to add it on top of an image. Is it possible directly in the SQL query?

(I can parse the SVG in JS or PHP and convert every point...)

Thanks in advance for your help.


texts:
I am trying to set up HTML text in a legend of this kind:
 O-1  - Kind of zone
, but the button render as HTML is only in the text properties.
How can I use HTML in a layout legend?


texts: Retrieve QgsPalLayerSettings object from a layer

texts:
A pyqgis stand-alone app labels layer features with QgsPalLayerSettings, in the same fashion as in Labelling line layer using QgsPalLayerSettings() QGIS failed, it works for me too.

How do I retrieve that pal layer object after it is "tied in" with a layer for further manipulation (enable/disable etc) or do I just define another pal layer object and writeTo that same layer again to overwrite it?

there is a fromLayer(layer) static function of the QgsPalLayerSettings class that seems to be able to get back a PalLayer object, but it looks like it is a cloned copy, I would like to just get the original pal layer object.


texts:
I'd like to use mapshaper to simplify shapefiles http://www.mapshaper.org/
The main reason is to reduce my file sizes in QGIS as the native tools are too rough / split countries apart.

I can import SHP or topojson OK to QGIS, however all the data tables from my shapefile are gone. How can I get mapshaper to retain country names and codes? I'm new to code.

Im using vectors from Natural Earth http://www.naturalearthdata.com/downloads/10m-cultural-vectors/


texts: I have a postal code file, I want to geocode to that polygon (not point) so that I can perform some analysis on each of the polygons.  Sometimes more than one point will share the same polygon so it can't be a 1 to 1 join.


texts:
I am looking for an overview map functionality for Leaflet. A participant at a conference was quite convinced that this is part of the Leaflet API - unfortunately I don't remember his name so I cannot ask him where to look for that.

The only overview map functionality for Leaflet I was able to find is
this code on github (also mentioned here on gis.stackexchange).

Assuming that an overview map is a rather commonly used feature and taking the low github popularity of 4 stars into account, I am wondering if there is a more popular solution.

Is an overview map funtionality included in the Leaflet API? Or is there another solution besides the aforementioned?


texts: Selecting features based on another feature attributes?

texts:
Here is my case: I have a layer of features (buildings from a city) that is populated with polygons. This polygons have attributes : street name and postal code

In the project I have another layer of features, that is populated with points. Each point has similar attributes : street name and postal code. For each polygon we must add a point that has the same values.

ex: we have a polygon whit attributes >
   - Street_name> Main Street
   - postal code> 96

after this polygon is drawn and it's attributes are updated I am required to create a point feature (the other layer) that has the same attributes
   - Street_name> Main Street
   - postal code> 96

What i need to do > sometimes my team forgets to add the point feature, and the project is quite big. Can I write some code in python or use some tools to detect if every polygon has an equivalent point feature?

ex: i created and updated attributes for the polygon
   - Street_name> Main Street
   - postal code> 96
but i forgot to add the point feature with the same attributes, or made an typo (main street becomes Main Stret)

the tool should scan all the database , all the polygons and all the point features, compare they'r values and finally display (or store in a new layer) the polygons that don't have an equivalent point feature.


texts: Converting Feature Collection of one element (consisting of Multipolygon geometries) into feature collection of many features

texts:
I took help from Converting Feature Collection of one element (Multi Point Geometry) into Feature Collection of different elements

I tried the code suggested in the above question but i get the following error: Error in map(ID=0):

  GeometryConstructors.Polygon, argument 'coordinates': Invalid type.
  Expected: List. Actual: Feature.

The feature collection has coordinates property which consists of List of geometries but its still giving the error. do i need to provide anything specific to apply the same code for Multipolygon as against Multipoint suggested in the question.

The FeatureCollection looks as follows when printed:

the code is as follows:

var mp = ee.FeatureCollection(result.map(function(p){
var poly = ee.Feature(ee.Geometry.Polygon(p), {})
return poly;}))

where result is my featureCollection.


texts: Hydrological conditioning using stream vector in GRASS?

texts:
I would like to make an hydrological conditioning of a DEM using an available stream network in vector form. I know r.hydrodem but it doesn't allow to use the stream vector. I need something to burn the stream network in the DEM to assure water flowing.

Is there a solution in GRASS?


texts: Finding highest path of DEM or contours file within area using ArcGIS Desktop?

texts:
I would like to create a polyline from one specified point to a target point. In between, I'd like to the vertices to fall on the highest point of elevation of that Y coordinate.

Essentially, this will be the highest path for that area - or highest profile for this extent.

Has anyone had luck doing this in ArcGIS Desktop?

Haven't been able to find a tool for this yet.


texts:
Can someone tell me how to add a ArcGIS FeatureServer to a QGIS plugin? I've found a post here on how to add maplayers which is simple enough but I can't find an description on how to do the same for FeatureServer layers. I'm trying to build a plugin to load a large number of various WMS, WFS, maplayer, featureserver datasets in QGIS without having to manually load them.

example:
https://services2.arcgis.com/chCSiGO4ORzXeSGk/arcgis/rest/services/Archeologie__Onderzoeken/FeatureServer

maplayer syntax that works:

iface.addRasterLayer("url='my_featuere_layer_url' layer='0'", "mylayer", "arcgismapserver")


texts: WMS layer not working in OpenLayers (Geoportal)

texts:
I tried to implement layers in WMS format in openlayers 6. This layer comes from the Geoportal of the Emilia Romagna region, but I cannot see it on the map. This is the code:
HTML:
<!doctype html>
<html lang="en">
  <head>`    <link rel="stylesheet" href="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.4.3/css/ol.css" type="text/css">

    <script src="Prova2.js"></script>
    <script src="https://unpkg.com/elm-pep"></script>
    <style>
      .map {
        height: 800px;
        width: 100%;
      }
    </style>
    <script src="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.4.3/build/ol.js"></script>
    <title>OpenLayers example</title>
  </head>

The editor does not let me insert the code with the div containing the map on this page <div id="map" class="map", but there is in the code.
JAVASCRIPT:
window.onload = init;

function init() {

  var layer2 = new ol.layer.Image({
    source: new ol.source.ImageWMS({
      url: 'http://servizigis.regione.emilia-romagna.it/wms/uso_del_suolo',
      params: {
        'LAYERS': '2017_uso_suolo_ed2020'
      },
      projection: "EPSG:4326",
      /* serverType: 'geoserver',*/
    }),
  });

  var map = new ol.Map({
    target: 'map',
    layers: [
      new ol.layer.Tile({
        source: new ol.source.OSM()
      }), layer2
    ],
    view: new ol.View({
      projection: "EPSG:4326",
      center: [12.42125, 41.838481],
      zoom: 7
    })
  })
}

I am a beginner.


texts: Overlap Analysis of two rasters (heat maps) in QGIS?

texts: I would like to analyse two different heat maps I have created from point data and to understand areas where hot spots overlap with each other. To add a back story I need to understand the overlap of areas of high level fishing activity and the spatial occurrence of marine animal species. Is there any way I can perform this analysis in QGIS?


texts:
I am making a document with descriptions of how to get to some points of interest. For each point I am making a document that contains two maps, one showing the general overview and the larger roads (Typically scale 1:50000 - 1:100000) and another map showing a more detailed map close to the point of interest and the smaller roads or paths to get there (typically scale > 1:10000)

As the exact area I want to use in each map depends on the way the road network is around the point, it i not possible to calculate the mapextent from the point, but I can draw two polygons for each point indiating the overview area and the local area for each point.

So my question is:

Is it possible to make an atlas where I am using two different data sets to define the area covered by two different maps on one page?

Or do I have to make two different atlas and combine the maps when I am making the final document?


texts:
I tried

>>>QgsVectorLayer.loadNamedStyle(r'C:\Some\Path\FOO.qlr')

Get this error:

Traceback (most recent call last):
  File "<input>", line 1, in <module>
TypeError: arguments did not match any overloaded call:
  QgsVectorLayer.loadNamedStyle(QString, bool) -> (QString, bool): first argument of unbound method must have type 'QgsVectorLayer'
  QgsVectorLayer.loadNamedStyle(QString) -> (QString, bool): first argument of unbound method must have type 'QgsVectorLayer'

Then i tried:

>>>vlayer = QgsVectorLayer
>>>vlayer = loadNamedStyle(r'C:\Some\Path\FOO.qlr')

and gets the same error.
The question is simple. How do i load a .qlr-file from path?

My research:
The developers archive,
The C++  documentation (which i suppose is good but I don't know C++),
Some one with a similar problem

I have a hard time finding Api for the py qgis.


texts: How to get shape with lines between points which have to lay on the same line entity?


texts: Not Able to See Feature Layer With Different Spatial reference On the API Base Map

texts:
I know ESRI is using following spatial reference for base maps

World_Street_Map (MapServer)
Spatial Reference: 102100  (3857)
102100  (3857)

Now I have a feature layer which is already published in this spatial reference

  Spatial Reference: 102671  (3435)

after adding the fc to map like

  var fl= new FeatureLayer("http://devesri01.domain.com:6080/arcgis/rest/services/ElectricVGA/FeatureServer/6");
  map.addLayer(fl);

but I am not seeing the layer in the map () No error Message either)? Do I have to publish the Future service with same spatial reference of the Base map?


texts: Creating a polygon around a "real estate" depending on my point "control" using QGIS

texts:
I am new to GIS and trying to learn as much as I can. Currently, I am not able to understand how one specific project works.
What I have:

Layer "Controls"
layer "Base Map of the city" (With polygons!)

What I want to do: I want the polygon of the real estate to change color depending on the entry "date" of my point (layer control).
Better: I want a new layer, which creates a polygon matching the real estate outlines where my "control" layer sets a point. The new polygon should change color depending on the entry "date" of my point (layer control)
I'll leave a screenshot of how I would like it to be:


texts: maskClouds vs filter clouds on Google Earth Engine

texts:
On Lab 6: Time series analysis from https://developers.google.com/earth-engine/edu, there is maskClouds:

var maskClouds = function(image) {
  var quality = image.select('BQA');
  var cloud01 = quality.eq(61440);
  var cloud02 = quality.eq(53248);
  var cloud03 = quality.eq(28672);
  var mask = cloud01.or(cloud02).or(cloud03).not();
  return image.updateMask(mask);
};

While in Lab 1: Intro to Remote Sensing and EE, clouds are filtered with a different approach:

// Sort the collection by a metadata property.
.sort('CLOUD_COVER')

// Get the first image out of this collection.
.first());

When should I use each of these approaches?


texts:
Is there a way to export my labels from AutoCAD into QGIS?

Right now i have to label everything again in QGIS. importing an dxf-file does not work. Only points, lines or polygons will be inserted.


texts: Plotting a shapefile directly from PostGIS using Python

texts:
I'm trying to plot a shapefile directly from PostGIS using Python. My shapefile is already stored in PostgreSQL. I know how to connect Python with PostgreSQL, but I just can't find anything that helps me plot my data into a map. I read somewhere that I should connect to PostgreSQL, query my shapefile table, select the geom attribute, store it into a geodataframe and then plot it.

Here's the code I'm using. Any ideas??

code example:

import psycopg2
import shapely
import osgeo.ogr
import shapely.wkt
import geopandas as gpd
import matplotlib.pyplot as plt

try:
    conn = psycopg2.connect("dbname='strokes' user='postgres'    host='localhost' password='****'")
except:
print "I am unable to connect to the database"
cur = conn.cursor()
cur.execute("CREATE INDEX bassin_index ON bassin USING GIST(geom)")
connection.commit()
cur.execute("SELECT st_astext(geom) AS wkt, fid_limite, codebassin FROM bassin")
rows = cur.fetchall()

rows_list=[]
for geom,fid_limite,codebassin in cursor:
    data={'codebassin':codeb,'fidb':fidlim,'geom':shapely.wkt.loads(geom)}
    rows_list.append(data)
    gdf=gpd.GeoDataFrame(rows_list).set_index('codebassin')
    gdf.head()
    gdf.plot(column='rows_list', scheme='QUANTILES', k=5, colormap='gray')
    plt.show()

conn.commit()
conn.close()


texts: Is there a way to upload non-spatial lookup tables in CartoDB?

texts:
I have a lookup table that only has 4 columns (bool, date, number, string) with 9,998,081 rows.

The CSV is relatively small (285mb), but it "blows up" when I upload it to CartoDB, I suspect, because of the other columns that CartoDB adds (the_geom, created_at, updated_at, etc).

It eats up the quota of my account and the table becomes 819.62 MB.

Would it be possible to tell CartoDB that this is just a lookup table that I'll use in joins, and I won't need these extra columns?


texts:
I am in trouble with the topology in a sde database.

I create a topology in a dataset from arccatalog, after managing to register it as aversion, I try to validate it from arccatalog, it tells me "The topology cannot be validated because it is registered as versioned.". If I don't register it as versioned I cannot start editing in arcmap. What should I do with the topology so that I can validate it and also start editing session.

I am using Arcmap 10.2 and sql server 2008 r2.


texts: Splitting of polygons and spatial join of lines to polygons in ArcGIS Desktop?

texts:
For the past couple of weeks I have been arduously creating flood depth rasters using flood zone polygons  and base flood elevation survey lines.

The issue is that some of the zones do not have a base flood elevation assigned. So far, I've just been using the lines to split the polygons and roughly shifting the lines a small distance diagonally, all at once, to perform a spatial join to the polygons. The lines are shifted because they are of course being used to split the polygons and are overlapping the splits. I then manually go through null values that resulted from a bad spatial join and manually change those.

So far, this is what I have been doing and everyone is satisfied with the result. It has proven accurate. The polygons are then used to create a raster.

Is there an easier way to do this? Not necessarily looking for a script.


texts: Merging hundreds of rasters in QGIS to create a basemap

texts:
I have a dataset of over 400 .tif-files, which I'm supposed to merge in QGIS into a single layer, representing a base map. These .tif files are pieces of a topography map (scale of 1 to 50000), and the file sizes are between 2 MB and 83 MB.

I'm relatively new to QGIS, and I have no idea what would be the smartest way to handle this (if there is any). Is it possible to create anything similar to ArcMap mosaic datasets in QGIS? I was specifically requested to find a way to do this in QGIS.

My research on the topic suggested that I could try to 1. build pyramids and then build a virtual raster or 2. use the merge function. Well, building pyramids for 400+ files seems painfully slow and manual.

Building a virtual raster without pyramids only gave me the error: "The process failed to start. Either the invoked program is missing, or you may have insufficient permissions to invoke the program". I have a write permission to the disk with the input and output locations, so that shouldn't be the problem.

Merge-function apparently only works for a few rasters at a time, not 400 of them.

Any ideas?


texts:
Is there any way to add Layer groups in ArcGIS Flex Viewer layer widget?
For example I have about 30 layers that are all other companies fiber networks. I would love for the user to be able to turn those all on or off at once. It is a burden for the end user to have to turn them all on and off one by one. The image below shows only some of the layers that need to be turned on individually, as you can see it would be cumbersome.


texts: Split polygons based on raster layer and defined target size

texts:
I have polygons with interpolated data like on the picture on right.

Now I want to split this polygon into smaller polygons with predefined similar area. Positions of these smaller polygons must be based on value from interpolated data below. So one smaller polygon should contains area with values for example 1-10.

Result should look like picture on the left.

Does anyone can help me or give idea how to do this in QGIS, python, R, javascript or any other GIS Software?


texts: Android devices and apps with support for Galileo positioning?

texts:
While the European's Galileo positioning system is only expected to be fully deployed by 2020, the first positioning made using already deployed Galileo satellites   was made back in 2013. Also, hardware devices supporting Galileo are already commercially available.

However, I am not aware of any Android smartphone supporting Galileo.

Are there any known, existing or planned, Android devices and android apps supporting Galileo positioning?

From what I understand, Galileo support might depend on support within the microprocessor itself. So the question regarding android devices supporting Galileo may boil down to a question about which microprocessor series support Galileo...


texts: Ingest 8-band satellite imagery into an array in C?

texts: I’ve been unable to locate a lightweight library implemented in C to ingest 8-band geotiff files.  GDAL can do this easily in python, though I’ve been unable to get the GDAL C libraries working; I’ve installed GDAL, but C still can’t find gdal.h.  Do you have any tips as to how to ingest 8-band images with another C library or tips on how to install and set the GDAL paths correctly?


texts: Editing different types of vocabulary in attribute table of shapefile in QGIS

texts:
I wanted to edit the name of the cities in the attribute table of a shapefile, for example for Turkey, there have cities like 'Eskişehir' and 'İzmir', How I can edit these names exactly the same in qgis?

I've copy-pasted the name for the city in Portugal 'A Coruña' was working (copy-pasted from Wikipedia), but for these cities of turkey weren't.

Is anyone knows how to edit them in qgis attribute table?


texts: ArcGIS Raster Calculator Con() set 0 values to null

texts:
In some vides appear a process separated and I would like to do this two steps in one step:

set all values greater or equal to a number to 0:

Con("MyRaster" >= 2800, 0 , 1)

the result is called "BandCloud.tif"

use created mask of 0 and 1 is to give 0 value to clouds and keep the rest of the data in "MyRater":

floar("MyRaster" * "BandClud")

the results is a ratser where are cloud values are equal to 0 and the rest is the data "MyRaster"

So, finally I get a raster of my data ("MyRaster") and 0 for clouds

But would liker to directly in the same first conditional (1) to set could values not to 0 but directly to null, I did it in a third step:

set 0 values to null:

SetNull("MyRaster" == 0, "MyRater")

How could this be done with raster calculator in ArcGIS?


texts:
Regarding my plugin pyqgis 2.2, i can load a vector layer succesfully. But when i try to zoom to the layer extent it does not work.

How to do in order to zoom to the layer extent?

I have just tried your proposition.I  pasted your piece of code into my plugin. But unfortunatly, when i load a vector layer. I try to zoom to the layer extent (thanks to right click -> zoom to layer extent). I have no error message but there 's no change zoom to the extent of vector layer.

You can read the python.py (github plugin link)  here :

https://github.com/picardie-nature/clicnat-qgis/blob/master/plugin.py

Any thoughts?


texts: ogr2ogr too many GCP points - exceeds maximum number of characters in the command line

texts:
I am trying to adjust a vector polygon dataset using ogr2ogr spline transformation with cmd commands like:

C:\\OSGeo4W64\\bin\\ogr2ogr.exe -f "ESRI Shapefile" C:\\path\\output.shp -order 3 -tps -gcp -13163.1627 271026.6107 -13274.2879 271201.236 -gcp 64464.2679538 -210398.439419 64390.5172491 -210240.678326 -gcp 98432.5310433 -181184.99898 98313.298804 -180985.232826 C:\\path\\input.shp

Also I try to do this inside Python:

import os
control_points = (
"-gcp -13163.1627 271026.6107 -13274.2879 271201.236 "
"-gcp -12706.1672 271461.5611 -12843.4862 271599.4093 "
"-gcp 98432.5310433 -181184.99898 98313.298804 -180985.232826 "
"-gcp 64464.2679538 -210398.439419 64390.5172491 -210240.678326 "
)
command = 'C:\\OSGeo4W64\\bin\\ogr2ogr.exe -f "ESRI Shapefile" C:\\path\\output.shp -order 3 -tps ' + control_points + 'C:\\path\\input.shp'
print command
os.system(command)

With a few points like this it works fine.
But actually I have much more control points, more than 1000. I need to use them all because the dataset covers the entire country. And when I write into the command all my GCPs, the Python script doesn't produce any output at all (no errors), and when I paste this command into the Command prompt it just cuts my command and then gives an error from ogr2ogr.

I found out that the maximum number of characters in the command prompt is 8191 but my command has 79792 characters.

So, how can I use ogr2ogr with hundreds of control points? Is there a way to store them in some batch file or variable? The GDAL page only says that -gcp option may be provided multiple times to provide a set of GCPs, directly in the command.


texts: How to query data from a vector map (shapefile)?

texts:
We want to obtain the specific place of the boundaries of a shapefile. We want to be able to get the points of each polygon corresponding to a specific place in order to use these data in an algorithm. GRASS is the GIS that we are using. We have the elevation map from SRTM and the vector map containing the boundaries from http://www.unsalb.org/. Can we query/retrieve the points contained in a specific boundary for all the polygons in the vector map?

Sorry first time to use a GIS


texts: How to copy features from a vector layer to another vector layer with only correspondent attributes in QGIS

texts:
I have two vector layers:

A first one with all cadastral parcels of a region. This layer has several attributes, it's only the geometry and 3 of the attributes I'm intrested in.

My second layer is a selection of parcels from the same region. This layer has 7 attributes. The first 3 attributes are identical to 3 attributes from the first layer with all cadastral parcels. The 4 other attributes are filled in afterwards.

How can I copy a parcel from the first layer to the second one, only with the geometry and matching attributes in the right field of the attribute table?


texts: Using Stamen toner-lines tiles as ArcMap base map?

texts:
Is it possible to use stamen toner-lines tiles as a base map in ArcMap? If so, how are they added as a base map?

Several Stamen maps are available through ArcGIS Online, though not the toner-lines tiles.

ArcBruTile can make external tiles usable as a base map in ArcMap, though the toner-lines tiles aren't one of the "supported maps".

In 2013, Kenneth Field wrote a blog post about using external tilesets in ArcGIS online maps, though the interface for uploading datasets has changed and some of the interfaces/options are no longer presented when following the steps.


texts: Force closest facilities take the same route

texts:
I have a routing problem that I intend to solve with ArcGIS network analyst. First, Closest facility works perfect for me to find the shortest paths to a number of demands. The results look like the left part of the figure:

For my analysis I would like to compare the left situation with the right one, when it is more cost-efficient to share the route. This far I have just painted the routes in the right part for clarity. In a normal Dijkstra I would just recursively adapt the weights of the links, setting the already chosen to a lower cost. However, I do not know how could I do this in ArcGIS closest facility.

Is there a way to modify the shortest path algorithm in Network analyst?


texts: converting vectors to raster (burning streams)

texts:
I converted my stream vector to a raster file (using polyline to raster conversion tool). I anticipated that I would get a No Data row in the raster file attribute table, so that I could use a conditional statement to burn the streams into the DEM:

con(stream_raster>0, (focalmin_dem - 2), original_dem)

However, since I did not get a No Data row, all this gives me is a a DEM img directly beneath the streams (not the remaining extent of the image).

What am I missing?


texts: OpenLayers.Layer.XYZ keep previous tile if zoom with no tile

texts:
I've got an OpenLayers web application with a custom XYZ layer which downloads tiles from the server itself:

OpenLayers.Layer.MOBAC_OSMDROID_ZIP = OpenLayers.Class(
          OpenLayers.Layer.XYZ, {
              name: "tiles base",
              transitionEffect: "resize",
              sphericalMercator: true,
              url: '../map/${z}/${x}/${y}.png',
              [...]
          });

The funny thing here is that I can use MOBAC to download selected parts of OpenStreetMap and serve them from my website. So far, so good.

The problem is: obviously I do not want to download the whole world from OpenStreetMap, but only, for example, zoom levels 0 and 1 for the world, 2-4 for Europe, 5-7 for Spain and so on. Then, when I am zooming in to Spain I get the right tiles. But if I move outside Spain, I just get the damn pink "no-tile" image. I can easily remove such no-image with the .olImageLoadError CSS class... but it still does not show a map.

What I would like is: if there is no tile, use the previous tile, scaled. This way you will probably have an overview of where you are when there are no specific tiles for such level.

Any idea? Is something like this included in OpenLayers? Any idea to implement it myself?

Since the "resize" effect does not remove the old tiles for a while, maybe keeping them "forever" would work...


texts: How can I convert a LineString GeoDataFrame to a Points GeoDataFrame while retaining the vertice order?

texts:
I have a LineString GeoDataframe that I am trying to convert into a Points GeoDataframe, but I want to retain the GroupBy and SortBy features inherent in a LineString (i.e., all the vertices that make up a line are grouped by some ID and sorted in a specific order).
A  similar question was asked here, but I don't understand from the answers (1) how to get my groupby/sortby requirement; and (2) why they use a one line function, it seems like there should be a cleaner way.
Below I have an example where I build a LineString from a Points GeoDataFrame, and I am basically trying to decompose it back to Points. In reality, I don't have the original Points GeoDataFrame, I just made one up here so that someone can have an easy copy/paste example to work with (per the questions guidelines).
Build Example LineString GeoDataFrame
%matplotlib inline
import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
from shapely.geometry import LineString
import pyproj
from pyproj import CRS

myid = [1, 1, 1, 2, 2]
myorder = [1, 2, 3, 1, 2]
lat = [36.42, 36.4, 36.32, 36.28, 36.17]
long = [-118.11, -118.12, -118.07, -117.95, -117.95]
df = pd.DataFrame(list(zip(myid, myorder, lat, long)), columns =['myid', 'myorder', 'lat', 'long'])
gdf_pt = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df['long'], df['lat']))
display(gdf_pt)
gdf_line = gdf_pt.sort_values(by=['myorder']).groupby(['myid'])['geometry'].apply(lambda x: LineString(x.tolist()))
gdf_line = gpd.GeoDataFrame(gdf_line, geometry='geometry')
gdf_line.crs = "EPSG:4326"
display(gdf_line)
ax = gdf_line.plot();
ax.set_aspect('equal')
ax.set_xticklabels(ax.get_xticklabels(), rotation=90);

Attempt
Below follows one of the answers from the linked question. It returns a Pandas Series, and I'm just not sure how to unpack it into a dataframe with GroupBy (based on "myid") and then create a SortBy based on the order.
mypoints = gdf_line.apply(lambda x: [y for y in x['geometry'].coords], axis=1)
print(mypoints)
print(type(mypoints))

System details:
Windows 10
conda 4.8.2
Python 3.8.3
shapely 1.7.0 py38hbf43935_3 conda-forge
pyproj 2.6.1.post1 py38h1dd9442_0 conda-forge


texts: QGIS - Split Vector Layer returns empty layer

texts:
I'm to split a vector layer by feature with "split vector layer" tool in qgis
I want to split the "switchboar" column

after I fill the tool

I receive the new empty layer

Does anybody know why I don't get different layers separated by the attribute?
I must say that if there is a NULL value I get only one layer with the NULL values and that's it.

I must say that I solved the issue by using SAGA software (split table tool)
but still it's less comfortable then using QGIS alone.


texts: Preserving nodata values when using gdal_translate -scale

texts:
I have four 12 - bit Bands in single images which I want to merge to one 8bit Geotiff using GDAL (GDAL 2.1.3)

First thing I do is building a VRT:

gdalbuildvrt -srcnodata 0 -separate bands_merge.vrt band1.jp2 band2.jp2 band3.jp2 band4.jp2

Then with gdal_translate I create the tif.

gdal_translate -of GTiff -ot Byte -scale 0 4095 0 255 bands_merge.vrt bands_merge.tif

That works just perfect. But as the nodata value of the source files is 0 and all values will be scaled, original values below ~8 will become 0.

My first thought on how to fix that issue was using -scale 1 4095 1 255 but that does not work because it assigns 1 to all previous 0 values.

How can I make sure that the -scale 0 4095 0 255 does not create false nodata values ?


texts: Openlayers draw multiple coordinates off KML data coming from the database


texts: How do I change an existing srs in QGis?

texts:
I am showing postgis data with srid 28992 (Amersfoort / RD New). But when I plot the geospatial data on a map (with openlayers plugin), the data is shown a hunderd meters off because the datum shift is not taken into account.

When I make a custom CRS with the correct datum shift and set my layer to that CRS is it shown correctly.

I want to update QGis' standard interpretation of 28992 so it takes into account this datum shift (+towgs84 parameter) correctly and when I open new data it is shown correctly immediately.

I found that QGis has a sqlite db called srs.db in /usr/share/qgis/resources/.  But even when I update the record for 28992 it still doesn't show the data correct on the map. When I inspect the layer in QGis I see it still has no +towgs84 parameter.

Does this mean that QGis does not actually use this sqlite db? Because the data I show is from a postgresql database I also edited spatial_ref_sys.sql to make sure that it also knows about the datum shift.

I also edited the epsg en esri files in /usr/share/proj/ just to be sure, but no effect.

At this moment I cannot think of another place where QGis gets its projection information from.


texts:
I'm getting this error when trying to merge 5 different files:
2020-08-25T15:52:53     WARNING    Traceback (most recent call last): File "C:/Users/.../AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\interference_zone\interference_zone.py", line 362, in interference_calculation
              processing.run("gdal:merge", parameters)
              File "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\tools\general.py", line 106, in run
              return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)
              File "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\core\Processing.py", line 137, in runAlgorithm
              raise QgsProcessingException(msg)
             _core.QgsProcessingException: Unable to execute algorithm
             Incorrect parameter value for INPUT

This is the code I am running:
parameters = {
        'INPUT': final_result_layers,
        'PCT' : True,
        'SPERATE' : True,
        'OUTPUT' : outputfile_path
    }
processing.run("gdal:merge", parameters)

The final_result_layers is a list of type QgsRasterLayer.


texts:
I want to obtain the quantity of precipitation for each municipality in Brazil on a given year. For that, I'm trying to use IDW interpolation on QGIS to interpolate the precipitation data from weather stations to municipality centroids. I have two shapefiles, one for the weather stations (with station long-lat and total precipitation data) and another with municipality polygons (from which I obtain the centroids). I'm new to QGIS, but I'm already able to perform an IDW interpolation for the entire map, but cannot find a way designate specific values to the centroids.

If that is possible on QGIS, how do I do it?


texts: Redistribute linestring vertices to specified distance

texts:
I use a groundwater model program called GMS, and it has a very useful function "Redistribute vertices" (more info). I use this function all the time to both simplify and provide regular spacings between points (vertices) in linestrings.

See this image from GMS where the black lines are original stream with irregular spacings, and the grey lines have redistributed vertices using 50.0 m

I haven't seen the equivalent function anywhere else in the GIS-processing world (Esri, JTS, PostGIS, etc.). Redistribute vertices requires a specified spacing distance, whereas simplify-type functions use a tolerance and don't use a specified spacing. Does anyone know if there is something equivalent in either ArcGIS or open source GIS projects?


texts: Trying execute lyr to kml conversion on a loop ArcPy - Only the first one works ok

texts:
I have a few folders, and inside each folder i have a few shapefiles. My task is duplicate these folder hierarchy, converting the shapefiles to kml. But only the first kml generated works ok. All the others doesn't work. I can generate all of than, on a loop, but they don't work, and no error is showing.
folder = r'D:\EXTERNO\20200721_Gustavo_Feitosa_AGRODATA\01_shp_para_KML'
folder_to_KML = r'D:\EXTERNO\20200721_Gustavo_Feitosa_AGRODATA\02_kml_FINAL'

import os

for mats in os.listdir(folder):
    list_info_folder = list(map(str, mats.replace('_', ',').split(',')))
    reg_code = list_info_folder[len(list_info_folder) - 1]
    for files in os.listdir(os.path.join(folder, mats)):
        if files.endswith('.shp'):
            if mats not in os.listdir(folder_to_KML):
                os.makedirs(os.path.join(folder_to_KML, mats))

            if "{}.kmz".format(files.replace('.shp', '').replace('.lyr', '')) not in os.listdir(os.path.join(folder_to_KML, mats)):
                import arcpy
                ###### criando_pastas_GIS.temp_composta --> it is a folder
                wrk_spc = criando_pastas_GIS.temp_composta

                arcpy.env.workspace = wrk_spc
                arcpy.env.overwriteOutput = True

                ### CAMINHOS ###
                shp_path = os.path.join(folder, mats, files)
                kml_path = os.path.join(folder_to_KML, mats,
                                        "{}.kmz".format(files.replace('.shp', '').replace('.lyr', '')))

                add_area = arcpy.AddGeometryAttributes_management(shp_path, 'AREA_GEODESIC', '#', 'HECTARES',
                                                                  coord_system)
                add_centroid = arcpy.AddGeometryAttributes_management(shp_path, 'CENTROID', '#', 'HECTARES',
                                                                      sirgas_2000)

                area = 0
                latitude = 0
                longitude = 0
                arquivo = ''
                with arcpy.da.SearchCursor(shp_path, ['AREA_GEO', 'CENTROID_Y', 'CENTROID_X']) as search:
                    for itens in search:
                        area += itens[0]
                        latitude = itens[1]
                        longitude = itens[2]
                    del itens
                del search

                # definindo projeção
                def_proj = arcpy.DefineProjection_management(shp_path, coord_system)

                # Lyr shp
                lyr = arcpy.mapping.Layer(shp_path)

                arcpy.MakeFeatureLayer_management(lyr, "{}.lyr".format(files.replace('.shp', '')))
                arcpy.SaveToLayerFile_management(lyr, "{}.lyr".format(files.replace('.shp', '')))
                f = open(arquivo_txt, 'a')

                f.write("{}, {}, {}, {}, {}, {}\n".format(files.replace('.shp', ''),
                                                          files.replace('.shp', '').replace('Lote_', ''), reg_code,
                                                          area, latitude, longitude))
                f.close()

                # Layer
                arcpy.LayerToKML_conversion("{}.lyr".format(files.replace('.shp', '')), kml_path, '0', 'NO_COMPOSITE', '#', '1024', '96', 'CLAMPED_TO_GROUND')


texts: Select from Attribute Table with step by 2

texts: I want to select from the Attribute table all items with step by 2 (select one item, I don't select the next one, etc). How can I do that? I'm using ArcMap 10.2.2.


texts: Using gdal_proximity inside conda environment

texts:
I want to batch gdal_proximity using python.
The following command is working inside OSGeo4W Shell:
C:\python ./gdal_proximity.py -srcband 1 -distunits PIXEL -values 1 -ot Float32 -of GTiff "[INPUT LAYER PATH]" "[OUTPUT PATH]"

But I get an error when using following script inside conda environment:
import gdal
import subprocess
subprocess.run(['python', './gdal_proximity.py', '-srcband', '1', '-distunits',
                'PIXEL', '-values', '1', '-ot', 'Float32', '-of', 'GTiff',
                '"[INPUT LAYER PATH]"', '"[OUTPUT PATH]"'])

The error:
ERROR 4: `"[INPUT LAYER PATH]"' does not exist in the file system,
and is not recognized as a supported dataset name.
Unable to open "[INPUT LAYER PATH]"


texts: Adding constant (0) image to masked image in Google Earth Engine?

texts:
I'm attempting to add one constant image iEmpty to a masked image iCTSurfaceWater, such that the value of the resulting image at a particular pixel is equal to iCTSurfaceWater, and regions where iCTSurfaceWater is mapped are equal to the constant value in iEmpty. I've tried adding, unmasking, and mosaicking, all without luck.

How can I accomplish this?

var findConnecticut /*:FeatureCollection */ = function() {

  var fcUSWatersheds = ee.FeatureCollection('ft:1IXfrLpTHX4dtdj1LcNXjJADBB-d93rkdJ9acSEWK');

  var fcLowerConnecticut = fcUSWatersheds.filterMetadata("name", "equals", "Lower Connecticut");
  var fcUpperConnecticut = fcUSWatersheds.filterMetadata("name", "equals", "Upper Connecticut");
  var fcConnecticut = fcLowerConnecticut.merge(fcUpperConnecticut);
  fcConnecticut = fcConnecticut.union();

  return fcConnecticut;
}

var fConnecticut = ee.Feature(findConnecticut().first());
Map.addLayer(fConnecticut, {}, "Connecticut River Basin");
print(fConnecticut);

var iCTSurfaceWater = iGlobalSurfaceWater.clip(fConnecticut);
print(iCTSurfaceWater);
Map.addLayer(iCTSurfaceWater, {}, "Water in CT");

var iEmpty = ee.Image(0).clip(fConnecticut);

print(iEmpty);

//var icCTAll = ee.ImageCollection([iEmpty, iCTSurfaceWater]);
var iCTAll = iCTSurfaceWater.add(iEmpty);

print(iCTAll);

Map.addLayer(iCTAll, {}, "All of CT");


texts: Two files with same coordinate systems should overlap but don't

texts:
I'm a complete beginner and I have been trying to solve this problem for hours now without success. I know there are other questions like this one and I've tried to apply the answers to my problem, but just don't seem to be working. So I have a TIFF file which is a map of my country and a shapefile which contains coordinates. Both files are in the same coordinate system - MGI 1901 Slovene National Grid, but they just don't overlap in ArcMap.
I tried Project and Define projection tools to try and change the coordinate system of the shapefile with my coordinates, but nothing works. If it helps - in QGIS they overlap perfectly without any kind of manipulation beforehand.
Does anyone have any idea of what else I can try?


texts: ArcGIS Field Calculator Not Working With Date Field

texts:
I am trying to calculate a field based on the date values in another field. Here is my code from the pre-logic script code box:

def getHalf(dfield):
  if dfield >= date '2017-08-27 00:00:00' AND dfield < date '2017-08-27
  12:00:00':
    return 1

The field calculator fails on line 2. The input field is formatted as a date field, its values look like this when viewed in the attribute table:

8/27/2017 10:05:00 PM

However they look like this when I perform a select by attributes on the field, so I am pretty sure I got the syntax right:

date '2017-08-27 22:05:00'

And then of course this is the line of code attempting to calculate the field:

getHalf(!timestamp_cst!)


texts: Creating shapefile with features in GISPro?

texts:
We cannot figure out how to create a shapefile that contains the features seen within the app GISPro (by Garafa) on an iPad.

I have been able to export a shapefile, but there are no features attached to it.
I have been able to create a KML/KMZ, and the features are there but buried within a ton of HTML and all crammed together.

If it helps, the original data was created in GISKit (also by Garafa) but we upgraded when we couldn't get useful exports hoping that the shapefiles (which we need to end up with) would auto create the features.

Any ideas?


texts: Random forest classification using Sentinel 2 data in Google earth engine

texts:
I am trying to do RF classification in GEE but i am getting an error
Classifier (Error)
Classifier.train, argument 'classProperty': Invalid type. Expected type: String. Actual type: List. Actual value: [name, id]
The code which i have entered is as follows:
/Loading Sentinel 2 data and filtering by date and bounds
var Sentinel2Collection=ee.ImageCollection("COPERNICUS/S2_SR")
    .filterDate('2019-01-01','2019-01-28')
    .filterBounds(geometry)
    .sort('CLOUD_COVER', true)

// Filter the collection by the IMAGE_QUALITY property.
var filtered = Sentinel2Collection
  .filterMetadata('IMAGE_QUALITY', 'equals', 9)

var Sentinel_med = Sentinel2Collection.median()
var Sentinel_crop = Sentinel_med.clip(geometry)

print(Sentinel_crop)

Map.addLayer(Sentinel_crop, {bands: ['B8', 'B4', 'B3'], gamma: 2}, 'S2 Image');

//Collecting training data
//Merging the three geometry layers into a single featurecollection
var newfc= forest.merge(urban).merge(agriculture);
print(newfc, 'newfc')

// Select the bands for training
var bands = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7','B8','B8A'];

// Sample the input imagery to get a FeatureCollection of training data.
var training = Sentinel_crop.select(bands).sampleRegions({
  collection: newfc,
  properties: ['name','id'],
  scale: 20
});
print(training, 'training')

// Making a Random Forest classifier and training it.
var classifier= ee.Classifier.smileRandomForest(10).train({
  features: training,
  classProperty: ['name','id'],
  inputProperties: bands
});
print(classifier)

//Classifying the input imagery
var classification= Sentinel_crop.select(bands).classify(classifier)

// Define a palette for the Land Use classification.
var palette = [
  'D3D3D3', // forest (0)  // red
  '0000FF', // urban (1)  // green
  '008000' //  agriculture (2) // blue
];
print(classification)
// Display the classification result and the input image.
Map.addLayer(classification, {min: 0, max: 2, palette: palette}, 'id');


texts: What is the most common way of displaying geodata from Postgis on Leaflet?

texts: For example, in my database I've some points with geographic coordinates. I want to display them on Leaflet JS layer. Which instruments should I use? Maybe there is some API...


texts:
I am trying to project a shapefile in QGIS that I downloaded from:
http://www.dartmouthatlas.org/downloads/geography/hrr_bdry.zip

As downloaded the shape files are not projected. The geographic coordinate system is North American Datum 1983. The problem I'm having is that the US is horizontally stretched, and from what I've read, this problem is because the shapefile is unprojected. This is my first experience working with GIS. I've tried Vector-->Data Management Tools-->Define Current Projection and changed Output CRS to North_America_Albers_Equal_Area_Conic, but this made the map disappear.


texts: GDAL project lat lon coordinates to UTM in C

texts:
I have dug up some old code I wrote 6 years ago that was doing the conversion from Lat Lon coordinates to UTM.

I have tried to run it today and the code returns NaN values.
I suspect there has been some significant changes with GDAL.
I am using GDAL 3.0 on Arch Linux.

Function latlon2utm takes X1, Y1 lat/lon coordinates and returns UTM coordinates X2, Y2

#include "ogr_api.h"
#include "ogr_srs_api.h"

void LatLon2UTM(double *X1, double *Y1, double *X2, double *Y2)
  {
  int zone;
  int err;

  printf("%f %f\n", *X1, *Y1);
  OGRGeometryH pt1;
  OGRSpatialReferenceH SREF;
  OGRSpatialReferenceH UTM;

  SREF = OSRNewSpatialReference(NULL);
  OSRImportFromEPSG(SREF, 4326);

  zone = (int)(((*X1 + 180.0) / 6.0) + 1.0);
  UTM = OSRNewSpatialReference(NULL);
  if(*Y1 < 0) OSRSetUTM(UTM,zone,FALSE);
  if(*Y1 > 0) OSRSetUTM(UTM,zone,TRUE);

  pt1 = OGR_G_CreateGeometry(wkbPoint);

  //printf("%f %f", *X1, *X2);
  OGR_G_SetPoint_2D(pt1, 0, *X1, *Y1);

  OGR_G_AssignSpatialReference(pt1, SREF);

  // Convert to UTM
  err = OGR_G_TransformTo(pt1,UTM);
  printf("%d", err);

  *X2 = OGR_G_GetX(pt1, 0);
  *Y2 = OGR_G_GetY(pt1, 0);
  printf("%f %f\n", *X2, *Y2);
  }

Is There anyone that could give me the proper way to do that?


texts: Converting Lines to Polygons in QGIS produces Sliver Polygons

texts:
While trying to polygonize the polyline layer in PyQGIS using

processing.runandload("qgis:linestopolygons",explode_path,polygon_path)

command in  python console, it does not create polygons corresponding to the line edges. Instead it creates sliver polygons along the curved edges. Why the polygons are not created along the line edges?

On polygonization, the black lines are supposed to be the edges of the created polygons. Instead only few slivers are created along the curved edges(blue polygons).


texts: Normalizing addresses using PostgreSQL and PostGIS?

texts:
I am trying to geocode addresses in Hawaii using the PostGres/PostGIS with the 13.1. Tiger Geocoder. In Hawaii, many street addresses take the following format: ##-####. This is problematic with the geocoder, which is setting the street address as 0 when using the geocode function.

For example, I run the following block:

SELECT g.rating,
    pprint_addy(g.addy),
    ST_X(g.geomout)::numeric(8,5) AS lon,
    ST_Y(g.geomout)::numeric(8,5) AS lat,
    g.geomout FROM geocode(normalize_address('92-419 Makakilo Dr, Kapolei, HI 96707')) AS g

And the output is:
0 Makakilo Dr, Makakilo, HI 96707

Can anyone share a workaround so that the house address is accounted for in the geocoder?


texts: Extracting number values from attribute string with mixed data

texts:
In my attribute table, I have a column "Name" which contains data about type number and name of streets. For example: “SR 18 DI FIUGGI”. What I need is to create another column in which there will be only numeric data from column “Name”, so just "18".

I know how to do that in Excel, and that I can export/import data to do it there, but my guess is that you can do it directly in QGIS.

Also, one thing that can be useful to me is dividing columns by spaces. So there different columns for "SR", "18", "DI FIUGGI". Is this possible?


texts: Different results extracting data from the same point from MODIS HDF with R

texts:
I am trying to extract NDVI data from MODIS HDF files using R. I have several files downloaded from EarthData.
The file I am using in this example is here to download.
These MODIS HDF files have a sinusoidal projection, so I have to reproject the coordinates of my points of interest from WGS84 to the projection of the .hdf file or viceversa. Trying both I get different results for the same point. How can I know which one is correct?
Here is my code:
library(sp) # SpatialPoints(), CRS()
library(MODIS)
library(rgdal) # readGDAL()

# defining the coordinates of interest
site = cbind(-71, -41)
latlon = CRS('+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0')
site_wgs84 = SpatialPoints(site, latlon)

# I have several files but here I open just one as an example
f <- "MOD13Q1.A2000049.h12v13.006.2015136104516.hdf"

# open the file
sds <- getSds(f)

# in this case I am interested in the "250m 16 days NDVI" dataset
ndvi_sinus <- raster(readGDAL(sds$SDS4gdal[1], as.is=TRUE))

# transform the site coordinates to the sinusoidal projection
sinus <- crs(ndvi_sinus)
site_sinus <- spTransform(site_wgs84, sinus)

# extract the data divided by the scale_factor
data_ndvi_sinus <- extract(ndvi_sinus, site_sinus) / 10000 # I get data_ndvi_sinus=0.1912

# Now I reproject the raster to WGS84 and then extract the data
ndvi_wgs84 <- projectRaster(ndvi_sinus, crs=latlon)
data_ndvi_wgs84 <- extract(ndvi_wgs84, site_wgs84) / 10000 # I get data_ndvi_wgs84=0.2162236


texts:
I have some data in a deprecated projection that might not have an EPSG code. I I have asked the provider to tell me which coordinate system is and the answer come with those file, probably to use FME.

The .tf is like this:

TRANSFORMATION  stj_lok_sweref991800    2009-03-04 10:14:40 Kempe
Kempe         2009-09-03 GPASS 3.51
Inpassning:   Stj_lok_090903_2.k -> Swe99_lok_090903_2.k
Plan:         KDIFF, 819  punkter.
+   Restfelsinterpolering     TRIAD 4.0   2009-09-03 16:12
+
The parameters are based on a least squares fit using a Transverse Mercator  projection.
Ellipsoid: GRS 80
TRANSVERSE MERCATOR PARAMETERS
Central meridian                18   2   9.348000 degr. min. sec.
Scale along central merdian        1.000025550000
False northing                      -6455066.2870 m
False easting                         133496.8020 m
+
COMPUTATIONAL MODEL
7-parameter transformation vector formula:
   X(to) = DeltaX + (1 + scale_correction) * R * X(from)
where the order of multiplication of the compound rotation matrix R is RZ*RY*RX. RX, RY and RZ being the matrices expressing the rotations a
The numerical values of the seven parameters are:
Translation along X-axis:      -40.0458223192 m
Translation along Y-axis:      128.1224164328 m
Translation along Z-axis:        -.3377596899 m
+
 Överräkning mellan SYSTEM:
 SWEREF 99 lat long-> SWEREF 99 18 00
/
FSYSTEM Södertälje_lokalt/
AREA              70000              80000             130000             130000 /
TSYSTEM SWEREF 99 18 00/
TRIAD   stj_korr.tr1/
SYSTEM Södertälje_lokalt/
ELLIPSOID GRS 80/
PROJ Gauss
G     18   2  9.3480000000 DEG
  -6455066.2870
    133496.8020
    1.000025550000 /
SYSTEM  Fictive lat long/
LATLONG   DEG/
ELLIPSOID GRS 80/
STACK H  /* Save height
Z/
HTRANS ELLIPSOID 0/ /* h=0
GEOCART/
FORMEL CART
BW 7
      -40.0458223192      128.1224164328        -.3377596899
      717.9984439633      228.9888458955     1259.7748095849
         .0159239170
/
ELLIPSOID GRS 80/
CARTGEO/
RETURN H
Z/
SYSTEM SWEREF 99 lat long/
LATLONG   DEG/
ELLIPSOID GRS 80/
PROJ    GAUSS
G 18 00 0 DEG
0
150000
1
/
STOP /

Which migh have the information I need, but I do not find any specification of this file.

Is there a way to use the info in this file to build a custom CRS using the proj.4 library?

BTW, I do not need too much precision, so any fair approximation to the meter is fine for my purposes


texts: "Invalid JSON primitive: geometry" when adding features

texts:
I'm trying to use the arcGIS Rest API to add features to a feature layer, but I get an "Invalid JSON primitive: geometry" error.  I get the same results if I leave out the spatialReference.

The feature layer was created on the arcGIS Developer Site, setting the geometry type to "Points". The fields are store_id (integer), address (string), items (string), and item_desc (string).

What am I doing wrong?

Python 2.7 code

import requests

def test_upload(token):
    test_locations = [
        {
            'store_id': 79,
            'address': 'address1',
            'latitude': 48.4680466419,
            'longitude': -94.4650303558,
            'item_id': '2 3 1',
            'desc': 'bar, baz, foo'
        },
        {
            'store_id': 89,
            'address': 'address2',
            'latitude': 48.7006190216,
            'longitude': -115.9411092720,
            'item_id': '2',
            'desc': 'bar'
        },
        {
            'store_id': 161,
            'address': 'address3',
            'latitude': 48.7425504958,
            'longitude': -111.5460764360,
            'item_id': '3 1',
            'desc': 'baz, foo'
        }
    ]

    features = []
    for location in test_locations:
        feature = {
            'attributes': {
                'store_id': location['store_id'],
                'address': location['address'],
                'items': location['item_id'],
                'item_desc': location['desc']
            },
            'geometry': {
                'x': location['longitude'],
                'y': location['latitude'],
                'spatialReference': {'wkid': 4326}
            }
        }
        features.append(feature)
    data = {
        'f': 'json',
        'features': features,
        'rollbackOnFailure': False,
        'token': token
    }

    print '############'
    print 'The JSON request looks like:'
    printable_data = data.copy()
    printable_data['token'] = '<withheld>'
    print json.dumps(printable_data, indent=2, separators=(',', ': '))
    print
    print '############'

    resp = requests.post(FEATURE_URL + '/addFeatures', data=data)
    print 'The response from the server is:'
    print 'HTTP Status Code:', resp.status_code
    print json.dumps(resp.json(), indent=2, separators=(',', ': '))

Output
When FEATURE_URL is defined and test_upload is passed a valid token, this gives the following output:

############
The JSON request looks like:
{
  "rollbackOnFailure": false,
  "features": [
    {
      "geometry": {
        "y": 48.4680466419,
        "x": -94.4650303558,
        "spatialReference": {
          "wkid": 4326
        }
      },
      "attributes": {
        "items": "2 3 1",
        "item_desc": "bar, baz, foo",
        "store_id": 79,
        "address": "address1"
      }
    },
    {
      "geometry": {
        "y": 48.7006190216,
        "x": -115.941109272,
        "spatialReference": {
          "wkid": 4326
        }
      },
      "attributes": {
        "items": "2",
        "item_desc": "bar",
        "store_id": 89,
        "address": "address2"
      }
    },
    {
      "geometry": {
        "y": 48.7425504958,
        "x": -111.546076436,
        "spatialReference": {
          "wkid": 4326
        }
      },
      "attributes": {
        "items": "3 1",
        "item_desc": "baz, foo",
        "store_id": 161,
        "address": "address3"
      }
    }
  ],
  "token": "<withheld>",
  "f": "json"
}

############
The response from the server is:
HTTP Status Code: 200
{
  "error": {
    "message": "Cannot perform operation. Invalid operation parameters.",
    "code": 400,
    "details": [
      "'features' parameter is invalid",
      "Invalid JSON primitive: geometry."
    ]
  }
}


texts: Error TypeError: 'NoneType' object is not iterable

texts:
I am trying to batch process KML polygons from a folder into a temp geodatabase.
From there, Iif this works and runs then I am trying to Append the polygons to result to a Server connection (Postgres withe append layer). Going to use here the ModelBuilder when the script runs.
But I am stuck on a error

  Traceback (most recent call last):
        File "C:\Users\Bjorn.tenBroeke\Documents\Bjorn\KML.py", line 20, in
          for kml in arcpy.ListFiles('.KM'):
      TypeError: 'NoneType' object is not iterable
      Failed to execute (Script).

# Name: BatchKML_to_GDB.py
# Description: Converts a directory of KMLs and copies the output into a single fGDB.
#              A 2 step process: first convert the KML files, and then copy the featureclases

# Import system models
import arcpy, os

# Set workspace (where all the KMLs are)
arcpy.env.workspace= (r"C:\DATA\KML")

# Set local variables and location for the consolidated file geodatabase
outLocation = "C:/WorkingData/fGDBs"
MasterGDB = 'AllKMLLayers.gdb'
MasterGDBLocation = os.path.join(outLocation, MasterGDB)

# Create the master FileGeodatabase
arcpy.CreateFileGDB_management(outLocation, MasterGDB)

# Convert all KMZ and KML files found in the current workspace
for kml in arcpy.ListFiles('*.KM*'):
  print "CONVERTING: " + os.path.join(arcpy.env.workspace,kml)
  arcpy.KMLToLayer_conversion(kml, outLocation)

# Change the workspace to fGDB location
arcpy.env.workspace = outLocation

# Loop through all the FileGeodatabases within the workspace
wks = arcpy.ListWorkspaces('*', 'FileGDB')
# Skip the Master GDB
wks.remove(MasterGDBLocation)

for fgdb in wks:

  # Change the workspace to the current FileGeodatabase
  arcpy.env.workspace = fgdb

  # For every Featureclass inside, copy it to the Master and use the name from the original fGDB
  featureClasses = arcpy.ListFeatureClasses('*', '', 'Placemarks')
  for fc in featureClasses:
    print "COPYING: " + fc + " FROM: " + fgdb
    fcCopy = fgdb + os.sep + 'Placemarks' + os.sep + fc
    arcpy.FeatureClassToFeatureClass_conversion(fcCopy, MasterGDBLocation, fgdb[fgdb.rfind(os.sep)+1:-4])

# Clean up
del kmz, wks, fc, featureClasses, fgdb


texts: Batch add column to layers and add layer prefix in PyQGIS

texts:
I have a number of layers each with a prefix like A_ B_ and so forth.

I'd like to:

Add a column (text) called DENSITY to each shp
Add the prefix letter from the layer name to each row in the DENSITY column for each shp file.

Can it be done?

I'm using QGIS 2.8.1

Update:

I think there was a typo here with a missing +

print "- layer: " child.layerName() + "  ID: " + child.layerId()

and I ran this:

# Iterate through the Legend to find all groups and layers
root = QgsProject.instance().layerTreeRoot()

for child in root.children():
    if isinstance(child, QgsLayerTreeGroup):
    print "- group: " + child.name()

# If we find a group, save the prefix name

    prefix = group.name().split('_')[0]
    elif isinstance(child, QgsLayerTreeLayer):
    print "- layer: " + child.layerName() + "  ID: " + child.layerId()
    layer = child.layer()

# Add code to add a field named density to the layer.
#######
# Add code to iterate through each feature in the layer and populate the field with prefix value
#######

but got these errors:

  line 4, in
      root = QgsProject.instance().layerTreeRoot()
  AttributeError: 'QgsProject' object has no attribute 'layerTreeRoot'


texts:
I have used GeoServer for rendering my own maps.
Can GeoServer be used to serve geocoding / reverse geocoding requests (i.e. you put a street address and it gives you coordinates and vice versa)?
I am using non-US addresses.


texts:
I have problem with python in QGIS. A module I have installed (AROSICS if it matters) uses a gdal_array.OpenNumPyArray function. It however crashes my QGIS when I try to run a script including functions from AROSICS. I have tracked down the error to the function and below is the (not)working example, which crashes QGIS when run from python console:

import numpy as np
from osgeo import gdal_array

array = np.array([[1,2,3],[4,5,6],[7,8,9]]) # a test array
ds = gdal_array.OpenNumPyArray(array) # command, which crashes QGIS

I'm using Windows 8.1 on Lenovo Y580 and QGIS 2.18.11 installed, alongside everything else (python, gdal, bindings...) with the newest version of OSGEO4W installer.
The problem persists after a reinstall of QGIS and on a completely different machine running Windows 10, so I assume it's something within GDAL. From a research on the net I got to know, that I should have a _gdal_array.pyd file somewhere in the installation and I do, however a _gdal_array.dll is missing. Is it important? What should I do to fix the problem?


texts: How to save a layer file with a relative path

texts:
In VBA, how do you save a layer file (.lyr) with a relative path?

I have some code that loads a shapefile into ArcMap but I am unsure of how I go about using this shapefile to then save a layer file to disk that has a relative path set.

The idea is after I have run my code and output all the layer files I can move the shapefiles and their associated layer files from a local drive which I use for processing my data to another location where other people can use it.


texts: Two rasters created from the same polygon have different extents?

texts:
My problem is that I am trying to export some rasters to Rockworks so that I can create a sand layer from them. Rockworks requires the grid files i use to have the exact same XY co-ordinates.

So i have a polygon of my sand layer - I create one raster of the surface (by using extract by mask from some AHN data I have) and the second one i create myself with a depth of -10.

When i check their layer properties however they are different! e.g. the top extent of the lower raster is 503002.235094748 and the top extent of the upper raster is 503000.235094748. The left and bottom extents are exactly the same..

They're both in the same co-ordinate system..?!


texts: Geoserver automatic histogram manipulation for WMS requests?

texts:
I'm having this weird thing when making two different WMs requests. The requests differ by the requested BBOX, meaning two different image-regions, yet the requested regions do overlap.
And now for the weird part; It seems that the resulting images have been through some kind of an histogram manipulation (equalization perhaps...) resulting the overlapped area between the two requests to be in different colors.

The problem worsens when making client-side applicative tiles via WMS, resulting in a bizarre transition between tiles.

Hence the assumption that GeoServer automaticly decides to make some-kind of manipulation for my WMS requests.

I've looked for a solution online and the only thing related was GeoServer's RasterSymbolizer which can be configured to make histogram manipulations if configured so for rasters in their SLD file.

Looking at my raster's SLD file shows it's only configured with the opacity option (and nothing about the histogram / normalize options).

Bottom line, Is my assumtion correct? and if so, what am I suppose to do to correct it (how to stop GeoServer from messing with my rasters?)


texts: Controlling number of gdal_translate cores when started via R's gdalUtils

texts:
I am processing satellite data in R and need to merge sat images of different UTM zones, which is where I ran into unexpected problems. I think I have solved them, but before moving to the production server I'd like to (have to) regain full control over the number of cores that are used. In addition to my own parallelizations, GDAL also runs commands in parallel / multiple times. Is it possible to set the number of workers that are being started?

In detail, I have the call

require(gdalUtils)
mosaic_rasters(gdalfile = list.of.two.files, dst_dataset = ftemp,
               gdalwarp_index = 2, of = "GTiff")

This results in all of my available cores starting an instance of gdal_translate. I want to limit this.

I could not find anything about this topic in the GDAL configuration files, and passing co = "NUM_THREADS=2" has no effect - rightfully so I think.


texts: Convert QgsVectorLayer.extent() to EPGS:4326

texts:
Many APIs such as overpass turbo require bounding box notation like (x_max, y_max, y_min, x_min) in EPSG:4326.
However the method QgsVectorLayer.extent() returns values in the source's CRS. A possible solution of course is to reproject the QgsRectangle as described here long time ago. But I wonder:

Is there a way to easily return these coordinates in EPSG:4326 without having to reproject all the data i want to derive extents from?

I mean Qgis nativliy derives these coordinates when using the 'Identify Features' tool and also by glancing at the sourcecode of the the plugin 'QickOSM' it seems capable of doing so.


texts: How can I run QGIS Flowmapper in batch mode?

texts:
I installed and used the Flowmapper plugin in QGIS. It works perfect. However, I have numerous matrices that I need to create desire lines. Doing it one by one is pretty tedious so I wanted to run them in batch. For this, I tried writing a python function which calls the function

shapefilemaker(typeofcalculation,CreateShpNodes,IncludeNodeNames,fulldirectorystring,outputfilename,outputfilenamenodes,fulldirODinput,fulldirPTinput,fulldirPTNamesinput,combotext)

inside flowpyv07.py However since it lacks proper documentation, I'm having some difficulty using it. Has anyone tried this?


texts: Combining PDFs in Folder using ArcPy for ArcGIS Pro?

texts:
This was able to combine two pdfs in a folder into one pdf, but unable to combine more than this.

import arcpy, os

# Set the workspace.
arcpy.env.workspace = r"C:\Output"

#Set file name and remove if it already exists
pdfPath = r"C:\Output\Map.pdf"
if os.path.exists(pdfPath):
    os.remove(pdfPath)

# List all of the pdf files in the output folder.
# pdfList is a Python List returned from the ListFiles function.
pdfList = arcpy.ListFiles("*.pdf")

#Create the file and append pages
pdfDoc = arcpy.mp.PDFDocumentCreate(pdfPath)

for f in pdfList:
    pdfDoc.appendPages(f)

pdfDoc.saveAndClose()

del pdfDoc


texts: visualising edge_graph from list pgr_dijkstra graph

texts:
I am creating a query that finds the shortest path from dangerous crossing nodes to hospitals (many-to-many). I have the proper output from pgr_dijkstra, with output as follows

seq | path_seq | start_vid | end_vid | node | edge | cost | agg_cost
-----+----------+-----------+---------+------+------+------+----------
   1 |        1 |         2 |       3 |    2 |    4 |    1 |        0
   2 |        2 |         2 |       3 |    5 |    8 |    1 |        1
   3 |        3 |         2 |       3 |    6 |    9 |    1 |        2
   4 |        4 |         2 |       3 |    9 |   16 |    1 |        3
   5 |        5 |         2 |       3 |    4 |    3 |    1 |        4
   6 |        6 |         2 |       3 |    3 |   -1 |    0 |        5
   7 |        1 |         2 |       5 |    2 |    4 |    1 |        0
   8 |        2 |         2 |       5 |    5 |   -1 |    0 |        1

I am trying to use st_makeline to create seperate lines for each individual node path, however thus far the function has been creating one big line
SELECT  st_transform(st_makeline(t1.the_geom),4326)
        FROM pgr_dijkstra(
                'SELECT osm_id as id, source, target, cost FROM spatial.melbourne_edges',
                array(select hospital_id from hospitals_in_vicinity),
                array(select crossing_id from dangerous_crossing_nodes), false) as route
                inner join spatial.melbourne_nodes t1 on route.node = t1.id --);

How am i able to make a column of all the different paths?


texts: Mapbox GL AddLayer: Where are the icon images coming from

texts:
A lot of the examples always load the icon image like this (For example here: https://www.mapbox.com/mapbox-gl-js/example/geojson-markers/):

map.addLayer({
    "id": "markers",
    "type": "symbol",
    "source": "markers",
    "layout": {
        "icon-image": "{marker-symbol}-15",
        "text-field": "{title}",
        "text-font": ["Open Sans Semibold", "Arial Unicode MS Bold"],
        "text-offset": [0, 0.6],
        "text-anchor": "top"
    }
});

Where is this icon stored and how can I link to a local png or svg? Or how can I put my own markersinto source? The example is not very well documented.


texts:
I have a world map shapefile and I can only seem to get two results, both of which are half right.

The first, using a projection like EPSG:3857, I get a map that doesn't cut the east tip of Russia off, but does include a weird bar streaming across the whole earth:

On another projection, I get this, but with the Russia getting cut up problem (which I think looks a bit odd):

This seems to replicate no matter what Shapefile I use, so I'm guessing it's the projection?


texts: Getting name of executing tool in ArcPy?

texts:
I created a script in a custom ToolBox and I defined a name (here ImportTopojis):

This name is displayed when you execute the tool :

I simply want to get this name because I want to use a single Python file for multiple tools and choose the proper class to use depending on this name.


texts: How do I force gpspipe to read from a USB Serial device rather than a server address with a port and everything

texts:
I'm desperately trying to grab the simplest piece of information I can out of my usb dongle from the command line. A program of mine needs to use it so it has to be a STDIN/STDOUT type solution.

I have a GlobalSat ND-100 which has a fix on my position (that Google Earth can and is currently using) which is very accurate. So I know the device is working. I have a mac.

I've tried these techniques so far to try to grab my current location:

-read directly from the usbserial (fails. binary? data)

  cat /dev/cu.usbserial
  #for some weird reason this doesn't return anything right now...it did before. anyway, it doesn't work.

-use gpspipe (installed along with gpsd) to grab the sentences

sudo gpspipe -rs /dev/cu.usbserial
#=> gpspipe: could not connect to gpsd localhost:2947, Unknown error: -6(-6)

So gpspipe is looking for a server, not a serial port. I need it to use the serial port. Or maybe I could host the usb serial port on a server. How might I do that?


texts: Search for geojson polygon in mongodb implementing a custom vector tileserver

texts:
I am building a vector tileserver in node.js using mongo 2.4.7 as database.

My database contains geojson data from multiple sources and they are in many layers. It is important to render the layers in the correct sort order.

Since I have little, or no control of the client rendering the individual tiles I have to clip the geojson data before returning it to the client so one “lower lying” polygon does not hide another higher lying tile in an adjacent tile.

I have created an simple rectangle ->polygon clipper that works for now - but this should be improved later on. If anyone knows an efficient clipper library that is usable from node.js please let me know. This is, however, not the issue.

My real challenge is the finding the correct polygons within the database.
I am currently using both $geoIntersects and $geoWithin in my query.

Let say that we have a polygon A that in zoom level x just fits within one tile bounding box (tile A). Here the Polygon will be returned by the $geoWithin clause and all is well.

Going several zoom levels in, the area from the tile A now fits within 9 tiles in a square. How do I find Polygon A when the server needs to render tile number 5 (the middle one)? And yes - in this situation I would like to render a square covering the complete tile due to the clipping.

I have created a drawing to illustrate my challenge.


texts: How to determine on which side of a line a polygon feature falls?

texts: I have parcel data intersecting line data.  In the parcel data there are some parcels that don't intersect the line.  How could I programmatically figure out if the non-intersecting parcel is on the Right side or the Left side of the line?  Thanks.


texts:
I've tried install Qgis 3.0 (and 2.18) using OSGEO. I'm using windows 10.

I've followed the steps I found in youtube (https://www.youtube.com/watch?v=G1WyUGYHC0I and https://www.youtube.com/watch?v=4NCtprMYyg4).

Summing up, OSGEO acomplish download, but the installation process always stops at 93% (sometimes 94%). In the window is showing:

Running
          qgis
          /etc/postinstall/qgis.bat

In Task Mannanger, for OSgeo4w-setup-x86 is indicated:

    1- Console Window Host
    2- NirCmd
    3- osgeo4w-setup-x86
    4- Windows comad processer (sorry for pour translate!)

I finished last task (4), and then this mensage was noticed:

Postinstall script errors These do not necessarily mean that affected packageswill fail to function properly but please check /var/log/setup.lo.full and report problems
    Package: gdal-python
             gdal-python

The installation was completed, Qgis is starded, but I can not open my projects.

Can someone help me solve this problem and install qgis?
I think this problem is related to .bat. Will be?


texts: Display XY coordinates (with different lengths) in ArcGIS does not work

texts:
I have a list of Excel lat long coordinates in a spreadsheet that I want to import as points in ArcMap by displaying the XY Data.

lat      lon
5235746   486124
5236799   487447
5238042   489453
523686    486487
5236528   49301

However, ArcMap does not seem to recognise/display my coordinate data as points. I think because there are some 6-numbered coordinates in there as well. How should I tackle this problem?


texts: Merging polygons that share one vertex (node neighbors) using ArcGIS Desktop?

texts:
I have a set of very large shapefiles and would like to merge polygons within these that share one vertex.

Currently they are all separate like so:

What tool can I use to merge all the polygons that share one vertex so in the above example the selected polygon would be part of the one to the left like so:


texts: Python script exports only first page of composer template

texts:
In a python Qgis plugin, with this code i create a pdf file, exported from a template generated with "print composer" tool:

    myMapRenderer = self.iface.mapCanvas().mapRenderer()
    # Load template from file
    myComposition = QgsComposition(myMapRenderer)
    myFile = "/...mypath.../pdfTemplate.qpt"
    myTemplateFile = file(myFile, 'rt')
    myTemplateContent = myTemplateFile.read()
    myTemplateFile.close()
    myDocument = QDomDocument()
    myDocument.setContent(myTemplateContent)
    myComposition.loadFromTemplate(myDocument)
    printer = QPrinter()
    printer.setOutputFormat(QPrinter.PdfFormat)
    printer.setOutputFileName("/...myoutputpath.../out.pdf")
    printer.setColorMode(QPrinter.Color)
    printer.setFromTo(0, 0) # convention for whole document
    pdfPainter = QPainter(printer)
    paperRectMM = printer.pageRect(QPrinter.Millimeter)
    paperRectPixel = printer.pageRect(QPrinter.DevicePixel)
    myComposition.setNumPages(3)
    myComposition.render(pdfPainter, paperRectPixel, paperRectMM)
    pdfPainter.end()

I do not understand why the template is composed of 3 pages but the generated pdf file contains only one, the first one.

odd thing: the number of pages of the composition is 3

  print myComposition.numPages()  --> output 3

Which command is missing to export all pages of the document?

Is it maybe a bug of the library?