texts: Update/change/commit an attribute/a field value "on the fly" in QGIS 3

texts:
I'm using QGIS 3.2.
I have 4 fields ("kk1"-"kk4") where at least one of the fields contains a number (real).
I want to calculate the mean of these 1-4 fields and save the result (the average) in another field called "mkviskrans".
I calculate the average (based on this answer) of these 1-4 fields using this function in the Function Editor window in the Layer Properties:

@qgsfunction(args=-1, group='Custom')
def mean_col(cols, feature, parent):
  vals = [feature[i] for i in cols if feature[i] is not None]
  return sum(vals)/float(len(vals))

and calling the function from the Expression window in the widget's Default value option in Layer Properties:

mean_col('kk1','kk2','kk3','kk4')

UPDATE: I tried return 2.3instead of return sum(vals)/float(len(vals)). In this case the attribute table is updated automatically with the value 2.3. This indicates that the values typed in by the user in "kk1"-"kk4" actually have to be updated/evaluated/"commited" before returning the mean value, i.e.the division-expression will always return "None"? I have tried to update the feature/field values of "kk1"-"kk4" before return, but without any luck so far.

However, the mean value does not automatically appear in the "mkviskrans" "Text Edit" widget (alias "Midlere kvistkransavstand (dm)"):

Nor in the attribute table:

When I click in a cell in the attribute table, then the mean value appears:

Is it possible to get/show/save the value "on the fly"? I have tried different variants of update and/or commit within the function:

@qgsfunction(args=-1, group='Custom')
def mean_col(cols, feature, parent):
  layer = qgis.utils.iface.activeLayer()
  fieldIndex = layer.fields().indexFromName( 'mkviskrans' )
  vals = [feature[i] for i in cols if feature[i] is not None]
  """return sum(vals)/float(len(vals))"""
  mv = sum(vals)/float(len(vals))
  layer.changeAttributeValue(feature.id(), fieldIndex, mv)
  layer.commitChanges()

If I save the edits without clicking in a cell in the attribute table, the average (the field "mkvistkrans") will not contain any mean values.

(I may of course calculate the mean values after the field work, but I want the user to see the mean value "on the fly" - if possible.)


texts: Accessing personal assets if connected through service account to Google Earth Engine

texts:
I'm connecting to the Earth Engine API through a service account :
service_account = 'test-sepal-ui@sepal-ui.iam.gserviceaccount.com'
credentials = ee.ServiceAccountCredentials(service_account, 'test.json')
ee.Initialize(credentials)

but I cannot access my assets. when I try to execute :
ee.data.getAssetRoots()[0]['id']

That should give me the root folder name of my assets. Instead I get the following error :

IndexError: list index out of range

because the list is empty...
Is it a normal behaviour ? How can I access my assets ?


texts: GDAL rasterizelayer function 's additive rasterizing option (add=True) is not working, what am I doing wrong?

texts:
I am trying to rasterize a set of polygons with integer features. Since I am very new using raster data as an initial step I am rasterizing the the polygon centroids first.
In the output array I want to see the aggregated values of the integer features. I found out in the documentation that setting "add=True" should work. But it is not working.
link to the documentation is here https://gdal.org/python/osgeo.gdal-module.html#RasterizeOptions
You can access my data here: https://drive.google.com/drive/folders/1uzQ6eZodQNN_Qzi1-G85-BaPvElU8KeS?usp=sharing
Code:
import numpy as np
import os
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point, Polygon, MultiPolygon, MultiPoint
from descartes import PolygonPatch
from osgeo import gdal
import ogr
import logging
from rasterio.plot import show

fn_ras = 'Data/Vector_Data/sanfrancisco.tif'

fn_vec ="Data/Vector_Data/vec_cbg.shp"

### Read the base tif file -- from
ras_ds = gdal.Open(fn_ras)

geot = ras_ds.GetGeoTransform()

vec_ds = ogr.Open(fn_vec,0)

lyr = vec_ds.GetLayer()

featureCount = lyr.GetFeatureCount()

print ("Number of features in {}: {}".format(os.path.basename(fn_vec),featureCount))

lyr.GetFeature(10).GetField('BlockGroup')

### CREATE THE TARGET RASTER

pixel_size = 0.01
source_srs = lyr.GetSpatialRef()

x_min, x_max, y_min, y_max = (-122.63749999999999, -121.39250000000001, 36.902499999999996, 38.1375)

# Create the destination data source
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)

### Change the geotransform
## the GT(2) and GT(4) coefficients are zero,
## and the GT(1) is pixel width, and GT(5) is pixel height.
## The (GT(0),GT(3)) position is the top left corner of the top left pixel of the raster.
# Create a new geomatrix for the image

geot= list(geot)
geot[0] = x_min
geot[3] = y_max
geot[1] = pixel_size
geot[5] = - pixel_size

# Create a new raster for the rasterized polygons
drv_tiff = gdal.GetDriverByName("GTiff")

target_ds = drv_tiff.Create('sf.tiff', x_res, y_res, 1, gdal.gdalconst.GDT_Float32)

target_ds.SetGeoTransform(geot)

target_ds.SetProjection(ras_ds.GetProjection())

##WHY IS THIS PRINTING OUT 0 ???

# Rasterize
gdal.RasterizeLayer(target_ds, [1], lyr , options=["attribute=BlockGroup", 'allTouched=True',"add=True"])

# Read as array
band = target_ds.GetRasterBand(1)
##this returns none
array = band.ReadAsArray()

print(array.shape)
print(array.sum()) ## something WRONG

geometry = gpd.read_file(fn_vec)
df = geometry[(geometry.bounds['minx'] >= x_min) &
         (geometry.bounds['maxx'] <= x_max) &
         (geometry.bounds['miny'] >= y_min) &
         (geometry.bounds['maxy'] <= y_max) ]
print(len(df))

print(df['BlockGroup'].sum())

print(df['BlockGroup'].max())

show(array)

I am double checking with the input sums and there is a huge difference. The values per geometry should be added to the corresponding cells. At least cells look ok by visual inspection.
I think add option is not working. But why?


texts: Manually transforming rotated lat/lon to regular lat/lon?

texts:
First I should clarify I don't have previous experience with the field, so I don't know the technical terminology. My question is as follows:

I have two weather datasets:

The first one has the regular coordinate system (I don't know if it has an specific name), ranging from -90 to 90 and -180 to 180, and the poles are at latitudes -90 and 90.
In the second one, although it should correspond to the same region, I noticed something different: latitude and longitude were not the same, as they have another reference point (in the description is called a rotated grid). Together with the lat/lon pairs, comes the following information: southern pole lat: -35.00, southern pole lon: -15.00, angle: 0.0.

I need to transform the second pair of lon/lat to the first one. It could be as simple as add 35 to the latitudes and 15 to the longitudes, since the angle is 0 and it seems a simple shifting, but I'm not sure.

Edit: The information I have about the coordinates is the following

http://rda.ucar.edu/docs/formats/grib/gribdoc/llgrid.html

Apparently, the second coordinate system is defined by a general rotation of the sphere

"One choice for these parameters is:

The geographic latitude in degrees of the southern pole of the coordinate system, thetap for example;
The geographic longitude in degrees of the southern pole of the coordinate system, lambdap for example;
The angle of rotation in degrees about the new polar axis (measured clockwise when looking from the southern to the northern pole) of the coordinate system, assuming the new axis to have been obtained by first rotating the sphere through lambdap degrees about the geographic polar axis, and then rotating through (90 + thetap) degrees so that the southern pole moved along the (previously rotated) Greenwich meridian."

but still I don't know how to convert this to the first one.


texts: MODIS MOD11A2 does not update frequently its LST images?

texts:
I am working with the MODIS satellite, product MOD11A2, band LST_Day_1km. To access data I use the MODIS REST SERVICE.

With the following url I get the data from date 2020050 (2020-02-19) to date 2020071 (2020-03-11):

https://modis.ornl.gov/rst/api/v1/MOD11A2/subset?latitude=10&longitude=10&band=LST_Day_1km&startDate=A2020050&endDate=A2020071&kmAboveBelow=1&kmLeftRight=1

The product has a periodicity of 8 days but I only get one image, the request above returns me only data from day 2020-02-26 when it should return me also data from date 2020-03-05 (8 days after)

I searched in the documentation information about how often the data source is updated with the new data but I can't find anything and it is strange than it takes more than 16 days to add a date.


texts:
How to transform Earth-Centered Inertial J2000 CRS to WGS84 CRS in QGIS? I searched J2000 in spatialreference.org,
however nothing found. I don't know the proj4 string of J2000.

Any advice would be thankful.


texts: How to visualize inundation using Sentinel 1 SAR Data?

texts:
I am trying to visualize flooding after hurricane Harvey that made landfall on August 25-26 2017.

Based on this video, I've been able to obtain the basics of Sentinel-1 SAR imagery. Sentinel 1 comes in many flavors I am not sure what is the best approach to visualize inundation.

Splitting up my question into smaller ones:

What is the best way to visualize inundation using SAR data?
1.1  What Acquisition Mode?  (I choose IW, based on ESA description)
1.2  What Polarization should I use?  (I randomly choose [VV,HV])
1.3  Should I handle ascending and descending imagery separately? (I choose to split it up)
1.4  How to go to create a false color map? (I've tried (r,g,b)(HH,HV,HH/HV) and (r,g,b)(HH,HV,HH-HV).

In addition it would be great to understand what min-max values to choose for each band.

I've created a little testing script in Google Earth Engine (link):

var ic = ee.ImageCollection("COPERNICUS/S1_GRD")
var ivp = {"opacity":1,"bands":["VV","VH","VV/VH"],"min":-20,"max":-5,"gamma":1};
var ivp2 = {"opacity":1,"bands":["VV","VH","VV-VH"],"min":-20,"max":0,"gamma":4.935};

// Harvey made landfall on August 25
var date_start = ee.Date("2017-08-25")
var date_end = ee.Date("2017-09-02")

var ic_vvvh = ee.ImageCollection('COPERNICUS/S1_GRD')
        .filterDate(date_start,date_end)
        .filterMetadata("transmitterReceiverPolarisation","equals",["VV","VH"])
        .filter(ee.Filter.eq('instrumentMode', 'IW'))
        .map(function(image) {
          var edge = image.lt(-50.0);
          var maskedImage = image.mask().and(edge.not());
          return image.updateMask(maskedImage);
        });

function add_ratio_band(image){
  var new_image = ee.Image(image)
  var i_ratio = image.select("VV").divide(image.select("VH"))
  i_ratio = i_ratio.rename("VV/VH")
  new_image = new_image.addBands(i_ratio)
  return new_image
}

function add_difference_band(image){
  var new_image = ee.Image(image)
  var i_ratio = image.select("VV").subtract(image.select("VH"))
  i_ratio = i_ratio.rename("VV-VH")
  new_image = new_image.addBands(i_ratio)
  return new_image
}

ic_vvvh = ic_vvvh.map(add_ratio_band)
ic_vvvh = ic_vvvh.map(add_difference_band)

var ic_vvvh_desc = ic_vvvh.filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'));
var ic_vvvh_asc = ic_vvvh.filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'));

Map.addLayer(ic_vvvh_asc,ivp,"S1 [VV,HV,VV/HV]")
Map.addLayer(ic_vvvh_asc,ivp2,"S1 [VV,HV,VV-HV]")
Map.setCenter(-95.10,29.844,14)

I am cross referencing with areal imagery:
https://storms.ngs.noaa.gov/storms/harvey/index.html#12/29.8395/-95.0695

And the result look promising:


texts: How to keep polygon from becoming polyline when converting polygon feature class to .dwg

texts:
Is it possible to convert a polygon feature class from a geodatabase to a CAD .dwg while keep the features as a polygon in AutoCAD? Currently the polygons appear only as polylines.

I am currently using ArcGIS 10.2 and AutoCAD Map 3D 2010.

The Conversion Tool - Export to CAD appears to produce the necessary point, polyline and polygon data in the ArcMap Table of Contents after the tool has run, but only the polyline appears when I try to open the .dwg file into AutoCAD Map 3D.

I have attempted to convert the polygon to a variety of CAD types including DWG_R2010 and DWG_2007. The DXF and DGN files do not appear to be recognized by AutoCAD when converted.

It has been several years since I have touched AutoCAD, and any hints (is it even possible?) would be greatly appreciated. Perhaps I am just not going about it in the correct manner?


texts: Replacing raster values with closest value from list using R?

texts:
The task: I'm creating a function to assess flood damage in R. The inputs are:

Dataframe with water depth values and corresponding damage values in percentage.
Raster representing water depth

The raster values must be approximated to the closest value in the first column of the df and then create a new raster with damage values using corresponding values from the second column).

library(raster)
damage <- cbind(c(0.0,0.2,0.4,0.6,0.8,1.0),c(0, 5, 10, 40, 70, 100))
r.wd <- raster(myraster.tif)

I've tried using which.min(abs(damage[,1]-r.wd)) but doesn't work.


texts: Error generating chart: The image collection is empty in GEE

texts:
I am building a Landsat image collection time series and then trying to plot a time series of values for my region of interest by mean reduction. However, upon charting, I get an error message stating the image collection is empty. Why would the collection be empty?

Here is my code:

var Akdarya = ee.FeatureCollection("users/xodjayevshovkat/Akdarya_Crop_Class_Wheat_dissolve");

// Import S2 image collection.
var S2 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
  // Filter collection.
  .filterDate('2018-01-01','2018-12-31')
  .filterBounds(Akdarya)
  .filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 20);

