texts: Query country border polygon from OpenStreetMap

texts:
Is it possible to download the polygons of countries from OpenStreetMap using Nominatim web API? If so, how should I query its points?

In the example responses from Nominatim on OpenStreetMap wiki there is a polygonpoints property (attribute) on place object (tag), but it's not always present in search results when querying for countries. Also, what happens when the country "consists" of more than one polygon - several islands, for instance, can I query all of them?

If it's impossible to do with Nominatim, is there any other freely available API with such functionality?


texts: CartoDB - Non-animated heatmap, but won't let me layer under

texts:
I have a simple heatmap layer that is not animated. I want to move a choropleth layer above that heatmap layer, but Carto will not let me, as it thinks my Heatmap layer is animated.

Is there any way I can get around this problem?


texts: Issue importing CSV with CSVT (fields format not applied - QGIS 3.18)

texts:
I'm trying to import a CSV file with some data and hoping that QGIS will properly format the fields type but so far no luck.
I have my data in a file called coords.csv and I made a single line text file for fields format called coords.csvt (both files located in same folder)
My CSV is like this:
picnum,east,south,stamp
1,145.214879,-15.473314,06/05/2021 3:53:12 PM
2,145.214854,-15.473339,06/05/2021 3:53:13 PM
3,145.214822,-15.473369,06/05/2021 3:53:14 PM
4,145.214786,-15.473408,06/05/2021 3:53:15 PM

My CSVT is like that:
"Integer","String(12)","String(12)","String(50)"

To import: Layer -> Add Layer -> Add Delimited Text Layer ;
I then select my CSV file, set the coordinates and CRS, and then add the layer.
But when checking fields format, I get:

When I am expecting to see this:

What am I doing wrong?
I tried using small caps for the CSVT file but no luck either. I'm using QGIS 3.18.2 on Windows 10


texts:
JTS Topology Suite TestBuilder supports WKT, WKB and GML inputs. What's a quick and easy way to get the geometry for a KML feature into one of these formats for input into JTS TestBuilder?

For clarification, I'd like to avoid using OGR, FME, etc. which will likely do their own validation and potentially mess with the original coordinates. I'm thinking I'm more looking for a dumb text-to-text transformation tool.


texts: Database connection fails when using db authentication, succeeds using OS

texts:
When i attempt to make a connection to a database on my SQL Server 2014 sever via Database connections in ArcCatalog using with db authentication i get

"Connection was attempted with an older version of SQL Server Client...."

error message. However this is incorrect as the client is up to date. I download and re-installed the 2014 Client.

Also when i use "OS Authentication" i can connect just fine.

ArcGIS 10.3, ArcGIS for Server 10.3, SQL Server 2014.

The account in question is not a account in our AD (active directory) just and account in the db
resolved.


texts: How to use processing *.py files manually?

texts:
I want to rasterize vector layer by using Python code. On QGIS 3, there is no rasterize tool. Others are not applicable for me. However I found rasterize.py under this folder:

C:\OSGeo4W64\apps\qgis-dev\python\plugins\processing\algs\gdal\rasterize.py

Is there any way using these py file manually to rasterize vector layers?

NOTE:
I have already tried to define these parameters and run:

processing.runAndLoadResults('gdal:rasterize', parameters)

But it does not work!


texts:
I am trying to change an plugin from QGIS 1.8 to 2.0.
In the 1.8 version I create a table of the field names like this:

layer = vlayer
provider = layer.dataProvider()
feat = QGSFeatrue()
allAttrs = provider.attributeIndexes()
provider.select(allAttrs)
fieldNames = []

while provider.nextFeature(feat)
    fields = layer.pendingFields()
    for (k, field) in fields.iteritems():
         fieldNames.append(field.name())

I tried this for the 2.0 version with

layer = vlayer
fieldNames = []

for feature in layer.getFeatures():
    fields = layer.pendingFields()
    for (k, field) in fields.iteritems():
         fieldNames.append(field.name())

But pendingFields() does not work for me. Is there now a new method to get the names of the fields?


texts: Creating Buffer zones that fill area of postcode

texts:
I created a Map with Geocoded Addresses from the UK. The Basis for the map is a UK districts shape-file, this constitutes my main layer. After Geocoding (via .cvs) 3 Addresses I have a new layer with the three Address-Points.
The goal is now to create a Buffer around the Address Point that spans or fills out the entire area of the postcode to which that Address belongs
In the picture I included an example with a normal Buffer (just a circle with radius = x) the goal would be to fill the whole postal area with the Buffer.


texts: Resources to integration with external web service at desktop level

texts:
I am using ESRI based GIS software, I looking a way to integrate ESRI solutions with external web service of other information system. For start I would like to join feature class to external tables and change there symbology dynamically by the on line data. later i would like to develop custom applications for the collection, conversion, maintenance, and analysis of data. Mainly visualization of external data in maps when data is streamed online.

1. is it possible to connect JSON web service to ArcGIS desktop without ArcGIS Server ? I'm using version 10.1 sp 1 Advance level

2. Where can i find resources to integrate GIS to external web service at desktop level ? I have no idea on where to start and to proceed further.

3. What ArcMap Tools and Commands host the trigger for calling the services ?


texts: OpenLayers doesn't show an overlay with a different SRS than the base layer

texts:
I need to add an overlay layer that uses EPSG:4326 to a map with a base layer that uses EPSG:32628.

Every time I add the overlay OpenLayers automatically switches the SRS param of the overlay to EPSG:32628 (SRS of baselayer); I can see the base layer, but the overlay isn't shown.

My code is this:

baselayer = new OpenLayers.Layer.WMS(
               baselayername, baselayerurl,
               {
                  width: varmap.size.w,
                  srs: 'EPSG:32628',
                  layers: baselayername,
                  height: varmap.size.h,
                  styles: style,
                  format: format
               }
               );

overlaylayer = new OpenLayers.Layer.WMS(
               overlaylayer name, overlaylayer url,
               {
                  width: varmap.size.w,
                  srs: 'EPSG:4326',
                  layers: layername,
                  height: varmap.size.h,
                  styles: style,
                  format: format,
                  transparent: true
               }
               );

varmap.addLayers([baselayer,overlaylayer]);

How can I do this?


texts:
I hope this question is not too general. I'm dealing with a lot of data from vehicles (1.000.000 GPS points per day) and I want to snap the data to a big city's street network.
The timely order is not so important, so I don't need corrective filtering using HMM like in other implementations. I just need to snap the point to the nearest road.

I want to use Spark (pySpark) for this Problem and have

gps_df: lons and lats
edges_df: road network as shapely LineStrings (each Road is a new line and loaded as WKT to the DataFrame)

My first try was using osmnx using the get_nearest_edge function, which is way too slow for the amount of GPS points.

Now I want to do this completely in Spark to make it fast and running parallel, but don't know where to start. Do I need kdtree? Which algorithm to chose (should I use shapely's nearest_point, which uses rtree according to another post, or kdtree which is faster, or completely build it myself)?


texts: Creating a raster surface from two vector shapefiles

texts:
I am preparing raster surfaces to do a multi criteria evaluation in ArcGIS 10. I need to create a raster surface that contains distance values for each cell from the electricity grid. However, my data for the grid is contained in two shapefiles - a point file (e.g. poles, pylons) and a line file (electricity lines). I can't merge the vector files because they are different types.

Is there a way of creating a raster surface where each cell value is the distance from any part of the grid; e.g. distance from a point or a line?


texts:
I want to transform a ee.Iamge that descripbe the ecozone of a specific country to a vector. I try to reproduce this example but I cannot get anything else than black vector...

//load the country
var country = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017').filter(ee.Filter.eq('country_co', 'CG'));

//load the ecozones
var gez_2010 = ee.Image('users/bornToBeAlive/gez_2010_wgs84')
var country_gez_2010 =  gez_2010.select('b1').clip(country)

country_gez_2010 has only one band b1 with int values from 0 to 21.
//construct the zones
var zones = country_gez_2010.eq(0)
var i
for (i=1; i<21;i++){
    zones = zones.add(country_gez_2010.eq(i))
}

//Convert the ecozones to vectors.
var vectors = zones.addBands(country_gez_2010).reduceToVectors({
    'geometry': country,
    'crs': country_gez_2010.projection(),
    'scale': 100,
    'geometryType': 'polygon',
    'eightConnected': false,
    'labelProperty': 'ecozone',
    'reducer': ee.Reducer.min(),
    'maxPixels': 1e9
})

//display it on the map
var display = ee.Image(0).updateMask(0).paint(vectors, '000000');
Map.addLayer(display, {palette: '000000'}, 'vectors');

I'm doubting on the zones definition and the reducer selection but maybe it's something else


texts: Creating python function to merge line segments if they are connected using ArcPy

texts:
I have a polyline shapefile consisting of several line segment features. I want to do a merge operation if the line segment features endpoints are connected, and if the features are disconnected do not merge it. Put it in separate.
So my Logic is -
fetch the co-ordinates of the vertices from the ROADS feature class & then using Polyline method and Insert Cursor create a new feature which will create a new road having all the co-ordinates; but I am getting an error.
Here is my code to fetch the co-ordinates & make new polyline
import arcpy

arcpy.env.workspace = r'C:\Users\Drona_Maps\Desktop\Office Work\Python task\Geodatabase\Road_Geodatabase.gdb\Road_featuredataset\ROAD'
arcpy.env.overwriteOutput = True
fc = r'C:\Users\Drona_Maps\Desktop\Office Work\Python task\Geodatabase\Road_Geodatabase.gdb\Road_featuredataset\ROAD'

with arcpy.da.SearchCursor(fc, ['OID@', 'SHAPE@']) as cursor:
    for row in cursor:
        array1 = row[1].getPart()
        for vertice in range(row[1].pointCount):
            pnt = array1.getObject(0).getObject(vertice)
            listOfCoordinates = [pnt.X,pnt.Y]
            print (listOfCoordinates)

array2 = arcpy.Array([arcpy.Point(listOfCoordinates)])

spatial_reference = arcpy.SpatialReference(4326)
singlepart_polyline = arcpy.Polyline(array2, spatial_reference)

cursor.insertRow([singlepart_polyline])

The Error I am facing is:
[-144.37499999999997, -63.020833332999985]
[-109.54166666699996, -26.812499999999943]
[-109.54166666699996, -26.812499999999943]
[-15.583333332999928, 26.354166667000072]
[-15.583333332999928, 26.354166667000072]
[104.50000000000006, 61.18750000000006]
[112.75000000000011, 64.85416666700007]
[163.16666666700007, 100.14583333300004]
Traceback (most recent call last):
  File "C:/Users/Drona_Maps/Desktop/Office Work/Python task/Initial try/2nd Method.py", line 16, in <module>
    array2 = arcpy.Array([arcpy.Point(listOfCoordinates)])
  File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\arcobjects\mixins.py", line 1082, in _init_
    setattr(self, attr, value)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\arcobjects\_base.py", line 89, in _set
    return setattr(self._arc_object, attr_name, cval(val))
RuntimeError: Point: Input value is not numeric

Process finished with exit code 1

I am doing this in Pycharm using Python version 2.7 and arcpy module and using ArcGIS Desktop 10.4.


texts: Issues importing feature class to SDE SQL Server 2012

texts:
I have received an error indicating there could be issues with my geometry when attempting to import a feature class into SDE. I ran the check geometry tool. No issues.
 I have removed all fields except for my lat/lng, which is a double type, my geometry type, which is point, and objectID.
When doing this, the feature class imports without error.

