texts:
I am using QGIS Version 2.14.3.
I have a map with polygons and a lot of different attributes. Some polygons have NULL values in some columns.
I've tried to color the polygons with properties->style->graduated. But in some columns half of the values are NULL and now they don't have any colour.

I know that i could color the NULL values when I chose "categorized" coloring, but I have way to many different entries to do that, because QGIS chooses a colour for every different entry.

I need to have intervals for coloring. Is there an easy way to do that?


texts: How do I zoom to the overall extent of two (or more) selected features from two layers?

texts:
I am creating an application where the user picks a value from a combo box, the application then selects that feature and then selects the corresponding features from a second layer, (both selections are in a different colour). What I want to do is as soon as the selection is made, the map will zoom to the envelope of the selected features. I know how to zoom to 2 individual layers (using pEnvlope.union...), but i need to do it on features.

Has anyone got any idea of the best way to do this?
(I'm using ArcGIS 9.3.1 and VBA).

Thanks


texts: Conducting linear correlation based on stacked raster in R?

texts:
I am analyzing average annual GPP during 2000-2014 of a area. I want to conduct linear correlation analysis between annual GPP and time, and I also want to get the slope and correlation coefficient values of the linear relationship of each cell and output the results as rasters.

Do you know how to calculate the slope and R^2 values in "raster" package in R?

Can calc() or lm () will make it?

I list the data as following:

library(raster)
GPP2000=raster("GPP2000.tif")
GPP2001=raster("GPP2001.tif")
GPP2002=raster("GPP2002.tif")
GPP2003=raster("GPP2003.tif")
GPP2004=raster("GPP2004.tif")
GPP2005=raster("GPP2005.tif")
GPP2006=raster("GPP2006.tif")
GPP2007=raster("GPP2007.tif")
GPP2008=raster("GPP2008.tif")
GPP2009=raster("GPP2009.tif")
GPP2010=raster("GPP2010.tif")
GPP2011=raster("GPP2011.tif")
GPP2012=raster("GPP2012.tif")
GPP2013=raster("GPP2013.tif")
GPP2014=raster("GPP2014.tif")
GPP=stack(GPP2000,GPP2001,GPP2002,GPP20003,GPP2004,GPP2005,GPP2006,GPP2007,GPP2008,GPP20009,GPP2010,GPP2011,GPP2012,GPP2013,GPP2014)


texts: Extracting data from four intersecting polygon layers in PostGIS?

texts:
I've been using the following code in PostgreSQL/PostGIS for the last year, as needed, which is about once a quarter. It works OK, it's just that it takes days to run. With over 25,000 parcels making up the county, maybe that's fine, but I'd like some expert input to see if there's a better way to do it.

An additional consideration, wetlandyn is a multi-polygon, the rest are single, though soils and luc could be converted to a multi-polygon.

Is querying against multi-polygons any more or less efficient than singles?

SELECT
    parcels.pin AS "PIN",
    soils.musym AS "Soils",
    parcels.camacls as "CAMACLS",
    luc.landcode as "Land Code",
    wetlandyn.inwetland as "In Wetland (y/n)",
    ROUND(((SUM(ST_AREA(ST_INTERSECTION((ST_INTERSECTION(parcels.geom,luc.geom)),
    (ST_INTERSECTION(wetlandyn.geom,soils.geom)))))/ST_AREA(parcels.geom))*parcels.acres)::numeric,4) AS "Legal Soils"

FROM
    soils,
    luc,
    parcels,
    wetlandyn

WHERE
    ST_INTERSECTS(parcels.geom, soils.geom)
    AND
    ST_INTERSECTS(parcels.geom, luc.geom)
    AND
    ST_INTERSECTS(parcels.geom, wetlandyn.geom)
    AND
    CAST(parcels.camacls AS int) >=100
    AND
    CAST(parcels.camacls AS int) <200
    AND
    parcels.acres > 0
    AND
    ((ST_AREA(ST_INTERSECTION(ST_INTERSECTION(ST_INTERSECTION(parcels.geom,luc.geom),
    ST_INTERSECTION(parcels.geom,soils.geom)),ST_INTERSECTION(parcels.geom,wetlandyn.geom)))/ST_AREA(parcels.geom))*parcels.acres)
    >= 0.00005

GROUP BY
    parcels.pin,
    soils.musym,
    parcels.camacls,
    luc.landcode,
    parcels.acres,
    parcels.geom,
    luc.geom,
    soils.geom,
    wetlandyn.geom,
    wetlandyn.inwetland

ORDER BY
    parcels.pin, soils.musym, luc.landcode


texts:
I use Quantum GIS 1.8.

I have a polyline layer. I can drag&drop polyline nodes by using "Node Tool".
Is it possible to drag&drop line segments of a polyline in order to preserve line segments' direction?


texts: Filtering out wastelands,scrubs,swamps, grasslands etc from LULC for India using ERDAS IMAGINE Model maker?

texts:
I have a LULC for India. i need to filter out non-agricultural, non-builtup, non-waterbody, non-forest. simply wastelands,scrubs,swamps, grasslands etc.

Please how to filter these using ERDAS IMAGINE Model maker?


texts: Problem with export vector layer by PyQGIS in QGIS

texts:
I found problem with exporting vector layer with "Replace all selected raw field values by displayed values" in PyQGIS script. Source layer contains value map for attribute "value_map":

PyQGIS script:
# -*- coding: utf-8 -*-
from qgis.PyQt.QtCore import QCoreApplication
from qgis.core import (QgsProcessing,
                       QgsFeatureSink,
                       QgsVectorFileWriter,
                       QgsProcessingException,
                       QgsProcessingAlgorithm,
                       QgsVectorFileWriter,
                       QgsProcessingParameterFeatureSource,
                       QgsProcessingParameterFeatureSink)
from qgis import processing
from qgis.utils import iface

class ExampleProcessingAlgorithm(QgsProcessingAlgorithm):
    INPUT = 'INPUT'

    def tr(self, string):
        return QCoreApplication.translate('Processing', string)

    def createInstance(self):
        return ExampleProcessingAlgorithm()

    def name(self):
        return 'myscript'

    def displayName(self):
        return self.tr('My Script')

    def group(self):
        return self.tr('Example scripts')

    def groupId(self):
        return 'examplescripts'

    def shortHelpString(self):
        return self.tr("Example algorithm short description")

    def initAlgorithm(self, config=None):
        # We add the input vector features source. It can have any kind of
        # geometry.
        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.INPUT,
                self.tr('Input layer'),
                [QgsProcessing.TypeVectorAnyGeometry]
            )
        )

    def processAlgorithm(self, parameters, context, feedback):
        layer = iface.activeLayer()
        feedback.pushInfo(str(layer))
        QgsVectorFileWriter.writeAsVectorFormat(
            layer=layer,
            fileName='C://Users//001-PC//Desktop//output_layer.json',
            fileEncoding="utf-8",
            driverName="GEOJSON",
            fieldValueConverter=QgsVectorFileWriter.FieldValueConverter()
        )
        return {}

How can i export layer values with values map? My script does not working with:

fieldValueConverter=QgsVectorFileWriter.FieldValueConverter()

Exported JSON (script result):
{
"type": "FeatureCollection",
"name": "output_layer",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::5514" } },
"features": [
{ "type": "Feature", "properties": { "fid": 1, "text": "aaa", "value_map": 1 }, "geometry": { "type": "Point", "coordinates": [ -553020.462518186424859, -1170227.009075149893761 ] } },
{ "type": "Feature", "properties": { "fid": 2, "text": "sadasd", "value_map": 2 }, "geometry": { "type": "Point", "coordinates": [ -553025.836183290113695, -1170234.612239605281502 ] } },
{ "type": "Feature", "properties": { "fid": 3, "text": "sdasdsad", "value_map": 3 }, "geometry": { "type": "Point", "coordinates": [ -553018.004352234653197, -1170230.038907602196559 ] } }
]
}

Source layer value map setting:


texts:
I have my environmental variables such as elevation, soil, vegetation and others as shape-files.

I want to use a raster command to upload them as tiff files for running a boosted regression tree in ecological niche modelling.

How do i convert my shape files to tiff format?


texts:
I'm trying to create a layer for each county (based on a layer of all counties) and join only those census tracts located in that county. However, the result I'm getting is a layer named for each county but has all counties and census tracks.

coFC = r"C:\TenCounties"
coField = "County_Name"
censusFC = r"C:CT.shp"

# Get a unique list of the counties
values = [row[0] for row in arcpy.da.SearchCursor(coFC, (coField))]
list_of_counties = set(values)

# Create a feature layer for the Census Tracts
arcpy.MakeFeatureLayer_management(censusFC, "lyr_Census")

for county in list_of_counties:
    print county
    # Create a feature layer that contains just the current county
    newout_name = arcpy.ValidateTableName(county)
    arcpy.MakeFeatureLayer_management(coFC, newout_name, "\"County_Name\" = '{0}'".format(county))
    # Select the census tracts that are located in the current county
    arcpy.SelectLayerByLocation_management("lyr_Census", "INTERSECT",    newout_name)
    #Create a layer for each County
    rcpy.CopyFeatures_management(coFC, "lyr_" + newout_name)


texts:
I have sentinel-1 dataset and I want to map Frost filter function on each image in the dataset.
The problem is that when I try to do this. using map, I get the error:

  Required argument (count) missing to function: List.repeat(value,
  count)

  Returns a new list containing value repeated count times.

  Args:

  value (Object): Undocumented.   count (Integer): Undocumented.

I have tried to insert those values, but it also didn't work and I keep get the same error.

this is my code:

var Sentinel1 = ee.ImageCollection('COPERNICUS/S1_GRD')
        .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
        .filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'))
        .filter(ee.Filter.eq('instrumentMode', 'IW'))
        .select('VH','VV','angle')
        .filterDate('2019-10-01', '2020-01-01')
        .filterBounds(geometry);

var clippedSAR= Sentinel1.map(function(im){
  return im.clip(geometry);
});

//Apply Frost filter for speckle correction
function Frost(image, frostDamp, ksize)
// Frost is really a simplified version of GammaMap
// Make sure to pass the frost damp factor (frostDamp) as a NEGATIVE number
{
  // ~~(ksize/2) does integer division in JavaScript
  var distance_kernel = ee.Kernel.euclidean(~~(ksize/2));

  // Square kernel, ksize should be odd (typically 3, 5 or 7)
  var weights = ee.List.repeat(ee.List.repeat(1,ksize),ksize);

  var kernel = ee.Kernel.fixed(ksize, ksize, weights, ~~(ksize/2), ~~(ksize/2), false);

  // Get mean and variance
  var mean = image.reduceNeighborhood(ee.Reducer.mean(), kernel);
  var variance = image.reduceNeighborhood(ee.Reducer.variance(), kernel);

  var B = variance.multiply(frostDamp).divide(mean.multiply(mean));

  var W = B.exp().reduceNeighborhood(ee.Reducer.mean(), distance_kernel);

  return image.multiply(W).reduceNeighborhood(ee.Reducer.sum(), kernel).divide(W.reduceNeighborhood(ee.Reducer.sum(), kernel));
}