// Define a cloud masking function.
var maskcloud1 = function(image) {
var QA60 = image.select(['QA60']);
  return image.updateMask(QA60.lt(1))//.clip(cotton);
};

// Define an index function (return only NDVI).
var NDVI = function(image) {
  return image.expression(
    '(NIR - RED) / (NIR + RED)',
    {
      'NIR': image.select('B8'),
      'RED': image.select('B4'),
    }).rename('NDVI').copyProperties(image, image.propertyNames());
};

// Apply the cloud mask function.
var NDVIcol = S2.map(maskcloud1)
 // Apply the NDVI transformation.
 .map(NDVI);

// Calculate the per-pixel median NDRE for all images in collection.
var NDVImed = NDVIcol.median();

// Define vis params.
var NDVIvis = {
  palette: ['#d73027', '#f46d43', '#fdae61', '#fee08b', '#d9ef8b', '#a6d96a',
  '#d73027', '#f46d43', '#fdae61', '#fee08b', '#d9ef8b', '#a6d96a',
  '#d73027', '#f46d43', '#fdae61', '#fee08b', '#d9ef8b', '#a6d96a',
  '#d73027', '#f46d43', '#fdae61', '#fee08b', '#d9ef8b', '#a6d96a',
  '#d73027', '#f46d43', '#fdae61', '#fee08b', '#d9ef8b', '#a6d96a',
  '#d73027', '#f46d43', '#fdae61', '#fee08b', '#d9ef8b', '#a6d96a'],
  min: -0.5,
  max: 0.9};

// Plot a time series of
var regions = Akdarya.toList(Akdarya.size())
var regionCount = Akdarya.size().evaluate(function (count) {
  for (var i = 0; i < count; i++) {
    var plotNDVI = ui.Chart.image.seriesByRegion({
      imageCollection:NDVIcol,
      regions: ee.FeatureCollection([ee.Feature(regions.get(i))]),
      reducer: ee.Reducer.mean(),  //median //mean
      xProperty: 'system:time_start',
      seriesProperty: 'Farmers'})
      .setChartType('LineChart').setOptions({
        title: 'NDVI Akdarya 2018',
        hAxis: {title: 'Date'},
        vAxis: {title: 'NDVI'},
    });
 print(plotNDVI);
  }
})
// Display the results.
//print(plotNDVI);
Map.addLayer(NDVImed.clipToCollection(Akdarya), NDVIvis, 'NDVI');


texts: Set the Cache-Control max-age of ArcGIS Server 10.1 tiles to a custom value

texts:
Using ArcGIS Server 10.1 I have a cached service. The service is set to allow clients to cache tiles and the images do contain max-age values:

Cache-Control   private, max-age=86400

However, I would like to change this value.

I could find the cacheControlMaxAge property in the ESRI documentation for 10.3.1, but this property does not exist in my 10.1 service. These are the available properties:

"properties": {
  "virtualCacheDir": "/rest/directories/arcgiscache",
  "maxImageHeight": "2048",
  "textAntialiasingMode": "Force",
  "maxImageWidth": "2048",
  "enableDynamicLayers": "false",
  "dynamicDataWorkspaces": "",
  "supportedImageReturnTypes": "URL",
  "disableIdentifyRelates": "false",
  "isCached": "true",
  "maxScale": "250",
  "maxBufferCount": "100",
  "maxRecordCount": "1000",
  "schemaLockingEnabled": "true",
  "filePath": "\\\\dir\\a.msd",
  "cacheOnDemand": "false",
  "useLocalCacheDir": "true",
  "virtualOutputDir": "/rest/directories/arcgisoutput",
  "outputDir": "\\\\dir\\arcgisserver\\directories\\arcgisoutput",
  "cacheExtent": "full",
  "minScale": "200000",
  "maxDomainCodeCount": "25000",
  "tilingScheme": "2",
  "ignoreCache": "false",
  "antialiasingMode": "None",
  "clientCachingAllowed": "true",
  "cacheDir": "\\\\dir\\arcgisserver\\directories\\arcgiscache"
 }

How can I set the max-age to be a custom value?


texts:
I'm trying to develop a script, something that can be run every hour/day/week etc that will look at some shapefiles (initially downloaded from a WFS) compare them to the WFS and download and append any changes.

I have absolutely no idea how to start developing something like this, i'm not even sure if python is the right thing. I have looked rather extensively online for help with this and the best I came up with was this website (http://geopython.github.io/OWSLib/#wfs) but there is nothing on WFS in it, which isn't helpful.

I found the documentation for WFS be that wasnt particularly helpful for me. My background is primarily Geography, geology and basic computer science. I've worked with JS, HTML, PHP, CSS, Python but am in no way a developer.

If anyone is able to walk me through the steps to beginning to work with WFS outside of ArcMap that would be amazing.


texts:
I'm a linguist, not a geographer, and a first-time poster here.

I have some data about the use of a particular phrase in England between the years of 1225 and 1500. I have geolocations for the places where various manuscripts (MSs) were produced. I have the year (roughly) of the MSs, and I have token counts of the phrase in each MS. There are about 100 MSs, and thousands of tokens, but only tens of locations.

I want to produce a map with a slider that controls the year. I want a dot to appear in the MS location in the year when the MS is produced, and the size of the dot should grow as the number of tokens from that location accumulates over time.

So, if there was one MS in London in 1250 with 10 tokens, a small dot should appear at London when the slider is at 1250. As the slider is moved to 1275, the dot should remain stable. In 1300, a new MS is produced with 20 more tokens, so the area of the dot should now grow.

I've produced this in CartoDB, but it's not quite right. The number for the slider should show the year, not just an integer from 1-64, and I can't vary the sizes of the dots.

Is there something I can fix here or is there a better tool for doing this? I'm not married to any one platform, but as I said, I'm not a geographer, and I'm also not a web designer.


texts: Adding minimum point Z of polygon using ArcGIS Field Calculator?

texts:
In the middle of a ModelBuilder tool I created, I have a Z enabled feature with multiple polygons.
I want to add on a specific field of this feature, the lowest Z value of each polygon. I know this is easy to do using calculate geometry, but this option is not available on ModelBuilder and the "Add Z Information" tool is giving me error 00245 and the "Add surface information" is taking to long to process.

So I want to do this using a Python expression on the calculate field tool.

If this is possible, what expression can I use?


texts: ogr2ogr keep geometry when converting to CSV?