What could be the issue?

My first guess is having to do with Nulls?

When an objectID doesn't allow null values, does this mean that the feature should be free from nulls in other fields also?

I am using ArcGIS 10.2 and sql server 2012 for my SDE connection.


texts: Project vector in epsg 3347 to wkt in Python using GDAL

texts:
I am new to working with geospatial data. I am currently using gdal and rasterio to work in Python.
I have a raster with the following crs definition:
GEOGCS["Coordinate System imported from GRIB file",DATUM["unnamed",SPHEROID["Sphere",6371229,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST]]

I believe this is a WKT projection definition.
I also have a vector (geojson) in EPSG:3347 that I would like to project into the same crs as my raster.
I modified the below code from the GDAL cookbook, however it is not correctly projecting my vector. I would like to reproject the vector to the same crs as my raster and save the projected vector to my disk.
from osgeo import ogr, osr
import os
import gdal

driver = ogr.GetDriverByName('GeoJSON')

# get the input layer
inDataSet = driver.Open("B:/Canada-d8-100m/Canada-d8-100m.geojson")
inLayer = inDataSet.GetLayer()

inSpatialRef = inLayer.GetSpatialRef()

# loading projection
sr = osr.SpatialReference(str(inSpatialRef))

# detecting EPSG/SRID
res = sr.AutoIdentifyEPSG()

srid = sr.GetAuthorityCode(None)

print(srid)

# input SpatialReference
inSpatialRef.ImportFromEPSG(int(srid))

#read in wkt projection from raster
ds = gdal.Open("B:/dd.weather.gc.ca/model_gem_global/15km/grib2/lat_lon/00/000/CMC_glb_TMP_TGL_2_latlon.15x.15_2020052500_P000_Copy.grib2")
inSRS_wkt = ds.GetProjection()
outSpatialRef = osr.SpatialReference()
outSpatialRef.ImportFromWkt(inSRS_wkt)

# create the CoordinateTransformation
coordTrans = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)

# create the output layer
outputShapefile = "B:/projected2.geojson"

if os.path.exists(outputShapefile):
    driver.DeleteDataSource(outputShapefile)

outDataSet = driver.CreateDataSource(outputShapefile)
outLayer = outDataSet.CreateLayer("ProjectedPoints2", geom_type=ogr.wkbPoint)

# add fields
inLayerDefn = inLayer.GetLayerDefn()

for i in range(0, inLayerDefn.GetFieldCount()):
    fieldDefn = inLayerDefn.GetFieldDefn(i)
    outLayer.CreateField(fieldDefn)

# get the output layer's feature definition
outLayerDefn = outLayer.GetLayerDefn()

# loop through the input features
inFeature = inLayer.GetNextFeature()

while inFeature:
    # get the input geometry
    geom = inFeature.GetGeometryRef()
    #print(geom)
    # reproject the geometry
    geom.Transform(coordTrans)
    #print(geom)
    # create a new feature
    outFeature = ogr.Feature(outLayerDefn)
    # set the geometry and attribute
    outFeature.SetGeometry(geom)
    for i in range(0, outLayerDefn.GetFieldCount()):
        outFeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))
    # add the feature to the shapefile
    outLayer.CreateFeature(outFeature)
    # destroy the features and get the next input feature
    outFeature.Destroy()
    inFeature.Destroy()
    inFeature = inLayer.GetNextFeature()

# close the shapefiles
inDataSet.Destroy()
outDataSet.Destroy()


texts: Mapserver - dynamic contrast stretching for custom COLORRANGE's

texts:
In MapServer, is there a way to use the PROCESSING "SCALE=AUTO" directive to get dynamic (current extent) contrast stretching with custom colorranges? When I test my mapfile (relevant layer part below) with the Openlayers viewer it gets the desired color ramps but does not do contrast stretching for the current extent. If I take the CLASS'es out the greyscale output is contrast stretched according to current extent. I have tried taking out the EXPRESSIONS and/or the DATARANGES but then it only takes the first COLORRANGE and does not stretch it.

# raster layer
LAYER
    NAME "mapservtest"
    TYPE RASTER
    PROCESSING "SCALE=AUTO"
    STATUS ON
    DEBUG ON
    PROJECTION
        "init=epsg:3857"
    END
    METADATA
        #WMS stuff here
    END
    TILEITEM "location" #location is a column in postgis table with varchar of the filepath to each image
    TILEINDEX "index"
    FILTER (`[datetime]` = `2015-01-01T00:00:00Z`)
    OFFSITE 0 0 0
    CLASS
        EXPRESSION ([pixel] >= 25 AND [pixel] < 30)
        STYLE
            COLORRANGE  "#00ff00"  "#0000ff"
            DATARANGE 25 30
        END # STYLE
    END # CLASS
    CLASS
        EXPRESSION ([pixel] < 25)
        STYLE
            COLORRANGE  "#ff0000"  "#00ff00"
            DATARANGE 20 25
        END # STYLE
    END # CLASS
END


texts: What does ArcGIS have that QGIS doesn't?

texts: I recently installed Ubuntu 15.10 and I'm comfortable with this linux OS. I installed latest QGIS and it works just fine and something that I'm curious and keen to get answers from you is, what can possibly be in ArcGIS that there is no way I can have it in QGIS? i.e Map layers available?


texts: Graduated Symbol classification not working with a decimal range

texts: I am classifying an .adf wind resource map which I vectorised using the Polygonise tool (which worked very well). I have classified the data using the Graduated symbol renderer using Natural Breaks (Jenks) mode.  Unfortunately although the Range column will take decimals, only columns with whole numbers will render.  Does anyone know if there is a way to do this. Am using QGIS 1.7.3.


texts:
I use this python code to calculate a greatcircle:
https://code.google.com/p/pyroms/source/browse/trunk/pyroms/greatcircle.py?r=39
After initializing the greatcircle I can calculate points on the line.

I would like to do the same with an arbitrary smallcircle. I would like to define the smallcircle using its center location (in latlon) and one location on the line (also latlon), so not using a diameter. Say as an input I would use a certain point spacing in meters, how would I then proceed to calculate the other locations on the smallcircle line?

Thanks!


texts: How I can do event with two clicks using python add-in wizard?

texts: How I can do event with two clicks using python add-in wizard, I have found a way to do event with one click by adding a New Tool in the wizard, but I need  "first click To copy the attributes from the features and in the second to update another feature"?  Can I use (def onMouseDownMap(self, x, y, button, shift)) twice on the same tool?


texts: Batch Cell-statistics from 2 different folders using ArcPy with ArcMap

texts:
I'm using these codes to generate Cell-Statistics, I have data in four different folders but right now I'm using two folders just for the testing but I couldn't get the result. I'm a beginner level of Python user.
Can you help me out with this code?
I'm getting this kind of error again and again:
import arcpy
from arcpy.sa import *
arcpy.CheckOutExtension ("Spatial")
u'CheckedOut'
arcpy.env.overwriteOutput = True
arcpy.env.workspace = r'I:\Test_Delete\AQUA_D'
inws = arcpy.env.workspace
inws1 = r'I:\Test_Delete\AQUA_N'
outws = r'I:\Test_Delete\OUTPUT'
rasters = arcpy.ListRasters()
for r in rasters:
    prefix = os.path.splitext(r)[0]
    suffix = os.path.splitext(r)[1]
    prefix1 = prefix.rstrip("Aqua_D")
    raster2 = os.path.join(inws1, prefix1 + "Aqua_N" + suffix)
    print raster2
    outraster = os.path.join(outws, prefix1 + "mean" + suffix)
    print outraster
    outCellStatistics = CellStatistics([arcpy.sa.Raster(r), arcpy.sa.Raster(raster2)], "MEAN")
    outCellStatistics.save(outraster)

I:\Test_Delete\AQUA_N\clip1_A2002185Aqua_N.tif
I:\Test_Delete\OUTPUT\clip1_A2002185mean.tif

Runtime error
Traceback (most recent call last):
  File "<string>", line 19, in <module>
RuntimeError: ERROR 000732: Input Raster: Dataset I:\Test_Delete\AQUA_N\clip1_A2002185Aqua_N.tif does not exist or is not supported


texts: Scale Bar Element moves when using Data Driven Pages?

texts:
I have made many mapbooks using data-driven pages. However I've encountered something strange and wanted to ask to see if anyone might know what's going on or has seen this before.

I have a scale bar in the title block and as I step through each page the scale bar moves slightly (0.0389 inches). I move it back but every time I try to export, the resulting pdf has the scale bar moved by 0.0389 on the first page, 0.0778 on the second, etc.

Seems like a bug, but maybe I'm missing something.

Page 2:

Page 21:


texts: Writing conditional (if/then) statements into ArcGIS Field Calculator using Python Parser?

texts:
I have gotten stuck on seemingly simple ArcGIS field calculator python expression.

I obtained a great answer to a another similar question a few years ago: Using Range in Python expression of ArcGIS Field Calculator?.
Now I need to populate a new field based on 6 ranges of values (bedrock depth).

My attempt at writing this expression is below:

I just need to populate the empty field "bedrockdp" using the ranges in the toc called brockdepmin.

screen shot

Pre-Logic Script Code:

def !GRIDCODE! (value):
    if value >0 and value < 50:
        return "0 - 50 cm from the surface"
    else:
        return "n/a"

bedrockdep =

GRIDCODE (!bedrockdep!)


texts: Combining MapInfo records with a duplicate value

texts:
I have a MapInfo table with a lot of points in it. The points are clustered together - they represent the first part of a post code, so they're in a similar position but not identical.

What I would like to do is combine the duplicate rows and attach a single point to each row that represents the rough centre of the other points.

In fact, making the centroid for a record is easy, so my real question is how do you combine rows that have a duplicate field?

There are about 1.5 million rows so doing anything by hand is a no-no!

Thanks in advance!


texts: How to create a shapefile or gml from a Oracle database using SDE.ST_Geometry?

texts:
I have a ORACLE database 10g(Oracle Database 10g Enterprise Edition Release 10.2.0.3.0) with a SDE.ST_Geometry field and i want to make a program in python so i can create several shapefiles based on this ST_Geometry field and use other fields as well so i can fill my shapes with attributes.
I only have database access and i have no idea which arcgis is related.

How can i do this? Is there any suggestion?


texts:
A colleague asked me to install ArcPad 10 on his Trimble JUNO 3B.  He runs an application with scripts that produce messages when the GPS is turned off.  The script works fine on his Trimble Nomad - but on the new Juno he gets the following error:

VBScript runtime error 800A0046 "permission denied: 'MsgBox("message")'

The error quotes line 4534 column 2 which has nothing to do with the msgbox.

Basically the message that should pop up is something along the lines of 'you can't do that, there is no GPS fix'

I have no skill whatsoever with VB and have no desire to learn.  Any one know if there is a quick fix aside from telling my colleague to figure it out on his own?  Could this be an issue with the Juno using a different OS than the Nomad?  Or perhaps a different version of ArcPad?


texts: Performing supervised classification on old aerial imagery taken in 1950s

texts:
I am working on a project to track the expansion of the city's built area in the last 70-80 years, and I have to deal with scanned images taken before satellites were first used, and I can't seem to find any resources on how I can do any sort of pixel classification on them.
I am not an expert and I'm still learning.


texts: Generating Pseudo-Absences in spatialEco?

texts:
I want to generate pseudo absences based on kernel density estimate from a linear feature spatial line layer, where pseudo absences have a higher probability of being sampled near a line.