//this is when I tried to insert values.
var frost_result=clippedSAR.map(Frost(-1,5));

//originally did that:
var frost_result=clippedSAR.map(Frost);

My end goal: To map this function on all the images in the image collection and then to create time series.

Edit:

I thought  maybe the error is because I have 2 bands, so I seperated them but still get the same error:

var vv=clippedSAR.select('VV');
var vh=clippedSAR.select('VH');

var frost_result=vv.map(Frost(-1,5));

//also tried:
var frost_result=vv.map(Frost);


texts: How to create custom Locator Styles in ArcGIS 10

texts: Does anyone know if there is any public documentation on how to create these?


texts:
I am creating my first map using OpenLayers and I am having trouble. The map was created with QGIS Desktop 2.6, and I am using Ubuntu 12.04. The WMS service is QGIS Service. This is my first try doing anything with javascript, so there's a lot about the code that I don't fully understand.

I first got the map into QGIS Server, and both GetCapabilities and GetMap are working the way they are supposed to. I was able to create this sample map without problems. I have the openlayers map I want to create saved as /usr/lib/cgi-bin/Indiana2/Indiana2.qgs and the projection is 4326. I used the code from here, and just changed the specifics to match my map. Here is the code I am using.

  <title>Indiana Map</title>

QGIS server seems to working fine. I can pull the map using GetMap and all of the layers are visible. I can't think of what the problem is,  but like I said, my javascript knowledge is neglible.

I am not getting any error message when I try to open Indiana2 in firefox. I just get a blank screen.

I am able to connect to QGIS Server using QGIS Desktop.


texts: Aggregate data from source value to mask layer by overlapping area

texts:
Let's assume I have a georeferenced data set seed_layer.geojson and a mask layer leech_layer.geojson.

test_leech.geojson:

{
"type": "FeatureCollection",
"features": [
{ "type": "Feature", "properties": { "content_leech": 1000 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 0.0, 1.0 ], [ 0.0, 3.0 ], [ 3.0, 3.0 ], [ 3.0, 1.0 ] ] ] } },
{ "type": "Feature", "properties": { "content_leech": 2000 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 4.0, 1.0 ], [ 4.0, 3.0 ], [ 6.0, 3.0 ], [ 6.0, 1.0 ]] ] } }
]
}

test_seed.geojson:

{
"type": "FeatureCollection",
"features": [
{ "type": "Feature", "properties": { "content_seed": 10 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 0.0, 0.0 ], [ 0.0, 2.0 ], [ 2.0, 2.0 ], [ 2.0, 0.0 ] ] ] } },
{ "type": "Feature", "properties": { "content_seed": 10 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 0.0, 2.0 ], [ 0.0, 4.0 ], [ 2.0, 4.0 ], [ 2.0, 2.0 ]] ] } },
{ "type": "Feature", "properties": { "content_seed": 10 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 2.0, 0.0 ], [ 2.0, 2.0 ], [ 4.0, 2.0 ], [ 4.0, 0.0 ] ] ] } },
{ "type": "Feature", "properties": { "content_seed": 200 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 2.0, 4.0 ], [ 4.0, 4.0 ], [ 4.0, 2.0 ], [ 2.0, 2.0 ] ] ] } },
{ "type": "Feature", "properties": { "content_seed": 10 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 4.0, 0.0 ], [ 6.0, 0.0 ], [ 6.0, 2.0 ], [ 4.0, 2.0 ] ] ] } },
{ "type": "Feature", "properties": { "content_seed": 10 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 6.0, 2.0 ], [ 6.0, 4 ], [ 4.0, 4.0 ], [ 4.0, 2.0 ] ] ] } }
]
}

I plot both layers ...

import geopandas as gpd

leech_layer = gpd.read_file('../temp/test_leech.geojson')

seed_layer = gpd.read_file('../temp/test_seed.geojson')

base = seed_layer.plot(figsize=(6,3), column='content_seed', legend=True, cmap='Accent', edgecolor='black', alpha=0.75)
base.text(1,1, '0 (10)')
base.text(1,3, '1 (10)')
base.text(3,1, '2 (10)')
base.text(3,3, '3 (200)')
base.text(5,1, '4 (10)')
base.text(5,3, '5 (10)')
leech_layer.plot(ax=base, color='red', alpha=0.5)

... and receive the following figure:

What's the best practice in GeoPandas to aggregate values in a column seed_layer.geojson (e.g. content_seed) to the geometries in leech_layer.geojson as the sumproduct of overlapping area and the value in the content_seed column?

Example:  The red shape on the left side would be assigned the following value:

(2*20+2*20+1*10+1*200)/6 = 41.67


texts: Alternative to Picasa for converting Geotagged photos to KML files

texts:
This is a follow-up question to another thread (Geotagging photos with smartphone camera?) I asked a few days ago.

Using my smartphone, I have figured out a way to geotag photos using an app called Geotag Photos. This app, once you follow all the steps, can create Geotagged jpeg photographs accessible from a Desktop PC.

My goal is to convert these to a Google Earth Photo Overlay (KML?) which places all the photographs in their appropriate geo-locations on Google Earth. I'm lucky enough to have an old version of Picasa installed which can do this very easily. It does this by creating a link with Google Earth and displaying all the photos in their appropriate locations. However, Picasa is being discontinued by Google, to be replaced by the much inferior Google Photos feature online.

I need an alternate way to do this since I will be using another computer which does not have Picasa installed, and also to find a future-proof way to do this. Also, why on earth does Google replace great tools with crappy ones? This reminds me of when they removed Google Street View from the MyMaps feature. I just don't get it...


texts: How to convert shapefile to KML file using ArcPy?

texts:
Python script that checks to see if a specific layer file exists in the current workspace
and, if it does, converts the layer file to a KML file. Does it work with a shapefile?

Code so far:

for file in arcpy.ListFiles("cities.shp"):
    desc = arcpy.Describe(file)
    #print "Name: " + desc.BaseName
    name = desc.baseName

    if name == "cities.shp":
         arcpy.LayerToKML_conversion("cities.shp", "cities.kmz")
         print "Converted"
    else:
        print ("There are no layers in ",env.workspace,".")


texts: How define new 'GeoExt.panel.Map' with OpenLayers-3.6.0 and Extjs-6.0?

texts:
I'm using OpenLayrs-3.6.0 and Extjs-6.0. I want to define GeoMap Panel class like this. In other words, I want to rewrite GeoExt.panel.Map in geoExt-2 library.
for this I define a class as follow:

Ext.define('MyApp2.view.geo.Map', {
    extend: 'Ext.panel.Panel',
    requires: [
    ],
    xtype: 'geoMap',
    layout: 'fit',
    map: null,
    initComponent: function() {
        var me = this;

        var html = '

';

        mm = me.map;
        me.callParent(arguments);
        Ext.apply(me, {html: html});

        me.map = new ol.Map({
            target: 'myId',
            layers: [
                new ol.layer.Tile({
                    title: "Global Imagery",
                    source: new ol.source.TileWMS({
                        url: 'http://localhost:8080/geoserver/world/wms?service=WMS',
                        params: {LAYERS: 'world:worldRaster', VERSION: '1.1.1'}
                    })
                }),
            ],
            view: new ol.View({
                projection: 'EPSG:4326',
                center: [35, 41],
                zoom: 4
            })
        });
    },
});

The strategy is that create a panel and set its html with:

. I then create a ol.map instance and set its target with: "map"(the "id" of div in html panel).

This solution does not work, How do I do?

or

Someone knows a better solution?


texts: Implementing Vehicle Routing Problem in ModelBuilder?

texts:
I am trying to solve the "Vehicle Routing Problem", using ModelBuilder.

For example, using Exercise 7: Servicing a set of orders with a fleet of vehicles from the Network Analyst Tutorial:

add a "Make Vehicle Routing Problem Layer" tool
add a "Add Locations" tool and set "stores" as "input location"
add a"Add Locations" tool and set "Distribution Center" as "input location"
add a "solve"tool to calculate.

That's all ... and I have a error "Routes class is empty."

I think I have to add that something connected with the route class, but I did not understand how to do it in ModelBuilder

Click here for higher resolution image

How to make the scheme correctly?

PS Sorry for my english


texts: gdaldem -alpha flag in Python bindings is missing

texts:
I've been trying to get the python bindings to work in order to keep the processes in memory, making it faster.

I have previously used:

subprocess.call('gdaldem color-relief file1.tif -alpha hramp.txt file2.tif')

which has worked the way I have wanted it to, however I want to move this command over to the python api. My problem is that the call to the color-relief in the python api doesnt add the alpha channel

This is my working code, with the test generations commented out

warp_options = gdal.WarpOptions(dstAlpha=True, srcSRS="EPSG:4326", dstSRS="EPSG:3857", format="MEM")

for file in files:
    try:
        original_file = 'downloads/' + file
        filename = 'combined/' + file[:-4]

        src_ds = gdal.Open(original_file)

        warp = gdal.Warp(destNameOrDestDS='',srcDSOrSrcDSTab=src_ds, options=warp_options)

        hillshade = gdal.DEMProcessing('', warp, 'hillshade', format = 'MEM', zFactor=3)
        hillshade_colored = gdal.DEMProcessing('', hillshade, "color-relief", colorFilename='hramp.txt', format='MEM')

        #gdal.Translate(filename + '_hill.tif', hillshade)

        #I want these to output the same file, however the inmemory hillshade doesnt have an alpha channel
        #gdal.Translate(filename + '_hillcolor1.tif', hillshade_colored, creationOptions=['compress=LZW', 'alpha=yes'])

        #This works the way I want it to
        #subprocess.call('gdaldem color-relief -q -co compress=lzw -of GTiff {0} -alpha hramp.txt  {1}'.format(filename + '_hill.tif', filename + '_hillcolor2.tif'))

        colors = gdal.DEMProcessing('', warp, 'color-relief', colorFilename='ramp.txt', format='MEM')

        overlaid = gdal.Warp('', [colors, hillshade_colored], format='MEM')
        gdal.Translate(filename + '.tif', overlaid, creationOptions=['compress=LZW', 'alpha=yes'])
        ...

I have checked the test page, and the python api page, but they don't seem to mention what the alpha flag in the console command actually does. It is mentioned here, but I'm not sure how to apply it to the python api calls.

So my question is, how do I implement the -alpha  flag in the gdal demprocessing python api calls?


texts: How to get a city road network as line segments?

texts:
I want to have road-network information of different cities(New York for example) in the format of line segments(i.e. every road is represented by segment(s)). I want this for experimental purposes in my college project. I have looked into using QGIS however I am unable to get what I want as I am totally new to the platform. To be precise the question is as follows:

Input : A city name

Output : File containing a list of segments which together approximates the road network of the input city.

I have installed the openlayers plugin in qgis, but as I am totally new I don't know how to use it and where to start. Any pointers will be appreciated.

UPDATE:

So far I was able to get upto a road network in QGIS. This was done as follows(e.g. for a city in India):

Downloaded India Shapefile from gadm.org
Extracted the city I wanted into a new shapefile(Select city from attribute table, and do a vector clip)
Got road shapefile for India from http://download.geofabrik.de/asia/india.html and again did a Geometric intersect with the shapefile obtained in step 2. Saved this as a new shapefile.