texts:
Using ogr2ogr, how to keep the geometry when converting to CSV?
( Nothing relevant here: http://www.gdal.org/drv_csv.html )

$ ogr2ogr -f "CSV" MyLayer.csv MyLayer.geojson

Here is the begin of the result file:

$ head -n4 MyLayer.csv
OBJECTID,ID,VAR1,VAR_TXT,VAR2,VAR_TXT
677406,145477,0,AB94,51,rocks � distribute
267075,145400,0,AB94,47,road acces � distribute
101700,145434,0,AB94,47,road acces � distribute

I would like to keep the geometry (X,Y) on the two first columns for example:

$ head -n4 MyLayer.csv
X,Y,OBJECTID,ID,VAR1,VAR_TXT,VAR2,VAR_TXT
152.282,635.743,677406,145477,0,AB94,51,rocks � distribute
150.624,632.052,267075,145400,0,AB94,47,road acces � distribute
151.194,631.049,101700,145434,0,AB94,47,road acces � distribute

And I encountered an error if I translate the same data but as a shapefile (or a *.mif Mapinfo file) to CSV:

$ ogr2ogr -f "CSV" MyLayer2.csv MyLayer.shp
FAILED: Layer MyLayer already exists, and -append not specified.
        Consider using -append, or -overwrite.
ERROR 1: Terminating translation prematurely after failed
translation of layer MyLayer (use -skipfailures to skip errors)

It seems it absolutely wants to keep using the previously named file (from the conversion from the GeoJSON file above).

Usage of ogr2ogr:

Usage: ogr2ogr [--help-general] [-skipfailures] [-append] [-update]
           [-select field_list] [-where restricted_where]
           [-progress] [-sql <sql statement>] [-dialect dialect]
           [-preserve_fid] [-fid FID]
           [-spat xmin ymin xmax ymax] [-geomfield field]
           [-a_srs srs_def] [-t_srs srs_def] [-s_srs srs_def]
           [-f format_name] [-overwrite] [[-dsco NAME=VALUE] ...]
           dst_datasource_name src_datasource_name
           [-lco NAME=VALUE] [-nln name] [-nlt type] [-dim 2|3|layer_dim] [layer [layer ...]]


texts: Invalid input syntax for type double precision with POST request to Maps API

texts:
I'm just trying to replicate this example but using a different account (libregis), using a raster table called foto_pnoa, created just dropping a tif file when creating a new dataset. And I'm getting this error

{"errors":["Postgis Plugin: ERROR:  invalid input syntax for type double precision: \"76.4293851819992
 ,10\"\nCONTEXT:  SQL function \"_raster_constraint_info_scale\" statement 1\n\nin executeQuery Full sql
  was: 'select r.r_table_schema sch, r.r_table_name tab, r.r_raster_column col, greatest(abs(r.scale_x
 ), abs(r.scale_y)) scl from raster_overviews o, raster_columns r where o.r_table_schema = 'public' and
  o.r_table_name = 'foto_pnoa' and o.r_raster_column = 'the_raster_webmercator' and r.r_table_schema
 = o.o_table_schema and r.r_table_name = o.o_table_name and r.r_raster_column = o.o_raster_column ORDER
  BY scl ASC'\n  encountered during parsing of layer 'layer0' in Layer"]}

when sending a POST request to http://libregis.cartodb.com/api/v1/map, using this as POST body (map config)

{"version":"1.3.1","layers":[{"type":"cartodb","options":{"sql":"select * from foto_pnoa","cartocss"
:"#foto_pnoa {raster-opacity: 0.5;}","cartocss_version":"2.3.0","geom_column":"the_raster_webmercator"
,"geom_type":"raster"}}]}


texts: GPS Device using non-proprietary software and connector

texts: Are there GPS devices that use non-proprietary software and connector? Ideally USB C for the connector. I've been using my Android phone with GPS Essentials and am frustrated by the constant dropping of tracking by the software.


texts: Import of GeoServer WFS into uDig not working

texts:
I'm trying to follow the example of how to integrate GeoServer with uDig that is on the uDig web page:

http://udig.refractions.net/files/docs/latest/user/getting_started/walkthrough2/WMSandWFSIntegration.html

I'm using a virgin install of both GeoServer 2.6.2 and uDig 1.5.0 RC1 on a Mac.

I don't seem to be able to add features using WFS from the GeoServer to uDig.

If I drag the WFS 1.0.0 link onto the map in uDig a 'processing drop event' popup is displayed temporarily, but nothing is added to the map. I can't see any errors in the GeoServer log, and I can't find a log in uDig - but no error messages are displayed.

I can confirm that the adding a WMS does work, which confirms that the GeoServer is working properly.

I've tried 'add data' in uDig - but when I select 'Web Feature Server' I get an error "An error occurred during page transition. This is most likely a programming error. Please report it". Is this related?

Can some one help shed some light on why the WMS isn't working?

It's possible that the documentation is out of date (it refers to an older version of geoserver), and it's also possible that the geoserver isn't configured out of the box the way that it was when the docs were written. I'm just learning these tools, and so I'm hoping that the answer is trivial.


texts: Using Cumulative Geometry Calculation to Add Field in ArcPy

texts:
I am using Desktop 10.6, Basic license, Spatial Analyst, Python 2.6, 2019 AD.

I have 4 overlapping features in a polygon feature class. I'd like to sort the features by the depth field, calculate the acreage for the largest depth, then move to the next depth, calculate its acreage and then add it to the previous depth, and so on. This code currently does not work past adding the cumulative field. In the attached photo, the cumulative depth is what I am trying to generate. The acres field is to show the progression.

Working code:

contours = arcpy.GetParameterAsText(0)

desc = arcpy.Describe(contours)
arcpy.env.workspace = desc.path

arcpy.AddField_management(contours, "Cumulative_AC", "DOUBLE")
arcpy.AddField_management(contours, "Acres", "DOUBLE")
CursorFieldNames = ["SHAPE@", "Cumulative_AC", "Acres"]
RunTl = 0

with arcpy.da.UpdateCursor(contours, CursorFieldNames, sql_clause=(None, 'ORDER BY DEPTH DESC')) as cursor:
    for row in cursor:
        AreaValue = row[0].getArea('GEODESIC', 'ACRES')  # Read area value as double
        RunTl += AreaValue  # Capture running total
        row[1] = RunTl  # Write total area value plus the previously added values to field
        row[2] = AreaValue
        cursor.updateRow(row)
del RunTl


texts:
As stated here: https://trac.osgeo.org/gdal/wiki/CMake
GDAL is not yet ready to be compiled with CMake:

  News

  2012-04-28:

  Pilot project of CMake configuration for GDAL is being developed internally by Kitware. It should be delivered soon,
  so stay tuned. Once first beta is released to the public, it will be
  available for testing, improvements

But this is pretty old.
Any idea on the roadmap if always planned?
(I'd like to know which flag I can switch to enable python bindings because now it cannot find osgeo within QGIS).

I'm using Ubuntu 18.04.


texts:
I am using ggmap package to generate routes. The resulting data frame will contain the routes information from the given source to the destination.

   library(ggmap)
    route_df <- route(from = "Cambridge, 02138, USA",
              to = "Massachusetts Ave, 02139, USA",
              mode = "driving", structure = "route", alternatives = TRUE)

I want to verify if a particular co-ordinate (-71.13857, 42.38039) fall in this route or not.
I also want to verify if a coordinate (-71.13857, 42.38039) is part (one of the waypoints) of the route "A" resulted from the "route_df" variable. Is there any package that I can use to do so?


texts: GeoWebCache: Is there a way to see cache size per zoom level?

texts:
I seeded a layer group (of data imported from OpenStreetMaps via osm2pgsql) from levels 0 to 15 (PNG; EPSG:900913), it completed successfully, and now I am seeding again only for levels 16 and 17. I would like to figure out how much disk space is being taken by each zoom level.

All I can see now is the total for the entire layer group -- 46 GB as I write this. Is there a relatively simple way to obtain the usage broken down by zoom level?


texts: Avoiding labeling features if overlapped by another layer in QGIS?

texts:
With QGIS 2.12.2, how can I set up layer labeling to avoid placing labels where features from another layer already exist?

For example, if I have a stream/river polyline layer that contains lake "centerlines", and I place a "lake" polygon layer above it in the drawing order, I don't want the river layer to place a label inside the lake.  Instead, I would rather have the river labeled outside of the lake (as needed).  That way, I can place labels from the lakes layer and I don't run into label collisions.

Here is an example, where (I have intentionally put the lines on top for visual purposes) what I am hoping to achieve is no river center-line labels shown inside the lake polygon:


texts:
I have a dataset containing buildings where each building has a certain value/score. What I want to do is make clusters of these buildings depending on their values/scores. What I did so far is; I created a buffer around each building, so it is possible to determine which and how many other building are intersected by a building. This buffer has the corresponding value/score of the building where the buffer is built around.

What I want to do now, is that I want to dissolve these buffers as one, when/if the buildings that the buffer intersects have the same or higher value/score than the building the buffer is being built around (origin building).

So I'm thinking that a 'simple' if or where statement should do the trick. However, I'm a beginner in PostGIS and therefore I'm not sure how and where to write this in PostGIS.

My SQL query so far is:

        WITH sum AS (
WITH buffertest AS (
    SELECT Building_score, Building_ID, ST_Buffer(geom, 30)::geometry(polygon,28992) AS buffer
     FROM postgis_export)

        SELECT b.Building_score AS bufferscores, b.Building_ID AS bufferID,  a.Building_score AS building, b.buffer
FROM postgis_export a, buffertest b WHERE ST_intersects(b.buffer, a.geom)  AND b.gebouw_sco - a.gebouw_sco <= 0 )

        SELECT bufferscores, bufferid, SUM(gebouw), COUNT(*), buffer
FROM sum
GROUP BY bufferscores, bufferid, buffer
ORDER BY bufferid


texts:
I'm really quite new to QGIS, I want to style polygons fill with raster imagws based on a attribute contained in the property table.

What do I have to write in the description in order to get for each polygon the appropriate raster image?


texts:
I am trying to extract country borders between two points, using Python. The border is a Shapely Polygon.exterior or LinearRing, the start and end points are Shapely Points.

My idea is the following:

Get the closest points or projections for both points along the border
These points split the border into two splines (or LineStrings). Select the shorter one, based on length (ideally reprojected to meters, but it might be overkill).
Decide if the first and/or the last point of the selected spline is needed, to avoid zig-zag like splines at the ends. The given start and end points must be included in the final spline.
Calculated border is: start_point + calculated_spline + end_point

So far I'm a bit stuck with the following concepts:

Shapely doesn't support projection on LinearRing, only LineString. How can I calculate the closest points efficiently on a LinearRing? I believe projection isn't needed, as I'll be selecting whether to keep the end points in step 3. Should I just loop over all points of the LinearRing and select the closest one?
How can I split a LinearRing into two LineStrings, based on two of their points?
Is it true that when comparing relative segment lengths of a LinearRing, I can just measure their length in WGS84?
When I have the LineString, should I use .project for deciding weather the endpoints should be included, to avoid zig-zag like splines?


texts:
I can't find OpenLayers Plugin to install in QGIS 3.6.

I want to open Google Earth in order to add coordinates and open satellite image to start mapping.


texts: Change the border colour and line width of a polygon with pyQGIS

texts:
I have worked out how to change the color of a polygon using pyQGIS but cannot change the stroke color or line with based on the questions detailed here and here.  My example code below that works to make the polygon transparent fill but not for the blue border and desired linewidth.
# Attach core modules
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
from qgis.core import *                        # attach main QGIS library
from qgis.utils import *                       # attach main python library
import os                                      # attach operating system library
import processing

# Set the working directory - recommended to use unix forward slash for path
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
wd = "C:/TEST"      # Set work directory
os.chdir(wd)                                              # Change the directory
os.getcwd()                                               # Confirm change

# Set a variable for the current project instance
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Prj = QgsProject.instance()                         # Object for current project

# Save the project to this file name
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
pnm = "200924_Test.qgs"                      # Project file name
pnm = wd + "/" + pnm                               # Concat. with path
Prj.write(pnm)                                     # Save the project

# Create an array [] object with the polygon vertices for MRE model area
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
vrtcs = []
vrtcs.append(QgsPointXY(541600,6504400))
vrtcs.append(QgsPointXY(541600,6505100))
vrtcs.append(QgsPointXY(542160,6505100))
vrtcs.append(QgsPointXY(542160,6504400))

# Create a polygon from the vertices
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
ply_01 = QgsGeometry.fromPolygonXY([vrtcs])

# Create a feature object then append the polygon into the feature
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
ftr = QgsFeature()
ftr.setGeometry(ply_01)
print(ftr.geometry())

# Create a  layer for the feature, in the desired CRS  28351  GDA94/MGA 51
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
lyr = QgsVectorLayer('Polygon?crs=epsg:28351', '200924_MRE_Area',"memory")
Prj.addMapLayers([lyr])

# Make the layer editable, add the feature with the polygon and save the layer
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
lyr.startEditing()
lyr.addFeature(ftr)
lyr.commitChanges()

# Save as a shapefile
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Fl_ou = 'Test'
Fl_ou = wd + '/' + Fl_ou

options = QgsVectorFileWriter.SaveVectorOptions()
options.driverName = "ESRI Shapefile"

# Change the fill, stroke colour and line width of the the polygon
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
rndr = lyr.renderer()                   # object to change rendering
rndr.dump()                             # report current rendering

sym = rndr.symbol()                     # object to change the symbology of layer
syms = sym[0]                           # use first layer in QGIS tree
syms.setColor(QColor("transparent"))                 # remove the fill colour
syms.symbolLayer(0).setStrokeColor(QColor("blue"))   # change the stroke colour (Fails)
syms.symbolLayer(0).setWidth(3)                      #
lyr.triggerRepaint()                                 # Refresh canvas

Additionally, the questions (linked) also suggest the following lines of code should be included at the end of this process - I would like to understand why they are required
iface.legendInterface().refreshLayerSymbology(layer)
    iface.mapCanvas().refresh()


texts: Mantel correlogram outside of the range [-1,1]

texts:
I want to represent the correlogram from the package ncf of the working residual my_glm$residual of logistic regression my_glm from glm function. When I represented it, I got a strange output that the Mantel similarity index sometimes is outside of the range [-1,1] (see the attached figure). I don't understand the problem here.
This is how I computed it :
knn <- spdep::knn2nb(spdep::knearneigh(xy_sp))
dist_class <- 2*(min(unlist(spdep::nbdists(knn, xy_sp)))) # double of the minimum distance between the sampled sites
RESIDU <- my_glm$residuals
rnd_residu <- sample(1:length(RESIDU))
spat.cor <- correlog(xy[rnd_residu,1], xy[rnd_residu, 2], RESIDU[rnd_residu], increment = dist_class, resamp = 10,latlon = FALSE)
par(mfcol=c(1,2))
plot(spat.cor, pch=16, col="firebrick3", ylab="Correlation", xlab="Distance", main="Spatial correlogram", font.lab=2)

xy_sp represent a spatialPoints :
Object of class SpatialPoints Coordinates: min       max x  467698.7  656887.7 y 5450967.1 5681736.1 Is projected: TRUE proj4string : [+proj=utm +zone=19 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0] Number of points: 196


texts: Applying symbology based on algorithm using ArcPy with ArcGIS Pro

texts:
I have 100 points and I want to apply a color palette to the points so that no two points within 25 miles have the same color. The issue I'm trying to prevent is having points close to each other being mistaken for similarity.
What I want this to do is, start with the first point, fill in a field with a value of "ONE" and then select a random point that's not within 25 miles of it. That next selected random point will also get a field value of "ONE". I want this to repeat until I have ten records (points) with a field value of "ONE" and that no two points in this sample is within 25 miles of any other point in the sample.
# Process: Select Layer By Attribute
arcpy.SelectLayerByAttribute_management('InputPoints100', 'NEW_SELECTION',"OBJECTID = 1")

# Process: Calculate Field
arcpy.management.CalculateField("InputPoints100", "Color", '"ONE"', "ARCADE", "ONE")

arcpy.SelectLayerByLocation_management('InputPoints100', 'WITHIN_A_DISTANCE', 'InputPoints100',"20 Miles", 'NEW_SELECTION', 'INVERT')

layer = "InputPoints100"
count = 1

def SelectRandomByCount (layer, count):
    import random
    layerCount = int (arcpy.GetCount_management (layer).getOutput (0))
    if layerCount < count:
        print("input count is greater than layer count")
        return
    oids = [oid for oid, in arcpy.da.SearchCursor (layer, "OID@")]
    oidFldName = arcpy.Describe (layer).OIDFieldName
    delimOidFld = arcpy.AddFieldDelimiters (layer, oidFldName)
    randOids = random.sample (oids, count)
    oidsStr = ", ".join (map (str, randOids))
    sql = "{0} IN ({1})".format (delimOidFld, oidsStr)
    arcpy.SelectLayerByAttribute_management (layer, "", sql)

SelectRandomByCount(layer, count)

arcpy.management.CalculateField("InputPoints100", "Color", '"ONE"', "ARCADE", "ONE")


texts: Naming output files where original file name retained without file extension using ModelBuilder iteration?

texts:
I am using iterate raster in ModelBuilder to create hillshade and slope simultaneously for many rasters within one folder.

My problem is when I name the output files.

Whenever I use, %Name%, the model fails because the resulting output names include the file extension.

If I use %n%, it only creates one file.

Is there another way of naming output files where the original file name is retained but without the file extension?


texts: Seeking USA congressional district to county relationship files?

texts:
This must be a common enough problem that other users have come across it. I've finagled my way around inconsistencies and dropped observations, but now I'm reaching out to see if anyone knows of a 'holy grail' type spatial dataset.

I am in search of USA congressional district to county relationship files going back to the 1980s. While many components of said files are available online via the US Census Bureau (https://www.census.gov/geo/maps-data/data/cd_national.html), they do not release them for every year.

I've tried combining different online sources (combining this CD dataset (http://cdmaps.polisci.ucla.edu/) with the atlas of historical county boundaries (https://publications.newberry.org/ahcbp/)) in a GIS environment; however, inconsistencies in the shapefiles (see here:

lead to too many (or too few) relationships.

Does anyone have such a file on hand or know where I could obtain it?


texts: How could I switch data connections in case of server failover?

texts: Here is the scenario.  I want to be able to switch the end users database connections automatically in case we ever have to switch to our backup SDE server.  My first thought was to set up sde connections to the backup server on the users machines, then write a batch file that would simply rename the connections to the users standard connections.  I'm hoping to find a more automated solution.  Ideally find some way to tell if a connection can't be made, then switch to the alternate connections.  Has anyone implemented something like this?  Any help much appreciated.


texts: QGIS 2.8.2 - Using spatial expressions in labelling

texts:
I'm creating some hard copy map output from QGIS.

I want to use an expression to suppress labels for some features on a layer when another feature from another layer intersects with those features.

E.g. I want to suppress labels for river (line) features when they intersect with lake (polygon) features.

From my understanding of the documentation the expression should look something like this in the Layer Properties label dialogue of the 'Rivers' layer:

CASE WHEN  intersects( ( $geometry , geometry(getfeature(Lakes_dams)) = 1 THEN  null  ElSE title( "NAME") END

This is invalid so what is the correct syntax?


texts: Trouble calculating areas in QGIS using field calculator

texts:
I am testing various projections to calculate areas using ICES statistical shapefiles and the QGIS field calculator.

I am using ICES areas downloadable from the following link: http://geo.ices.dk/download.php?dataset=ices_ref:ices_areas. (Please note that the calculations described below made by both ICES and myself were made using a densified version of the data available on the website).

I have been in contact with ICES and they are getting very different results when they calculate their areas using the same projections as me. Can anyone see what I'm doing wrong with my projection or process? To change the projection of my layer I used the 'save as' option to convert it and then imported the new file into my project.

Please see a copy of the layers 'properties' found under the metadata tab all the way at the bottom of this question.

This is what I entered into field calculator: I am certain this is fine.

The results: the left column titled 'Area_km2' represents the correct ICES calculations and the the column on the right titled '3035area' represents my calculations which seem very wrong. Both calculations used CRS EPSG 3035 - ETRS89/LAEA.

Meta data --> Properties:

General
Storage type of this layer
ESRI Shapefile

Description of this provider

OGR data provider (compiled against GDAL/OGR library version 1.11.0, running against GDAL/OGR library version 1.11.0)

Source for this layer
/Users/aninahenggeler/Documents/Humboldt/Thesis /High seas proposal/GIS & Google Earth/ICES/NewICES densified/ICESdensifiedareasclippedat899EPSG3035.shp

Geometry type of the features in this layer
Polygon

The number of features in this layer
65

Editing capabilities of this layer
Add Features, Delete Features, Change Attribute Values, Add Attributes, Delete Attributes, Create Spatial Index, Fast Access to Features at ID, Change Geometries, Simplify Geometries, Simplify Geometries with topological validation

Extents
In layer spatial reference system units
xMin,yMin -6150.61,1579195.62 : xMax,yMax 6185275.19,7365876.42

Layer Spatial Reference System
+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs


texts: How to get a bbox selection of huge shapefile as GeoJSON?

texts:
I have a 4.29 GB Shapefile from which I want to extract features as GeoJSON, defined by a bounding box.

This should be as fast as possible, because I want to use it in a web application.

My first attempt at this was to convert to SQL using shp2psql and import it into PostGIS. I created a spatial index (gist) and did a full ANALYZE VACUUM on the table. The VM has 4 GB of memory.

The following query works, but took 623 seconds(!):

SELECT
    ST_AsGeoJSON(geom)AS geometry,
FROM
    shapes
WHERE
    geom && ST_GeogFromText(
        'POLYGON((5.117705762386322 52.09243663877818,5.118895322084427 52.09243663877818,5.118895322084427 52.09314276719956,5.117705762386322 52.09314276719956,5.117705762386322 52.09243663877818))'
    )
AND ST_Intersects(
    geom,
    ST_GeogFromText(
        'POLYGON((5.117705762386322 52.09243663877818,5.118895322084427 52.09243663877818,5.118895322084427 52.09314276719956,5.117705762386322 52.09314276719956,5.117705762386322 52.09243663877818))'
    )
)

An equivalent using ST_DWithin took even longer.

I decided to try an Amazon RDS PostgreSQL 9.3.1 instance with 68 GB of memory and 26 CPU's. Query time dropped to around 30 seconds. Faster but still not fast enough.

Preprocessing will take forever, even on Amazon RDS. And it will cost a lot of money.

I also tried Fiona and 'ogr2ogr -spat', but they're about as slow as PostGIS.

What other options exist for quickly getting GeoJSON out of a large Shapefile for a certain spatial extent?

Thanks!!


texts: Looking for Linux tools to visualize GPS+compass+focal-length EXIF data?

texts:
I have JPG images tagged with GPS and compass Exif-data. I took these images in a remote country especially to use it for mapping OSM. Now I am looking for a tool that would help me to use these images effectively:

I would like to see

an areal background image,
a dot using the GPS coordinates,
an arrow pointing to where the camera was looking (compass),
and a circular sector into compass direction that represents zoom/focal
length.

I have only ever seen these features in GeoSetter beta on Windows. But using Geosetter on Linux is not easy to run (background map missing) and also it is not really actively developed nor is it open-source.

Does anyone know a tool or set of tools to accomplish this on Linux?

Any recommendation for a QGis plugin that can do some of the features (so I might actually extend it myself)?


texts: Calculating polygon area using 2 methods in QGIS

texts:
I am trying to calculate the polygon area using 2 different methods I saw on previous posts here.
The first method: Layer > Open attribute table > Field Calculator/Ctrl+I or right-click shapefile > Open attribute table > Field Calculator/Ctrl+I). There is a function $area that will calculate the area of each row in the table.
Second method: Vector | Geometry Tools | Add/export geometry columns, which creates a new shapefile with area and perimeter (or length) columns added.
But I still get different results. I would like to know why it gives different results, even though, it's for the same area?
The problem I am trying to solve is to find the percentage of the population who live within 10 minutes walk from parks, and I need to be as precise as possible but I am super confused with these 2 different results.


texts: QGIS expressions don't evaluate as expected

texts:
I am using QGIS with a vector layer and the attribute table stored at this link:
https://docs.google.com/spreadsheets/d/1K7rRdlcSheclJ6MqzpBNwSZei911WDV1Bl2zBGV38-s/edit?usp=sharing

I then have the following code:

CASE

/*Kieran Winner*/
WHEN "r_Kieran" > "r_Wynne"
        THEN ramp_color('PCMarg', scale_linear(("r_Kieran"/"r_TOTAL"-"r_Possian"/"r_TOTAL"),  0,  1, 0, 1))
ELSE color_rgb(64, 224, 208)
END

Here's the problem:
The code evaluates incorrectly.

Take row 15 of the attribute table for example. (Where PD_Number=14)
You would expect that the condition would evaluate TRUE since
"r_Kieran"=397 and "r_Wynne"=225.

But, it evaluates FALSE in QGIS, as indicated by the colour that is displayed on the map!

I've spent hours simplifying code, and pulling my hair out over this.


texts: Identifying small polygons that are inside larger bounding polygons within same feature class using ArcPy?

texts:
I've tried digging around for an answer to this question, but I've only been able to find answers that are either too broad or too different from my case, so here goes!

Context:

I have a layer of city parcels, and I am trying to find all the smaller parcels that lie within bigger container parcels (i.e. I am trying to identify condos, trailers, or any other type of housing that has its own unique address, owner name, and ID that is different from the container parcel). The image below illustrates my task: the red parcel is the container and the blue circled parcel is one of several that I am trying to programmatically identify.

My goal is to create a dictionary where the key of each entry is the smaller parcel's ID# and the value is the container parcel's ID# (basically a relational table of bounded parcels).

I have tried working with the Polygon Neighbors tool, but it's redundant in my case for two reasons:

At the city-scale, there are ~44,000 parcels, and each of those parcels has anywhere from one to tens of neighbors, and the list would be frighteningly endless. You would think that I could just run the tool and then query out all the entries that only have one neighbor (meaning they only have one enclosing neighbor), but...
There are small parcels within larger parcels that are not completely enclosed AND have more than one neighbor (see picture below), so even if I ran the tool and weeded out table keys with repeating IDs, I would risk losing some info.

How I think the problem can be solved:

My fundamental idea has been to compare extents between and among parcels. Essentially, I'm thinking that if the bounding box of the smaller parcel is completely within the bounding box of the bigger parcel, I can have my script write those IDs to a dictionary. As far as coding that workflow, it seems that I should make use of two Search Cursors on the same feature class and compare geometries line-by-line.

My code:

parcels = path_to_testing_fc

parcelDict = {}
with arcpy.da.SearchCursor(parcels, ['SHAPE@', 'PARCEL_ID']) as sCurs1, sCurs2:
    for inside in sCurs1: #small parcels' loop
        for container in sCurs2: #container parcels' loop
            if inside[0].extent.within(container[0], 'PROPER'): #"proper" ensures my completely within rule
                parcelDict[str(inside[1])] = 'Contains ' + str(container[1])

When I run this, I get no errors, but the "parcelDict" dictionary is also empty, so I'm kind of in limbo.

Is the double search cursor a viable, efficient way to get this task done?In the same vein, is there a better, more efficient way to accomplish this goal?

Getting the snippet to work on my small test feature class is entirely different from getting the code to run through 44,000 parcels in a timely manner. Perhaps I'll look at the code for Polygon Neighbors and see if I can't glean something.

I opened a new console and my code from above suddenly started throwing errors (I probably forgot to delete cursors from other attempts at solving this problem), so here's an update on what my code looks like (still not giving errors and still not populating my dictionary...)

parcels = path_to_test_fc

parcelDict = {}
with arcpy.da.SearchCursor(parcels, ['SHAPE@', 'PARCEL_ID']) as sCurs1:
    for inside in sCurs1: #small parcels' loop
        with arcpy.da.SearchCursor(parcels, ['SHAPE@', 'PARCEL_ID']) as sCurs2:
            for container in sCurs2: #container parcels' loop
                if inside[0].extent.within(container[0], 'PROPER'):
                    parcelDict[str(inside[1])] = 'Contains ' + str(container[1])


texts: How to find the intersection areas of overlapping buffer zones in single shapefile?

texts:
Find the intersections of the buffer zones

I found an example where they tried finding the intersection of two different shapefiles -Intersect Shapefiles using Shapely. I've tried to adapt my code to do the same intersection, but I am only working with one shapefile with buffer features inside. The program is able to write the shapefiles, but the intersection aren't occurring, thus no data written to my output shapefile Anyone willing to help me understand where the problem lies or has an alternate way of finding intersection in a single shapefile?

# creation of the new shapefile with the intersection
for i in threshold:
    with fiona.open("threshold_buffer_shape_"+str(i)+".shp","r") as layer1:
        with fiona.open("threshold_buffer_shape_"+str(i)+".shp","r") as layer2:
            # We copy schema and add the  new property for the new resulting shp
            schema=layer2.schema.copy()
            schema['properties']['uGID'] = 'int:10'
            # create an empty spatial index object
            with fiona.open("intersection_region_"+str(i)+".shp","w","ESRI Shapefile",schema) as layer3:
                index = rtree.index.Index()
                # populate the spatial index
                for feat1 in layer1:
                    fid = int(feat1['properties']['GID'])
                    geom1 = shape(feat1['geometry'])
                    index.insert(fid,geom1.bounds)

            # get list of fids where bounding boxes intersect
                for feat2 in layer2:
                    geom2=shape(feat2['geometry'])

                    for fid in list(index.intersection(geom2.bounds)):
                        if fid!= int(feat2['properties']['GID']):
                            feat1 = layer1[fid]
                            geom1 = shape(feat1['geometry'])
                            if geom1.intersects(geom2):
                                props = feat2['properties']
                                props['uGID'] = feat1['properties']['uGID']
                                layer3.write({
                                    'properties':props,
                                    'geometry':mapping(geom1.intersection(geom2))
                                    })

  Overlapping buffer regions


texts: How to use OpenLayers Geolocate with the OpenGeo Suite Client SDK?

texts:
I’d like to have GPS functionality added to my OpenGeo Suite Client SDK application.

From what I understand it should be possible to use the Geolocate control from OpenLayers.   Preferably there would be a checkbox (placed on map.tbar) with the option of continuous tracking. If that’s not possible, a button that gets the current position would do as well (like this example).

How would I go about incorporating Geolocate?

Thanks.


texts: Extract .mxd file + data from map package .mpk?

texts:
I have two applications - a pc app which uses .mpk and a mobile app which uses .tpk
I can generate both, given the .mxd file - but I'd like an easy way to go from a .mpk to .tpk
The idea being that the .mpk is provided to the PC app which at a later time generates .tpk given an area of interest.

Given the .mpk I don't see an easy way of generating a .tpk without going back to an .mxd
I don't see anything in ArcPy (ArgGIS 10.2) which can:
a) open a .mpk and share as a .tpk (arcpy.CreateMapTilePackage_management wants .mxd input)
b) open a .mpk turn in back into an .mxd so I can use function highlighted in a)

I've noted that a .mpk is only an archive which can be opened and probably the .mxd extracted along with all relevant data, still if I didn't have to re-invent the wheel, that would be nice...


texts: Is there a limit for how many layers can be put into a group?

texts:
I am Running QGIS Lisboa 1.8.0 on a MacBook Pro. I have several vector layers with similar information which clutter the Layers panel.  I created a group to put them in but when I drag and drop the layers into the new group I get just so many into the new group and then the layers panel becomes unstable. I can start to drag another layer and the new group disappears. QGIS then exits generating a large error report.
I am using the method given in the manual for QGIS 1.7.0 and it appears to work up to the point where it all goes wrong.

Is this a known problem with this version of QGIS? I havent been able to find a reference to this particular problem.

Nigel


texts: Coloration of orthophoto different between QGIS and ArcGIS

texts:
I have an orthophoto created with Agisoft Photoscan that I would like to use for a map.

Why is there a difference in coloration between QGIS and ArcGIS?

For the record, the color showing in QGIS is identical to how it looks in Photoscan, whereas in ArcGIS it appears too red.


texts:
I want to draw a perfect rectangle. How can I do that with QGIS v.3+?

Here is a possible way to draw a rectangle in old QGIS versions.

Is there a similar solution for the latest QGIS version?


texts:
In previous versions of Carto it was easy to set the zoom and center-point of a map embed via URL parameters

http://<username>.carto.com/viz/<viz_id>/embed_map?zoom=3&center_lat=0&center_lon=0

via source

This was a very handy feature as it enabled you to save out a set of maps with different zoom and center-points from the same base map.

Does anyone know how to do this with Carto Builder?

These same query parameters do not work.


texts: Is there an inverse operation to ST_Shift_Longitude in PostGIS < 2.3?

texts:
ST_Shift_Longitude(geom) (renamed to ST_Shift_Longitude in PostGIS 2.2) adds 360 to longitudes below 0, thereby (for lon-lat coordinates like WGS84) moving the part of the world "left" (west) of the 0-meridian to the "right" (east) of the "right" (eastern) antimeridian. So this turns a flat -180°…+180° map "split" around the international date line into a flat 0°…+360° map "split" at Greenwich.

Is there any function that does the reverse (i.e. subtracting 360 from longitudes above 180) to get back from a flat 0°…+360° map "split" at Greenwich to a flat -180°…+180° map "split" at the antimeridian?

In PostGIS ≥ 2.3 I could use ST_WrapX(shifted_geom, 180, -360), I guess, but in PostGIS 2.1 and 2.2 ST_WrapX isn't available, yet.

I tried ST_MakeValid(shifted_geom) and ST_Transform(shifted_geom, 4326), but both return the shifted_geom unchanged.


texts: Adding a control for geolocation in Leaflet

texts:
In my Leaflet app, I am trying to add a simple geolocator control on my map. I am following the instructions from https://github.com/domoritz/leaflet-locatecontrol. I am already using Leaflet control (L.Control.locate). However, it does not show up on the map. I've tried adding it in different places in my app.js file. If I simply add it as var locate = L.control.locate().addTo(map); it keeps my point icons from showing up. Currently I have it as below. My point icons show up and the Layer controller and Zoom controller show up, but the Locate controller still does not show up.
Any suggestions? Side-Note: I know that the geolocation might not work at this point with the current setview, bounds, and max/min zoom settings (if the user is outside them). I just want the control to show up in the map and then I will try to add an alert if the user tries launching it outside the bounds/zoom limits.

var map = L.map('map').setView([35.7787, -78.6397], 14.50);
map.options.minZoom = 14;
map.options.maxZoom = 18;
map.setMaxBounds([[35.82389, -78.40], [35.704567, -78.80]]);

var layers = L.tileLayer('http://{s}.tile.openstreetmap.se/hydda/full/{z}/{x}/{y}.png', {attribution: '&copy; OpenStreetMapcontributors'}).addTo(map);

//icon code goes here...
var types = ['beer','bar','restaurant','music','culture','cafe', 'occult shop'];

var layerControl = L.control.layers().addTo(map);

types.forEach(function(type) {
    var layer = L.geoJson(dtr_points, {
        filter: function(feature, layer) {
            return feature.properties.Type == type;
        },
        onEachFeature: function (feature, layer) {
            var link_url = "" + feature.properties.Name + "";
            layer.bindPopup(link_url);

            if(feature.properties.Type == "beer") {
            layer.setIcon(beer);
            };
           //More icon code...
        }
    }).addTo(map);

    layerControl.addOverlay(layer, type);
});

var locator = L.control.locate({
    position: 'topright',
    strings: {
        title: "Show me where I am, yo!"
    }
}).addTo(map);


texts: Finding crossing point between geodesic line and parallel?

texts:
Given a geodesic line (starting point and initial azimuth) is there an algorithm to find the longitude of the [first] crossing of a given parallel?

I have already checked with Clairaut's constant that the geodesic will reach said parallel, so the crossing exists.

Or should I iterate the direct problem with different lengths until I get the endpoint on the parallel?


texts: Spatio-temporal block kriging with R package gstat?

texts: In the documentation of sp I saw that the coords argument of the SpatialPolygons function only accepts 2 columns. Is it not possible to create 3-dimensional polygons? I am doing spatio-temporal kriging using package gstat by adding time as a third spatial dimension of a SpatialPointsDataFrame. I would like to extend that to block kriging with polygons as blocks. Is there any way to do that?


texts: Is it possible to change the plugin folder location?

texts:
I am teaching a GIS course in a local collage using QGIS (qgis4w).
The way our lab pc's work is that they delete all added material every start up from the C:\ drive. They also have a D:\ drive where data can be saved.
My question is, is there a way to enable plugin installation in this configuration. I thought of something like moving the plugin definition file and the plugin folder to D:\


texts: Resetting system index of feature collection in Google Earth Engine?

texts:
I am using the sample to repeatedly sample from an image and because I am sampling many times and merging the collections, the ids are coming like 1_1_1_1_1_2_8. I am pretty new to EE so I am not sure if I am doing it the optimal way.

How can I proceed?

SAMPLE = ee.FeatureCollection(PSU_CP.iterate(function(cur,coll){
      img = image.clip(cur.geometry());
      fc = img.sample({numPixels: sSample , scale: 30, seed:seed_key, geometries:true})
    return ee.FeatureCollection(coll).merge(fc)
    }, ee.FeatureCollection([first]) ));


texts: Survey123 info file not reflecting changes

texts:
I am using Survey123 for collecting data. I am using Survey123 Connect for designing form. Now i have requirement to use one of our own map in survey123 for collecting geopoint.

After doing some research i came to know that custom maps can be added in tile map form. And to do so you have to modify .info file.

I did that as asked the problem is when i do a change and press save nothing happens of survey123 connect rather the info files gets back to its original state.


texts:
I have a query in Oracle to select polygons from my_polygon_table that interact or reside inside of a single polygon residing in border_table:

Select t1.*, t2.*, sdo_geom.sdo_intersection(t1.geoloc, t2.geom, 2 ) as genericgeo From      ( sdo_join ('my_border_table', 'geoloc', 'my_polygon_table', 'geom', 'mask=ANYINTERACT' ) ) j,
 my_border_table t1,my_polygon_table t2
 Where
 (j.rowId1 = t1.rowId and j.rowId2 = t2.rowId) and (t1.geoloc.sdo_gtype = 2003 or t1.geoloc.sdo_gtype = 2007) and (t2.geom.sdo_gtype = 2003 or t2.geom.sdo_gtype = 2007)

So far I've installed geometry with CREATE EXTENSION postgis; and copied tables. When I run a query in pgAdmin I get

ERROR:  syntax error at or near "table"
LINE 2: ...ntersection(t1.geoloc, t2.geom, 2 ) as genericgeo From table
                                                                  ^
SQL state: 42601
Character: 129

How to do the same thing with PostgreSQL?


texts: Adding offset to overlapping lines in qgis in same layer?

texts:
I have ONE layer of linestrings and some of them are overlapping, I need to display the overlapping lines with an offset.

I have seen the solution using the offset in style for the whole layer but it applies it to all the lines.

Is there a way to do it either with qgis or programatically using pyqgis ?


texts: GRASS GIS - V.VOL.RST avoid segmentation (paneling)

texts:
In order to reduce the segmentation effect where the interpolation appears to have different square regions, is this due to segment that are too small (segmax?).

Is the solution to avoid this to increase the segment size (which increases the run time)?

What should be the ratio of npmin and mpmax vs segmax as these two parameters need to be greater than segmax?

The default setting for the module is 50 for segmax, 200 for npmin, and 700 for mpmax.


texts: How to use gdal_calc.py to process multiple bands

texts:
I have 2 rasters (raster1 and raster2) with 5 bands for each one.
I want to do some calculation using all bands and write the output to a file with 5 bands.

I try this command :

gdal_calc.py -A raster1 -B ratser2 --outfile=result.tiff --calc="A+B"

but it generates a wrong output file with only 1 band.


texts: How to download the MODIS data from R using MODISTools library

texts:
I want to download MODIS product MOD13Q1 to use the Scientific Data Sets (SDS) layers of "250m_16_days_EVI", "250m_16_days_pixel_reliability". MODISTools library in R software provides the capability to download the MODIS data from LP DAAC data center.

I have read the MODISTools pdf guide which provides the basic understanding of how to use the library to subset the MOD13Q1 data layers by using the points location.
The following subset function which I want to use:

library(`MODISTools`)
MODISSubsets(LoadDat = modis.subset, Products = "MOD13Q1",
Bands = c("250m_16_days_EVI", "250m_16_days_pixel_reliability"),
Size = c(1,1))

Here LoadDat = modis.subset is the data frame of latitude and longitude points in decimal format and according to Size = c(1,1) parameter that is used to extract the 1 Km left and 1 Km right from the focal point as extent. But I want to used the extent of whole study area like

lat <- c(23.6947, 37.0894)
long <- c(60.8786,  79.3074)
period <- data.frame(lat=lat,long=long, start.date=2013,end.date=2014,id=1)
MODISSubsets(LoadDat = period, Products = "MOD13Q1", Bands = "250m_16_days_EVI", "250m_16_days_pixel_reliability", Size = c(0,0), TimeSeriesLength = 1)

But here I am confused about using Size parameter that is used to extract the point values. How can I use it to extract all the available data in this extent?
And second one is can I use this point information in ASCII file to display a raster in R software? I don't want to download the full tile which contains 250MB and in my case this will be 250GB data if I download all the available data starting from 2000 to 2015.


texts: How to extract EVI from MOD13Q1 with MRTBatch?

texts:
I want to extract EVI from Modis MOD13Q1 from some scene of images. I use ModisTool and MRTBatch. I have selected ONLY 250m_16_days_EVI as selected bands. All is correct until I save the parameter file.

The problem is: MRTBatch extract NDVI band, instead of EVI, as I selected.

Anybody knows the problem?

Here is the parameter file:

#The "INPUT_FILENAMES" field would be commented. If you want to load multiple input files please uncomment the "INPUT_FILENAMES" field and comment the"INPUT_FILENAME" field.
#Also the "ORIG_SPECTRAL_SUBSET" field needs to be uncommented and changed to "SPECTRAL_SUBSET". The initial "SPECTRAL_SUBSET" field should be deleted.

#INPUT_FILENAMES = ( E:\MAP_DATABASE\RASTER\MODIS\MOD13Q1_H30V08_H30V09_H31V09_2010-2014\MOD13Q1.A2010001.h30v08.006.2015198101220.hdf E:\MAP_DATABASE\RASTER\MODIS\MOD13Q1_H30V08_H30V09_H31V09_2010-2014\MOD13Q1.A2010001.h30v09.006.2015198101350.hdf E:\MAP_DATABASE\RASTER\MODIS\MOD13Q1_H30V08_H30V09_H31V09_2010-2014\MOD13Q1.A2010001.h31v08.006.2015198101343.hdf E:\MAP_DATABASE\RASTER\MODIS\MOD13Q1_H30V08_H30V09_H31V09_2010-2014\MOD13Q1.A2010001.h31v09.006.2015198100925.hdf )

INPUT_FILENAME = E:\EVI-batch\TmpMosaic.hdf

SPECTRAL_SUBSET = ( 1 )
#ORIG_SPECTRAL_SUBSET = ( 0 1 0 0 0 0 0 0 0 0 0 0 )

SPATIAL_SUBSET_TYPE = INPUT_LAT_LONG

SPATIAL_SUBSET_UL_CORNER = ( 9.999999999 121.851193415 )
SPATIAL_SUBSET_LR_CORNER = ( -9.999999999 142.159725651 )

OUTPUT_FILENAME = E:\EVI-batch\MOD13Q1.A2010001.h30v08.006.2015198101220.tif

RESAMPLING_TYPE = NEAREST_NEIGHBOR

OUTPUT_PROJECTION_TYPE = GEO

OUTPUT_PROJECTION_PARAMETERS = (
 0.0 0.0 0.0
 0.0 0.0 0.0
 0.0 0.0 0.0
 0.0 0.0 0.0
 0.0 0.0 0.0 )

DATUM = WGS84


texts: ArcGIS - OD matrix: Quality of results with a fragmented rail network

texts:
I have a question about the results of an OD-matrix. I have data of harbours (destination and origins) and want to determine the distance between them. I will use a road network. It seems that not all railroads are connected. However I still get a distance between these specific points.

Does anyone know how they are measured and if this is correct. The distance does not correspond with euclidean distance, so that is not what is happening if the tool can not find a rail connection.


texts:
I have a shapefile that show this description under Geometry type of the features in this layer:

Point (WKB type: "PointZ")

I need to use is it in some plugins that seems to work only with points:

Point (WKB type: "Point")

I know how to access the Z value,  z($geometry ), but how can I convert the geometry type?


texts: What is the target .NET Framework of ESRI ArcObjects SDK 10.2

texts:
I had ArcMap 10 Add-In project built with ArcObjects SDK for .NET 10.1 where I could not reference external COM DLL build on .NET 4.0 due to a .NET forward compatibility issue. (ArcObjects SDK 10.1 uses .NET 3.5 SP1)

I wonder if the newer ArcObjects SDK 10.2 release targets the .NET 4.0 framework.
That will solve my problem.

Does anyone know?
Thanks


texts: Arcpy batch raster reclassify script unable to save

texts:
I have a folder with  multiple rasters that need to be reclassified. Using spatial analyst reclassify tool, I have the following script:

import arcpy
import pandas as pd
arcpy.CheckOutExtension("Spatial")
from arcpy import env
from arcpy.sa import *
arcpy.env.overwriteOutput = True

# Set environment settings
env.workspace = r"C:\Users\\Documents\ND_files\NASS_CDL"

nrasterList = arcpy.ListDatasets("*","Raster")
     for targetRaster in nrasterList:
     inRaster = targetRaster
     remap = RemapValue(remap1)
     #remap1 is a list of lists ex. [[0,1],[1,3],[44,3]]
     outReclassify = Reclassify(inRaster, "Value", remap, "NODATA")
     # Value is the name of the column to be reclassified
     outReclassify.save(r"C:\Users\Documents\ND_files\NASS_CDL\reclass")

the error I receive:

RuntimeError: ERROR 010240: Could not save raster dataset to C:\Users\Documents\ND_files\NASS_CDL\reclass with output format GRID.

I read somewhere else that people were having issues saving the NASS Crop Data  Layer and working with  GRID format.


texts: In what real-world examples do users still digitize features from a map?

texts:
Is digitizing of map features completely dead?

For example, when and why do you still digitize features from a map?

It can be heads digitizing up from a scanned paper map (or image, or whatever) or using a traditional digitizing tablet. It can also just be digital and heads-up (i.e. you never had the paper map).

Do you ever still utilize manual feature digitization to create points, lines, or polygons?

I need this information for a statistical technique I am developing.

UPDATE:

Any VGIS examples?
Any Archaeology examples?


texts:
The python buffer by percentage plugin produces error everytime I try to run it in QGIS 2.16.1 on Windows OS. Here is the full error description:

An error has occurred while executing Python code:

IOError: [Errno 9] Bad file descriptor
Traceback (most recent call last):
File "C:/Users/GenotypePhenotype/.qgis2/python/plugins\BufferByPercentage\bufferbypercentage.py", line 184, in run
    self.target_factor, self.segments)
  File "C:/Users/GenPhen/.qgis2/python/plugins\BufferByPercentage\bufferbypercentage.py", line 57, in find_buffer_length
    area_unscaled, target_factor)
  File "C:/Users/Genotype-Phenotype/.qgis2/python/plugins\BufferByPercentage\bufferbypercentage.py", line 78, in secant
    oldf, f = func(oldx, *args), func(x, *args)
  File "C:/Users/Genotype-Phenotype/.qgis2/python/plugins\BufferByPercentage\bufferbypercentage.py", line 67, in calculateError
    print area_scaled