I have tried to do this using the pseudo.absence function in the spatialEco package; however, the input pattern has to be a point pattern and I don't see how to specify a bandwidth parameter for the KDE, which I would like to do.

I tried using spatstat to create a density.psp, which I did successfully and am happy with that output, but I cannot figure out how to generate pseudo-absence points based on the KDE.

My bootleg solution was to use the spsample function in sp to create points along my line features and then use that as the point pattern input to the pseudo.absence function. This seems to work; however, I want the sample of points to be based on the inverse of the KDE output from the function.

Is there a way to do this within spatialEco?


texts:
I'm using fleviewer 3.0 with wmts.
I've got no result.. [white screen]
from server-side i'm using both mapproxy.1.5a & geoserver2.2

my config.xml for flexviewer:

 <layer label="Gwc-wmts" type="wmts" visible="true" imageformat="png"
    url="http://xxx.xxx.xxx.xxx:8080/geoserver/gwc/service/wmts"
                servicemode="KVP" layerid="cite:my_layer"
                tileMatrixSetId="EPSG:4326"

                />

        <layer label="WMtS"
                type="wmts"
                url="http://xxx.xxx.xxx.xxx/mapproxy/service"
                layerid="my_layer"
                servicemode="KVP"
                imageformat="png"
                visible="true"
                tileMatrixSetId="grids_demo"
                />

I've succesfully display with openlayers.
how can i use wmts with flexviewer?


texts: How to use raster extract/clip raster in Qgis?

texts:
In arcgis, raster like a shp file, could extract another raster.But in Qgis,Raster->Extraction-clip, the mask file must be shp file.
My mask raster are compound of  many adjacent cells. If converted to shp file,
it look like some adjacent cells would combine as polygon. It will reduce extraction accuracy.


texts: Selecting only points which are intersecting line layer using QGIS?

texts:
I have large number of points(As a point layer) and few line(Roads basically), how do I remove points which are not falling on line layer.

I tried several methods like-

Join attribute by location,
Intersect
Clipping etc,

Some how QGIS gets crashed while running algorithm.

I am using QGIS 2.18.14 and QGIS 2.18.14


texts: Calculating lacunae of MCP to calculate home range?

texts:
I'm using QGIS 2.18.26 and I am trying to calculate the home range of howler monkeys I was tracking. Using the AniMove plug in I have managed to produce Minimum Convex Polygons (MCPs). However, I would like to calculate the 'Lacunae' (Contiguous areas within the MCP that were not entered by the monkeys) and am struggling to find information on how to do this in QGIS.

I am attempting to replicate the method used by Ostro et al (1999) in 'A Geographic Information System Method for Estimating Home Range Size' which can be found here--> https://pdfs.semanticscholar.org/72aa/0a32838490a106a3526c9bd017e978134492.pdf


texts:
I want to connect to this REST Service, from QGIS Essen, using its URI:

http://maps.six.nsw.gov.au/arcgis/rest/services/public/NSW_Property/MapServer/2?f=pjson

Any tips please?

I tried adding it as a Vector layer, protocal GeoJSON. I also tried the OGR method, constructed as follows; ogr2ogr -f GeoJSON NSW_Property_20160314.json "maps.six.nsw.gov.au/arcgis/rest/services/public/NSW_Property/…*&f=json" OGRGeoJSON


texts: Use field calculator QGIS: if value in one column, then value in new column

texts:
I am having troubles using the field calculator to get what I want.
One column in my attribute table is named 'Taxa' and has a set of species in it. Species are separated by #.

I want to make a new column using the field calculator that returns a value, based on the species in the 'Taxa' column.

An example of the column I want to use to calculate a value in the new column:

Below i have some examples of values per species:

Calluna_vulgaris =1
Baeomyces_rufus = 1
Porpidia_crustulata =2
Rhizocarpon reductum =2
Micarea_erratica = 2
Hypnum_jutlandicum = 1/3

To summarize: If the Taxa column says 'Calluna_vulgaris#Porpidia_crustulata#Micarea_erratica'

I want the new column to have the number 5.

I understand this is not the most convenient way to group species in a column, but I used the TomBio addon which does this automatically.


texts:
I want to use rgdal to plot a shapefile graduating the colour of features according to the values of one field (density). I couldn't find online detailed indications on how to do it and the package documentation doesn't explain the options for plot().

This is what I composed so far

library(rgdal)

shapefile = "shapefile.shp"
map <- readOGR(dsn = shapefile, layer = "shapefile", verbose=FALSE)

plot(map, xlim = c(6.70, 18.32), ylim = c(35.2, 47.6), border=NA, add=TRUE)

This is the style in QGIS I am trying to replicate:


texts:
I'm having a crash on startup that appears to've come from nowhere—I'd not used QGIS in a while and then tried to launch it and it crashed. This happened on both my work and home computers, so I suspect it's related to a system update (running MacOS 10.12.4); QGIS 2.18.2/2.14.3 (depending on computer) installed from KyngChaos along with python lib.

Here's the crash report I'm getting:

Process:               QGIS [1403]
Path:                  /Applications/QGIS.app/Contents/MacOS/QGIS
Identifier:            org.qgis.qgis2
Version:               2.18.2 (2.18.2 [exported])
Code Type:             X86-64 (Native)
Parent Process:        ??? [1]
Responsible:           QGIS [1403]
User ID:               502

Date/Time:             2017-04-24 11:01:35.517 +0100
OS Version:            Mac OS X 10.12.4 (16E195)
Report Version:        12
Anonymous UUID:        CE6880BA-8529-B952-7BBC-A72CAED27607

Time Awake Since Boot: 540 seconds

System Integrity Protection: enabled

Crashed Thread:        0  Dispatch queue: com.apple.main-thread

Exception Type:        EXC_BAD_INSTRUCTION (SIGILL)
Exception Codes:       0x0000000000000001, 0x0000000000000000
Exception Note:        EXC_CORPSE_NOTIFY

Termination Signal:    Illegal instruction: 4
Termination Reason:    Namespace SIGNAL, Code 0x4
Terminating Process:   exc handler [0]

Application Specific Information:
dyld: in dlopen()
BUG IN LIBTRACE: MH not found

What's going on here? How can I get QGIS to launch?


texts:
I have problems with the code, I do not know where to place the part with the popup. This does not work if I want to use the data attribute in the popup window:

var okruh1= L.geoJSON(trate, {
    color: "#3399FF",
    weight: 5 ,
    filter: function(feature, layer) {
         return (feature.properties.OBJECTID == "1" );
    },
});

okruh1.bindPopup(layer.properties.OBJECTID);


texts: How to get Drainage Divide and Streams from DTM

texts:
Does anyone know how you can get a map like this attached, having the digital terrain model?

I'm sorry, but I do not know what is the exact technical term in English for them.
I hope you understand from the images.
(maybe drainage divide, water divide, divide, ridgeline, watershed, water parting)
My professor called it compluvium and impluvium. But another professor colled it
articulation of the morphological type of the soil.
(in the Italian language: Carta dei compluvi e displuvi / Articolazione morfologica terreno )

Thanks in advance for the answers


texts: Debugging Indentation error on Python Parser of ArcGIS Field Calculator?

texts:
Would anyone please help me find out why I get an indentation error in line 3

def calc(f1,f2,f3):
remove_list = ['01ST', '02ND', '03RD', '04TH', '05TH', '06TH', '07TH', '08TH', '09TH']
    if f2 in remove_list:
        f2 = f2.replace('0','')
    else:
        f2 = f2
f1 = f1.strip()
f2 = f2.strip()
f3 = f3.strip()
f1 = str(f1).replace('.0','')
x = f1+' '+f2+' '+f3
return x

THe function is simply supposed to evaluate the if statement then concatenate fields f1, f2 and f3. Do I need to end the if statement somehow?


texts:
Is it possible to perform a cookie cut in R?

Firstly, I have a polygon shapefile of many polygon some of which overlap. I then have a second polygon shapefile with contiguous non-overlapping polygons. I want to establish the parts of the first (potentially overlapping) polygons that fall within each of the second (non-overlapping) polygons. The first shapefile is a number of layers of cookie dough to cut and the second is the cutter to apply. Alternatively I can use QGIS if people think that's the best way to go.


texts:
This is similar to this question regarding shifts in raster output, but the difference I'm experiencing is fairly large-- the output raster is shifted 15 KM to the south. Here's an image:

When viewed in ArcGIS, the original raster (green) claims to be in spatial reference GCS_Greenwich. When the properties of this coordinate system (Angular Unit, Prime Meridian, Semimajor Axis, Semiminor Axis, Inverse Flattening) are compared to WGS84, they are identical.

As such, I used the EPSG for wgs84 when to project the raster to web mercator:

gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3857 input.tif output.tif

I'm using GDAL 2.1.0 dev. Is the coordinate system the issue?

Update: using esri's CopyRaster_management function and selecting WGS_1984_Web_Mercator_Auxiliary_Sphere as the output coordinate system worked perfectly. I'd prefer to stick to GDAL however, because this will be running on a Linux server without arcpy.


texts: Getting different results using Select by Location with ArcPy vs ArcMap?

texts:
I have two different layers. One layer contains features that are supposed to sit directly on top of features in the other layer.

With ArcPy I'm doing a SelectLayerByLocation with "ARE_IDENTICAL_TO" and then a "SWITCH_SELECTION" in order to find which features from layer one aren't directly on top of layer two, and then updating the geometry accordingly.

However, ArcPy is giving me different (incorrect) results than using the Select by Location in ArcMap.

Granted the feature is only 0.01 ft away from the other, but it's still not identical.

When using ArcPy, the feature above with the S is selected, but when I use Select by Location in ArcMap it's not.

layer1_fl = arcpy.MakeFeatureLayer_management(layer1, layer1 + "_fl")
layer2_fl = arcpy.MakeFeatureLayer_management(layer2, layer2 + "_fl")

arcpy.SelectLayerByLocation_management(layer1_fl, "ARE_IDENTICAL_TO", layer2_fl)

Any idea why it's acting this way?

As suggested by @FelixIP

I added arcpy.env.outputCoordinateSystem = arcpy.SpatialReference(102726) and all worked as expected.


texts: OpenLayers behaves problematic after certain zoom level in Google Maps Satellite view

texts:
I have some polygon data in .shp format to show in Google Maps Satellite view. I can show them well until certain zoom level.

But after certain level (zoom level 17 or 18), OpenLayers layer and Google Maps Satellite layer displace each other like this.

We know Google Maps Satellite View changes it camera angle after certain zoom level and put it to 45 degrees. How can i resolve this with OpenLayers.


texts: Using DATA_DIR relative path in GeoServer?

texts:
I am using REST API with curl to populate geoserver with layers on the server. I am not setting GEOSERVER_HOME so Geoserver itself decides what the actual path is. However, when I download the populated geoserver to my laptop and run it with startup.sh script, layers are not published as the datastore location still aims to the server location (/var/www/user/...).

How do I make geoserver independent on the environment it is run in? I would like the datastores to use relative paths, so instead of /var/www/user/... they should be using data_dir/datastore_dir.

Is this possible?


texts: QGIS - Create Garmin Custom Map plugin - Python Error