Now there is only one hurdle remaining. I have a road network for a city in my hands as a shapefile. How do I convert the road network into co-ordinates for line segments which approximate the network(i.e. each straight segment in the shapefile is represented by two co-ordinates, I want all these co-ordinates).
Any pointers on this would be really helpful.


texts: Joining Attributes by Location to nearest point?

texts:
I have a point feature and I want to perform a spatial join with the nearest point of another point feature (maybe even within a certain buffer distance). I know QGIS has a nice spatial join tool "Join Attributes by Location" and I figured out what the geometric predicates do from Documentation "Join Attributes by Location" QGIS 2.18.

However, there isn't an option nearest, although the QGIS API Documentation for the QgsGeometry class suggests that there is a class nearestPoint().

Does anyone know how to perform this kind of analysis in QGIS?


texts: Divide a polygon into multiple small polygons using a MultiLineString

texts:
I'm trying to divide a polygon grid into multiple smaller polygons basis their intersection with a MultiLineString using Shapely. Following code works good for LineString but not for intersecting MultiLineString.

from shapely.geometry import Polygon, LineString, MultiLineString
from shapely.ops import linemerge, unary_union, polygonize

poly = Polygon([(0,0), (4,0), (4,4), (0,4)])
lines = [LineString([(0,1), (1,1), (2,1), (5,1)]), LineString([(1,1), (2,1), (3,3), (5,3)])]

multi_line = MultiLineString(lines)
merged = linemerge([poly.boundary, multi_line])
# borders = unary_union(merged)
# polygons = polygonize(borders)

While performing linemerge operation, it throws an error -

NotImplementedError: Multi-part geometries do not provide a coordinate sequence

One possible reason could be that I can't do merged_line.coords for intersecting lines like I do for normal linestring. But still I am not able to solve it.

Any thoughts on solving or any other alternative to divide a Polygon basis intersection with MultiLineString?


texts: Splitting large amount of features using arcpy

texts:
We are publishing a geoprocessing service that will split a large amount of features in a file geodatabase.  The script uses the arcpy.defense.SplitFeatures() method, a ResourceSrf layer or a user passed in shapefile grid to split the features, so any features that fall within the grid or ResourceSrf should be split by the grid/ResourceSrf cell size.

The issue we are facing is that some features are split, while others are not, particularly polygon features, the Physiography surfaces are giving trouble as well as some line features.  I am wondering if there is any reason why features would not split when called by this method? Of course the documentation for this method is sparse at best.  Below is the code that splits the features.

def splitter(gdb_to_split, use_alt_surface):

gdb_exchange_folder = os.path.join(r'F:\Exchange\Rework_Temp', gdb_to_split[:-4])
gdb_to_split = os.path.join(gdb_exchange_folder, gdb_to_split)
gdb_to_split_tds = os.path.join(gdb_to_split, 'TDS')

if use_alt_surface:
    gdb_resource_srf = os.path.join(gdb_exchange_folder, 'clip_shape.shp')
else:
    gdb_resource_srf = os.path.join(gdb_to_split, 'TDS', 'ResourceSrf')

physio_srf = os.path.join(gdb_to_split_tds, 'PhysiographySrf')

target_features = r''
arcpy.env.workspace = gdb_to_split_tds
fcl = arcpy.ListFeatureClasses()
for fclass in fcl:
    desc = arcpy.Describe(fclass)
    if desc.shapeType in ('Polygon', 'Polyline') and desc.name not in ('ResourceSrf', 'MetadataSrf', 'PhysiographySrf'):
        target_features += os.path.join(gdb_to_split_tds, fclass) + ';'
target_features = target_features[:-1]

arcpy.management.MakeFeatureLayer(gdb_resource_srf, "split_shape")
arcpy.management.SelectLayerByLocation("split_shape", "ARE_IDENTICAL_TO", gdb_resource_srf)

arcpy.CheckOutExtension('defense')
arcpy.defense.SplitFeatures("split_shape", target_features)

arcpy.Delete_management("split_shape")


texts:
I have georeferenced some aerial images, loaded in QGIS and saved them (files .tiff, aux.xml and .points). These images overlap, however I found that some details in the borders did not match exactly after georeferencing (perhaps due to the lack of GCPs - I am using TPS transformation). Since I am building a mosaic, I wanted to add a few more GCPs to improve the precision and match the elements in the borders of these images.

In order to do so, I reopened Georeferencer, loaded the unprocessed image (.tiff file) as well as the GCPs previously created (.points file) over that same image. The red GCPs showed on the Georeferencer window, but not on QGIS. When added extra points, they showed on both windows, but I noticed that the calculated residual values disappeared (became zeros) in the GCPs table. From what I understand, no residuals means no georeferencing process. When trying to load the raster into QGIS, there are no error messages, only a grey rectangle on a strange scale, different from the projects CRS.

Is this procedure (adding/editing previously saved GCPs) even possible? What I am doing wrong?


texts: Is it possible to put multiple pages in an InfoWindow - ArcGIS API 3.27

texts: I noticed that when you use IdentifyTask to identify features at a given point on the map, it shows an InfoWindow with multiple pages. I was wondering if it was possible to add pages to an InfoWindow without explicitly having to select an area with multiple features.


texts: Applying stretch value to raster output in ArcMap ModelBuilder

texts:
I've been dabbling in ModelBuilder for a couple of days now on ArcMap 10.8 trying to automate a process for work. We make custom maps for the public and for topographic maps we make a mix of a hillshade and a TIN grid for the background. I managed to build the model and it works great but the output hillshade at the end gets a stretch value of standard deviations, and I want a Min-Max stretch.
Is it possible to specify the stretch value for the hillshade in ModelBuilder?
It is not a big deal to just use Image Analysis and change it with one click, but if I can automate the whole process that would be great.


texts: QGIS: Expression for previous value in a child layer

texts: Is it possible to show not the newest but the previous value in a child layer with an expression? (Many child layers are joined to one parent.)


texts:
I have one layer of roads(polylines) and another one that marks the intersections(points) of those roads.

The task is to create laps(one circuit of a racecourse or track) of 2000m, 5000m and so on. They have to be rounds so the startng point is the same as the endpoint.

I have built a network in ArcMap and I did a spatial join between these two layers, but I don't know how to go from there or if I'm even on the right track.

I'd like a solution for ArcMap 10.2.2.


texts:
I have a GRIB file from NOAA that contains several variables for a specific subregion in Norht-East Atlantic.

I'm using ECMWF grip-api Python library to parse the file. While I'm able to get the values for each {latitude, longitude} pair for each GRIB message with this code:

def iterate_latitude_longitude():
f = open(INPUT)

while 1:
    gid = grib_new_from_file(f)
    if gid is None: break

    iterid = grib_iterator_new(gid,0)

    #Get the missingValue for this GRIB message (e.g. NaN)
    missingValue = grib_get_double(gid,"missingValue")

    i=0
    while 1:
        #Get the value from the geo iterator obtained previously, value is a tuple with [latitude, longitude, value]
        result = grib_iterator_next(iterid)
        if not result: break
        [lat,lon,value] = result

        sys.stdout.write("-Point # %d - lat=%.6f lon=%.6f value=" % (i,lat,lon))
        #Check if missing value is present or not
        if value == missingValue:
            print "missing"
        else:
            print "%.6f" % value
        #increase iterator
        i += 1

    grib_iterator_delete(iterid)
    grib_release(gid)

f.close()

While this is great, this file contains forecast data, it was extracted from this place, more precisely, this file:

multi_1.nww3.t00z.grib2

Which contains world data. There are other files for different times: 00-60-12-18-00, thus, every 6 hours.

What I'm missing is the relation between the value obtained for each latitude,longitude and the time dimension. Does anybody know how to extract this info from a grib file, whether using Python or Java?

I know I can use things like wgrib2 or pygrib to convert to CSV and do the other way around...but it just looks so inefficient


texts: How to get layers in x number of subgroups?

texts:
I have a project which has dozens of groups and subgroups. E.g.:

Main_Group > Sub_Group_1 > Sub_Group_1A > Layer_1A
                                        > Layer_2A

How can I access those layers which are contained within Sub_Group_1A?

I currently have this:

root = QgsProject.instance().layerTreeRoot()
for child in root.children():
  if isinstance(child, QgsLayerTreeGroup):
    if child.name() == "Main_Group":
        for subChild in child.children():
            if isinstance(subChild, QgsLayerTreeGroup):
                if subChild.name() == 'Sub_Group_1':
                    # What next?

I have read through the following sources and can get layers in Sub_Group_1 but not in Sub_Group_1A:

http://www.lutraconsulting.co.uk/blog/2014/07/06/qgis-layer-tree-api-part-1/
pyqgis traversing subgroups and group layers


texts: Run a python script every time the page is changed in Data Driven Pages? (ArcGIS 10)

texts: I'm setting up a Data Driven Pages mapbook mxd for a series of parking locations. I have dynamic text that updates according to the current map sheet, but I also have text that needs to update using a script. I need the dynamic text to return the LotType attribute, but it is coded using a domain, so it returns an integer rather than the string associated with the integer (such as "Special Permit", "Fac/Staff", "Metered", etc). I have a VBscript from an old MPS Atlas MXD that returns the associated string, but I don't know of any way to use a script with non-MPS Atlas dynamic text. Is it possible to embed a python script into the MXD that runs whenever the user changes pages in the map book? Or is there a way to include scripts within the core Arc10 dynamic text?


texts:
I am trying to write a script where you preemptively select features and then go to the tool which calls for the feature layer and it populates fields based on certain criteria which I will design of the selected features.

I am simply trying to populate fields right now using the Update Cursor, but the field still come up blank.

here is my code so far:

input= arcpy.GetParameterAsText(0)

fields= ['Conf','Com']
with arcpy.da.UpdateCursor(input,fields) as cursor:
    for row in cursor:
        row[0] == '2'
        cursor.updateRow(row)
        del row
    del cursor

Why are the selected features, or any feature for that matter not updated?


texts:
I´m making an app using opengeo suite SDK and I added a LayerManager to show the layer list and the corresponding leyends. Layers load fine but now I want some of them off by default, added the "selected:false" option on those layers but it doesn´t work and the app loads with all layers enabled.

Any idea how can I load the layers but disabled ?


texts: Batch exporting suites of similar maps in ArcGIS 10.1 / 10.2

texts:
I work for a flood risk consultancy where one of the most common outputs from our projects is a suite of maps created in ArcMap 10.1 / 10.2 showing the results of our hydraulic model simulations. We open a template mxd layout with our logo, a title box, legend etc and drop in background mapping, the model boundary (polygon shapefile), the development site boundary (another polygon shapefile) and then a set of grids (ASCII format) showing results such as maximum flood depth, velocity etc. We end up with a single mxd that can be quickly updated to contain the latest simulation results.