IOError: [Errno 9] Bad file descriptor

Python version: 2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]
QGIS version: 2.16.1-Nødebo Nødebo, 8545b3b

Python Path:
C:/PROGRA~1/QGIS21~1.1/apps/qgis/./python/plugins\processing
C:/Users/Genotype-Phenotype/.qgis2/python/plugins\LecoS
C:/PROGRA~1/QGIS21~1.1/apps/qgis/./python
C:/Users/Genotype-Phenotype/.qgis2/python
C:/Users/Genotype-Phenotype/.qgis2/python/plugins
C:/PROGRA~1/QGIS21~1.1/apps/qgis/./python/plugins
C:\PROGRA~1\QGIS21~1.1\apps\Python27\lib\site-packages\matplotlib-1.3.1-py2.7-win-amd64.egg
C:\PROGRA~1\QGIS21~1.1\apps\Python27\lib\site-packages\nose-1.3.3-py2.7.egg
C:\PROGRA~1\QGIS21~1.1\apps\Python27\lib\site-packages\tornado-4.0.1-py2.7-win-amd64.egg
C:\PROGRA~1\QGIS21~1.1\apps\Python27\lib\site-packages\backports.ssl_match_hostname-3.4.0.2-py2.7.egg
C:\PROGRA~1\QGIS21~1.1\apps\Python27\lib\site-packages\certifi-14.05.14-py2.7.egg
C:\PROGRA~1\QGIS21~1.1\bin\python27.zip
C:\PROGRA~1\QGIS21~1.1\apps\Python27\DLLs
C:\PROGRA~1\QGIS21~1.1\apps\Python27\lib
C:\PROGRA~1\QGIS21~1.1\apps\Python27\lib\plat-win
C:\PROGRA~1\QGIS21~1.1\apps\Python27\lib\lib-tk
C:\PROGRA~1\QGIS21~1.1\bin
C:\PROGRA~1\QGIS21~1.1\apps\Python27
C:\PROGRA~1\QGIS21~1.1\apps\Python27\lib\site-packages
C:\PROGRA~1\QGIS21~1.1\apps\Python27\lib\site-packages\PIL
C:\PROGRA~1\QGIS21~1.1\apps\Python27\lib\site-packages\jinja2-2.7.2-py2.7.egg
C:\PROGRA~1\QGIS21~1.1\apps\Python27\lib\site-packages\markupsafe-0.23-py2.7-win-amd64.egg
C:\PROGRA~1\QGIS21~1.1\apps\Python27\lib\site-packages\pytz-2012j-py2.7.egg
C:\PROGRA~1\QGIS21~1.1\apps\Python27\lib\site-packages\win32
C:\PROGRA~1\QGIS21~1.1\apps\Python27\lib\site-packages\win32\lib
C:\PROGRA~1\QGIS21~1.1\apps\Python27\lib\site-packages\Pythonwin
C:\PROGRA~1\QGIS21~1.1\apps\Python27\lib\site-packages\Shapely-1.2.18-py2.7-win-amd64.egg
C:\PROGRA~1\QGIS21~1.1\apps\Python27\lib\site-packages\wx-2.8-msw-unicode
C:\PROGRA~1\QGIS21~1.1\apps\Python27\lib\site-packages\xlrd-0.9.2-py2.7.egg
C:\PROGRA~1\QGIS21~1.1\apps\Python27\lib\site-packages\xlwt-0.7.5-py2.7.egg
C:/Users/Genotype-Phenotype/.qgis2//python