texts:
I am trying to create a KMZ file using the Garmin Custom Map plugin and at times it fails to execute (and I haven't been able to establish a pattern). I get an Python Error, and the map will not be created.
I will be the first to admit, I don't fully understand how CRS works, and would love for someone to point me to a good tutorial or info on understanding CRS.

I create the map I want to convert to a KMZ and then use the plugin
I get the message (Warning - Project CRS differs from WGS84 (EPSG:4326))
Sometimes it executes and produces the map, and other times I get the below message
WARNING    Python error : An error has occurred while executing Python code: See message log (Python Error) for more details.
When the above happens the process of producing tiles stalls


texts:
I have a NetCDF dataset which I normally call using the xarray function. I have recently been doing some transect analysis on the dataset, which looks something like this (the area between the horizontal blue lines are the transect for analysis):

For now, I am only able to extract the data, which either falls within the vertical or horizontal transect (as they are easy to extract.

But now I want to select a certain shape within the data-set, on which I would like to analyse. And these shapes might not be symmetrical or in a straight line. Some examples are shown below:

Basically, these shapes can be irregular (with known co-ordinates of each point).

Is it possible to extract the dataset/values (with co-ordinates), specifically only for the region of interest?

These steps are easy to perform just using ArcMap's Clip function or Google Earth Engine clip function. But I am not able to use this in Python (as I want solely use Python for all steps). Can someone provide some suggestions on how to this? If anyone knows any package, that has great integration with xarray, that would be great.

Additional Info (after import rioxarray):

Properties of output raster:

import xarray as xr
import rioxarray as rx

Treecover  = xr.open_rasterio('/home/chandra/data/Treecover_MOD44B_2000_250m_AMAZON.tif')
[Output]:
<xarray.DataArray (band: 1, y: 32093, x: 20818)>
[668112074 values with dtype=float64]
Coordinates:
  * band     (band) int64 1
  * y        (y) float64 13.71 13.71 13.71 13.71 ... -58.35 -58.35 -58.36 -58.36
  * x        (x) float64 -81.38 -81.37 -81.37 -81.37 ... -34.63 -34.63 -34.62
Attributes:
    transform:   (0.002245788210298804, 0.0, -81.37613580017715, 0.0, -0.0022...
    crs:         +init=epsg:4326
    res:         (0.002245788210298804, 0.002245788210298804)
    is_tiled:    0
    nodatavals:  (nan,)

geometries = [
    {
        'type': 'Polygon',
        'coordinates': [[
            [-46.23140155225633, -21.53505449239459],
          [-44.91304217725633, -20.221175092759253],
          [-70.22554217725633, 1.5816072875439455],
          [-71.36812030225633, 0.5271132528460204]
        ]]
    }
]

Treecover_clipped = Treecover.rio.clip(geometries, Treecover.rio.crs)
[Output]:
<xarray.DataArray (band: 1, y: 10293, x: 11779)>
array([[[nan, nan, ..., nan, nan],
        [nan, nan, ..., nan, nan],
        ...,
        [nan, nan, ..., nan, nan],
        [nan, nan, ..., nan, nan]]])
Coordinates:
  * band         (band) int64 1
  * y            (y) float64 1.58 1.578 1.575 1.573 ... -21.53 -21.53 -21.53
  * x            (x) float64 -71.37 -71.36 -71.36 ... -44.92 -44.92 -44.91
    spatial_ref  int64 0
Attributes:
    transform:     (0.0022457882102988043, 0.0, -71.36665774687539, 0.0, -0.0...
    _FillValue:    nan
    grid_mapping:  spatial_ref

Treecover_clipped.plot()


texts: Converting point data into gridded dataframe for histogram analysis using R?

texts:
I am very new at using GIS data and only modestly experienced with R.  I've been reading about how to analyze spatial data using the spatial-analyst.net PDF book, so I'm not completely lost, but I thought I could describe my problem and people might suggest ideas.

I have a dataset with about 2000 measurements at different lat/long coordinates, although I will probably subdivide this dataset as the data was collected over 3 years and conditions changed over time.  Let's call the variable being measured "IP."

I want to create a map of IP in the full area under question using Kriging or some other interpolation method on the sample data.  Then I want to create a histogram measuring the amount of land in various IP buckets.  I will also need to create a histogram which shows the number of samples in each bucket (note a sample could have a higher or lower actual IP than what kriging predicts for its land).

I follow how to load the data into a SpatialPointsDataFrame and run a kriging analysis, where I'm having trouble is how to convert that data into a gridded dataframe so I can do the histogram analysis.

Any suggestions for converting points into grids?


texts: How to create proper PostGIS points out of GRIB2 file?

texts:
I'm trying to load NCEP NOAA weather data from GRIB2 format into PostGIS and then visualize it using QGIS etc.

I'm using wgrib2 library to convert the GRIB2 data to csv first.

wgrib2 -match "(:TMP:surface:anl:|:TMP:2 m above ground:anl:)" gfs.t00z.sfluxgrbf00.grib2 -csv fcast.csv

Then I copy CSV to PostGIS and create a geometry column:

COPY tmp (begindate, enddate, name, name2, value, lat, lon)
FROM
'/Users/nl/Work/ncep/fcast.csv'
WITH DELIMITER ',' CSV;
SELECT addGeometryColumn ('public', 'tmp', 'the_geom', 4326, 'POINT', 2);
UPDATE tmp SET lon = lon-360 WHERE lon > 180;
UPDATE tmp SET the_geom=ST_SetSRID(ST_MakePoint(lat,lon), 4326);

When looking at the data in QGIS, I see a strange result:

When I expect to see the regular grid - something like this:

I see, that someone has done it already here:
How to process GRIB2 Data for Visualization?

Can anyone give some details?


texts:
I need to vectorize the cell boundaries of a raster map and I have forgotten how to do this.

For instance, if I have a raster map of a region that is just map=1 with a cell size of 500m, I would get a map that shows a grid that is 500m on each side.


texts: Offsetting Overlapping Symbols in ArcMap

texts:
I am working on a recreational atlas in ArcMap 10.5 that highlights points along the coast including marinas, mooring spots, kayaking and diving points etc. The problem I am facing is that I have so many points that it is just way to crowded with the labels and symbols and the points overlap and it is hard to see what is actually in that spot. I have an example here.

I was wondering if there was any tool or option to offset the symbols like in labeling or if I will have to offset the layers symbols individually.


texts: Run processing only in layers with selection

texts:
First of all, I am really begginer with PyQGis...

I am trying to run a random point algorithm in different vector layers of a folder. But I dont want to do it for every polygon, only for selected ones. The selection is done with a geometric filter and with an attribute filter.

When I run the code it is working perfectly for the layers in which the filters selects different polygons, but in those in which the selection gives no answer (because they are out of the geographical filter) the algorithm runs on every polygon... and I need exactly the oposite, not to run in any polygon.

I have tried a lot of things but nothing seems to work... Here is the code:

import os
import processing

#Selection of layers
path="/home/basi/Escritorio/GEA Forestal/SATree/MFE/Todos/Seleccion/"
capas = os.listdir(path)
for capa in capas:
    if os.path.splitext(capa)[1]==".shp":
        layer=iface.addVectorLayer(path+capa,"_","ogr")
#Filters
        Area=QgsRectangle(-4,39,-4.2,44)
        exp = QgsExpression('FORM_ARB_D ILIKE \'%pino albar%\'')
        request = QgsFeatureRequest(exp)
        request.setFilterRect(Area)
        features = layer.getFeatures(request)
#Selection
        ids = [i.id() for i in features]
        layer.setSelectedFeatures( ids )
#Random point algorithm
        Output_layer="/home/basi/Escritorio/GEA Forestal/SATree/Eleccion poligonos prueba Solo Sentinel/Programacion seleccion/parcelas_longitud_PS_"+os.path.splitext(capa)[0]+".shp"
        processing.runalg('qgis:randompointsinsidepolygonsfixed',layer,0,3,0.0006,Output_layer)


texts:
I am looking for a way to create a raster grid from a vector map.

I have a map consisting of shapefiles which I want to export as a raster with UTM coordinates. I was thinking of something using the identify tool in QGIS but not really sure how to apply it.

I want a file which has x,y coordinates and feature id as the variables. Does anyone have any suggestions?


texts: How do I calculate acreage grouped by owner, masked by a poly?

texts: I have statewide land ownership data and have created a mask to project ownership for only one county.  I need to calculate the acreage of each owner within that county.  Is there a quick way to achieve this?  I'm using ArcGIS 10.


texts:
I am trying to export a specific part of a dataset, based on an attribute query.

My dataset has a field REGION

and I would like to select all features where the region is equal to 2.

Here is my command:

ogr2ogr -f "esri shapefile" -clipdstwhere "REGION = 2" africa.shp world.shp

Although an output is created it is an exact copy of the original. In other words, the query does not keep the command from running but it is not correctly executed nevertheless.

I have tried different variations (putting REGION in double quotes and escaping them etc., but I am unable to get it to work).


texts: OL3 Geoserver WFS error OWS:ExceptionReport

texts:
I am attempting to load and draw a WFS layer from geoserver. It is installed on a Windows 2012 Machine. I have enabled CORS and I am able to draw and query WMS layers without fail. Below is my syntax for the layer

var vectorSource = new ol.source.ServerVector({
  format: new ol.format.GeoJSON(),
  loader: function(extent, resolution, projection) {
    var url = 'http://gis.halker.com:8080/geoserver/sandbox/wfs?service=WFS&version=1.1.0&request=GetMap&layers=sandbox:junk__poly&styles=&bbox='
     + extent.join(',') + ',EPSG:4269';
    $.ajax({
      url: url,
      dataType: 'jsonp'
    });
  },
  strategy: ol.loadingstrategy.createTile(new ol.tilegrid.XYZ({
    maxZoom: 19
  })),
  projection: 'EPSG:3857'
});

var loadFeatures = function(response) {
  vectorSource.addFeatures(vectorSource.readFeatures(response));
};

// Vector layer
var vectorLayerFixed = new ol.layer.Vector({
  source: vectorSource,
  style: new ol.style.Style({
    stroke: new ol.style.Stroke({
      color: 'green',
      width: 2
    })
  }),
  name: 'vectorLayerFixed'
});

I believe this is correct because in firebug I am getting errors that lead me to believe something is being requested. the error is as follows.

<ows:ExceptionReport xmlns:ows="http://www.opengis.net/ows"

I am really up a creek here because this is well beyond my understanding. I have looked over stack exchange and can't find anything I can make sense of.


texts: How to use ArcView GIS extension in QGIS or how to convert one to do so?

texts: I really would like to know whether is a way to can work with ArcView GIS extensions at QGIS? or at least a way to be converted to plugins at QGIS?


texts: Virtual layer with spatial query crashes in QGIS

texts:
Working with QGIS 3.10. Details:

This should be really simple, but I am unable to solve it:
I've got a point layer with two numeric fields "max" and "min", and a polygon layer. I want to create a Virtual layer with the polygons and the sums of "max" and "min" of the points that are within each polygon. All my tries around the following code result in QGIS crash
SELECT polygon.name,
       polygon.geometry,
       SUM(point.max) AS maxim,
       SUM(point.min) AS minim
FROM polygon
JOIN point ON st_within (point.geometry, polygon.geometry) = 1
GROUP BY polygon.name, polygon.geometry

My layers have shapefiles as their data source.


texts: Installing 3rd party python libraries for QGIS on Windows

texts:
How can I use 3rd party libraries on QGIS plugins on Windows?

I've developed a plugin that uses rasterio and numpy for a customer, but he's having problems installing rasterio and numpy.

Actually rasterio and numpy were installed in it's main system Python (C:\Python27), but I need QGIS Python to recognize it.


texts:
I have issue while exporting PDF, when i export PDF including Parcel layer its working fine but when add other layers in same file its not editable on PDF editor.
1-Pic: editable PDF with parcel layer
2-Pic : PDF along with other layers is not editable


texts: ListFields not working from Python Toolbox

texts:
I created a small arcpy script tool which works fine when I run it from ArcMap.  All it does is list the fields from a specified layer.

import arcpy
myLayer = arcpy.GetParameter(0)
fieldList = arcpy.ListFields(myLayer)
for field in fieldList:
    arcpy.AddMessage(field.name)

However if I create a Python Toolbox with this same script, it fails on fieldList = arcpy.ListFields(myLayer) with this error:

Traceback (most recent call last):
  File "<string>", line 37, in execute
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\__init__.py", line 1131, in ListFields
    return gp.listFields(dataset, wild_card, field_type)
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\geoprocessing\_base.py", line 344, in listFields
    self._gp.ListFields(*gp_fixargs(args, True)))
IOError: "" does not exist

Failed to execute (TestPYTListFields).

If I replace fieldList = arcpy.ListFields(myLayer) with fieldList = arcpy.Describe(myLayer).fields it works as expected.

Why does my ListFields() fail from a Python Toolbox, and what can I put in place to make it work?

I suspect it may be something to do with the Feature Layer object, as I'm also having Search Cursors fail when referencing the Feature layer objects as well, although I haven't tested that to the same extent as the ListFields.

Below is my entire PYT:

import arcpy

class Toolbox(object):
    def __init__(self):
        self.label = u'ListFields'
        self.alias = ''
        self.tools = [TestPYTListFields]

class TestPYTListFields(object):
    def __init__(self):
        self.label = u'TestPYTListFields'
        self.description = u''
        self.canRunInBackground = False

    def getParameterInfo(self):
        param_1 = arcpy.Parameter()
        param_1.name = u'Input_Layer'
        param_1.displayName = u'Input Layer'
        param_1.parameterType = 'Required'
        param_1.direction = 'Input'
        param_1.datatype = u'Feature Layer'
        return [param_1]

    def isLicensed(self):
        return True

    def updateParameters(self, parameters):
        return

    def updateMessages(self, parameters):
        return

    def execute(self, parameters, messages):
        import arcpy
        myLayer = parameters[0]

        fieldList = arcpy.ListFields(myLayer)         # ListFields doesn't work from PYT
        #fieldList = arcpy.Describe(myLayer).fields   # Describe().fields does work from PYT

        for field in fieldList:
            messages.AddMessage(field.name)
        return


texts: How to load data programmatically from shapefile into feature class (geodatabase)?

texts: Can some Genuis here help me programmatically (using python and comtypes) load shapefile data in to an existing geodatabase (feature class). The attributes of the shapefile has to match with the feature class fields. For instance i have field "Height" in shapefile which has to go into "Height" field in the feature dataset.


texts:
Postgresql11 - postgis2.4
Here is a table mytbl with points and an attribute var :

geom | var
pt      a
pt      b
pt1     a
pt1     a
pt1     b
pt2     c
...

I want to "rotate" the table like this :

geom count(a) count(b) count(c)
pt      1       1
pt1     2       1
pt2                      1
...

For now I have with...

select geom,
 (case when var like 'a' then count(a) end) ,
 (case when var like 'b' then count(b) end) ,
 (case when var like 'c' then count(c) end)
from mytbl
group by geom, var;

...the result :

geom|count(a)|count(b)|count(c)
pt      1
pt                  1
pt1     2
pt1                 1
pt2                         1
....

How can I "group again" the features with the same geom to get just one line per geom with all count ?


texts: Using Legend properties in ArcGIS Desktop?

texts:
I am trying to input acress right next to values in the Legend; however, I was not able to . See here

Even though I know how to do the Legend Properties. See here

Can you point me where can I point to make them move wider ?


texts:
I need to do integrated viewshed analysis of a target with several observer locations. I am going to use r.viewshed analysis function in GRASS. With this function, it is possible to perform analysis based on single observer location. My aim is to add output of each single observer location into integrated viewshed output. Each output, i am going to flag with -b where o denotes non-visible while 1 denotes visible. Through a loop of observer locations, output of viewshed of each location can be integrated with r.mapcalc function. After integration, values = 0 in the rater map dentoes non-visible areas and value >= 1 denotes visible areas. Python code in GRASS for r.viewshed is ok. Under the loop of observer locations how to use r.mapcalc function with each viewshed adding? After getting the integrated viewshed raster, i need to extract its boundary depicting visible area in vector format. How to do that in GRASS? Any idea for this integrated analysis is highly appreciated.

If the integrated output is intoutput, can i use a code like this.
in the loop if i=0 -> r.mapcalc 'intoutput = viewshedoutput'
if i > 1 -> r.mapcalc 'intoutput = viewshedoutput  + intoutput'


texts:
Is there a way to add metadata to a shapefile in QGIS?

I have the shapefile and an XML metadata file and I want to import in shapefile.


texts: How to use Planet Lab tiles in Jupyter Notebook?

texts:
I just want to use Planet tiles as a basemap for a Jupyter Notebook application, using ipyleaflet widget.

In my environment, I imported these packages:

from ipyleaflet import Map, basemaps, basemap_to_tiles, TileLayer

Now, I am trying to initialize a map and streaming Planet tiles using my API.

planet_tiles_url=("https://tiles{0}.planet.com/basemaps/v1/planet-tiles/global_monthly_2016_05_mosaic/gmap/{x}/{y}/{z}.png?api_key=abc123")
tile_layer=TileLayer(url=planet_tiles_url)

m = Map(
    default_tiles=tile_layer,
    center=center,
    zoom=14,
    scroll_wheel_zoom=True)

m.add_layer(tile_layer)
m.add_control(LayersControl())
m

No error pops up. But no Planet basemap is displayed and not even the layer controls. See image..

What am I doing wrong?


texts: Calculate dynamic radius based on population and map resolution

texts:
Background

Create a radius around the centre of a city that corresponds to the city's population, taking into consideration the map's current resolution. The populations range from 39 to 1090936, with an average population of 2127. (Most of the populations are an estimated 150 people: hamlets.)

Google Fu returns people looking to calculate the population within a given radius -- this question is the exact opposite: calculate a radius for a given population.

City Data

The following code calculates a radius for a city based on the city's population:

function logn( base, n ) {
  return Math.log( n ) / Math.log( base );
}

function addCity( cities, lon, lat, population ) {
  cities.push(
    new OpenLayers.Feature.Vector(
      getPoint( lon, lat ),
      {
        tooltip: population,
        radius: Math.round( logn( 2, population ) ) * 2
      }
    )
  );
}

function getCities() {
  var cities = [];

  addCity( cities, -116.383062683871, 56.4680558322581, 39 );
  addCity( cities, -111.07205136129, 55.6315222967742, 208 );
  addCity( cities, -116.019763045161, 58.3878539741935, 739 );
  addCity( cities, -113.330203110508, 53.823233143442, 2848 );
  addCity( cities, -113.400259352921, 49.7204543294903, 3072 );
  addCity( cities, -114.088890487523, 52.3022833957834, 11115 );
  addCity( cities, -113.538803745389, 53.2670385040641, 24139 );
  addCity( cities, -113.797128210771, 52.281512598323, 65733 );
  addCity( cities, -114.052214116579, 51.0350223014551, 1090936 );

  return cities;
}

Radius Style

The radius is styled as follows:

var cityStyle = new OpenLayers.Style(
  {
    fillColor: '#0065BD',
    fillOpacity: 0.1,
    strokeColor: '#0065BD',
    strokeOpacity: 0.25,
    strokeWidth: 1.5,
    pointRadius: '${radius}'
  }
);

var cityStyleMap = new OpenLayers.StyleMap( cityStyle );

cityLayer = new OpenLayers.Layer.Vector( 'Cities', {
  styleMap: cityStyleMap
});

Problem

The radius is fixed, regardless of zoom level. This leads to incorrectly sized circles, depending on the zoom (i.e., map.getResolution()).

Zoomed In

Zoomed Out

Questions

Is log2(n) a good idea to calculate the relative radii (as opposed to a linear calculation, such as simple division)? (This is not overly important.)
How would you adjust the radius based on the map's resolution?

Note that the circumscription need not align with the city's centroid; a simple circle with a radius that is relative to the population and zoom level, centred on the city's location is sufficient.

Example - Zoomed Out

Example - Zoomed In

Thank you!

Related Links

https://stackoverflow.com/questions/10080358/how-can-i-calculate-dynamic-default-radius-for-map-search-based-on-population-de
Resizing and Reshaping features in Openlayers
http://users8.jabry.com/captdragon/ol/modify1.html


texts: How to create a multiple buffer using geopandas

texts:
I want to create buffers from any point in a multiple point file.

I can do it with a constant buffer using geopandas.buffer()

track.geometry = track['geometry'].buffer(0.05)
track_line = track.iloc[2:].head()

but the buffer is applied to all points.

I would like to apply different buffer values. For example, using values in another column of the data frame or a new csv column value.

For example,

 ID,latitude,longitude,buffer_value
 0,45.4,-156.0,0.3
 1,45.53,-156.98,0.4
 2,45.66,-156.97,0.34
 3,45.79,-156.957,0.23
 4,45.93,-156.942,0.45

How can I do it?


texts: Extracting tree crown areas from remote sensing data (visual images and LiDAR)

texts:
I am looking for a method to process a remote sensing image and extract the crown areas of the individual trees from the image.

I have both visual wavelength areal imagery, and lidar data from the area. The location in question is a desert area, so the tree cover isn't as dense as a forest area. The resolution of the aerial imagery is 0.5 feet by 0.5 feet. The lidar resolution is approximately 1 x 1 feet. Both the visual data and the lidar come from a Pima County, Arizona dataset. A sample of the type of aerial imagery I have is at the end of this post.

This question Single Tree detection in ArcMap? seems to be the same issue, but there does not seem to be a good answer there.

I can obtain a reasonable classification of the vegetation types (and information about the overall percent cover) in the area by using the Iso Cluster classification in Arcmap, but this provides little information on individual trees. The closest I have to what I want is the results of passing the output of the isocluster classification through the Raster to Polygon feature in Arcmap. The problem is that this method merges near by trees into a single polygon.

Edit: I probably should have included some more detail about what I have. The raw datasets I have are:

Full las data, and a tiff raster generated from it.
Visual imagery (like the sample image shown, but covering a much wider area)
Manual direct measurements of a subset of the trees in the area.

From these I have generated:

The ground/vegetation classifications.
The DEM/DSM rasters.


texts:
I have Apache Tomcat 8.0 and GeoServer 2.6 on Windows 7.

I have placed GeoExplorer folder into webbapp (The files web.xml and weblogic are included) as described, but it does not work!

Apache Tomcat works properly, GeoServer too, but GeoExplorer does not.

Please could you help me?


texts: Converting GeoDataFrame to DXF using GeoPandas

texts:
I have ten geodataframe;
I would like to export each one of them to DXF:
geodataframes= [gdf1,gfd2,gdf2]
for geodataframe in geodataframes:
    #geodataframe to dxf

it would be possible?


texts: Labeling features and converting them to annotations with ArcPy?

texts:
how could I do the following steps with arcpy:

right-click on the layer with the features
Label Features
Convert Labels to Annotations ...

The annotions should be saved in a new shapefile.

How could I do this in python?


texts: German Census Demographic Data @ zip-code level

texts:
The US Census has an API to pull Census/ACS data at the ZTCA (zip-code tabulation area) level. I'm looking for the same functionality with German data. So far I've identified several websites such as https://www-genesis.destatis.de, https://www.destatis.de/EN/Homepage.html and http://www.forschungsdatenzentrum.de/ for German demographic data, but they seem to go down to the "Gemeinde" and "Verbandsgemeinden" level, not the "Postleitzahl" level.

Also, data seems packaged in Tables, rather than extractable via APIs. I'm still going through the documentation.

Any pointers?

(BTW, http://www.programmableweb.com/api/genesis-online points to the urls above)


texts: Merge/Dissolve Polygon overlapping within

texts: I have a polygon that has been broken down to an intersecting analysis point. Each segment has its own acreage within an acre field on the attribute table. I want the intersecting lines to be removed and the polygon to be one whole with the acreage of each small segment all summed up together. I have tried dissolving by the acre field, and unchecking multi-part lines but the lines and small segments remain.


texts: Any way to make QGIS display band names of a stack properly?

texts:
I'm working with several layer stacks that have 100+ bands in QGIS.
The stacks are ENVI files, all band names are properly stored in the headers. Unfortunately, QGIS will not diplay them when choosing the style or when extracting info by clicking on the map. However, QGIS will show all the names Properties->Metadata.

Isn't there any way to get the band names to show properly?


texts: Disabling the editing option in arcgis.com pop-ups?

texts:
I am the owner/administrator of an arcgis.com web map containing feature points with pop-ups.  How do I disable the editing option from the pop-up window?  (see image below)

I unchecked attribute editing in the Configure pane (below), but no success.

I couldn't find an answer in the ArcGIS Online help pages...

Edit:

I finally found ESRI help at https://doc.arcgis.com/en/arcgis-online/create-maps/disable-editing.htm

Here are the pertinent instructions from that web page:

Verify that you are signed in and open the map in Map Viewer.
Click Details and click Contents.
Browse to the editable feature layer, click the More Options button More Options, and click Disable Editing.

But when I follow these steps, there's no Disable Editing option:

I'm still unable to disable the editing option...

EDIT:

Following Colin Lang's kind input (below)...

Apparently I don't have any Feature Layers.  What I have are Web Maps containing layers that I downloaded from my local computer (as zipped shapefiles) to the web map I was creating in Arcgis.com, as shown here:

When I filter for Feature Layers, as Colin suggests, nothing is returned:

In summary, I'm still stuck with the original problem: How to disable the edit option in my web maps and web applications?


texts:
I'm am trying to use a searchcursor to identify a list of ID's from a dataset. Using the list, I am now trying to delete those rows with the same ID's in another dataset.

I'm not sure where I'm going wrong.

rtc_id_values = []
with arcpy.da.SearchCursor(rtc_shp, "ID") as scursor:
    for srow in scursor:
        rtc_id_values.append(srow)

with arcpy.da.UpdateCursor(corporate_rtc_dataset, "ID") as ucursor:
    for urow in ucursor:
        if urow in rtc_id_values:
            row.deleteRow()


texts:
I am using catalog_apply() to normalize las points using lasnormalize() and segment individual trees using lastrees(). This is modified from the catalog_apply() lidR documentation (p.10):

library(lidR)

my_tree_detection_method <- function(cluster, ws)
{
  las <- readLAS(cluster)
  if (is.empty(las)) return(NULL)
  las_n = lasnormalize(las, tin())
  ttops <- lastrees(las_n, li2012(R = 3, speed_up = 5))
  # ttops is a SpatialPointsDataFrame that contains the tree tops in our region of interest
  # plus the trees tops in the buffered area. We need to remove the buffer otherwise we will get
  # some trees more than once.
  # bbox <- raster::extent(cluster)
  # ttops <- raster::crop(ttops, bbox)
  return(ttops)
}

ws <- "/path/to/my/las/data"
ctg <- readLAScatalog(ws)
lidR:::catalog_laxindex(ctg) # Build index

opt_chunk_buffer(ctg) <- 10
opt_chunk_size(ctg) <- 1000

opt <- list(need_buffer = TRUE) # catalog_apply will throw an error if buffer = 0
output <- catalog_apply(ctg, my_tree_detection_method, ws = 5, .options = opt)

The example uses the following to remove the processing buffer:

# bbox <- raster::extent(cluster)
# ttops <- raster::crop(ttops, bbox)

However, my modified example using lastrees() is not a SpatialPointsDataFrame and, I believe, cannot utilize extent() and crop() to remove the processing buffer.

What is the preferred way to remove the buffer used for processing in catalog_apply()?


texts: Exporting a group of raster layer bounds from QGIS to a GeoJSON

texts:
I have a QGIS project with a lot of raster layers; and I want to build a Mapbox map based on this project; that would display those images on a website's map.

I need a way to export a group of raster layers from QGIS to a GeoJSON collection of polygons; that I could load in Mapbox.
I already had a look on GeoTIFF files and PGW files; but still, exporting every layer one by one will be painful - as I'll have to update the web project from QGIS quite often.
So I'm looking for a solution that would automate the job.
I would like to get a GeoJSON file that would look like this:
{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": {
        "layer_name": "AND17",
        "image": "images/png/world.png",
      },
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [
              4.3813338557822,
              50.842573058675
            ],
            [
              4.418685805334,
              50.842573058675
            ],
            [
              4.418685805334,
              50.825893099756
            ],
            [
              4.3813338557822,
              50.825893099756
            ],
            [
              4.3813338557822,
              50.842573058675
            ]
          ]
        ]
      }
    }
  }