When it comes to outputting the maps, we currently have to manually switch on/off each ASCII grid and then change the text within the title / description boxes (to state the visible dataset) before exporting as a pdf. Doing this once or twice is okay, but we quite regularly have 30+ datasets to map per job, and the client often comes back asking us to make modifications to the model meaning re-pdfing all 30+ maps again - not only time-consuming but also very prone to human error due to boredom!

Is there a tool / plugin available to do this or is this something I might have to write a new python script to do?


texts: Set picture marker symbol for a layer with arcpy - ArcGIS Pro 2.2.4 & ArcGIS 10.6

texts:
A little UP about the question in the title. I found yet this question on a topic : Set picture marker symbol for a layer with arcpy (ArcGIS)

And my question is the same.

My purpose is about how to batch import many picture to picture marker symbol on .style file. The batch should work on ArcPY as possible.

I know that .style could be considered as a .mdb file. I opened previous .style I created on ArcCatalog to watch the file structure.

Is it a right to considered to import the pictures as on a database with ArcPy.

And how scripting this one ?


texts:
My geoscientists like to work with MapInfo and ArcGIS. I require a common 1 band raster format that can hold obscene amounts of data without damaging the data (compression). I have tried .tif and found it to be rather unhappy once it's starts to gain a few GB. I'm now inclined to move to .ers Question 1: What universal analytical raster is best?

In this particular case we have CSIRO file geodatabase rasters (around 7GB uncompressed). I can not see them directly with: ERDAS ER Mapper, or ERDAS Imagine, nor MapInfo Professional 12.5+ (with the fgdb api). Therefore I need to bring them out to an intermediate raster so I can export them to .ers or ecw with ERDAS.  Question 2: What intermediate raster is best?

So far ArcGIS exports to massive .tif has failed as a universal analytic raster and intermediate raster. I'm now testing the Esri grid. Does anyone have good advice?


texts: Extracting site location from original MODIS sinusoidal projection from Google Earth Engine?

texts:
I'm currently trying to find a way to extract site location data from Google Earth Engine for several MODIS products in its original projection (sinusoidal). Yet I have not managed to receive the exact pixel values as provided from MODIS itself.

So far I have used the commands like this.

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

def _data_extract(site_id, sensor, lon, lat, radius, product, bands, res, start, end):

    geometry = ee.Geometry.Point(lon, lat).buffer(radius)

    col = ee.ImageCollection(product).\
                    filterDate(start, end);

projection = col.first().pixelLonLat().projection()

data2 = col.getRegion(geometry, int(res), projection.crs()).getInfo()

df2 = pd.DataFrame.from_records(data2[1:])


texts: How to diagnose an "invalid data" error message using the r.sun GRASS tool?

texts:
I'm using the r.sun GRASS tool to map solar inputs and shade or lack of it in Pacific NW salmon streams. I clipped out the input raster data (a canopy height model, aspect, slope) for a 15 hectares test area. Everything worked fine, after I got everything into EPSG 102748, and converted my canopy height data into meters. Now, trying the same analysis for a larger (580km2) watershed I get no outputs except for log messages that the "Data source is invalid". Same source LiDAR data to derive the CHM, aspect, slope. Same data prep steps. I generated a latitude and a longitude raster map for this run. But the model doesn't like any of the data.

Are there a couple of steps I can take to track down the issue?

Is there an extent limit for the r.run tool, so the data is "invalid" because its extent is too large?


texts: Creating Stream Order with using thresholding

texts: While creating stream network, I implement a thresholding. All cells's value is "1". Later, while creating stream segments, Do I use the second thresholding for separating link to segments ?


texts: Limit Zoom Levels in ArcGIS Online WebMap

texts:
I am adding a feature class to a WebMap that I am then incorporating into a WebApp using WebApp Builder for ArcGIS. I have set the visibility scale on my feature class as I would like it. When accessing the WebApp, and zooming to a feature, it zooms to the closest possible level.

Is there anyway to limit how far the zoom-in goes? Does this have to be done on the WebMap side, some configuration I can set? Or can I edit some .json config files from the WebApp files to achieve this.


texts: Access QGIS rendering map context within QGIS Python function

texts:
I maintain the maps for my local Nordic ski trails. Each winter during the 'Tour de Soup' skiers travel from one soup station to another with the aid of a special event map.

Because the map is already fairly busy I want to offset the soup station icons and link to their true locations with a label line. I am aware of the EasyCustomLabeling plugin but it is not suitable for my needs. EasyCustomLabeling only supports label positioning with a coordinate. That coordinate may be suitable for display at one scale, but zoom in or out and the label will disappear off-screen or start to overlap the feature. The Tour de Soup map includes an inset map, so in the same PDF there are two very different scales. My goal is to offset an icon using display units, e.g. millimetres, so that it displays consistently at multiple scales.

I have an approach that respects the scale of the map in the main QGIS window, using a custom function (below). If I zoom in or out in the main QGIS window the icons correctly maintain the same display offset distance. However this approach is not sensitive to the layout context.

If the main QGIS window's scale is similar to the layout's main map scale then the layout's main map displays correctly. However, because the layout's main and inset maps have different scales I cannot have both display correctly.

EDIT: David Galt's answer below shows how to determine the scale of the different maps in the layout, but this only solves half of the problem. I also need to know which map is currently rendering the feature so that I can choose the appropriate scale from those available. I find the API documentation difficult to navigate and cannot tell if what I want is possible.

Custom Function:

from qgis.core import *
from qgis.gui import *
from qgis.utils import iface

# displayOffsetStr: "100,-55" is 100 units X offset, -55 units Y offset
# conversionFactor: ratio of offset unit to coordinate system's distance unit, so 1000 gives mm to m conversion
@qgsfunction(args='auto', group='Custom')
def offset_point_by_display_distance(displayOffsetStr, conversionFactor, feature, parent):
    # *** this is where the main QGIS window's scale comes from
    scale = iface.mapCanvas().scale()
    locationPoint = (QgsPoint)(feature.geometry().asPoint())
    if displayOffsetStr != NULL:
        offsetParts = displayOffsetStr.split(",")
        offsetX = float(offsetParts[0])
        offsetY = float(offsetParts[1])
        labelPoint = locationPoint.project(offsetX / conversionFactor * scale, 90)
        labelPoint = labelPoint.project(offsetY / conversionFactor * scale, 180)
        return (QgsGeometry)(labelPoint)
    else:
        return (QgsGeometry)(locationPoint)

This function is referenced by a Geometry Generator symbol rule with

make_line($geometry, offset_point_by_display_distance("DISP_OFF", 1000))

The problem comes from iface.mapCanvas().scale() - that is what gives me the main QGIS window's scale instead of the scale of the map that is rendering the feature. The question is "how can I tell which map is currently rendering the feature for which my custom Python function is being called?"

This screenshot shows the current behaviour. The arrows identify the same feature in both the main and inset maps. As you can see the inset map's scale is incorrect.


texts:
I want to use the LDOPE tool unpack_sds_bits to extract the bit information from the State QA layer of the MYD09GA Aqua Surface Reflectance product.

However, I have never run a command line before and have no experience with LDOPE. I have it installed. How do I go about writing a simple line to extract the bits?


texts:
Leaflet provides the possibility of creating non-geographical maps by using CRS.Simple where

  [...] one horizontal map unit is mapped to one horizontal pixel, and idem with vertical.

For the sake of simplicity it can be said that this allows creating maps using pixels as units of measurement. However, I do not use Leaflet, I use OpenLayers (which works pretty well).

Is there any way to implement non-geographical maps in OpenLayers? Looking at its API or examples I couldn't find anything useful so far.


texts: How to use ogr2ogr to only add new features to an existing layer?

texts:
It seems like the correct usage is as follows:

ogr2ogr -update -append -lco SCHEMA=mySchema -f MSSQLSpatial "MSSQL:server=myServer;database=myDb;trusted_connection=yes;" "myShapeFile.shp"

But it is attempting to create a new layer, as if I didn't set the -update and -append flags:

ERROR 1: Layer Addresses already exists, CreateLayer failed.
Use the layer creation option OVERWRITE=YES to replace it.

What is the proper way to -append newly created items?  I'm fine with updating the database with any detected changes to the shapefile if the tool allows it (-update).


texts: How does fTools' Distance Matrix compute the distance and in what units?

texts: I am new to GIS in general, and QGIS in particular. I am trying to find information about fTools' Distance Matrix function in QGIS, but am having trouble finding documentation…I managed to make a matrix, but I would like to know how the distance is computed, and in what units it is measured.


texts: .nc file flipped and rotated in qgis, no CRS, need to be able to project properly to make a masl

texts:
I have two separate files, one is my Land Surface Temperature and the other is the Grid with the lat/lon both are .nc
I wanted to add my grid, which comes in as an attribute table, to my Land Surface Temperature.

Is there a way to merge/join/combine them?


texts: Create spatial view that references multiple databases

texts:
Based on this link, it sounds like I can create a view from separate spatial databases.

  Cross-database queries are permitted and ArcSDE views can reference
  data from multiple databases provided the CROSS_DB_QUERY_FILTER
  parameter in the SDE_dbtune table is set to 0.

But they don't show how. Does anyone know how I can go about doing this?

I'm using MSSQL and Oracle databases.


texts: Can input parameter be optional and Python script tool work when input is not provided?

texts:
I'm building a script tool that processes shapefiles ( GPS data collected on weekly basis).  The tool is set up to take an input for each shapefile ( in this case Sewer, Storm, Sign).  On some days, data for one of these assests may not be collected, so there will be no dataset to provide to the tool.  I thought O.K., i'll just set the input parameter to optional, and if no dataset is provided, it'll just continue on to the next dataset.  Wrong answer, play again.  Is it possible to make this work?

Here's a snippet of the code.

#SEWER#
sewer = arcpy.GetParameterAsText(0)
arcpy.AddField_management(sewer, "UID", "TEXT", 10, "", "", "", "NON_NULLABLE", "", "")
arcpy.AddField_management(sewer, "Phot1_Link", "TEXT", 240, "", "", "", "NON_NULLABLE", "", "")

#Calculate UID field
fieldName = "UID"
expression = "id_gen()"
codeblock = """def id_gen(size = 9, chars = string.ascii_uppercase + string.digits):
return ''.join(random.choice(chars) for x in range (size))"""
arcpy.CalculateField_management(sewer, fieldName, expression, "PYTHON", codeblock)


texts: Recurring SQL Expression Parsing Error on Query Builder

texts:
I am using Query Builder in QGIS 2.14.2-Essen to select one country (Peru) among 237 for which to show two layers (states_provinces, boundary lines) from the Natural Earth database. I am paying close attention to corresponding attribute tables. However, a recurring error message is generated: "OGR[3]: SQL Expression Parsing Error; syntax error, unexpected end of string." The screenshot of the message is below.

I have been unsuccesful at resolving this issue, despite widely searching the web.


texts: Limit access to WMS service after number of requests

texts:
Is it possible to limit access to WMS service for users who made more than X number of requests?
Is there any way this can be done?


texts: Converting an image into a font to label point features

texts: How do you convert an image (ex: company logo) into a font and store it within arcmap so that it can be called upon to label point features with?