texts: Geoserver cql_filter on wms json/vector requests

texts:
I'm trying to make a cql filtered wms request for a non image based format like geojson or mapbox vector but the filters aren't applying. If i change the format to image/png the filters behave as expected.

Is this functionality not supported by geoserver? We're running version 2.11.

Working request (filtered)

{{host}}/ows?bbox=16115371.547420282,-4546474.442402283,16115983.043646563,-4545862.946176002&format=image/png&service=WMS&version=1.1.1&request=GetMap&srs=EPSG:3857&transparent=true&width=512&height=512&layers=LINE_NB&CQL_FILTER=SEASON_KEY='BASE_2018_SUMMER'

Failing request (not filtered)

{{host}}/ows?bbox=16115371.547420282,-4546474.442402283,16115983.043646563,-4545862.946176002&format=application/x-protobuf;type=mapbox-vector&service=WMS&version=1.1.1&request=GetMap&srs=EPSG:3857&transparent=true&width=512&height=512&layers=LINE_NB&CQL_FILTER=SEASON_KEY='BASE_2018_SUMMER'

URL encoding the query string makes no difference to the result.


texts: What's the way to run a PostGIS function in PyQgis

texts:
I'm trying to use PostGIS a function in a PyQgis script.

My code is:

qsql="(select gid, id, St_Buffer(the_geom,200) as geom from archi)"
uri = QgsDataSourceURI()
uri.setConnection("ip", "port", "x", "x", "x")
uri.setDataSource("",qsql,  "the_geom","","gid")
layer = QgsVectorLayer(uri.uri(), "test_buffer", "postgres")
if not layer.isValid:
     print 'error'