As you can see, it contains the georeferencing informations of the image, and its filename.
If I can get that, i'm able to load and place the image on Mapbox.
Have you got any ideas of how I could do that; maybe using a plugin?
I've also seen that you can export a group to a *.qlr file.
It seems that it contains some kind of geospatial datas :
<extent>
    <xmin>475442.96307834662729874</xmin>
    <ymin>6588302.91859502345323563</ymin>
    <xmax>483907.46308698022039607</xmax>
    <ymax>6594287.91860113199800253</ymax>
</extent>

Any tips ?


texts: How to specify file format of Split Vector Layer in QGIS

texts:
I'm attempting to use the Split Vector Layer tool in QGIS (see screenshot below) in order to separate a larger shapefile into its constituent components by object ID. The tool accomplishes this perfectly, but by default exports everything as a geopackage file.

I instead require these files as shapefiles, not geopackage files. Is there any way to specify within QGIS that I want to export these files as such?


texts: Changing line width by classification in QGIS

texts:
I am using QGIS 2.18.18

I have a map of rivers and I want to style them based on their width. I have width as an attribute of the rivers layer.

Just as we can style colors by category or graduated styling, I want to do the same with width.
Currently I use an expression editor for data driven override (circled in the picture) after which I have to manually classify everything.
If classification could be done by  quantiles, jenks or other breaks automatically, then it would be easier for me to see which classification has what effect on my map. As I said, that is possible with color, but I cannot figure out how to do it with width