texts:
I'm currently running QGIS on Win7 64bit. My installation says it was compiled against and running against PyQt 4.7.1, although according to the Osgeo4W installer and PyQt4.QtCore.PYQT_VERSION_STR, I actually have PyQt 4.8.3 installed under the Osgeo4w directory.

I'm running into a bug in my PyQt code (related to QNetworkAccessManager) that, according to the bug tracking log, was fixed in PyQt 4.8.0. It seems, then, that the PyQt 4.7.1 module is being used, and the 4.8.3 module is being ignored. I'd like to resolve this issue.

I've seen a few posts on various forums that suggest that the generic distribution of QGIS 1.8 is compiled against PyQt 4.8. I've tried reinstalling QGIS using the Osgeo4w installer after removing all references to PyQt from my system path, but I still end up with QGIS compiled and running against PyQt 4.7.1.

Could anyone shed some light on this issue for me?

Should QGIS 1.8 installed through Osgeo indeed be a distribution compiled against PyQt 4.8+?

That is, am I doing something wrong in the installation process?

Or will I need to compile QGIS manually if I want it to be compiled against PyQt 4.8?


texts: Efficiently querying a region-partitioned PostGIS table

texts:
I have two PostGIS tables: one of streets, and one of regions, which, apart from spanning the same spatial extent, are a priori unrelated. The streets queries I perform can be restricted to a given region, by supplying a known id. To exploit this, I added a region_id column to the streets table (using a st_intersects-based update query, which was quite expensive, but performed only once).

My question: once this is in place, will the use of an index on the streets.region_id column be enough to make queries of this type maximally efficient:

select * from streets where region_id = <id> and ...

Or is there an alternative mechanism that would perform better?


texts: RGB composite of Landsat 8 bands not plotting correctly (Python, GDAL)

texts:
I am trying to create a natural color (RGB) image from Landsat 8 data (11 bands). Bands 4, 3, 2, correspond to red, green, blue respectively. When I stack them, however, the images i am getting are not correct. I suspect the issue is with my astype() but I am not sure what 'type' I need.

import numpy as np
from osgeo import gdal, gdal_array

red = gdal.Open(r'/Users/myfiles/Satelite/LC08_L1TP/LC08_L1TP_039037_20130616_20170309_01_T1/LC08_L1TP_039037_20130616_20170309_01_T1_B4.TIF')
green = gdal.Open(r'/Users/myfiles/Satelite/LC08_L1TP/LC08_L1TP_039037_20130616_20170309_01_T1/LC08_L1TP_039037_20130616_20170309_01_T1_B3.TIF')
blue = gdal.Open(r'/Users/myfiles/Satelite/LC08_L1TP/LC08_L1TP_039037_20130616_20170309_01_T1/LC08_L1TP_039037_20130616_20170309_01_T1_B2.TIF')

destination = r'/Users/myfiles/Satelite/rgb.tif'

red_band = red.ReadAsArray()
green_band = green.ReadAsArray()
blue_band = blue.ReadAsArray()

rgb = np.array([red_band, green_band, blue_band])
gdal_array.SaveArray(rgb.astype('uint8'), destination, 'GTiff', red)

The above code with astype('uint8'), which I have seen recommended for others with this issue, produces the following image:

using astype(np.byte) produces a very similar image. Omitting the astype() produces this image:

which, frankly, is bizarre. I have tried creating an empty .tif and writing to it:

gtiff_driver = gdal.GetDriverByName('GTiff')

rgb = gtiff_driver.Create(destination, red.RasterXSize, red.RasterYSize, 3, gdal.GDT_Byte)

rgb.SetProjection(red.GetProjection())
rgb.SetGeoTransform(red.GetGeoTransform())

rgb.GetRasterBand(1).WriteArray(red_band)
rgb.GetRasterBand(2).WriteArray(green_band)
rgb.GetRasterBand(3).WriteArray(blue_band)

but that returns an all white raster. And using np.dstack() returns the exact same image as the first example.

I am not really sure where to go from here, or where I am going wrong. As far as I can tell I am selecting the correct RGB bands. I need the image to maintain its georeference/ would like to automate it for a number of images otherwise I would just do this in photoshop (see last image).


texts: QGIS plugin (or other method) to create multi-layer buffer simultaneousely based on parameter input?

texts:
I'm looking for a way to create buffers around objects from more than one layer at the same time, with specific buffer distances, and possibly rendering live the global buffer (resulting of all the merged buffers).

An input box would ask which layer needs to be considered (selected objects or not...), and for each, the buffer required.

A result layer could then be saved for any parameter combination chosen.

The goal is to obtain an instant view of what land is available with a given set of constraints from different layers.


texts: What Spatial Analyst tool to use when comparing competitive species distributions?

texts:
A spatial analysis question:

I'm trying to compare species distributions and specifically overlap between two species (Water Vole and European Rabbit) that are competitive species. I also want to look how the common shared predator the Minx causes this distribution.

Can anyone suggest what ArcTools I can use?


texts: How to see the 'example provider' in processing toolbox of QGIS?

texts:
I would like to create my own provider with the processing toolbox. I see the example provider in my folder C:\QgisWien\bin...\algs\ExampleProvider.

It seems to be a good start but unfortunately, I don't see the example provider in the processing toolbox.

I see this post : Is there a way to integrate my own custom tools into the Processing toolbox in QGIS?

I try to copy the folder ExampleProvider to my qgis plugins path but nothing.

I have Windows se7en and qgis wien installed.


texts: Placing data driven dynamic text in data view?

texts:
I am using data driven text in some of my map label elements, but I would like those text boxes to move with the data as opposed to being placed in relation to the layout view.  In this way if I have two areas to show on two different maps but each have unique labeling elements around them, I do not have to shift around the dynamic text boxes when I output one map or the other.

Is this possible?  If so, how do you do it?


texts:
I have a choropleth map which displays crime percentage change. It is plotted from a geodataframe and the choropleth intervals are determined in the background by Python.

Unfortunately, the first label in the legend (for white polygons) shows "0.00 - -50.00" when it is in actual fact starting from -286, not 0. Similarly, the last label shows "30.18 - nan".

How do I change these to "Less than -50.0" and "Greater than 30.18"?

Please see my code below:

f, ax = plt.subplots(1, figsize=(16, 16))

PoliceDataMap.plot(column='PercChange', scheme='quantile', k=5, legend=True, axes=ax, colormap='Reds', linewidth=0.2, alpha=0.8)

ax.set_axis_off()

ax.set_title("Crime Percentage Change between Jan 2011 & Jan 2015", **TitleFont)

plt.axis('equal')
plt.show()


texts: Resample raster to a higher resolution with ArcPy

texts:
I am looking for a way to upsample a raster to a higher resolution without interpolating in ArcGIS and arcpy (because I have geoprocessing script which treats a lot of rasters).

For example if I have a 10m resolution raster such as in the schematic attached below on the left I would like to achieve a higher resolution of 2m but with all the same values in the single cells such as on the right side of the schematic below. All the resampling techniques (cubic and bilinear) with arcpy.Resample_management I tried so far are interpolating between the values, however I would like to have a higher resolution with same values per 10m cell (basically what I am looking for is a inversion of the Aggregate Tool in ArcGIS).

I found a solution which works for a single raster by using ArcGIS Export Raster Data, but I am looking for a solution in arcpy because I have a lot of rasters to process and a manual export is not an option. I am looking for a solution specifically in ArcPy/ArcGIS as I have no admin rights on my computer at work and cannot install additional packages/software (I am aware that there is a GDAL command that does this, e.g. here).


texts: Export data driven pages only in changed mxd, with arcpy script

texts:
I have an arcpy script below to export a map book to pdf. The index layer are forty acre boundaries (township and range quarter sections). The data in the mxd is on an oracle database and updated nightly.
I want to use task manager to run the script nightly but I would like to add code so that only mxd's with edits/changes to the data in its geographic extent from the previous days version will export.
The reason for this is there are approximately 12,000 of these boundaries and it takes about 21 hours to run. I am hoping to dramatically reduce that. Here is the code:

import arcpy
arcpy.env.overwriteOutput = True
mxd = arcpy.mapping.MapDocument(r"C:\SamplePath\SampleMXD.mxd")
for pageNum in range(1, mxd.dataDrivenPages.pageCount + 1):
  mxd.dataDrivenPages.currentPageID = pageNum
  pageName = mxd.dataDrivenPages.pageRow.FORTYAC
  arcpy.mapping.ExportToPDF(mxd, r"C:\SamplePath\SampleFileName_" + str(pageName) + ".pdf")
del mxd

Any suggestions are much appreciated!


texts: Normalising image bands in google earth engine before classification

texts:
I have an image that has four bands. Three of them are derived from Sentinel-2 to give NDVI, NDBI and SAVI. The fourth band is derived from a Sentinel 1 SAR image collection. The Sentinel 2 bands have pixel values in the range -1 to +1, but the Sentinel 1 band is in the range +2 to +17.
I want to take these four bands, together with some training samples, to generate a landcover classification.
Do I need to normalise all the pixel values (say 0 to 1) across all four bands, prior to running a classifier such as ee.Classifier.smileRandomForest or can the smile classifier handle regressors that are on a different scale?


texts: Preventing labels from overlapping in ArcMap?

texts:
Im trying to prevent a mxd that has cluster of lables from overlapping these labels. What options do i need to select to ensure this?

Labels have a prefered offset with a 500% Max offset.
I have set the label buffer to 50% as a hard constraint.
Label position is set to best. (I have tried different location settings)
I need every label to be visible.
Turning labels to annotation is not an option as the mxd is used to view results of ground investigation.

I have also tried different variation of the settings listed above.


texts:
I am working with the 10m DEMs and see some round depressions in Texas. I found this in two areas of about 12 sq km but do not have the tiles for the areas in between, so the area where I see this could be much larger. It is more or less at lat 33-34 and long 110-111.

The DEMs look like swiss cheese and the depressions have diameters between 500 and 3000m and the difference in depth between the edge of the depression and the deepest point is about 10-20m.

Is this a problem of these DEMs or do structures like that exist in Texas?


texts:
We have been testing the effectiveness of the k-NN query in MongoDB and Postgis on New Year City's taxi dataset on two aspects: run-time and accuracy. The staggering outcome is that MongoDB detects the k-NN of a given point instantly regardless the value of k (even for k=1 million).

In addition, we have relied on the Haversine and Vincenty functions to evaluate the accuracy and the detected neighbours. The neighbours detected by MongoDB are also more closer to the query point. In other words, Postgis missed some of the closer points.

How could MongoDB be super-fast (i.e. O(1)) and and also be accurate?

I seek scientific evidence.

Code is available openly avaiable on Github.

MongoDB Query

db.nyc2015.find( {
geometry_pk: {
$nearSphere: {
$geometry: {
type: "Point",
coordinates: [point_lon, point_lat] } } } }
).limit(k)

Postgres Query

select id
from trips
order by l_pickup <->
(select l_pickup
from trips
where id = point_ID)
limit k


texts: Creating Bathymetry using .DAT file format in ArcGIS Desktop?

texts: I have the bathymetry data in .DAT format and I want to create the bathymetry map in ArcGIS but don't know how to do that. I know how to create DEM from elevation data points but have never used .DAT format data, neither the bathymetry.