else:
     QgsMapLayerRegistry.instance().addMapLayer(layer)

The vector layer is loaded but there is NO buffer.

What's the way to run PostGIS function in PyQgis?


texts:
I was trying the following line of codes in python (that I learnt from a tutorial),

import folium
import pandas

data=pandas.read_csv("Volcanoes.txt")
lat = list(data["LAT"])
lon = list(data["LON"])
elev = list(data["ELEV"])

def color_producer(elevation):
    if elevation < 1000:
        return 'green'
    elif 1000 <= elevation < 3000:
        return 'orange'
    else:
        return 'red'
map = folium.Map(location=[38.58, -99.09], zoom_start=6, tiles="Mapbox Bright")
fg = folium.FeatureGroup(name="My Map")
for lt, ln, el in zip(lat, lon, elev):
    fg.add_child(folium.CircleMarker(location=[lt, ln], radius = 6, popup=str(el)+" m",
    fill=True, fill_color=color_producer(el), color = 'grey', fill_opacity=0.7))

fg.add_child(folium.GeoJson(data=(open('world.json', 'r', encoding='utf-8-sig'))))

map.add_child(fg)
map.save("Map1.html")

I receive an error saying

Traceback (most recent call last):
  File "map1.py", line 23, in <module>
    fg.add_child(folium.GeoJson(data=(open('world.json', 'r', encoding='
utf-8-sig'))))
  File "C:\Users\TOSHIBA\AppData\Local\Programs\Python\Python36-32\lib\s
ite-packages\folium\features.py", line 495, in __init__
    raise ValueError('Unhandled object {!r}.'.format(data))
ValueError: Unhandled object <_io.TextIOWrapper name='world.json' mode='
r' encoding='utf-8-sig'>.

I typed exactly the code what I learnt from tutorial course, but seems the code not working for me. Please help me identify the problem ! I was trying it on Atom 1.21.0 or is there any changes made in folium ?


texts: Openlayers 3 attribution - source/ credit

texts:
I can't get attribution (source/ credit) to show up for WMS layers added to and OL3 map.
The attribution for base layers (OSM etc) works fine but when I try and use the attribution parameter to add attribution to my WMS layers no attribution shows up - any ideas?

        }),
        title:"Special Protection Areas",
        opacity:"0.4",
        attribution: "where it came from"
      })