If there is a way to do it without code that would be great, as I haven't learnt that yet. If there isn't, then I'll implement it once I learn code.


texts:
I have a shapefile called 'Roads' with on the fly CRS which is ok but not for what I want to achieve. I already know the correct CRS value I want to essential duplicate the Road shapefile with the correct CRS and then continue editing the new shapefile. I have currently duplicated the road shapefile changed the CRS value and then tried to re-save it but once I do that and then try to open it again QGIS 3 comes up with an error message

I need a step by step guide as I am clearly missing something.


texts: Invalid Polygons from ST_MakeSolid(ST_ConvexHull(...))

texts:
The aim: a TINZ has been derived from at least two lidar scans with the results stitched together into a single mesh. The original scans are unavailable, and the stitched mesh is rejected by sfcgal functions as having an invalid structure.

I'm exploring the approach of reducing the mesh to points again and trying to reconstruct the mesh using postGis functions. The following test example uses a simple 3D cube as the POLYHEDRALSURFACEZ source:

WITH getpoints AS
(SELECT (ST_DumpPoints(ST_MakeSolid('POLYHEDRALSURFACEZ(
                            ((0 0 0, 0 0 1, 0 1 1, 0 1 0, 0 0 0)),
                            ((0 0 0, 0 1 0, 1 1 0, 1 0 0, 0 0 0)),
                            ((0 0 0, 1 0 0, 1 0 1, 0 0 1, 0 0 0)),
                            ((1 1 0, 1 1 1, 1 0 1, 1 0 0, 1 1 0)),
                            ((0 1 0, 0 1 1, 1 1 1, 1 1 0, 0 1 0)),
                            ((0 0 1, 1 0 1, 1 1 1, 0 1 1, 0 0 1)) )'))).geom)

SELECT ST_AsText(ST_Tesselate(ST_MakeSolid(ST_ConvexHull(ST_Collect(getpoints.geom)))))

FROM getpoints;

However, this fails with the error:

ERROR: Polygon is invalid : points don't lie in the same plane : POLYGON((0/1 0/1 0/1,0/1 1/1 1/1,1/1 1/1 0/1,1/1 0/1 0/1,0/1 0/1 0/1)) SQL state: XX000

I've checked the ST_ConvexHull and ST_MakeSolid functions using ST_IsValid and they both return true. Have also tried ST_MakeValid on ST_ConvexHull and ST_MakeSolid, but this makes no difference.

If I tesselate the test cube directly like this, it works fine:

SELECT ST_AsText(ST_Tesselate('POLYHEDRALSURFACEZ(
                            ((0 0 0, 0 0 1, 0 1 1, 0 1 0, 0 0 0)),
                            ((0 0 0, 0 1 0, 1 1 0, 1 0 0, 0 0 0)),
                            ((0 0 0, 1 0 0, 1 0 1, 0 0 1, 0 0 0)),
                            ((1 1 0, 1 1 1, 1 0 1, 1 0 0, 1 1 0)),
                            ((0 1 0, 0 1 1, 1 1 1, 1 1 0, 0 1 0)),
                            ((0 0 1, 1 0 1, 1 1 1, 0 1 1, 0 0 1)) )'));

So why are the outputs of the other functions not valid for ST_Tesselate???


texts: How to manage different SRID data in same database table?

texts:
My data is in two different cities and three different EPSG code area. (Database is PostgreSQL/PostGIS.)

So I think, I can not set SRID values by data. Making 3 different table is difficult. How can I do it?
Data will  be collected with high accuracy GPS devices. It is also difficult to change coordinates in city borders.

If I use global EPSG:3857 for all data, the calculation of length will be different from actual collected data data. How can I solve this problem?


texts: How to update a DATE field when moving PostGIS geometry in QGIS?

texts:
I've been correcting our address point file using various methods to identify points that need to be corrected, with the final step being a manual move the point (or by the Numerical Vertex Edit) to its corrected location using QGIS.

What I'd like to do is update the MODIFIED_DATE field automatically when the geometry is moved.

The address points live in PostGIS.

I can see that by setting the Field Edit Widget type to DATE/TIME, and then identifying the point and enabling the 'Auto open form' option, the date will be updated.

But can the date value be updated when the geometry is moved, by some type of dynamic field, or other method?


texts: Spectral Signature Extraction using landsat 7 and MODIS

texts: I am working on differentiating cocoa plants from other vegetation types using 4 epoch bwt 1980-2010. I need to extract the spectral signature for cocoa plants and other known features. I plan using landsat or MODIS imagery. How can I do this using either ArcGIS, Erdas, Envi or Idrisi?


texts: How to calculate raster statistics for polygons?

texts:
how can I calculate sums, averages etc of raster-points (multi-band) per polygon of a vector-layer. I was told that this is called "zonal statistics". I tried that with QGIS first.

There is a way to do it but that is much too slow (convert raster to vector, intersect with second vector-layer, calculate geometry, export numbers, calculate statistics with spreadsheet or other program, re-import the results, takes forever for my 350.000 raster-points).

I was also given the hint to use saga-gis. That has "zonal statistics" but those are based on categories from a raster-band, not on polygons from a vector-layer. So to use this I would have to convert my vector-layer to raster and then calculate the statistics.

This seems to be the wrong way to solve this. There would be no way to account for raster-points that belong to 2 or more polygons because they are intersected by the polygon-boundary. I assume that polygon-based statistics should be able to handle this so I also assume that I haven't found the correct module yet.

Saga-gis has really many modules. Please let me know which one is the right one for this application.


texts:
I just started working with ArcMap and I have a question:

How can I make the add in I made in C# to open a map automatically in ArcMap? Just want to know what function to use or if it's even possible.


texts: How can several vector layers have different controls, using Openlayers?

texts:
Using openlayers I have two vector layers. One layer is called lineLayer consisting of vector lines, and the other is called pointLayer consisting of points. When I hover on the lines I want to be able to highlight the line, and identify the vector object when I click on it. When I hover on the point I want to highlight the point and be able to drag it. My code works separately for lines and points, but not when I put it together.

This is how I manage the line layer:

// handle interactions with lines
var highlightCtrl = new OpenLayers.Control.SelectFeature(lineLayer, {
    hover: true,
    highlightOnly: true,
    renderIntent: "temporary"
});

map.addControl(highlightCtrl);
highlightCtrl.activate();

var selectCtrl = new OpenLayers.Control.SelectFeature(lineLayer, {
    clickout: true,
    onSelect: function(feature){console.log(feature)}
});

map.addControl(selectCtrl);
selectCtrl.activate();

And this is how I manage the points:

// handle interactions with points
var highlightpoint = new OpenLayers.Control.SelectFeature(pointLayer, {
    hover: true,
    highlightOnly: true,
    renderIntent: "temporary",
});

map.addControl(highlightpoint);
highlightpoint.activate();

var dragStart = new OpenLayers.Control.DragFeature(pointLayer, {
    autoActivate: true
});

map.addControl(dragStart);
dragStart.activate();

From reading other examples I have noticed that the order seems to matter when a control is added. However, I've tried all the possible combinations of layer adding, and it does not seem to make a difference.

Is there some way to do this that I am overlooking?


texts: I'm using ArcMap 10.1. When I make an attribute based buffer, the result had shown as overlapped. How can I solve this issues?


texts: How to pull locations/data from an embedded Google map on another website?

texts: Is there a way to download the point locations from an embedded Google map? I want to compare the locations of these businesses to a rail network.


texts: QGIS 3.4.4 Attribute table expression not working anymore

texts:
After uptdating Qgis to 3.4.4 the following expression doesn't working anymore:

"art" LIKE '%'  || attribute ( @atlas_feature , 'artere')  || '%'

We're using it in the printcomposer, specifically on an attribute table.
It was used to display lines with art value containing the same value as the layer we're using to generate the atlas, which always has the column artere.

On previous version, it was working. Not anymore.

Did we do something wrong?


texts: Math related sum distance problem with TSP content

texts:
I hava an interesting problem related with GIS / MATH / TSP.

If I calculate the distance between A and B on a x/y-style coordinate system it is obvious to use the following formula:

sqrt((Ax-Bx)^2+(Ay-By)^2) = Distance

If I compare now 2 distances AB and CD it's obvious to do this with the formula on top:

sqrt((Ax-Bx)^2+(Ay-By)^2) = Distance
AND
sqrt((Cx-Dx)^2+(Cy-Dy)^2) = Distance

Now if you need to do a lot of distance calculation, each SQRT() operation eats up valuable CPU time. And as we know 2<3 and 2^2<3^2 you can leave the sqrt() command away - just to enhance performance! SQRT is very expensive!

And now where it starts to get complicated for me. WHY doesn't this work with ABCDEF paths?
Here is a classical TSP 3-opt optimization on these points.(example ABCDEF, ABCDFE, ABCFED, etc. permutations...). It's obvious that:

sqrt((Ax-Bx)^2+(Ay-By)^2) + sqrt((Cx-Dx)^2+(Cy-Dy)^2)

is not the same

(Ax-Bx)^2+(Ay-By)^2 + (Cx-Dx)^2+(Cy-Dy)^2

is obvious, but it should be proportional, or is this a wrong assumption?

I have the following log which shows that my theory is wrong!

Distanz without SQRT(): 2423932.0 - Distanz SQRT(): 24131.530843013228
Distanz without SQRT(): 2421100.0 - Distanz SQRT(): 24124.90131901531
Distanz without SQRT(): 2419674.0 - Distanz SQRT(): 24126.12515189262
Distanz without SQRT(): 2414714.0 - Distanz SQRT(): 24117.38148735822
Distanz without SQRT(): 2414304.0 - Distanz SQRT(): 24112.91523442742
Distanz without SQRT(): 2409924.0 - Distanz SQRT(): 24117.550376378276
Distanz without SQRT(): 2403676.0 - Distanz SQRT(): 24114.96478009344


texts:
I am using GeoServer with OpenLayers.

I would like to return all the feature type details in a specific layer exactly as GetFeatureInfo does for a specified location. I tried using GetFeatureInfo with the x and y parameters but it insists on a location being specified. Basically, I need the details of all the feature types in the layer, including the geometry.  There may be a better way to do this, however, after hours of searching, I have been unable to find one.

Can anyone help?


texts: How to cut a raster into smaller, equal area parts?

texts:
I have a large raster file that I need to cut into smaller parts.
I already have a polygon shapefile that has the smaller parts that I want to use as an index. The reason for this is that end users will only have to load the tiles they need. The files will be stored on disk.

I would like to cut the raster and have the output files have a suffix of the index. i.e. R1C2

Qgis 2.8 on W7


texts: Creating a code that computes the distance between two properties given northings and eastings

texts:
I was given a file with different properties in different sheets containing coordinates for each corners. I need to create a program that computes the distance between two lots. I'm not sure how to approach making it, do I just need to compute for the centroid of each property and then use the distance formula? Or do I need to use each point corners?
This is my code so far.
import pandas as pd
class Lot:
    def __init__(self,name,vertices):
        self.name = name #name of the lot (in the sheets it's A1, A2 ..)
        self.vertices = vertices #list of tuples (eastings,northings)
        self.centroid = self.__getCentroid()

    def __getCentroid(self):
        '''A method that computes the centroid of the lot and returns a tuple containing the Eastings and Northings of the centroid.
        Used to determine the centroid attribute.'''
        self.e = 0
        self.n = 0
        count = len(self.vertices)
        for e in self.vertices:
            self.e += e[0]
        for n in self.vertices:
            self.n += n[1]
        centroid = (self.e/count , self.n/count)
        return centroid

    def distanceTo(self,other_Lot): #float ouput
       '''A method that computes the distance between the lot and the lot given'''
       pass

    #def nearest(lotList:List):
    #    '''Determines the nearest lot from a list of lots'''


texts: Labeling shapefile address numbers in AutoCAD?

texts:
I am trying to figure out how to label 'ADDR' field in autocad only showing the address number and not the street.

I try using LTRIM and RTRIM expressions and nothing changes the way it shows

How it shows : 123 main st

how i want it: 123


texts: Unable to load datasets from datacube using x,y, lat,lon or geopolygon

texts:
I get an empty array if I load the Landsat 8 dataset using x,y, lat,lon or geopolygon.
ds = dc.load(product='ls8_usgs_sr_scene',
          time="2014-01-26",
          geopolygon=geom,
          output_crs="EPSG:4326",
          resolution=(-0.00027,0.00027))

However, I'm able to load the dataset if those parameters are omitted, but that means I had to load the entire scene.
Is it something missing from the product definition or the metadata?


texts: Discontinued Sequence with pgr_ksp() result

texts:
I'm having some troubles with the result I got with pgr_ksp, which returns me a result with missing data.

Here is a preview of the answer, where you can see the sequence gets discontinued, and data missing:

Sequence goes from 3 to 11, and also skipping 13 and 15.

Here is the example of the query:

SELECT seq, path_id AS route, path_seq, node, edge,
    cost as duration,
    e.distance as distance,
    e.route as route,
    ST_AsGeoJSON(e.the_geom)::json->'coordinates' as polyline,
    agg_cost
FROM pgr_ksp('SELECT e.id AS id,
        source::integer,
        target::integer,
        duration::double precision AS cost,
        reverse_cost::double precision as reverse_cost
    FROM edges e
    GROUP BY e.id, source, target, distance, cost, reverse_cost
    ORDER BY e.id',
    2967, 17670, 6, true, false)  AS ksp
LEFT JOIN edges AS e ON e.id = ksp.edge

Would it be because the edge table got too many data?


texts: 2 legend objects with arcgis javascript api

texts:
Can i create 2 legend objects in one script?

I'm currently working on a web-application in which i already do have a legend/TOC running and it looks like this:

var legendLayers=[];
    legendLayers.push({layer:raum, title:'Raumkategorien', slider:true});
        map.on("layers-add-result", function(){
        var layerInfo = new TOC({
                map: map,
                layerInfos: legendLayers,
        }, "legendDiv");
        layerInfo.startup();
    });
    map.addLayers([raum]);

Now I added some geoprocessing functionality and i wanted to add a new TOC when the result map is loaded in the web-app.
I tried this by doing the following:

function gpJobComplete(jobinfo){
        var mapurl = mapServiceURL + "/" + jobinfo.jobId;
        var energy = new ArcGISDynamicMapServiceLayer(mapurl,{
            id:"Beste Energieformen",
            "opacity": 0.85
        });
        map.addLayers([energy]);

        var outputLayers = [];
        outputLayers.push({layer:energy, title:"Beste Energieformen"});
        map.on("layers-add-result", function(){
            domUtils.show(dom.byId('outputDiv'));
            var outputInfo = new TOC({
                map: map,
                layerInfos: outputLayers
            }, "outputDiv");
        outputInfo.startup();
        });
    }

I get a result map and it works just nice in the web app, but there's no new TOC.
Am i trying it the wrong way?


texts: Getting incorrect slope values from an ASTER DEM in QGIS?

texts:
I have a DEM layer from ASTER satellite and I want to get the slope information from it.

I use terrain analysis > slope, but the output is so wrong. Why is that?

Here is the DEM:

And here is the "Slope" that is totally wrong. It should stretch the values from 0 - 89.9097 but it does not.


texts: I need to get the start and the end point of each line in a shapefile. In my case these lines represents pipelines. I need to get the start and the end points of each pipeline in an Excel sheet and to refer to the number of the line that is connected between any two points.


texts:
I am trying to extract the dates and mean NDVI values of a processed Google Earth Engine image collection thru Python. All goes well until I try to use getInfo(), where I always encounter this error:

  ee.ee_exception.EEException: String.slice, argument 'string': Invalid type. Expected: String. Actual: Float.

But when I run the same code within GEE, it doesn't generate that error. Here is a sample code:

// Initialize //
var roi = ee.Geometry.Polygon(
        [[[122.10367346625048, 7.389991097840974],
          [122.10670974995855, 7.389948538640855],
          [122.10667756321027, 7.392938312477742],
          [122.1036841951634, 7.392938312477729]]]);
var collection = ee.ImageCollection("LANDSAT/LE07/C01/T1_SR")
  .filterDate('1999-01-01', '2018-12-31')
  .filterBounds(roi)
  .select(['B3','B4','pixel_qa']);
print('Collection size (Initial): ', collection.size());
var cfilt = 0; // Threshold for cloud cover %

// Cloud Filtering //
var ctROI = ee.Number(collection.first().select('B3').reduceRegion({ // # of ROI pixels
  reducer:ee.Reducer.count(),
  geometry: roi
  }).get('B3'));

var cmaskL457sr = function(image){
    var qa = image.select('pixel_qa');
    var cloud = qa.bitwiseAnd(1 << 5).and(qa.bitwiseAnd(1 << 7)).or(qa.bitwiseAnd(1 << 3)).or(qa.bitwiseAnd(1 << 4));
    var mask2 = image.mask().reduce(ee.Reducer.min());
    return image.updateMask(cloud.not()).updateMask(mask2);
};

var cloudScore = function(image){
    var ctCov = image.select('pixel_qa').reduceRegion({
      reducer: ee.Reducer.count(),
      geometry: roi
    });
    var ctpct = ee.Number(100).subtract(ee.Number(ctCov.get('pixel_qa')).divide(ctROI).multiply(100));
    return image.set('cloud_cover_percentage', ctpct);
};

var cMasked = collection.map(cmaskL457sr); // Apply cloud mask
var dcFiltered = cMasked.map(cloudScore).filter(ee.Filter.lte('cloud_cover_percentage', cfilt)); // Apply cloud scoring and filter to threshold
print('Collection size (After cloud filtering): ',dcFiltered.size());

// NDVI Calc //
var addNDVI = function(image){
    var ndvi = image.normalizedDifference(['B4', 'B3']);
    var meanNDVI = ndvi.reduceRegion(ee.Reducer.mean(), roi);
    return image.addBands(ndvi).set(meanNDVI); // Attach mean NDVI values
};

var final = dcFiltered.map(addNDVI).select('nd'); // Final image collection

// Fetch Data //
var dates_list = ee.List(final.aggregate_array('system:time_start')).getInfo();
var ndmeans_list = ee.List(final.aggregate_array('nd')).getInfo();

Any ideas on why the getInfo() fails? Why does the error say that it expects a string, when getInfo() is used to fetch different kinds of data. Also, it's weird to me that the exact code for GEE works, but not when I use the Python version (Note that the Python code works fine except for that last getInfo() part). Is there a special casting to be done when working with Python GEE code? I have tried various castings but none have worked so far.

Lastly, in the future if I want to generate a seriesbyRegion for a massive (50K+) FeatureCollection, what would be a better way to do that than looping my current code across all regions inside?


texts:
How to make a spatial join with ArcGIS with these HUGE datasets?

Target: almost 3 million squares 500x500m size (the attribute table is basically empty)
Join feature: 12000 polylines
Operation: one to many
Match: intersect

All my attempts have failed:

From the Spatial Join tool in the Toolbox: after 15 min the output is a file where no join has been made (Join_Count field is 0 in all rows)
From the dialog box: It seems to work but my estimate is that it takes 4 days to complete - without any guarantee the output will be correct.

It works when I do a Spatial Join of a small area.

Any ideas on how to do this in a proper way?

I am using file geodatabase, ArcGIS 10.2, background processing, 32 GB RAM

(The idea is to overlap those polylines to the grid (the 500x500m squares). I will then make a dissolve to find out how many lines are in each square. Finally, I will make a density map.)


texts: Merge Ordnance Survey VectorMap District from different folders

texts:
I have four different folders with Ordnance Survey-Shapefiles from different regions.

I want to merge all these files theme by theme (for example all buildings to just one shapefile). The problem is that these are loads of shapefiles and all have different names like SS-Buildings.shp, ST_Buildings.shp and so on.

Has anyone else before managed to merge these files from different regions automatically (ogr2ogr / processing-toolbox?)?

Any ideas are very welcome.

What I found so far is this question:
Batch merging shapefiles from different folders, different subfolders

The problem with the OS-files is that the different shapefiles from one topic like "Buildings" don't have the same name:

(Here I put together all "buildings"-shapefiles in one folder manually )

So I have to find a way to merge the files ignoring the first three characters of the filename and merging them together from different subfolders... seems to be pretty challenging


texts: Finding list of all shapefiles that overlap with single shapefile using QGIS?

texts:
I have a list of 34 shapefiles. These shapefiles are range maps for different species. I would like to take 4 of these shapefiles and document a list of those shapefiles (within the 30 shapefiles), whose ranges overlap with these 4 shapefiles. I wouldn't mind comparing them one at a time ie. Shapefile 1 with the other 30 shapefiles; Shapefile2 with the other 30 shapefile..

Is there a way to do this in QGIS/PyQGIS?

I tried using the Intersect tool, but that would only give me the features that are common to both the files chosen.

I basically need a table that gives me : shapefile1 and all the names of the shapefiles that overlap with shapefile1

EDIT: Following @xunik's suggestion:

ISSUE I am having right now is that the script is not picking up every possible combination(even if they are slightly touching each other). Attached is an image that shows the script being written to a csv file.

EDIT 2

The script suggested by @xunilk below is not working and is actually showing incorrect information. Any suggestions?

See Image below. The output is showing that these two features intersect, when in reality, they dont. Please Help!


texts:
I'm having problems creating a heatmap. I want to show the cumulative increase of points through time manager with heatmap symbology in QGIS 2.18.1 (through layer symbology, not the raster heatmap plugin). The points are meant to display from dark to bright as point density increases. What seems to be happening is as points increase the relative brightness of points where density is low decreases. Ideally I'd like absolute cumulative density to display, not a changing relative density. How do I achieve this?

This behavior is also reflected when I try regular (non time-manager) heatmap. As I pan across the map brightness of areas is increased or decreased depending on the relative density of other areas displayed in the window.


texts: Web Application and ArcGIS server are two different machines and unable to load layers using rest services

texts:
ArcGIS 10.2 server is installed in one machine which doesn't have public ip. My web application is hosted in another machine and  IIS7 is configured and it has public ip. Both the machines are in LAN connection.

When we are calling the ArcGIS rest layer service from web application within LAN connection by using public/local ip, the Layers are loaded and no issues.

But when we are trying the same from outside network by using the public ip link of web application, ArcGIS rest layer services are not loaded. Because ArcGIS service url (http://xxx.xxx.x.xx:6080/arcgis/rest/) has local ip address. However the request is coming through web application which has public ip.

Is there any possibility to set the proxy to access ArcGIS server layers, or does the ArcGIS server also require public ip access?

This is my first web application with arcgis server.


texts:
I tried to add a border to my osm file to make is smaller.

The command line is as follows:

osmconvert china.osm -b=22.30,114.11,22.38,114.22 --out-osm -o=hk.osm

an error occurred like this:

osmconvert Error: use border format:  -b="x1,y1,x2,y2"

I then tried following the border format by adding the two colons but still same error


texts: Using GeoServer API how can I create multiple layers from same datastore?

texts:
The datastore is a connection to a solr core. I can create this just fine. I find when I create a new layer by posting a new featuretype xml, I am forced to name the file the same name as the data source (solr core). This only allows me to create one layer, with the same name as the solr core.

Is it possible to create multiple layers from the same datastore using the GeoServer API?


texts: Using delimited text layer to keep only those polygons and records of a shapefile which match a column in the delimited text layer?

texts: What I want to do is trim the shapefile so that only those records and associated shapes that match the same values in a column in the delimited text layer (which in my case is a CSV) are kept in a new shapefile file.


texts: I have a .dxf file to convert to .shp.  Using ArcMap I've added the .dxf and tried exporting to .shp (various types, polyline, point, polygon) and also tried exporting to a gdb.  For the most part it works but some of the original polyline (or polygon) files come into ArcMap as points.  Has anyone come across this problem before and know of a solution?


texts: Enable editing of only one feature in Polygon Layer using PyQGIS?

texts:
I was thinking of locking neighboring polygons (parcels) by disabling editing on them, while editing would be possible on the polygon at the center. This means other polygons can also be edited in case when the coordinate or border of that editable polygon moves.

Is that possible to do in PyQGIS?

I want this to lock editing on the other features.


texts: Delineating coastline with panchromatic imagery

texts:
I have a panchromatic satellite image with a large coastline area and my goal is to delineate this coastline with a maximum accuracy and minimum amount of manual work. What is the best option to accomplish this?

I do not have either DEM nor other spectral bands, so these solutions do not work:

Delineating large scale watershed?
Watershed delineation with meandering river?
Extracting coastline as vector from landsat 8 imagery?

Here is a small area of the original 2-meter satellite image I have:


texts: Creating random location within a network from a given location and a given distance

texts:
For a habitat selection study, I am trying to identify availability. My species lives along waterbodies and I could assess distances moved between two consecutive locations using the network analyst (Calculating distances between consecutive points along rivers (polygon) using ArcGIS Desktop?).

I would like to get a potential location the animal could have reached within the same distance.

For better understanding:
I have location A and the distance d from location A to B
I have a network of rivers and other waterbodies

I now want to get the location C using A as a starting location and distance d as the distance along the network.

Is that possible to do and if so, how?


texts: ArcPy - Dividing every raster cell value by the maximum cell value

texts:
I want to divide every cell of a raster by the raster's maximum cell value using ArcPy. Here is what I have done so far.

densityRaster = LineDensity(shapefileName, None, "", "", "")
maxVal = arcpy.GetRasterProperties_management(densityRaster ,"MAXIMUM")
normalRaster = densityRaster/maxVal #maxValue i get is 5.66

This gives me an error described below. Anybody has clues to where I get it wrong?

  Runtime error : ERROR 000732: Input
  Raster: Dataset 5.666 does not exist or is not supported

If I enter normalRaster = densityRaster/5.666, I get the correct result.


texts: I saved a .qgs file with home server postgis layers. Had to work remotely, but will lose unhandled layers upon opening the file. I can connect to these through my GeoServer, but I'm going to lose my styles and queries and have to literally rebuild my map from scratch. Is there a way to keep QGIS from dropping these unhandled layers so I can refer to them when rebuilding the map remotely?


texts: Problems opening vector layers after shifting the location of the saved layer

texts:
At work, we have recently changed from MapInfo to QGIS and our operating system from Windows to Ubuntu. We run 4 computers each with its own copy of QGIS, all connected to a server. The following are some of our problems:
1.  We cannot save straight from QGIS to the server. From our reading, we expect that this is caused by having separate copies of QGIS and that if we had one copy of QGIS on the server, we could save directly to the server.

Is this true please?

Although tedious, this problem is easy enough to circumvent: Save to the desktop and then shift it to the server.
However, after creating a new vector layer, I’ll save all the associated files and project and then carefully shift all the files to the server. Returning to the same project later one of three things happen with no consistency!!??

(i) the vector layer will opens perfectly.
(ii)    the layer ‘opens’ but the attribute table is empty and thus nothing exists in the layer.
(iii)   the vector layer will open but we are unable to edit it.
Any ideas please?


texts: Adding Raster Layer without lyr file using ArcPy?

texts:
i wonder how i can add a raster layer to a remote (not the current) mxd file with arcpy.
I need to not use a lyr file for that.

For vector Layer i have the following lines working fine :

mxdPath = "D:\FormationSpecialisee\IGN_IFI\doc1.mxd"
mxd = arcpy.mapping.MapDocument(mxdPath)
layer = arcpy.mapping.Layer(fcADD)
arcpy.mapping.AddLayer(df, layer , "AUTO_ARRANGE")

but it seems to not work for Raster Layer.