texts: Using ArcPy to copy/paste text boxes from template?

texts:
I am in the process of generating hundreds of maps in ArcMap 10.3 from a few master templates. I recently found out that I have to add a small text box in the bottom-right hand corner of all the maps containing additional information (all the boxes will contain the same text and be placed in the exact same location).

Is there a way I can add the text box to the master templates and use Arcpy to copy & paste the text box on all the maps?

Trying to find a way not to do this manually.


texts: Draw svg-circle on a coordinate with polymaps

texts:
I'm making a tracker system which uses polymaps to display a satellite image. The code is copied from the bing-example on the website of polymaps.

I now want to periodically call an update function which retrieves the current location of my trackers in geoJSON format:

{
"type":"FeatureCollection",
"features":[{
            "type":"Feature",
            "id":"T_xxxxxxx",
            "geometry":{
                          "type":"Point",
                          "coordinates":["5.xx","51.xx"]  // 51N, 5E; Netherlands
                        },
            "properties": {
                          "snelheid":"5.00",
                          "richting":"0.00",
                          "satAantal":"5",
                          "devMessage":"",
                          "devBatVolt":"4.11",
                          "devBatStatus":"F",
                          "devBatLaden":"0",
                          "dt":"2012-06-28 07:07:48","dif":"-6667"
                          }
            }]
}

The update function:

function update(){
  $.getJSON('/_xml/trck/current.php', function(data) {
    $.each(data.features, function(key, val) {
      //update the tracker information in the list
          [....]
      //display an svg circle
            /*
             *  But how?
             */
    }
   });
};

So far so good, the information does get updated with the information from the geoJSON.

Can anybody point me in the right direction how to add an svg-circle on my polymap?


texts:
I want to merge these two features into one (N.B: One polygon shapefile, lots of features). But an error occurred.

Any possible solution please?

QGIS version: 2.14.20


texts:
I'm trying to figure out how I can preserve my GCP's after a gdalwarp. I need the georeferenced GCP's data, so I can calculate the Root Mean Squared Error.

Is there a way to preserve this data, and somehow use it, or are there other ways to calculate the RMSE of a georeferenced raster file?


texts: What purpose do service definition packages serve?

texts:
I'm a software developer. In my world, when you bundle up a set of scripts or binaries to be run on another machine, your goal is usually to publish it to multiple environments. Sometimes the intention is that someone else will use your scripts or binaries inside of an application (libraries), and sometimes you're actually trying to publish an application. In the former case, you generally do all you can to make your tool as environment agnostic as possible, allowing users to pass in information about what resources to use (files, databases, web services). In the latter, you give the user some ability to configure the application so that it can connect to the correct resources. Making an application configurable like this still has uses even when you're only giving it to one client. In particular, it makes it simpler to have multiple testing copies of the application before you actually push to production. (E.g., one for developers to experiment on without bothering anybody, one for internal QA, and of course production.)

When I first heard about service definition files (.sd files) for geoprocessing services, that's what I thought they were supposed to do: let you bundle your code up so you could deploy it to different environments. However, I have found that when you use them, there are zero options for configuring anything at publish time. (See Specify geoprocessing service connection at publish and Prevent parameter defaults in service definition file). This means that a service definition is locked to a single environment, and you have to create a new one for each environment. As a result, I cannot see any way in which service definitions are useful for geoprocessing services since you must recreate them in each environment you want to deploy to.

What is their intended purpose?

What problem do they solve?

When should I be using them?


texts: Which items do not support editing metadata in ArcGIS Pro?

texts:
According to the Best Practices for Editing Metadata,

  You can create metadata to describe most items in ArcGIS Pro.

Further, it states:

  The geoprocessing tools that support managing metadata in
  ArcMap are not yet available in ArcGIS Pro. Importing, exporting, and
  validating metadata must be managed using other applications in the
  ArcGIS platform.

This is confirmed by the Tools that are not available in ArcGIS Pro page.

However, I can't seem to find any resources that show which items do not support metadata creation.

For example, I'm attempting to edit metadata for a map in my Project and the metadata I enter won't save. Even after saving the project.

Step 1: Edit Metadata for one of the maps in my project.

Step 2: Apply Changes

Here I apply the metadata changes to my map.

Note the Save icon is disabled. When I hover over the save icon, it says: "This command is currently disabled. Available only when editing metadata for items stored outside the project such as feature classes and tools."

You can see in the image below, that it appears that my metadata has been applied.

However, when I click on a different map and back on the first one, the metadata is gone. Even if I save the project.

My suspicion is that items within an ArcGIS Project do not support metadata. But I can't find anything to support this. I did successfully saved the metadata for a feature class in a file geodatabase. But not for a map in the project.

What items do not support metadata in ArcGIS Pro (2.1.2)?

It's worth noting that: **Esri's roadmap** lists Metadata as a near-term update (Presumably at 2.2).

  Metadata Import, Export and Sharing - Support majority of metadata
  functionality including the ability to import and export using
  standard metadata formats and sharing metadata across items in ArcGIS
  Enterprise and ArcGIS Online.

Although the word "majority" in the above statement is problematic.


texts: Problem: Matching Data (Attributes) with the Intersection-Algorithm of QGIS

texts:
I'm fairly new to QGIS and struggle with a problem for a while now: I have two Shapefiles that are nearly identical, but have a few different Shapes in them. They also contain a few attributes. I want to merge them together.
So I merged the two Shapfiles together with the Intersection-Algorithm of QGIS. The Intersected Shapefile now contains the attributes i want, but also a lot of double attributes, like in the picture below the attribute in row 1 and 2. The one in row 1 is correct, and the other one incorrect (it doesn't fit to the attribute wbz_2)
I want a list with only the correct matched attributes. I tried to dissolve the list after the attribute "wbz" but that produces a lot of false matches. How can i only get the right matches? (right matches = where wbz equals wbz_2)


texts: How would one get the end points of a polyline?

texts:
I am attempting to write a script to get the end points of a series polylines, and then use those end points to create a polygon.   I'm not sure which tools to begin to use on this project.

I would prefer to use QGIS or even better a python script using open source tools.


texts: QGIS GRASS using v.drape FileNotFoundError: No such file or directory: '/usr/local/Cellar/osgeo-grass/7.6.1_5/grass76/etc/lock':

texts:
I'm running QGIS 3.8.3 and am trying to use v.drape. When I run it from the processing toolbox, I get a FileNotFoundError that's pointing to a strange /usr/local/Cellar/osgeo-grass.... folder.

I'm running OSX 10.14.6

I'm not sure where it's grabbing that user directory. I've been a constant upgrader over the years, so not sure if there is some older dependency or something that needs to be update... or even where to start.

Full log is:

         Starting GRASS GIS...
         Traceback (most recent call last):
         File "/Applications/QGIS3.8.app/Contents/Resources/grass7/bin/_grass76", line 2207, in
         main()
         File "/Applications/QGIS3.8.app/Contents/Resources/grass7/bin/_grass76", line 2143, in main
         grass_gui=grass_gui)
         File "/Applications/QGIS3.8.app/Contents/Resources/grass7/bin/_grass76", line 1390, in lock_mapset
         ret = call([gpath("etc", "lock"), lockfile, "%d" % os.getpid()])
         File "/Applications/QGIS3.8.app/Contents/Resources/grass7/bin/_grass76", line 247, in call
         return subprocess.call(cmd, **kwargs)
         File "/Applications/QGIS3.8.app/Contents/Frameworks/Python.framework/Versions/Current/lib/python3.7/subprocess.py", line 323, in call
         with Popen(*popenargs, **kwargs) as p:
         File "/Applications/QGIS3.8.app/Contents/Frameworks/Python.framework/Versions/Current/lib/python3.7/subprocess.py", line 775, in __init__
         restore_signals, start_new_session)
         File "/Applications/QGIS3.8.app/Contents/Frameworks/Python.framework/Versions/Current/lib/python3.7/subprocess.py", line 1522, in _execute_child
         raise child_exception_type(errno_num, err_msg, err_filename)
         FileNotFoundError: [Errno 2] No such file or directory: '/usr/local/Cellar/osgeo-grass/7.6.1_5/grass76/etc/lock': '/usr/local/Cellar/osgeo-grass/7.6.1_5/grass76/etc/lock'
         /Applications/QGIS3.8.app/Contents/Resources/grass7/grass76.sh: line 20: /private/var/folders/p4/212f8sh94g5fyf4ys8bw39yc0000gn/T/processing_7c31a02da3a34505a56e56b476bb1792/grassdata/temp_location/PERMANENT: is a directory


texts: 'QgsLayerTreeLayer' is not defined when running from python plugin but not when run in the console in QGIS 3

texts:
When I run the following in a plugin I get

  name 'QgsLayerTreeLayer' is not defined

print('Loading raster layer...')
#QgsProject.instance().removeMapLayer(tLayer)
print("%s called %s" % (filePath, fname))
tLayer = QgsRasterLayer(filePath, fname, 'gdal')
QgsProject.instance().addMapLayer(tLayer, False)
layerTree = iface.layerTreeCanvasBridge().rootGroup()
layerTree.insertChildNode(-1, QgsLayerTreeLayer(tLayer))

However in python console it works

filePath=r'c:\resources\landcover.img'
tLayer = QgsRasterLayer(filePath, 'imagery', 'gdal')
layerTree.insertChildNode(-1, QgsLayerTreeLayer(tLayer))

Note: When I try the plugin without
layerTree = iface.layerTreeCanvasBridge().rootGroup()
I get

  name 'layerTree' is not defined

The plugin imports

from qgis.core import QgsApplication, QgsVectorLayer, QgsRasterLayer, QgsProject
from qgis.utils import iface

Should I import something else? This is based on https://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/loadlayer.html

Purpose is to get the raster to load under all layers and also to stay on the bottom if other vectors are loaded after it.


texts: Distance to points tool giving errror ('Error in .pointsToMatrix(xy) : A points matrix should have 2 columns') in R

texts:
I am using the DistanceToPoints tool to make a raster where cell values contain the distance to the nearest point.

r <- raster('./Raster.tif')

s <- readOGR('./Shapefile.shp')

dfp <- distanceFromPoints(r, s)

It is giving me an error: Error in .pointsToMatrix(xy) : A points matrix should have 2 columns

Raster and Point layer are in the same projection.

> summary(s)
Object of class SpatialPointsDataFrame
Coordinates:
              min         max
coords.x1 4352987 5334474.056
coords.x2 2180706 2887232.448
coords.x3       0    1100.361
Is projected: TRUE
proj4string :
[+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]
Number of points: 792668
Data attributes:
    OBJECTID            id
 595    :416281   Min.   :   0
 18     : 92323   1st Qu.:1486
 261    : 19468   Median :2972
 463    : 17326   Mean   :2972
 472    : 12768   3rd Qu.:4457
 596    : 11914   Max.   :5943
 (Other):222588   NA's   :786724

As you can see in the picture, the raster is not one whole 'thing' but patchy. Do you think that could be the problem? The tool was working fine with another whole,'normal' raster.


texts: How do I properly model a transit map loop in QGIS?