texts: Looping over 5D masked Numpy array Python

texts:
I am a beginner in Python and I faced a problem in maintaining the same size of the last 2 dimensions in the 5D array that I have while looping over it to do some calculations, as they represent the latitude and longitude of the map that I will plot as result from this calculation. The calculation that I want to apply is in a function called (line), as I call this function inside the loop to do the calculation on the (wind_data) 5D array. But I have an error (ValueError: setting an array element with a sequence.)
import numpy as np

Swept_Area=5026
Air_Density= 1.225

def line(v):
    prob=np.empty(shape=v.shape)
    idx_low=np.where(v<=4)
    idx_mid=np.where((v >= 4) & (v <= 15))
    idx_mid2=np.where((v > 15) & (v <= 25))
    idx_high=np.where(v>25)

    m=  0.5 * Swept_Area * Air_Density *0.29 *  v**3   * 10**-3
    prob[idx_low]=0
    prob[idx_mid]=m[(v >= 4) & (v <= 15)]
    prob[idx_mid2]=3000
    prob[idx_high]=0

    return prob

prob=np.empty(shape=(wind_data.shape[3],wind_data.shape[4]))

for lat_index in range(wind_data.shape[3]):
   for long_index in range(wind_data.shape[4]):
        prob[lat_index,long_index]=line(wind_data[:,:,:,lat_index,long_index])


texts: Running PyQGIS Scripts in Las Palmas from Models created in Dufour?

texts:
I am programming with PyQQIS using QGIS DuFour(old version) and QGIS Las Palmas(latest version) and Python 2.7 on my system.

On QGIS DuFour I can export my model to Python script, on QGIS Las Palmas I can't.

I have synchronize my Python system with PyQgis and QGIS Las Palmas so I can create Python QGIS script standalone OK with that.

Now I want to export all my model from QGIS DuFour to Python scripts and run these scripts with  PyQgis from QGIS Las Palmas.

How can I do that?


texts: Domain created, applied to feature class, but pick list not showing in table?

texts:
Domain was created, applied to a field in a feature class, but the pick list not showing in the table when I edit the feature class. In ArcGIS Desktop 10.2.

I am wondering what I have missed?


texts: Wrong shape display with North_Pole_Azimuthal_Equidistant projection

texts:
I need to create a map that displays the breeding and wintering ground of a bird species that has a circular arctic breeding distribution and winters in Africa. The map should lool pretty much like this:

However, the graticules are not correctly displayed. Depending on scale and frame one or two are missing.

I use the following projection: North_Pole_Azimuthal_Equidistant

Basic layers are graticules and land shapes from naturalearthdata.com. Is there maybe something wrong with the basic shapes and is there another free alternative?

I created grid-lines myself with Vector > Research tools > gridlines too and also tried the grid function from the MMQGIS package but but doesn´t help and gives the same problem.


texts:
I have a stand-alone PyQGIS application that can load vectors/rasters etc. The problem I am facing is when the application is launched, it has an empty QgsMapCanvas as map display, when it is empty, I mean has no layer features on the canvas, the canvas appears to have empty extent. If I do the following:

canvas.extent().isEmpty() (this is true)
canvas.extent().width() (this is 0.0)
canvas.extent().height() (this is 0.0)

At this point, if I put a point on the canvas, and call canvas.zoomToFullExtent(), then the canvasMoveEvent will stop working (stop capturing mouse move's coordinates)

So, my plan is to get the canvas's visible center and manually set its extent to the visible view port that is displayed in the application when it is first launched. But I couldn't find a way to do that or even if this is possible. Basically, something like canvas.extent().center() or something...


texts: FME Segment Line with To/From Distance Values from Excel?

texts:
I have a large road network with some basic attributes like ROAD NAME, ROAD TYPE, LENGTH. Each road has been segmented but these segments are wrong for the new schema. I have successfully (mostly!) joined these line segments in preparation for creating the new segments with new attributes/values.

What I really need help on is: how to segment these roads individually based on the To/From distance values I have in an excel file, which is the basis of the new schema for whole road network. Whole bunch of new attributes and their values here too.

Ive read up on the SNIPPER and MEASURE SETTER but unsure of how to proceed. Should I be SCHEMA MAPPING or ATTRIBUTE CREATING first prior to creating segments? I would have thought I would need the segments to be made first so the road segments could be matched to the schema where the segment length values are coming from?


texts: ST_Intersection: Intersection of all geometries in a table

texts:
Must be a fairly trivial task, but PostGIS gives me errors..

I have a set of geometries geom in a single table t1, and I would like to find intersection area of all of the rows foo where value is bar.

What I am doing is:

SELECT ST_INTERSECTION(geom) FROM t1 WHERE foo = 'bar' and ST_INTERSECTS(geom)

I understand examples where I have two geometries or two tables, but I can't figure out how to do that in a single table across rows. I have several thousands records to filter, so specifying geometries explicitly is not an option.

And yet I am not getting the results I expect from St_intersection (yeah, I assume not the method's fault):

If I apply the method as described above I will get an intersection between each and every row in the table. It's easier to see what I mean, check the picture below.

Instead of having a single geometry (polygon) which represents the shared area between all geometries I get multiple geometries. I am still wondering how to get the area shared by all geometries, not a set of intersections for all pairs in the table.


texts: Change Control MousePosition displayProjection

texts:
I added mouse control on the map with Lat/Long as display projection.
I added a dropbox on the map for the user to choose a coordinate system (3 choices).
Depending on the selection, I would like to change the displayed coordinates. What I'm missing is the syntax to do that.

This is my initial coordinate system:

map.addControl(new OpenLayers.Control.MousePosition({displayProjection: 'EPSG:4326'}));

If the choice is EPSG:32187, I tried to change it:

displayProjection: new OpenLayers.Projection("EPSG:32187") });

OpenLayers.Control.MousePosition({displayProjection: 'EPSG:32187'});

map.displayProjection = new OpenLayers.Projection("EPSG:32187");

map.MousePosition = new OpenLayers.Projection("EPSG:32187");

Nothing works.

The only thing it worked is if I add another MousePosition control:

map.addControl(new OpenLayers.Control.MousePosition({displayProjection: 'EPSG:32187'}));

But that they are on on top of the other and I would like to reuse the existing Mouse control.


texts: Crayfish 3.0.0 not working in QGIS 3.4.3

texts:
I tried to install Crayfish 3.0 in QGIS 3.4 and although I didn't get any errors during the installation of the plugin I don't have access to all Crayfish features. I can only use the 'plot' and 'export animations' features but I don't have the option to load xmdf files or .dat or use any other crayfish tools.

Any ideas how to solve this?


texts: How to modify the automatic GetLegendGraphic in QGIS when using QgsRasterLayer in the python interface to load data from a ncWMS server

texts:
tl;dr: How to add an extra parameter (using python) to the GetLegendGraphic request made by QGIS when a WMS layer is added using QgsRasterLayer?

UK Reading University developed a map server called ncWMS. It is basically a regular WMS server with some extend capabilities to serve meteorological and oceanographic data stored as NetCDF files (a format widely used in the metocean community). The extended capabilities are due to having to deal with the special characteristics of metocean data: time and depth/height dimensions, variable color scale, etc.

It's possible to see the server in action here. I managed to load a layer from this server in QGIS using the python console:

title = 'Test'
server = 'http://godiva.rdg.ac.uk/ncWMS2/wms'
layer = 'foam/U:V-group'
style = 'default-vector/default'
time = '2011-01-12T00:00:00.000Z'

finalurl = ('crs=EPSG:4326'
            '&srs=EPSG:4326'
            '&format=image/png'
            '&service=WMS'
            f'&layers={layer}'
            f'&styles={style}'
            f'&url={server}')

resultLayer = QgsRasterLayer(finalurl, f'{title} @ {time}', 'wms')

# HACK taken from Anita Graser's Time Manager:
# https://github.com/anitagraser/TimeManager/blob/master/raster/wmstlayer.py
resultLayer.dataProvider().setDataSourceUri('IgnoreGetMapUrl=1&'
    + resultLayer.dataProvider().dataSourceUri()
    + (f'?time={time}'
       '%26colorScaleRange=0.0%2C1.0'
       '%26numColorBands=250'
       '%26logScale=false'
       '%26bbox=-180,-90,180,90'))

print('url query: {}'.format(resultLayer.dataProvider().dataSourceUri()))

if resultLayer.isValid():
    QgsProject.instance().addMapLayer(resultLayer)
else:
    raise ValueError('Not a valid layer.')

resulting in

The only problem is that the extra parameter colorScaleRange that defines the minimum and maximum values for the color scale (in the example above 0 to 1 m/s for ocean currents speed) is respected by the REQUEST=GetMap but not by the REQUEST=GetLegendGraphic issued by QGIS when the layer is added. You can see that the legend goes from 0 to 2 m/s.

Checking the Tomcat log for a ncWMS server that I have installed on my computer (using different data and behind a proxy, so not public like the UK Reading University used in the example above)

127.0.0.1 - - [30/Jan/2020:11:58:36 -0200] "GET /ncWMS2/wms?time=2020-01-10T00:00:00.000Z&colorScaleRange=16.0,32.0&numColorBands=250&logScale=false&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&BBOX=-38.00499999999999545,-59.00499999999999545,9.004999999999999005,-24.99500000000000099&CRS=EPSG:4326&WIDTH=670&HEIGHT=925&LAYERS=nasa_mur_sst/sst&STYLES=default-scalar/default&FORMAT=image/png&DPI=96&MAP_RESOLUTION=96&FORMAT_OPTIONS=dpi:96&TRANSPARENT=TRUE HTTP/1.1" 200 160372
127.0.0.1 - - [30/Jan/2020:11:58:37 -0200] "GET /ncWMS2/wms?REQUEST=GetLegendGraphic&PALETTE=default&LAYERS=nasa_mur_sst/sst&STYLES=default-scalar/default&SERVICE=WMS&VERSION=1.3.0&SLD_VERSION=1.1.0&FORMAT=image/png&LAYER=nasa_mur_sst/sst&STYLE=default-scalar/default&TRANSPARENT=true HTTP/1.1" 200 2916

it is possible to see that in the GetMap request the parameter colorScaleRange=min,max is present, but in the GetLegendGraphic it is not.

I am sure ncWMS is capable of dealing with this parameter in the GetLegendGraphic as shown here.

So my question is just how (using Python) to insert this extra parameter in the legend request before QGIS makes the GetLegendGraphic request.


texts: Change default layout in GeoExplorer from MapQuest to something else

texts:
By default if you open a view or a table in geoxplorer from geoserver it loads your layer on top of the mapQuest layout. But mapQuest is not available anymore.

I tried looking in every config options and file i found i cannot find how to load another layout by default. I can change it manually no problem but not by default.

Can someone help ?


texts: How can I calculate the coordinates of a point on the intersection of two parallel lines with a bisecting line with pyQGIS?

texts:
I am building some sort of homemade pyQGIS-Buffer-tool as the existing flat-end-buffer-tools are not working for me ( see v.distance.buffer does not create proper buffers with flat endcaps)

Therefore I am looking for a solution how to get the polygon-vertices of my buffer-polygon where the linesegments have an angle <>180° to each other (see image below).

For start- and endpoint I just make a linear interpolation with 90° angle to the line direction.

Now i wonder how I could get the two points that are marked red in the picture?
( The orange line is my source line I would like to "buffer".)
Note: I am not looking for a flat buffer tool but just for a pythonic solution how to get the vertexpoints that are marked in the picture.
My first idea was to make parallel lines to my source line, elongate them if nesccessary and find the intersection but this seems a bit cumbersome.

Any ideas how to solve this with pyQGIS or some math.xxx magic?

Here a picture how it looks like at the moment if I just interpolate every vertex 5m with -90° to the line direction:

and here a picture how it should look like:


texts: GDALRasterizeGeometries doesn't burn correct Values

texts:
What I have is a vector shapefile with a string field. What I've done is written a CSV file with the field ID, field_value and now I would like to rasterize the vector layer as a byte raster with the raster pixel value corresponding to the field ID.

I'm not sure what's going wrong but I get a blank raster every single time.

Here's how I'm trying to do it:

OGRFeature * ogrFeat;
std::vector<OGRGeometryH> ogrBurnGeometries;
std::vector<double> dBurnValues;

OGRGeometryH hPolygon;

// Create a list of burn-in values
poLayer->ResetReading();
while( (ogrFeat = poLayer->GetNextFeature() ) != NULL ){

    OGRGeometry * ogrGeom;
    // No geometry found. Move along.
    if( ogrFeat->GetGeometryRef() == NULL )
    {
        delete ogrFeat;
        continue;
    }
    ogrGeom = ogrFeat->GetGeometryRef();
    ogrBurnGeometries.push_back( (OGRGeometryH) ogrGeom->clone() );

    if (fieldType == OFTString){
        dBurnValues.push_back( ogrFeat->GetFID() );
    }
    else {
        dBurnValues.push_back( ogrFeat->GetFieldAsDouble(psFieldName) );
    }
    delete ogrFeat;
}
int band = 1;
char **papszRasterizeOptions = NULL;
papszRasterizeOptions =
    CSLSetNameValue( papszRasterizeOptions, "ALL_TOUCHED", "TRUE" );

CPLErr err = GDALRasterizeGeometries( pDSOutput, 1, &band,
                                      1,
                                      &(ogrBurnGeometries[0]),
                                      NULL, NULL,
                                      &(dBurnValues[0]),
                                      NULL,
                                      NULL, NULL );
ogrBurnGeometries.clear();
dBurnValues.clear();

EDIT: I've fixed a memory problem as suggested by Kyle in the comments. I'm now getting an empty raster with the right extents but the geometries are still not being written.

The CPLErr is empty and the ogrBurnGeometries and dBurnValues each contain the right values. I am confused.


texts: Does anyone know of an server implementation of the WFS 2.0.2 standard

texts:
I'm looking for either a public server than implements the WFS 2.0.2 standard or a reference implementation that I could run on my own copy of ArcGIS server.

I've googled for a long time and searched data.gov but found nothing but 2.0 implementations.


texts:
I need to create a monthly NDVI time series from GIMMS in GEE (i.e. from January 1981 to December 2013) from the sum of fortnight values (or mean).
I've been trying to create a function to do so, but can't figure it out.

I try this code, but generates annual sums, not monthly sums...

var monthly =
mageCollection.filter(ee.Filter.calendarRange(1981,2013,'year'))
.filter(ee.Filter.calendarRange(1,12,'month')).sum();


texts:
I have Spline plugin installed. It is activated. When I try to draw a curved line it doesn't change into a curve. The lines remain the way I have drawn them. I am using QIS 2.10.1 with Windows 8 or 7.
I also get this Warning, but I don't understand what does it mean and how to deal with it:

2015-10-20T16:22:15 1   warning:C:/Users/katy/.qgis2/python/plugins\spline\spline.py:116: DeprecationWarning: QgsMapCanvas.mapRenderer() is deprecated
              if self.canvas.mapRenderer().hasCrsTransformEnabled() and layer.crs() != self.canvas.mapRenderer().destinationCrs():

        traceback:  File "C:/Users/katy/.qgis2/python/plugins\spline\spline.py", line 96, in canvasPressEvent
            self.createFeature()

Please. Help me with Spline or suggest another way to digitize curved lines.


texts:
I'm looking for a good way to learn more ArcGIS techniques than I came across in my undergraduate Geography degree to really use the power of the software.

Browsing around and asking questions on this site has helped quite a bit e.g. I did not realised that you could link dataframes together, or even that you could have multiple data frames ([see How to provide automatic overview map in ArcGIS Desktop?).

Similarly, I suspect there are lots of ArcGIS features that could really make my life easier, but that I just don't know about e.g. geodatabases.

What I'm really looking for are tutorial recommendations, in any medium (books, videos, websites, whatever). Obviously free or cheap are better, but I do have a budget for books as part of my PhD. As I probably hinted at above, I'm more interested in general usage of ArcGIS rather than specific tools (like DEM processing or route-finding), and something that is relatively structured would be good, so that I can work through it in a sensible order.

I know this is a bit of a tall order, but any ideas?


texts: Cropping an image while maintaining the coordinates in Python

texts:
I have a 5120 X 5120 TIFF image, and I want to crop it to 25 small images of 1024 X 1024. How can I do it while maintaining the image coordinates of 5120 X 5120 TIFF image into my small cropped images using Python?
How can I incorporate the coordinate info of the original image into my cropped images, so I get the same coordinates when checking them in QGIS?


texts: How to add images of buildings inside a map?

texts:
I'm going to do a trip and I want to have a custom map on my website. My idea is to have a map similar to this one for example.

So I need to add images (svg or jpg) to the base map. I want a clean base map, so I don't like using Google Maps for this, that's why I'm trying a more custom approach (apart from learning new things).

I know nothing about GIS yet, but as I developer I'm eager to learn. I've seen two options:

With a gis editing program, import a base layer for the area of my trip, edit the layer and add those svgs images. Is it possible to do this with QGIS and how? I fear that it can take too much work.
Use something like OpenLayers and add those building images as overlays. I wonder if I can position it well to hide what's behind (I think with the position and positioning parameters it could work).

Which one do you think is a better alternative? Any other option? I can spend some time, but I want to go on the right direction.


texts: How can I dynamically add features to a GeoJSON file?

texts:
I have tons of pictures which I am showing in markers on a LeafletJs map and currently have a fixed GeoJSON file. So I want to know is there any easy way to dynamically add new features to the GeoJSON file?

I was thinking a script(Php/JS) which scans a folder for new images, then gets the coordinates from a GPX file based on the image name and then add those to the GeoJSON feature collection. Is this kind of solution possible?


texts: Versioning in PostGIS on PostgreSQL database

texts:
I have a cadastral database which requires multi-user editing/updating on daily basis.

I am trying to migrate from Esri enterprise geodatabase to PostGIS/PostgreSQL. I am not sure PostGIS offers versioning option similar to the Esri versioning.


texts: Select by Attributes query ArcMap 10.5.1

texts: I am trying to find capitals with a population between 200,000 and 300,000 and I do not know how to build the query correctly. Here is what I have, but it doesn't work:


texts: Getting Waypoints between two coordinates from CloudMade

texts:
Suppose I have A (source point) and B (destination point), then I want to get a series of route point or waypoints between A to B from CloudMade. The result that I need is something like this:
(LonA, LatA) (Lon1,Lat1) (Lon2,Lat2) (Lon3, Lat3)................(Lon_n, Lat_n) (LonB, LatB)

Are there any available methods or tools to do that?


texts: Returning percentage of area of polygon intersecting another polygon using shapely

texts: Using the shapely library, suppose I have a shp file V and another shp file Z. For each polygon v in V, and for each polygon z in Z, I would like to know what percentage of the area of z falls within (intersects) v.


texts:
I search a plugin or another good way to generate different maps, with the same layout. With this plugin a user should be able to create a map from his data. I set the parameters for the map, so every map should have the same layout, the same position and style for the legend etc.

Is there maybe any plugin with this functions or an other possibilty?