texts:
I am currently trying to model the transit network of my city in QGIS, and I have come across a problem: I cannot model the following loop right:

I have given both lines an offset to make them appear parallel to each other, but I cannot get them to connect to the main line properly. Unfortunately, I cannot currently post more images to describe my problem as I lack the reputation to do so.


texts: Google Earth Engine: label objects on map

texts:
In my Google Earth Engine script, the user can click on a point on the map, and then a time series history chart of that pixel value is displayed in the console. However, after the user has clicked on a few points, it's not clear anymore which chart belongs to which point.

That's why I'd like to add labels (consecutive numbers probably) to all the points on the map and to the chart. However, I have no clue how to add text at a coordinate location on the map.

I know that printing the point coordinates is a workaround but that's not user-friendly at all.

Here's what I have (without the chart because that's not really part of the problem):

Map.style().set({cursor: 'crosshair'});

var i=0;
var print_point = function(coords, map) {
  i++;
  var coord_array = Object.keys(coords).map(function (key) { return coords[key]; });
  var point = ee.Geometry.Point(coord_array);
  print('point ' + i, point);
  print('placeholder for time series chart');
  Map.addLayer(point);
  // put i somewhere near that point on the map
};

Map.onClick(print_point);

How do I add labels to the points on the map? (I mean label in the generic sense, it doesn't have to be a ui.Label.

EDIT: Any answer that makes it easy for the user to relate each chart to a point will be fine. It doesn't need to be a label on the map. For instance a solution where the chart pops up when a point is clicked is good too. However, it's important that GEE doesn't have to recalculate the chart every time if the same point is clicked multiple times.


texts: Setting style during table creation in postgis

texts: Is there any way to set a style from public.layer_styles table to a postgis table while creation?


texts: Placing GeoServer WMS layer within OpenLayers 3

texts:
I have the following basic OpenLayers 3 code that opens a OSM background mapping and this works fine:

<!doctype html>

        <title>Hello PMS</title>
        <link rel="stylesheet" href="../assets/ol3/css/ol.css" type="text/css" />
        <link rel="stylesheet" href="../assets/css/samples.css" type="text/css" />

I am now looking to bring in a WMS layer from GeoServer with the code below, where would it be placed and what needs to be changed?

        var pms = new ol.layer.WMS("PMS",
        "http://10.0.1.127:8080/geoserver/PMS/wms",
        {layers: 'PMS:QGIS'} );
        map.addLayer(pms);


texts: Dynamically changing spatial reference using ArcGIS API for JavaScript?

texts:
We have a web application that has tiled layers published in GCS.

We want to publish a replica of the layers in a different projection system such that :

When a user is at a specific zoom level, GCS based tiles are rendered.
  However when the user zooms in, the spatial reference of tha map is changed to a different projection system and tiled published for that projection system are rendered.

We are using ArcGIS JavaScript API 3.25.


texts:
Wondering if anyone knows how to get multiple geometries of what a user has drawn on an openlayers 3 map.

At the moment i am using

        var coords = feature.getGeometry().getCoordinates();

Is there a openlayers method to get the geometries of all the polygons/points etc. drawn?

hope that makes sense


texts: I need a DEM file for Scotland, UK. Do you know where I can get it from for free? I need to create a watershed for an area of study.


texts: How do I get layer feature type in OpenLayers?

texts:
I am new to OpenLayers and I wanted to add a modify feature function in my application, for this I need to know the feature type of the layer (i.e. point, polygon or line).  Because as per the feature type I have to give permission to the user to draw a line, point or polygon.

So is there a way that I can get the layer feature type?


texts: clipping ESRI shapefiles along with shx and dbf

texts:
I've seen other answers detailing how one can clip a .shp file with ogr2ogr but this just returns a clipped set of geographic features.  I also want to clip the .shx and the .dbf in the same manner.  Is this possible?

For instance, if I have all-of-rhode-island{.shp, .shx, .dbf} in my current directory, I could try:

ogr2ogr -clipsrc 41.31 -71.87 41.50 -71.53 smaller-ri all-of-rhode-island

But ogr2ogr can't figure out which driver to use with that input file.

My actual goal is to take a TIGER Census file and reduce its size drastically.  I only need a few features, but I do need "matching" .dbf and .shx files.


texts:
I have a polygon layer which looks like this .

I used the QGIS multiparts to singleparts tool and also polygon to lines. In both cases i  got this type of result

.

I want to extract polygon into individual lines like this


texts: Open Source GIS tools for using LiDAR to estimate storage behind dams/culvert

texts: Can you recommend an open source GIS tools for modeling the water impounded on the upstream side of a dam using LiDAR data?


texts: Cannot view WMS of GeoServer WFS unless logged in?

texts:
I'll use the import extension to import a KML file. When I go to Layer Preview and use the OpenLayers option (basically, a WMS GetMap url), it displays. When I try to use the url in a browser and I'm not logged in, however, I get "Could not find layer". Same result when I use WFS with OpenLayers -- the layer shows up if I'm logged in to GeoServer, doesn't show up otherwise.

I've installed this through the geonode package from the geonode/stable ppa onto Ubuntu 14.04. All of the settings are default. I can view layers and vectors uploaded through GeoNode without being logged in to GeoServer, but using GeoServer's import extension is what leads to the issue. Beyond that, I don't know enough about GeoServer's security to know what to change. I want the WFS and WMS in particular to be available to everyone, or better yet available to a select URL.


texts: Creating (x,y) coordinates in .nc from .grib2

texts:
I have to use a .nc file as input in a model and such model requires projection_x_coordinate and projection_y_coordinate as variables but I don't have them. Is it possible to create them using pyproj, for instance? (file available here). My projparams are:

{'a': 6371229,
 'b': 6371229,
 'proj': 'stere',
 'lat_ts': 60.0,
 'lat_0': 90.0,
 'lon_0': 264.0}

My first idea was to use the information of the number of points (Nx in the file), the horizontal resolution (Dx) and the first points (e.g. latitudeOfFirstGridPointInDegrees) to create the desired variables but I don't know if it is the right thing to do.


texts: Using Calculate Field to determine if any segment is greater than a certain length in ModelBuilder

texts:
I have a model that isolates a place's road network based on my input of the place's name. What I need to determine is if any of the road segments in the network are greater than 500m. To do so I've written the following into a Calculate Value tab in ModelBuilder:

def isLongerThanAllowed():
import arcpy
fc = "C:\Users\USER1\Documents\ArcGIS\Scratch.gdb\placeNetwork"
fields = "LENGTH_GEO"
with arcpy.da.SearchCursor (fc, fields) as cursor:
    #iterate rows
    for row in cursor:
        if row > 500:
            return "TRUE"
        else:
            return "FALSE"

For some reason it is returning True even when all road lengths are less than 500m.

The end goal is to use another model I created that splits the road features greater than a certain distance in two so that no segments are greater than 500m. So the logic of the end result will be "while the road network has roads greater than 500m, run the split model".


texts: Automating update of network dataset in ArcGIS Server

texts:
I am trying to set up an automated process that does a nightly extract of the network dataset from SDE and publish to:

ArcGIS Server 10.3 as a network service
fGDB as network dataset

It seems you can't dissolve a network dataset in SDE so you have to move it to a fGDB. Ideally I would have liked to copy and paste from SDE directly into a fGDB but can't find a way to automate this. I have done the following as a work around:

Arcpy script that imports data to fGDB, updates FCID
Manually create new network dataset in fGDB
Arcpy script that updates alternate ID and dissolves network
Manually push network dataset to ArcGIS Server (will script this)

I would be interested to hear in other peoples experiences in


texts: ArcGIS Calculate Fields gives Warning 000405 No records within table?

texts:
Situation:

I have feature layer sitting outside processing extent in Environment settings.
If I open layer's table and use field calculator - Calculate Geometry on some numeric field:

it works fine.

However if I am trying to populate the same or other field with say number of 1, I'll get this error:

Doesn't sound like a big deal, it is enough to fix extent. However a few times when extent was changed long ago I had a hard time to figure out what is wrong with my calculations or maybe layer itself.

Question:
What is a rationale behind blocking simple calculations?


texts:
I am getting the following  error message when trying to utilise an osm pbf from the Geofabrik (latest Europe file):

org.xml.sax.SAXParseException; lineNumber: 1; columnNumber: 1; Content is not allowed in prolog.

The information on it seem to be more Java related than OSM, but I might be mistaken and it has something to do with white spaces in the osm file.

Can somebody help here?


texts: Gravity related height from ellipsoidal height using pyproj

texts:
I'm using the pyproj library (http://pyproj4.github.io/pyproj/stable/index.html) to get translate ellipsoidal heights into gravity related height.

Code:

from pyproj import CRS
from pyproj import Transformer

crs7416 = CRS.from_epsg(7416)  # dk1 dvr90
crs7417 = CRS.from_epsg(7417)  # dk2 dvr90
crs4326 = CRS.from_epsg(4326)  # ellipsoid 2d 1984
crs4979 = CRS.from_epsg(4979)  # ellipsoid 3d 1984
crs7912 = CRS.from_epsg(7912)  # ellipsoid 3d 1980
crs7789 = CRS.from_epsg(7789)  # geocentric reference

geocentrictrans = Transformer.from_crs(7789, 7417, always_xy=True)
newgeocentrictrans = geocentrictrans.transform(12.568337, 55.676098, 100)
print(newgeocentrictrans)

transformer = Transformer.from_crs(4979, 7417, always_xy=True)
newcoord = transformer.transform(12.568337, 55.676098, 100)
print(newcoord)

Output:

(9451442.134520395, 0.0, -6378079.9229382025)
(347093.2566329163, 6172711.4819890065, 100.0)

The transformation seem to work with the geocentric transformation but for some reason the height isn't transformed when i go from ellipsoidal- to gravity related height. As far as i can understand i am using the correct EPSG-references, but correct me if i'm wrong.

It seems like the pyproj library should be able to handle this but i'm missing some piece of the puzzle.


texts: Convert MultiLineString to LineString using QGIS

texts:
In QGIS, I've tried the following options to convert my GeoJSON file from MultiLineString to one or multiple LineStrings (translating the Danish version):

Vector -> Geometry Tools -> Polygons to Lines

Vector -> Geometry Tools -> Multipart to Single parts

Vector -> Geoprocess Tools -> Dissolve

Each option produces just the same MultiLineString.

How do I convert the MultiLineString to LineString(s) using either QGIS or another GIS software?

http://bl.ocks.org/d/880c32ce5acfd18a289564c39f7b4f16


texts: Is there a python equivalent of gdaltindex?

texts:
I am using the following gdal shell command (gdaltindex) to create a raster tile index from a directory of raster (tif) files.

find /path/to/raster/directory -name "*.tif" -exec gdaltindex tile-index.shp {} \;

Is there an open source Python approach that can produce the same raster tile index shapefile?


texts: Layer styling: custom text as a polygon fill in QGIS

texts:
How would one fill a polygon with a custom text? I know you can use this approach to create a label-wise polygon filling (= labeling polygon with their attribute's values), but is there an option to just specify a simple text (not a simple character!) as an filling attribute?

This is the desired result for the polygons with the the example text "Example":


texts:
I'm trying to join 8000+ daily rain tables to one climatic stations shape. The goal is to get 8000+ daily rain point features. I used to work with the next arcgis model, but now i'm trying to put it in a script.

As I'm new to python, I worked my way to this script:

            # Import arcpy module
import arcpy

# Definir espacio de trabajo:
arcpy.env.workspace = r'E:\CLD_WEAP\ppt.gdb'
arcpy.env.overwriteOutput = True
# Definir variables locales:
layerName = r'E:\CLD_WEAP\ppt.gdb\estaciones'
in_field = "code"
join_table = arcpy.ListTables ("*tr")
join_field = "Transposed"

print join_table

for table in join_table:
    output = table + "_j"
    arcpy.AddJoin_management(layerName, in_field, join_table, join_field, "KEEP_ALL")
    arcpy.CopyFeatures_management(layerName, output)
    arcpy.RemoveJoin_management (layerName)
    print "output"
print "P R O C E S O   T E R M I N A D O"

I'm getting this error, and don't know what to do...

Traceback (most recent call last):
  File "e:\CLD_WEAP\scripts\join_1.py", line 22, in <module>
    arcpy.AddJoin_management(layerName, in_field, join_table, join_field, "KEEP_
ALL")
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\management.py", line 6120, in AddJoin raise e RuntimeError: Object: Error in executing tool

Can anyone help me with this?


texts: Is there a reason I should use gdal's Python library instead of just using os.system(gdal command)?

texts:
From my understanding, there is gdal installation that you can just access through the shell then there is the Python library version for it with 'Python bindings' that map to the original.
My current code is a mix of working in the GDAL Python library with the functions and classes in that library, and just calling os.system and running the gdal commands directly through that when that feels more convenient.
Is there a reason why that might not be good to do? This is shared code so I think there could theoretically be an issue with finding gdal/ogr if ran through shell with os.system, but by default it seems to work, and you would need to install gdal/ogr anyway to use the Python library.
I guess I'm only asking because it feels like this might be something I shouldn't be doing for some reason or would be frowned upon.


texts: How to put GeoServer layers behind authentication?

texts:
I want to publish map layers using OpenLayers and GeoServer. While this is easy enough to set up, I'm a bit lost when it comes to authentication possibilities. What I want is to create different users that will, by logging in with their own user name and password, be directed to their own map window. How can this be securely done so that none of my users will be able to access each other's layers?

I have understood that it is quite easy to "hack" the setup by, for example, making queries to the GeoServer where one simply asks it to list all the available layers.

What are my options? I suppose apache could be of help when setting it up in my ubuntu server? Any tutorials on this subject?


texts: Which projection is best when implementing regional datasets in a global one?

texts:
I have a question regarding the projection. I would like to set up an interactive atlas where the base layer is a global overview of e.g. peatland. In a next step, I would like to set into the global overview point layers which enables to zoom in and to have a look on a regional map of peatland. For the global layer, I would like to use WGS84; for the regional maps I would like to avoid using several different projections, i.e. I would like to use max 2 different projections for regional maps worldwide.
Do you have suggestions which projections are best for the above-mentioned purpose?


texts: What causes the striations in this slope map?

texts:
I have created this slope map using a 10m DEM from alabamaview.org and I just noticed it has small striations (zoom to see):

What could cause this, and is it a big deal, or am I over-analyzing the map?


texts: GRASS GIS: Interpolate DEM from irregularly spaced coastal survey data

texts:
I am attempting to interpolate the elevations of data sets of cross-shore survey transects. The data set size range from 50k to 300k points with ~100m between the transects and 1-5m between the points within a transect. Find an image here. Apologies, I'm still working on quickly drafting figures in GRASS.

I interpolated the data to a 2m resolution raster with v.surf.rst. I used the default parameters and set a mask based on the data's concave hull. During the computations I received the following messages:

Ignoring xxx points (too dense)

and

WARNING: Taking too long to find points for interpolation - please change the region to area where your points are. Continuing calculations...

The RMSD which GRASS output with the rasters averaged about 0.35m absolute and 3% of the range. Can anyone tell me if these errors are typical?

I also tried to set scalex and theta (other params were default), but during the computation I received this warning:

WARNING: Overshoot - increase in tension suggested. Overshoot occurs at (1340,1023) cell. Z-value -779923520635822222878965760.000000, zmin -8.433820, zmax 3.060190.

And one of the cells had that large negative value.

I don't really know how to tune the parameters. I'll probably script an optimization using the -c flag, if it seems warranted after a quick lit review. Any other suggestions on changing the default params?


texts: How to create a table with one polygon and one point?

texts:
using postgresql and postgis, i would like to create a table with one polygon and one point. then use st_contains to check if the point is inside the polygon or not.

so first, I should create a table :

CREATE TABLE my_schema.my_table (id serial);

then add

SELECT AddGeometryColumn( 'my_table', 'geo_polygon', 4326, 'POLYGON', 2);
UPDATE my_table
 SET geo_polygon=GeomFromText('POLYGON((0 0,4 0,4 4,0 4,0 0))',4326)
 WHERE id='p1';

SELECT AddGeometryColumn( 'my_table', 'geo_point', 4326, 'POINT', 2);
UPDATE my_table
 SET geo_point=GeomFromText('POINT((2 2))',4326)
 WHERE id='p2';

Finally, to check containing i should use :

SELECT ST_Contains(geo_polygon, geo_ponts, 4326))
    FROM my_table
    WHERE ???

I am not sure about the condition here. do I have to use the ids ? and how ?


texts: OGC: Retrieve all layers with given column name

texts:
I have a given column name city and a list of layers from a GetCapabilities on an OGC server request. Is there a way to get all of the layers that have a column called city without sending a request for each layer?

My ultimate goal would be to provide something like a query box for a user where they can enter city = springfield and come back with all of the results from the layers that have a column called city.


texts: Resample/Increase the resolution of Digital Elevation Model by QGIS

texts:
I have a DEM of 5m pixel. I wonder is it okay to convert it to 2m pixel? How to do that using QGIS?

And what would happen to the original data, I mean the accuracy? Would it be reliable?


texts: python code multirasters zonal statistics

texts:
env.workspace = "C:\Users\\Documents\pythonMODIS\extract"

# Set local variables
inZoneData = "mask.shp"
zoneField = "Classes"
inValueRaster = "mi1511171_clip.tif"
outTable = "zonalstattblout02.dbf"

# Execute ZonalStatisticsAsTable
outZSaT = ZonalStatisticsAsTable(inZoneData, zoneField, inValueRaster,
                                 outTable, "NODATA", "MEAN")

I have several modis images and one layer mask for 3 regional field classes. I want to iterate zonal statistics in all rasters data set,  get a summary dbf table but keeping each raster name, mean, min, max by region class.

I had run once for each image but I want to iterate all the images together. I wonder if I need to define a rasterlist in the workspace and to include that in outZsat

I need some help with the code.


texts: Finding position of character within string in Spatialite?

texts:
I have a bunch of data generated as a .dbf by ArcMap that I imported in Spatialite for processing.

One of the columns (text data) is of critical importance since it contains the ID of a point and a second ID of another point separated by a dash: 54 - 145 for instance.

What I need is to be able to extract the number that is before the "-". In order to do this, I found about the instr(X, Y) function wihch should be able to locate the position of the "-":

  The instr(X,Y) function finds the first occurrence of string Y within string X and returns the number of prior characters plus 1, or 0 if Y
  is nowhere found within X.

After this, substr(X,Y,Z) should be enough to get only the first number. It should be noted that the length of the first number ranges from one to six digits, extracting the substring based on a constant length is this not possible.

  The substr(X,Y,Z) function returns a substring of input string X that begins with the Y-th character and which is Z characters long.

After several tests, it seems that instr(X,Y) does not exist in Spatialite. Does anyone know about a function similar to instr(X,Y) in Spatialite?

(I am using Spatialite-GUI 1.6.0 with Spatialite 4.0.0).


texts: Extracting individual trees from DSM in QGIS

texts: I am trying to extract individual trees from a high resolution DSM in QGIS. I need to identify the trees and extract their individual height from the raster. Assuming I have identified the areas that I know contain trees, what tool can I use to get this information? Perhaps using some kind of local maxima?


texts:
i'm updating a SOE from pre 10.1 arcgis server to arcgis server 10.2.2.

The orginial code included the following part, which is not working anymore:

IMapServer3 mapServer = serverObjectHelper.ServerObject as IMapServer3;
IMapServerObjects3 mapServerObjects = mapServer as IMapServerObjects3;
IMap map = mapServerObjects.get_Map(mapServer.DefaultMapName);

I found that pre 10.1, this code would only work for MXD and not for MSD documents. ESRI says this:

  All map services now use the fast drawing engine that was associated
  with MSDs in versions 10.0 and previous.

So it seems logic this piece of code would not work anymore. Anyone has an idea on how I get my IMap variable in 10.2?

Thanks


texts: About grass algorithm differenece between v.to.rast.value and v.to.rast.attribute in QGIS

texts:
When i use the v.to.raster.value algorithm in QGIS to convert vector polygon shapefile to raster, the algorithm runs but at later stage, during completion, it shows an error message that the output layer could not be loaded. Why is this happening?

I am using qgis v 2.8.2 32 bit in windows 8.1 64 bit os.

And my another question is what is the difference between two algorithms of grass v.to.rast.value and v.to.rast.attribute?

Can they be used selectively or they have some specific purposes?


texts: raster::rasterize set wrong values with enclosed polygons

texts:
Similarly to this post

rasterize (R package raster) fail to rasterize island polygons?

the rasterize function in the raster package is giving me problems with enclosed polygons. In my case, some enclosed polygons are rasterized with cell values equal to the enclosing polygons, as they were somehow 'merged'.

Here's my sample code

require(sp)
require(rgdal)
require(rgeos)
require(raster)

## set wd
wd <- 'M:/foo'
setwd(wd)

## open shapefile
sample.poly <- readOGR(dsn=wd, layer='sample')

## rasterize
resol <- 10
r <- raster(res=resol, ext=extent(sample.poly))
sample.r <- rasterize(sample.poly, r, 'Code')

## export files
writeRaster(x=sample.r, filename='sample.tif', format='GTiff')
writeOGR(obj=sample.poly, dsn=wd, layer='sample', driver='ESRI Shapefile')

Here's the structure of table of attributes

str(sample.poly@data)
'data.frame':   11 obs. of  2 variables:
$ OBJECTID: num  6249 14593 14614 15434 15683 ...
$ Code    : int  33200 32110 33300 32410 32210 31210 33300 32110 32210 31210 ..

data$Code are CORINE LCM-like codes. I already tried to set them as characters instead of integers, but the same error occurred.

Here's how the sample shapefile looks like (in QGIS)

and here's the result of rasterize

As you can see, the enclosed polygon marked with yellow arrow is properly rasterized, while the one marked with the red arrow is 'dissolved' in the enclosing polygon.

Any suggestion? Is it a bug? This is just a sample taken from a larger shapefile, with dozens of such errors.