texts: What do the warnings generated from `raster::crs` in `sp::SpatialPointsDataFrame` mean? [R]

texts:
Take the latest versions as described here:
install.packages("rgdal")
# also installing the dependency ‘sp’
packageVersion("sp")
# [1] ‘1.4.2’
packageVersion("rgdal")
# [1] ‘1.5.16’

In this case, the over function returns many warnings which we don't understand:
# Warnung in showSRID(uprojargs, format = "PROJ", multiline = "NO")
# Discarded datum WGS_1984 in CRS definition,
# but +towgs84= values preserved

Details:
class(spdf)
# [1] "SpatialPointsDataFrame"
# attr(,"package")
# [1] "sp"
class(departement_boundaries)
# [1] "SpatialPolygonsDataFrame"
# attr(,"package")
# [1] "sp"
raster::crs(spdf)
# CRS arguments:
#  +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0
raster::crs(departement_boundaries)
# CRS arguments:
# +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0

So everything seems ok to me, but:
df_in <- sp::over(spdf , departement_boundaries , fn = NULL)
# Warnung in showSRID(uprojargs, format = "PROJ", multiline = "NO")
# Discarded datum WGS_1984 in CRS definition,
# but +towgs84= values preserved
# Warnung in showSRID(uprojargs, format = "PROJ", multiline = "NO")
# Discarded datum WGS_1984 in CRS definition,
# but +towgs84= values preserved
# Warnung in showSRID(uprojargs, format = "PROJ", multiline = "NO")
# Discarded datum WGS_1984 in CRS definition,
# but +towgs84= values preserved
# Warnung in showSRID(uprojargs, format = "PROJ", multiline = "NO")
# Discarded datum WGS_1984 in CRS definition,
# but +towgs84= values preserved

Comments:

I have problems to make a small reproducible example - the departement_boundaries from Zurich are open source, but its a huge string from dput.
The docs ?sp::over discuss the data.frame cases x = "SpatialPointsDataFrame", y = "SpatialPolygons" and x = "SpatialPoints", y = "SpatialPolygonsDataFrame" but this code ran with previous versions of sp and rgdal (see link above in this question) without errors or warnings.


texts: Merging multiple fields with same parcel id yet different attribute values, while keeping all fields?

texts:
I am using ArcGIS 10.3 for Desktop.

I have a property parcel layer that when Intersected, each parcel is broken up into hectares of different land covers. My question is can this be merged (or something similar) back into one attribute with each field populated with the different hectares? Might be a very vague question but it looks similar to this.

parcelID     field1      field2      field3      field4       field5

[1]            0           60          0       =[field2] 60     0

[1]            0           0           42          0        =[field1]+[field3] 67

[1]            25          0           0           0        =[field1]+[field3] 67
_________________________________________________________________________________

parcelID     field1      field2      field3      field4       field5

[1]            25          60          42      =[field2] 60 =[field1]+[field3] 67


texts:
My question is about the calculator in the QGIS graphical modeler. I explain my problem issue first:

I want to build a modeler that creates buffers around layers and makes a selection on a second layer if the polygons are inside the buffers. You see a picture of that modeler here under.

In fact, my modeler is working well until I add a calculator. My input layers are in WGS84 projection, that's why the buffer is measured in decimal degrees. I want to add a calculator that divides the distance by 111 (approximately) so that I can seize the distance in kilometers in the input field.

The calculator is a modeler function that specially exists for that issue of modifying numerical values (entered by the user or generated by some other algorithm). The documentation only said to use the modelers value as a variable. So, I tried differents syntaxes:

a / 111
Distance1 / 111
Distance1 / 111
"Distance1" / 111

etc.
But I always get the same error message:

According to some discussions, it means that the assignment of the variable is not understood what seems logical to me. Can anybody share experiences of the use of the Calculator in the graphical modeler with me?

The documentation about longer function series in the modeler is not so good, that's why I didn't find any tutorial that uses the calculator.

How should I define my variables well?
I am using QGIS Essen as I changed my computer recently, can this be a source of problems and I should go back to QGIS Vienna stable version?


texts: What to consider when hiring aerial LiDAR survey?

texts:
My organization is considering hiring a company to capture new LiDAR data to help with current stormwater issues that we are currently experiencing. I am a GIS analyst that has been assigned with this task.

Our current quote is for roughly 70 square km and the project deliverables will include:

1m resolution orthophoto TIFF imagery.
50cm resolution orthophoto TIFF imagery.
Tiled 11cm resolution orthophoto TIFF imagery - each image
delivered in ATS sections.
LiDAR data collected with approximate 7 points per sq meter.
Raster bare earth DEM produced from the LiDAR data. Two files will be delivered:

1 - 50cm grid spacing.
1 - 1m grid spacing.

6 new survey checkpoints.

We have been quoted roughly $30000 CAD (canadian dollars).

Is cost of new LiDAR data reasonable for deliverables provided?
What else should be considered when hiring aerial LiDAR survey?

Any feedback would be great. I have never had the luxury of being included in data purchases so this is all new to me.


texts: Calculating a coordinate based on approximate distances from known locations

texts:
My goal is to calculate a position based on approximate distances from known nearby locations.

The distance observations include inherent error.  The number of known locations is often 12 or more.

Is there a approach for calculating a position that decreases the impact of the distance errors?


texts: Accessing GeoJSON from ArcGIS online rest API?

texts: I have some features published on ArcGIS online. I'm accessing it as a feature service. I need to use that feature service as a GeoJSON and display it on the map. Is there any sandbox example on how to display GeoJSON on to the map or can any one share an example code on it using ArcGIS javascript API?


texts: OpenLayers 3 with GeoWebCache - WMTS vs. WMS

texts:
I am using OpenLayers 3 with GeoWebCache (GeoServer). Are there any performance differences between WMS and WMTS when using GeoWebCache?

I'm currently using url http://localhost:8080/geoserver/gwc/service/wms? and I would like to know if it will be faster to use WMTS instead


texts: Autofill field parameter based on parent layer in Processing Provider

texts:
I'm building a custom Processing Provider and would like my field dropbox to be dynamically filled with the first valid field of the associated parent layer. Currently, the dropbox list is updated if I change the parent layer (great), but it doesn't select any field (keeps the box empty). Is is possible to query this same list and autoselect the first field in the list? I should add that I don't want to hard code the name of the field as it may change.
class customAlgorithm(QgsProcessingAlgorithm):
    INPUT_vector = 'INPUT_vector'

    def initAlgorithm(self, config):
        self.addParameter(QgsProcessingParameterVectorLayer(self.INPUT_vector,
                                                            self.tr('Vector Layer'),
                                                            [QgsProcessing.TypeVectorPolygon]))

        self.addParameter(QgsProcessingParameterField(self.INPUT_field,
                                                      self.tr('Selected Field'),
                                                      defaultValue = None,
                                                      parentLayerParameterName = 'INPUT_vector',
                                                      type = QgsProcessingParameterField.Numeric))

There is this other similar question, but it didn't really got a full answer (and was with QGIS 2.x).


texts: Cannot find Grass7 folder in QGIS 2.18 content folder

texts:
I have freshly re-installed QGIS 2.18.28 in the prescribed order.
I was under the impression Grass7 in included in one of the install packages.
But when trying to point QGIS to the grass7 folder under "processing options" in QGIS, I cannot find that folder.
See screenshot attached.


texts: How to store a password during one QGIS work session?

texts:
I'd like to get into qgis-plugin development, but am rather new to both python and QGIS development.

I intend to write a plugin which connects data from a mysql database to shapefiles. As the data is related to persons, I need a kind of login before accessing the data. To avoid the user from typing his password every time he invokes the script, I wonder if it is possible to retain values between calls of the script. In the python cookbook, I found global project settings - but that is not quite what I'm looking for as the values should be forgotten on closing QGIS.

I hope I made my intentions clear in my clumpy English and somebody could give me a hint with which way to proceed.


texts: gdalwarp UTM Zone 1 to WGS84 - fails to transform

texts:
I have a global dataset consisting of tiles in UTM projection that I want to mosaic and therefore need to reproject to a common projection.

Can I solve my reprojection error with  options to gdalwarp or do I need to manipulate the input files?

Reprojection is done with gdalwarp and the following command:

gdalwarp -t_srs EPSG:4326 -of GTiff -co "COMPRESS=LZW" N01_50.TIF N01_50_wgs84.TIF

Tiles in UTM Zone01 produce the following error and the result is visibly distorted (spans around the globe in E-W direction). All other tiles work fine.

ERROR 1: Too many points (441 out of 441) failed to transform,
unable to compute output bounds.
Warning 1: Unable to compute source region for output window 6006,0,6007,340,

Here is the output of gdalsrsinfo and the download link of the input file (7.2 Mb):

PROJCS["WGS 84 / UTM zone 1N",
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-177],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["meters",1],
AUTHORITY["EPSG","32601"]]


texts: Can I simulate the user interaction described with a series of scripts/tools/buttons?

texts:
I’m trying to find the best solution for adding/update a new/old record in a feature class in ArcMap 10.3 environment using a GUI. I’ve done some research for the last couple of days trying to identify the solution that will best fit my needs but because I haven’t got experience with developing for ArcMap I’ve ended up on this forum. My experience in programming is related to using VBA with Excel and AutoCAD and some Java and Python codding.
The main structure of the input is described in the picture bellow.

Basically it should collect some data based on a user input, get some additional input based on the initial input and add a new record in a polygon feature class.
Considering what I’ve understood so far it looks like I can use ArcObjects or a combination between Custom Script Tool and Python with a possibility of including PyQT as well . I’ve spent some time in trying to see if ArcObjects for .Net is the proper tool and looks like I will be able to get what I need but I was wondering if there is a better approach or if I can do it without ArcObjects. I am most worried about the select/create feature from the drawing with snapping capabilities and branching functionalities.


texts: QGIS 3.14 - How to make OSM layer Permanent as .shp with too many field and/or fields containing semi-colon : characters?

texts:
I imported an .osm file with QuickOSM and layers were created in QGIS (in duplicate but it does not matter). Now, I would like to make these layers Permanent with a rigth-click. But choosing .shp as format, it gives an an error:
Export to vector file failed.
Error: Creation of field ref:lukr:breytandi failed (OGR error: Can't create field ref_lukr15 in Shape DBF file, reason unknown.)
It is either probably caused by too many fields and/or field containing semi-colon : special character
How can I export to .shp?
Is there any way to count the total number of fields in a layer?
Is it possible to batch replace semi-colon : by an underscore _ in field names?
Thanks


texts: Given a bunch of pushpins in Bing Maps, is there a quick and easy way to zoom in on those pushpins such that they are all visible, and fairly spread out? For example, say I have hundreds of pins all over Australia. Then I would want to zoom in to Australia as a continent, but seeing other continents would be too big, yet seeing only Sydney would be too small.


texts: Join two SpatialPolygons keeping any polygons which overlap

texts:
The best way I can describe my problem is that I’m after a left inner join for two spatial objects???

As you can see below, I have two spatial objects: (1) State of California [in blue] and (2) Grid polygons [in red]. I am trying merge both objects but only for the grid polygons which "touch" the California object. E.g. keep the blue grid polygons and any red grid polygons which have blue in them.

I tried "raster::intersect", but this clips the coast/boundary.  Here's my code ....

library(rgdal)
library(raster)
library(rgeos)

# Download US shape information from data.gov
URL <- "http://www2.census.gov/geo/tiger/GENZ2016/shp/cb_2016_us_state_500k.zip"
td <- tempdir()
setwd(td)
temp <- tempfile(fileext = ".zip")
download.file(URL, temp)
unzip(temp)

# Get US shape file and filter states for just California
us_states <- raster::shapefile("cb_2016_us_state_500k.shp")
ca <- subset(us_states, STUSPS == "CA")

# Create a grid
grid <- raster(extent(ca))
res(grid) <- c(1)
proj4string(grid) <- proj4string(ca)
gridpolygon <- rasterToPolygons(grid)

# Merge objects together
merged <- raster::intersect(ca, gridpolygon)

plot(gridpolygon, add=TRUE, col='red')
plot(merged, add=TRUE, col='blue')  # This clips the grids


texts:
Just getting into QGIS and coming from MapInfo. I have some custom projections that I would like transferred. I've been trying to figure out the parameters today without any luck. I have the following from MapInfo.

> 8 = Transverse Mercator 62 = North American 1927 8 = US Survey Feet
> -115.5833333333002 = latitude Origin
> 34.75 = Longitude Origin
> 0.9999 = Scale
> 499999.999998 = false easting
> -999999.9999960001 = false northing

Thought this would work in QGIS but it doesn't.

+proj=tmerc
+datum=NAD27
+units=ft
+centerLat=-115.5833333333002
+centerLong=34.75
+k_0=0.9999
+x_0=499999.999998
+y_0=-999999.9999960001

What am I doing wrong?


texts: OpenLayers - when zooming out OSM layer "slips" East

texts:
I have some code that displays an OpenLayers map composed of various layers (vectors, markers and features) on top of an OpenStreetMaps basemap.  Everything worked fine until some months back.  Without any change in code, now, when i zoom my map, all superimposed layers lose their position relative to basemap/baselayer.

The OpenStreets basemap "slips" east as one zooms out.

Checking latlon pointer/mouse position reveals 4326 projection is working properly (as i zoom out, pointer is shown in a more westerly position, but retains same latlon pair value).

It is as if the basemap projected "wrong".   All 6 layers are in sync with one another: if one layer "slips" west, all layers, except basemap, "slip" along with it.

I have tested this using various browsers (chromium, chrome and firefox) and the behavior is identical regardless of the browser.

var proj4326 = new OpenLayers.Projection("EPSG:4326");
var markerslayer = new OpenLayers.Layer.Markers( "New Markers" , {
    styleMap: new OpenLayers.StyleMap({
        "default": {
            graphicYOffset: 16,
            graphicYOffset: -6
        }
    })
});

var layerOSM = new OpenLayers.Layer.OSM("Street Map");
var map = new OpenLayers.Map({
 div: "map",
 layers: [
    layerOSM,
    markerslayer,
    //layerKMLClient,
    //layerKMLTarget,
    //layerKMLTrack,
    //layerKMLStops
]
});

map.events.register("click", map, function(e) {
 var position = this.events.getMousePosition(e);
 var lonlat = map.getLonLatFromPixel(position);
 var lonlatTransf = lonlat.transform(map.getProjectionObject(), proj4326);
 var latTemp = lonlatTransf.lat;
 var lonTemp = lonlatTransf.lon;
 alert ('lonlatTrans=> lat='+lonlatTransf.lat+' lon='+lonlatTransf.lon+'\nlonlat=>'+lonlat+'\nposition=>'+position);
}


texts: QGIS 2.14 2.5d render appears different in composer and map views

texts:
I'm experimenting with 2.5d rendering for the first time.  It seems to me that the objects are being rendered in a different order between the main working environment and the map composer.

I'm wondering if others are having this issue and if it's a bug, or simply user error?


texts: Using custom icons for Leaflet Layer groups?

texts:
I have searched posts and tutorials and could not find the answer to my specific application of this feature

Almost all the examples use L.marker to simply place a marker on the map using various custom icons. This would be fine except I am not just placing markers but rather adding all the geoJson data to layers, then adding those layergroups to the control and then to the map. - I.E. I want all the layer routeca1 to use a red icon. I think I have narrowed this to needing to use the point to layer function:

pointToLayer: function (feature, latlng) {
            return L.marker(latlng, {icon: redIcon});
    }

I am not using the onEachFeature function, as I have seen an example on how to do this using that function but cant make that work for what I am doing here.

The map works great without trying to do a different color icon. I have all my points and polygons appearing and controlable by the layercontrol, but when I try to add the point to layer I just get a white map or the same map with no changes depending on where I try to add it in.

I added the icon sections below based off of the Leaflet Tutorial -

 var routeca1 = L.layerGroup();

 var geojsonFeature1 = [{***random geojson data points here****}  ];

    L.geoJSON(geojsonFeature1)
.bindPopup(function (layer) {
    return ''+'Location ID: '+ layer.feature.properties.locationid +'' + '
' + '
'+ layer.feature.properties.Description + '
' + layer.feature.properties.model + '
' + '
'+ layer.feature.properties.addresscomponents;
}).addTo(routeca1);

    var streets  = L.tileLayer('https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {attribution: '&copy; OpenStreetMap', maxZoom: 17, minZoom: 6});

    var map = L.map('map', {
        center: [34.068, -118.248],
        zoom: 10,
        layers: [streets, routeca1]
    });

var baseLayers = {
        "Streets": streets
    };

    var overlays = {
        "CA1": routeca1
    };

var MarkerIcon = L.Icon.extend({
    options: {
        shadowUrl: './leaflet/images/marker-shadow.png',
        iconSize:     [38, 95],
        shadowSize:   [50, 64],
        iconAnchor:   [22, 94],
        shadowAnchor: [4, 62],
        popupAnchor:  [-3, -76]
    }
});

var greenIcon = new MarkerIcon({iconUrl: './leaflet/images/marker-icon-green.png'}),
    redIcon = new MarkerIcon({iconUrl: './leaflet/images/marker-icon-red.png'}),
    blueIcon = new MarkerIcon({iconUrl: './leaflet/images/marker-blue.png'});

    L.icon = function (options) {
    return new L.Icon(options);
};

    L.control.layers(baseLayers, overlays).addTo(map);


texts: Show selected point in ol3 with maxscale layer

texts:
I have to display a layer with 14000 object. I use maxscale to limit zoom level displaying features. I want to select item on init but in the top level no item are loaded (maxscale).
Do you have a solution to show selected object on top level?


texts: Auto populate a field in the attribute table with unique IDs

texts:
I am a GIS specialist for an engineering company and I am setting up Arc Collector on ipads for field workers. The field workers go to residential properties and collect soil samples. Each property has an ID that starts with "PIA0" followed by a unique four digit number, such as "PIA01123". They will be using the software to create a point for each sample they collect. Each point will be labeled like so: start with the ID # and finish with a number based on the number of samples they collected. For example, the first point they create on a property with the ID "PIA00555" will have the label, "PIA00555-01", the next point they create will be "PIA00555-02", and so on. Usually they have 10 samples per property and thus would create 10 points in Arc Collector.

I am trying to make the field crews spend as little time as possible using the keyboard on the ipad for efficiency, so I don't want them to manually enter the ID every time they create a point. Is there a way to autopopulate the unique ID every time the field workers create a point? I have a shapefile that represents parcels of every property we are sampling that includes the address along with other information, but not the ID I mentioned earlier. My only solution to this problem I came up with is to enter the ID for every property and create a field that == the ID plus the sequential number, but since there are over 10,000 properties, I don't have time to do that. I feel like there is a python script that creates something like this, but my knowledge of python is limited.


texts:
I have a simple d3 map displaying data from Canada in topojson format, but the entire map is rendering sideways.  I have been playing with the rotation parameters but I can't seem to get it just right.  I'm wondering if anybody else has had this problem because it seems to me like there should be a very simple swap of an x and a y somewhere that would fix this, but I cannot find it since I am rather new to d3.

This is the code:

window.onload = initialize(); //start script once HTML is loaded

function initialize(){ //the first function called once the html is loaded
    setMap();
};

//set choropleth map parameters
function setMap(){

    //map frame dimensions
    var width = 1000;
    var height = 750;

    //create a new svg element with the above dimensions
    var map  = d3.select('body')
        .append('svg')
        .attr('width', width)
        .attr('height', height);

    //create Canada Alberts equal area conic projection, centred on Canada
    var projection = d3.geo.albers()
        .center([-56, 75])
        .rotate([0, 0, 0])
        .parallels([50, 70])
        .scale(750)
        .translate([width / 2, height / 2]);

    //create svg path generator using the projection
    var path = d3.geo.path()
        .projection(projection);

    //use queue.js to parallelize asynchronous data loading
    queue()
        .defer(d3.json, "data/canada_districts_4326.topojson") //load geometry
        .await(callback); //trigger callback function once data is loaded

    function callback(error, ridingData) {

        //add electoral boundaries geometry to map
        var ridings = map.append('path')
            .datum(topojson.feature(ridingData, ridingData.objects.canada_districts_4326))
            .attr('class', 'ridings')
            .attr('d', path)

        };
}

And this is the result:


texts: ArcMap: Apply symbology from layer, remove zero counts from values

texts:
We have to generate 100's of pdf's (from mxd) for 100's of shapefiles based on a column value "dummy". We have a template lyr file with predefined symbols for each category type in the "dummy" column. But when we use

  arcpy.ApplySymbologyFromLayer_management(input_layer,source_layer)

code, all the category type are getting listed in the ToC. As in the image, though "Borefield, Coreyard, etc..." has zero counts - still they are getting listed in the ToC and subsequently in the legend.


texts: Is there an export feature to shapefile with send to FTP model out there?

texts:
Is there a model already built for this simple process? Can ModelBuilder do this?

Export feature to shapefile
Delete fields
Rename shapefile
Post to FTP site
Run this process every Friday

It would be way easier if something like this was already built.


texts: GeoDjango 'NoneType' error when using 'distance' function

texts:
I have a site that's working on our staging server but after copying it to the live server I am getting the following error when trying to perform a distance() call (code cut down to minimum)

>>> Buddy.objects.all().distance(pnt[0].geom, field_name='postalcode__geom')
Traceback (most recent call last):
  File "<console>", line 1, in <module>
  File "/home/buddyup/sites/buddyup/lib/python2.6/site-packages/django/contrib/gis/db/models/query.py", line 100, in distance
    return self._distance_attribute('distance', geom, **kwargs)
  File "/home/buddyup/sites/buddyup/lib/python2.6/site-packages/django/contrib/gis/db/models/query.py", line 582, in _distance_attribute
    geodetic = geo_field.geodetic(connection)
  File "/home/buddyup/sites/buddyup/lib/python2.6/site-packages/django/contrib/gis/db/models/fields.py", line 131, in geodetic
    return self.units_name(connection) in self.geodetic_units
  File "/home/buddyup/sites/buddyup/lib/python2.6/site-packages/django/contrib/gis/db/models/fields.py", line 122, in units_name
    self._get_srid_info(connection)
  File "/home/buddyup/sites/buddyup/lib/python2.6/site-packages/django/contrib/gis/db/models/fields.py", line 108, in _get_srid_info
    self._units, self._units_name, self._spheroid = get_srid_info(self.srid, connection)
  File "/home/buddyup/sites/buddyup/lib/python2.6/site-packages/django/contrib/gis/db/models/fields.py", line 35, in get_srid_info
    units, units_name = sr.units
  File "/home/buddyup/sites/buddyup/lib/python2.6/site-packages/django/contrib/gis/db/backends/base.py", line 281, in units
    return (self.linear_units, self.linear_name)
  File "/home/buddyup/sites/buddyup/lib/python2.6/site-packages/django/contrib/gis/db/backends/base.py", line 253, in linear_units
    return m.group('unit')
AttributeError: 'NoneType' object has no attribute 'group'

>>> pnt[0].geom
<Point object at 0x1fb0370>

I've been unable to find any information elsewhere this error. Any help as to what might be causing this would be very welcome.

EDIT

Not having a clue what I'm doing, I've come up with possibly the best patch in the world /s

--- /home/buddyup/base.py       2012-01-18 22:17:04.000000000 -0500
+++ /home/buddyup/sites/buddyup/lib/python2.6/site-packages/django/contrib/gis/db/backends/base.py      2012-01-18 22:16:09.000000000 -0500
@@ -239,7 +239,8 @@
             return None
         else:
             m = self.units_regex.match(self.wkt)
-            return m.group('unit_name')
+            return 'meter'
+            #return m.group('unit_name')

     @property
     def linear_units(self):
@@ -250,7 +251,8 @@
             return None
         else:
             m = self.units_regex.match(self.wkt)
-            return m.group('unit')
+            return 'm'
+            #return m.group('unit')

     @property
     def angular_name(self):

I have no idea why it is unable to find the units or the units name, but for now that works for me. I would love to hear what the actual problem could be and find a valid solution.

EDIT 2

yum install gdal

After talking to pizzaparty in the #django IRC channel he pointed out that installed gdal should stop Django trying to get the units from a Point object that doesn't have any. That seems to be a much better solution.


texts: Setting color of symbol from attribute of graphic from ArcGIS API for Flex?

texts:
I have a feature layer used in a client flex application.
The feature layer is used to publish a feature class that has many attributes. Among those attributes there is one named "color". As its name indicates, this attribute contains the color of the feature. This attribute (color) is meant to be modified by the user.

The symbology of the feature layer is set by a renderer created on the MXML file.
I'd like to be able to set the color of each graphic of the feature layer, from its color attribute.

Is it possible to tell the mxml renderer to set the color of its symbol to the value contained in the color attribute?
I'd like to do that on the mxml side without using ActionScript.


texts:
I have an operation:

var warpTif = Gdal.wrapper_GDALWarpDestName(ecwTif, 1, dss, warpOptions,
                    new Gdal.GDALProgressFuncDelegate(GdalProgressFunc), $"Warping {vrtName}:");

that uses the warp option -wm 6000 that I'm expecting to throw an integer overflow error and stop execution as it would via command line utility. However, the operation continues with warptif being null therefore throwing null reference exceptions later down the line. How does one properly handle gdal errors in C#? I have tried this but nothing happens...

PS - Obviously -wm 6000 is wrong and should be something more appropriate but I would still like to know how to capture this error and log it's message.


texts:
I would like to use SQL to create a virtual layer in QGIS by clipping information from a large data set using a bounding box around my study area.

I have two questions.  Is this an efficient way to work with large data sets (i.e., more efficient than just loading a large shapefile)?

How do I write the query? for it?  I don't really know what I am doing, so far my query looks like this:

Select F* Where ST_Intersects(F.geometry,C.geometry)

At the moment, I keep crashing my computer with my queries.  I will post back if I figure this out.


texts:
I have 3D map in QGIS2threejs where I have vector line and for Object type I use wall and for style I set Other side Z expression to -2000 which represents depth of -2km.
This is now vertical but this should be in angle so is there are a way to manipulate data visible in attribute table or to change some settings that this is in angle like blue line on the image I draw?


texts:
I'm new to GRASS.
We want to delineate a watershed using the LIDAR data we have our hands on.
I tried to follow several procedures to do so, but despite all my efforts, I can only get it to select a single cell out of the DEM.

I'm using GRASS 7.0.5

What I've done:

Created the DEM using r.patch and the necessary tiles
Defined my work region (g.region)
r.patch the DEM to fix the bridges
Used r.watershed to create the Accumulation and Drainage direction layers (tried with 0, 5k and 50k thresholds)
Used the Drainage layer to r.water.outlet

While I eventually get a result, it is only a single pixel, although I made sure to be very (VERY) precise (5 digits precise...) on the stream coordinates. Picture of the output.

I believe part of the problem is that my DEM is huge... but I find it hard to believe that I'm the only one using this kind of data.

Any thoughts?


texts: I have a line .kml which is needed to converted to shape. And I need to make a relation between "Folder" and "Placemark". I have a common field (kml_id in "Folder" and kml_parent in Placemark). How can I join them ?


texts: QGIS:joinattributesbylocation gives empty shapefiles and QGIS:joinbylocation does not give summary of all fields

texts:
I am trying to make a spatial join to calculate the mean value of many variables (100+ columns) of census tracts by city. In order to do that, I would, in QGIS Desktop, run a spatial join (join attributes by location). I did it and it worked. However, I need to automate the process and use the command line "processing" commands ("processing.runalg...").

Even tough it works in QGIS and should work in command line with the following line:
processing.runalg("qgis:joinattributesbylocation", neighborsPath, censusPath, "['intersects']",1,"mean", 0, layerOutput)
The process works, but gives an empty shapefile (1kb size).

As suggested in many posts (Why does "Join attributes by location" fail to compute sums when used in a model?) and (QGIS Join Attributes by Location not working from Python script), I downgraded processing to v.2.2.0, and used the joinbylocation, the result only gives me stats for the first column of my joined shapefile and not for the 100+ columns I have.
This is the line that I am using:
processing.runalg("qgis:joinbylocation", neighborsPath, censusPath, 1, 'mean', 0,1, layerOutput)

I don't know if one of you encountered this problem and successfully fixed it.
(p.s. I tried all versions of Processing)
I hope someone can help me! (This is my first question on stackexchange, I hope it's well formulated enough)

Also, qgis:extractbylocation also gives empty shapefiles, for an unknown reason... So the cause might be the same in both cases.

The files are all in WGS84 CRS

Here is the code I am using:

import os
import sys
from qgis.core import *
from PyQt4.QtGui import *
from os import walk
app = QApplication([])
QgsApplication.setPrefixPath("C:/OSGEO4~1/apps/qgis", True)
QgsApplication.initQgis()
sys.path.append(r'C:\OSGeo4W64\apps\qgis\python\plugins')
ogr2ogr = r"C:\OSGeo4W64\bin\ogr2ogr.exe"
import processing
from processing.core.Processing import Processing
Processing.initialize()
Processing.updateAlgsList()

neighborsPath = r"C:\Users\JoeNavut\Desktop\tests\wgs84\corrected_geom\calgary.shp" #PATH TO NEIGHBOR FILE
censusPath =  r"C:\Users\JoeNavut\Desktop\tests\wgs84\corrected_geom\censusData.shp" #PATH TO CENSUS FILE

layerOutput = neighborsPath.replace(".shp", "_Criteria.shp")
layerInput = QgsVectorLayer(neighborsPath, "file_to_calculate", "ogr")
censusInput = QgsVectorLayer(censusPath, "file_for_calc", "ogr")
print layerInput.isValid()
print censusInput.isValid()
try:
    processing.runalg("qgis:joinattributesbylocation", neighborsPath, censusPath, "['intersects']",1,"mean", 0, layerOutput)
    print "Criteria calculations complete"
except:
   print "Criteria calculations failed"
QgsApplication.exitQgis()
QApplication.exit()

And in QGIS, we see that the table has been created, but is empty...


texts:
I have just created some tiles on my project on QGIS with the plugin Qtiles; now I get several tiles in .png!
My question is, how could put these tiles .png on QGIS Server? Is it possible?


texts: Duplicating existing polygons based upon number in attribute field using ArcPy gives AttributeError?

texts:
We are working with parcel data that includes a single polygon for each individual tax parcel. We need to develop a list that includes an attribute for each residential address. Some tax parcels could have multiple residential addresses, e.g. condos, apartments, etc. Although this could be done in a spreadsheet, it would likely be time consuming. The though is to copy each polygon x number of times dependent on the number of units in order to begin this list, this would allow for all the attribute data to be copied as well. In researching a way to do this is we came across Exploding features based on numerical field in QGIS?. I provides a Python script and minimal directions to run the script.The script is as follows:

layer.startEditing()
for f in layer.getFeatures():
    # Create X-1 new features since we already have the original
    for i in range(f['X']-1):
        layer.addFeature(f)
layer.commitChanges()

I am actually trying to replicate this in ArcGIS. I have translated this to the following:

import arcpy
mxd = arcpy.mapping.MapDocument(r"G:\LUCA_2020\WorkingMap_4.19.mxd")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
lyr = arcpy.mapping.ListLayers(mxd, "newburgh_city_parcels_master_april2012_joinRPS", df)[0]
lyrfile = arcpy.mapping.layer(r"G:\GIS_Chad\Scratch\Delete\test.lyr")
arcpy.mapping.UpdateLayer(df, lyr, lyrfile, True)
layer.startEditing('newburgh_city_parcels_master_april2012_joinRPS')
for f in layer.getFeatures('total_units'):
    # Create X-1 new features since we already have the original
    for i in range(f['X']-1):
        layer.addFeature(f)
layer.commitChanges('newburgh_city_parcels_master_april2012_joinRPS')

When I run the script in the Python window, i get the following error:

Runtime error
Traceback (most recent call last):
  File "<string>", line 5, in <module>
AttributeError: 'module' object has no attribute 'layer'

I assume there may be other, put I am not sure. I am looking for guidance on correcting and running the script in ArcGIS 10.3.


texts: How to get the number of objects at a pixel position?

texts: Is it possible to get the total number of objects (Lines,Points) at a given pixel position in the map?


texts: Getting parent layer of feature in QGIS PyQGIS Custom Function?

texts:
I am writing a custom function in PyQGIS to use in expression editor: the aim of the function is to retrieve the first feature that contains the current feature from a polygon layer.

def get_containing_feature(value1, feature, parent):
    cont_feature = None
    cont_layer = QgsMapLayerRegistry.instance().mapLayersByName(value1)[0]
    crsSrc = QgsCoordinateReferenceSystem(4326)    # WGS 84
    crsDest = cont_layer.crs()
    xform = QgsCoordinateTransform(crsSrc, crsDest)
    curr_geom = feature.geometry()
    curr_geom2 = xform.transform(curr_geom.asPoint())
    features = cont_layer.getFeatures()
    for feat in features:
        if feat.geometry().contains(curr_geom2):
            cont_feature = feat
            break
    return cont_feature

The layers aren't in the same crs so I need a transformation
The polygon layer is accessed by name and I can get its crs using .crs() method, but I don't know how to get the feature's layer and crs (currently EPSG: 4326)


texts:
With QGIS, how can we analyse an electrical network to find loops?
Using Python is not a problem, however there is not a solution in the QGIS Python API.

The QgsGraphAnalyzer has dijkstra and shortestTree functions, but no loop finding function (QgsGraphAnalyzer Class Reference).

Example:
The image below includes 2 vector layers with wires that connect electrically. The loop is an error - it can't exist in a radial distribution network.
Is it possible to find these loops using QGIS and Python?


texts: Using GeoWebCache to cache Google maps tiles

texts: I'm new to GIS and I would like to know if anybody could give me a quick example of using GeoWebCache with Google maps tiles.


texts: Joining and Updating PostGIS data using QGIS

texts:
I have a PostGIS database with which I am currently using QGIS to manage the layers. Say I have a layer in PostGIS that I want to do a join to to add new data in. Is it possible to do this join in-situ without having to download the layer locally, do the join and re-upload the layer, overwriting the previous one?
If this isn't possible, is the best alternative to write a script to handle the downloading, joining and re-uploading all together?


texts:
I want to create a multibuffer view based on a point view in Postgis AND view the result in QGIS.

The points are in a view called "location"

To create the multibuffer I create a sequence in Postgis "multibuffer_vid_seq". This will be used to have unique ID.
Creating the multibuffer

 SELECT nextval('multibuffer_vid_seq'::regclass) AS bid,
  loc.id,
  st_buffer(loc.geom, 100::double precision) AS geom
 FROM location loc
 GROUP BY loc.id, loc.geom
UNION ALL
 SELECT nextval('multibuffer_vid_seq'::regclass) AS bid,
  loc.id,
  st_difference(st_union(st_buffer(loc.geom, 500::double precision)),
  st_buffer(loc.geom, 100::double precision)) AS geom
 FROM location loc
 GROUP BY loc.id, loc.geom
UNION ALL
 SELECT nextval('multibuffer_vid_seq'::regclass) AS bid,
  loc.id,
  st_difference(st_union(st_buffer(loc.geom, 1000::double precision)),
  st_buffer(loc.geom, 500::double precision)) AS geom
 FROM location loc
 GROUP BY loc.id, loc.geom;

This seems to work fine: I end up with a view  having a unique ID, the location ID and the multibuffer geom named geom.

QGIS sees the layer and I can bring it in the layer panel. However no geometry show up. It is empty.

I don't understand why so as doing a select in Postgis does return some proper records with a geom.


texts:
I'm trying to get the extent of a point CSV file, but I can't do it.
Before in my script making, I did it with shapefile with :

reproj_1 = processing.runalg('qgis:reprojectlayer', point, 'EPSG:2154', None, progress=None)

enveloppe = processing.runalg('qgis:concavehull', reproj_1['OUTPUT'], 0.25, False, False, None, progress=None)

layer = processing.getObject(enveloppe['OUTPUT'])

ext = "%f,%f,%f,%f" %(layer.extent().xMinimum(), layer.extent().xMaximum(), layer.extent().yMinimum(), layer.extent().yMaximum())

And It work.
But now i'm using csv and it returns me this classical error : AttributeError: 'NoneType' object has no attribute 'extent'

Somebody have an idea ?


texts: PostGIS raster query "out of range raster coordinates"

texts:
I have a large raster (national scale) that I have loaded into PostGIS using raster2pgsql and a 100x100 tile size: raster2pgsql -I -C -s 27700 my.tif -t 100x100.  I am now trying to query a large number of points to find the raster values at each location (values are held in a separate decode table keyed by the pixel value), my query is as follows:

INSERT INTO raster_values
SELECT p.id, d.*
FROM my_point_table p
JOIN my_raster r ON ST_Intersects(p.the_geom, r.rast)
JOIN my_raster_decode d ON (ST_Value(r.rast, p.the_geom) = d.value);

However when I run the query every so often I get a warning from PostGIS:

NOTICE:  Attempting to get pixel value with out of range raster coordinates: (28, 100)
NOTICE:  Attempting to get pixel value with out of range raster coordinates: (64, 100)
NOTICE:  Attempting to get pixel value with out of range raster coordinates: (100, 1)
NOTICE:  Attempting to get pixel value with out of range raster coordinates: (100, 17)

The presence of a 100 in every error would indicate there is an issue with the bounds I am using but I am struggling to work out what I have done wrong.  I have checked the SRID of both the point and raster table and confirmed that these are consistent (27700):

"enforce_srid_rast" CHECK (st_srid(rast) = 27700)
the_geom         | geometry(Geometry,27700)


texts: A graduated-categorized hybrid styling of features

texts:
I am trying to figure out how to do a styling of features on a mapi qgis, that is a sort of hybrid between the categorized and graduated styling.

I am doing a map of road trees in a city, and each tree is tagged with the species of the tree. There are about 40 different species among the 17000 trees in the attribute table.

What I would like to do is, instead of doing a simple color-coded categorized map, i would like to style the trees so that the least common species in the attribute table is the lightest in color, and the most popular species is the darkest in color.

Is this achievable?


texts: Is it possible to use QGIS to convert a vector that contains points/polygons into a csv file and, if yes, in which way?


texts: Counting points within 100km Buffer of country borders using ArcGIS Desktop?

texts:
I have a dataset composed of points and projected on a shapefile describing European countries.
My goal is to create a new shapefile, with zones specified 100 kilometers around the borders in-between all the States, since I want to see if there are more points close to the borders, or not.
I thought that Buffer would be a nice way to do it but did not manage to find out how to practically do this. I converted the polygon shapefile once to line, and then used the function Buffer, but there is nothing appearing. I did the same without converting before.


texts:
I am unable to create custom CRSes in QGIS. The 'save' button doesn't seem to have any effect. My best guess is that I don't have permission to save in the folder or database where custom CRSes go. I have looked at the screenshots of other users who are having trouble with this, and I usually see "1 of 1" in the custom CRS browser. On my system, I see "* of 0" instead

Any assistance appreciated.

QGIS version: 1.8.0

Ubuntu 12.04 LTS


texts: Deleting ArcGIS Project created using ArcGIS Pro?

texts:
I have deleted a few ArcGIS Projects, created using ArcGIS Pro, by using Windows Explorer to delete the folder that contains them.

However, I thought that, as long as the project (not the current project) was not already open somewhere else, then it may have been possible to use the Project pane to delete it.

I added a folder connection in the Project pane to a folder above the project that I wanted to delete but when I right-clicked on the project there is no Delete option offered (and hitting the Delete key while it is highlighted does nothing either).

Is the recommended way to delete an ArcGIS Project (created using ArcGIS Pro) simply to use Windows Explorer to delete its folder?

Some research and results subsequent to asking this question appears below.

When I tried to use the Delete tool from the Geoprocessing pane to delete a project (not the one open, nor one that had been used in the session, but one that could be easily deleted from Windows Explorer after closing ArcGIS Pro) it resulted in an error:

  ERROR 000464: Cannot get exclusive schema lock.  Either being edited
  or in use by another application. Failed to execute (Delete).

The other complication is that using Windows Explorer works fine for deleting a project for which a folder was created but, when that checkbox is not chosen at the time of project creation, the project parts can be placed in an existing folder amongst many other files and subfolders. I thought there might be something that would use the *.aprx to identify and delete its project parts.

I don't yet understand the use case for when you would want to leave that checkbox unticked and potentially have multiple projects "free-floating" in a single folder. Until/unless I do, I'll be recommending to always leave the checkbox ticked and to delete project folders using Windows Explorer.


texts:
I am trying to understand the scale parameter in earth engine better. My goal is to calculate monthly sums of precipitation across a geometry. Here is the python pseudo code:

area = ee.Geometry.Polygon(
        [[[-120.29277343749999, 37.319672371647655],
          [-118.51298828124999, 35.13985393638227],
          [-117.10673828124999, 35.552070101618895],
          [-118.46904296874999, 37.73787899128162]]])

col = ee.ImageCollection('OREGONSTATE/PRISM/AN81m')
t = col.filter(ee.Filter.calendarRange(2014,2014,'year'))
    .filter(ee.Filter.calendarRange(4,4,'month'))
    .filterBounds(area).select("ppt")

result = t.getRegion(area, scale).getInfo()

However this gives very different sums depending on the scale parameter value that is supplied. Here is a plot for the image sum as a function of scale parameter:

Is the best practice to reduce as sum before calling getRegion and getInfo? I have read the documentation on scale in earth engine, and the native resolution suggested by the dataset is 2.5 arc minutes or approximately 4.5 km. I am curious what the best practices are for supplying a scale parameter, especially when working with multiple data sources over the same area.


texts: Adding raster image to Raster field with Python?

texts:
I have a feature class with with field that is of Raster data type.

I want to add an image to each record in the feature class.

Is it possible to use Python to add the image to the Raster field for each record?

I've been searching but can't find anything.

Thanks.


texts: Discrepancy in the measurement of distance(measure line) in QGIS

texts: I am and my fellow colleague both of us using same version QGIS i.e., QGIS 2.8.3 and same CRS - NAD27 EPSG:4267, the problem is when we measure certain length of the feature(Using measure line option in QGIS), in her QGIS it displays say for ex: 300m and in mine the same length will be displayed as 450m. I crosschecked the results with some manual measurement the one displayed as 300m is correct. Also, we figured that this discrepancy is only for x- co-ordinate( Horizontally). Why such discrepancy even though we are using same version and same CRS. What settings can help us to get same and correct resutls?


texts: Using submodel within model in ArcGIS Desktop

texts:
I have two databases 175 shapefiles each I want to clip each one of them using the corresponding shapefile from the second database .. as I can't use two iterators I have to do a nested submodel within a model.

How can I do that?


texts:
I am running PostGIS on an Ubuntu server. I am able to add the Topology extension via CREATE EXTENSION, but when I try to use it as anything other than an admin I am denied permission.
I ran GRANT USAGE ON SCHEMA topology TO public, as suggested in this chain. When I run CreateTopology() I get a permissions error like the following:

ERROR:  permission denied for sequence topology_id_seq
CONTEXT:  PL/pgSQL function topology.createtopology(character varying,integer,double precision,boolean) line 17 at FOR over SELECT rows

I'm guessing this means that I need to monkey with more permissions on the topology schema. However,

That seems an odd design flaw. Shouldn't the Topology extension automatically take care of permissions when it is created, just like the PostGIS extension does for the spatial_ref_sys table?
I'm not sure how far down the permissions rabbit hole I need to go in order to make the extension usable.

Has anyone else encountered this issue?
Here's the output of SELECT postgis_full_version() from the regular user:
"POSTGIS="2.1.5 r13152" GEOS="3.4.2-CAPI-1.8.2 r3921" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.11.2, released 2015/02/10" LIBXML="2.9.1" LIBJSON="UNKNOWN" TOPOLOGY RASTER"


texts: QGIS Server Capabilities Response doesn't contain any Layers

texts:
This question is similar as this one QGIS Server not working (problem with fast CGI), but with more details:

I installed the latest QGIS from trunk on Linux (Ubuntu 11.04) in the non-default location '/usr/local/qgis/'. I installed apache2 from the repository.

Next,I followed the steps outlined in this post on Linfiniti's blog (http://linfiniti.com/2010/08/qgis-mapserver-a-wms-server-for-the-masses/), but probably I did not understand the part  about telling apache about the library path for my CGI; where do I add the ScriptAlias code to?

I created a qgis project with two vector (shapefiles). Under WMS server tab of the project properties I checked 'service capabilities' and gave the title 'trial'. I copied the qgs project file to the /usr/lib/cgi-bin/VECEA directory.

Trying http://localhost/cgi-bin/VECEA/qgis_mapserv.fcgi?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities in the browser gets me the 'capability' xml file. It contains the user information, but there is no information on the data layers (the qgis project has two shapefile vector layers) and the bounding box is not defined (all set at 0), is that right (I did define an 'advertised extent' in the qgis project under the WMS tab)

In QGIS I can connect to the server. However, it just gives me one row with the name of the project (ID 0), whereas I would expect it to show the two vector layers that were part of the project. When adding the layer, it is empty.

Edit: and now it stopped working altogether. When trying the getCapabilities, I am getting a 500 Internal Server error. Looking in the apache2 error log I am finding the following:

[Tue Sep 27 20:53:03 2011] [error] [client 127.0.0.1] Premature end of script headers: qgis_mapserv.fcgi
/usr/lib/cgi-bin/VECEA/qgis_mapserv.fcgi: error while loading shared libraries: libqgis_core.so.1.8.0: cannot open shared object file: No such file or directory

Edit: Got the 'getCapabilities' part to work again. In http://lists.osgeo.org/pipermail/qgis-user/2011-March/011549.html, it was suggested to add LD_LIBRARY_PATH to /etc/apache2/mods-> enabled/fcgi.conf, adding the line "DefaultInitEnv LD_LIBRARY_PATH /usr/local/qgis_git/lib". Only, I had to edit '/etc/apache2/mods-available/fcgid.conf'

Still stuck with no layers showing when connecting to the server in QGIS


texts: What is the difference between DEM data sources?

texts:
I have started experimenting creating basic 3D elevation models in QGIS recently. It turned out for me that there are several free sources of DEMs, like USGS's EarthExplorer, CGIAR-CSI's Digital Elevation Database and also EU-DEM from the Copernicus Program.

As far as I could have found out, the best resolutions of these data are the following:

SRTM elevation data from USGS's EarthExplorer site is available in 30 meter resolution - 1 arc second
SRTM elevation data from CGIAR-CSI site is available in 90 meter resolution - 3 arc second
SRTM elevation data from EU-DEM site is available in 25 meter resolution - 1 arc second

Is there any significant difference between these data sources and their accuracy?
Especially between USGS EarthExplorer and EU-DEM? I would like to use data of Europe and european cities so when I say difference, I mean the resolution, accuracy or quality and not the geographic coverage.

(Basically, I would like to find out which source is the best choice for a 3D model and I cannot see the difference between the American and the European DEM)


texts: Is D8 flow accumulation algorithm still suitable for watershed computation if terrain contains lots of flat area?

texts:
My question is inspired ( and similar) with my previous question. But this question is more concerned about the suitability of D8 ( or in general non-dispersive) flow accumulation algorithm for terrain with a lot of flat area, such as highly developed township.

For watershed analysis, the post here says that D8, as a non-dispersive algorithm, is better over other dispersive algorithms, such as D-infinity algorithm:

  Non-dispersive flow algorithms are particularly suited to watershed
  delineation for this reason. With a dispersive flow algorithm, grid
  cells near the boundary of a watershed can actually belong to both
  watersheds, i.e. two of the multiple flowpaths that issue from a grid
  cell situated near a watershed boundary may actually terminate at
  different outlets entirely. In reality of course watershed boundaries
  are somewhat fuzzy and ambiguously defined, but we generally want
  delineated watersheds that are mutually exclusive. Since the fuzziness
  caused by the flow dispersion inherent in the flow algorithm isn't
  necessarily related to the physical processes that actually govern the
  ambiguity of flow near watershed boundaries, we have no basis for
  interpreting the dual membership of bordering grid cells.
  That is, we can't say, 'this grid cell belongs to this watershed 90%
  of the time and its neighbouring watershed 10% of the time' solely
  based on the artificial nature of the dispersion imposed by the
  algorithm design. Dispersion is particularly problematic therefore.

But this is in direct contradiction to the answer I received:

  D8 is very sensitive to errors in the input digital elevation model.
  All the accumulated flow is directed to one neighbor, and sometimes
  the pixel is erroneous, which mean all the flow accumulation can be
  directed to a wrong pixel. Other algorithms, which distribute flow to
  more than one neighbor, are more robust. For example consider an area
  of 10 by 10 m (10 m resolution), quite often the water is not flowing
  to one direction, but to two or even more. I'd prefer multiple flow
  direction algorithm (by Quinn, Freeman), or the Tarboton's D infinity.

So which is correct? is D8 ( non-dispersive) flow accumulation algorithm  still suitable for watershed computation, especially if terrain contains lots of flat area?


texts: Iterating through a user-specified field, extracting unique values to separate feature classes

texts:
I'm trying to iterate through a user-specified field, extracting the unique values of that field into a list. I then want to loop through that unique list, doing a select by attributes and then extracting that group into separate feature classes.

My current script reads as follows:

import arcpy, os
from arcpy import env

env.overwriteOutput = True

#set input parameters
inputFeatures = arcpy.GetParameterAsText(0) # input feature layer
fieldName = arcpy.GetParameterAsText(1) # selecting features field name
outputWorkspace = arcpy.GetParameterAsText(2) #output .gdb workspace
outputFilePrefix = arcpy.GetParameterAsText(3) # string prefix for outputs

env.workspace = outputWorkspace

values= []

with arcpy.da.SearchCursor(inputFeatures, fieldName) as cursor:
    for row in cursor:
        values.append(row[0])
        uniqueValues = set(values)
        valueList = list(uniqueValues)

inputLayer = arcpy.MakeFeatureLayer_management(inputFeatures, "inputLayer")

arcpy.AddMessage("...Extracting unique values: ")
with arcpy.da.SearchCursor(inputLayer, "*") as cursor:
    for selection in valueList:
        arcpy.AddMessage(selection)
        arcpy.SelectLayerByAttribute_management(inputLayer, "NEW_SELECTION", "\'fieldName\' = \'%s\'" %(selection))
        arcpy.CopyFeatures_management(inputLayer, outputFilePrefix + str(selection))

The script executes without any traceback errors, however nothing is actually in the output feature classes, they are empty. Could anyone suggest why this might be the case?


texts: GeoTiff.lib: TiePoints only without ModelPixelScale do not get georeferenced in QGis

texts:
I have an image that has been georeferenced (geotiff.lib) by more than 3 tie points (x y 0 long lat 0 ...) and added following geotiff keys

    GTIFKeySet(gtif, GTModelTypeGeoKey, TYPE_SHORT, 1, ModelTypeGeographic);
    GTIFKeySet(gtif, GTRasterTypeGeoKey, TYPE_SHORT, 1, RasterPixelIsArea);
    GTIFKeySet(gtif, GTCitationGeoKey, TYPE_ASCII, 9, "GEODETIC");
    GTIFKeySet(gtif, GeographicTypeGeoKey, TYPE_SHORT, 1, GCS_WGS_84);
    GTIFKeySet(gtif, GeogCitationGeoKey, TYPE_ASCII, 7, "WGS 84");
    GTIFKeySet(gtif, GeogAngularUnitsGeoKey, TYPE_SHORT, 1, Angular_Degree);

without using the TIFFTAG_GEOPIXELSCALE tag entry.

If I open it in Quantum GIS the scale is wrong and it does not get georeferenced. There is no difference if I take the corner points of the image as tie points or randomized points of interest.

If I open the same image in e.g. GlobalMapper the image gets the correct position and the scale is correct.

As per advice from GeoTiff Specification it should hopefully be a correct GeoTiff:

 Case 1: The model-location of a raster point (x,y) is known, but not
      the scale or orientations:

        Use the ModelTiepointTag to define the (X,Y,Z) coordinates
        of the known raster point.

 Case 2: The location of three non-collinear raster points are known
      exactly, but the linearity of the transformation is not known.

      Use the ModelTiepointTag to define the (X,Y,Z) coordinates
      of all three known raster points. Do not compute or define the
      ModelPixelScale or ModelTransformation tag.

Is there any missing entry for QGis? How would you suggest to georeference images by using geotiff.lib without knowing the transformation parameters?

Thanks for helping!


texts: osm2pgsql: Extracting full buildings polygons + addresses text from OSM

texts:
I'm trying to extract all the polygons for the buildings in my country, however, it seems osm2pgsql is extracting very little information compared to what I need.

My idea was that I will rely on the empty.style file and I've set everything to have the "delete" flag, and keep only the tag=[building=yes], importing the building footprints should work. A brief extract from the style file:

# OsmType   Tag                     Type    Flags
# Insert your own columns here, or change polygon,nocolumn to polygon below
....
node,way    amenity                 text    delete
node,way    building=yes            text    polygon,nocolumn
....

The problem with this is I am seeing more unwanted polygons (administrative delimitation in Luxembourg) and less of the needed polygons (all the "mapped" buildings footprint polygons).

This is how my file looks like atm: http://pastebin.com/hwtYth22

This is how things look where the tag building=yes, notice the administrative areas polygons making it into the database, but also the lack of all the known polygons:

If I revert from "building=yes", and I leave it only "building", I get more polygons, but still, not all the buildings, yet more compared to using the "building=yes" tag.

Any advice on what I can do to:

prevent the administrative boundaries from showing up in postgres
have all the building footprints in postgres so I can filter them out by Office/amenity/etc.

TY,
Mark


texts:
I'm developing a QGIS plugin, and under certain circumstances, it causes QGIS itself to crash. How can I debug this, given that I am on Windows 7? I've tried using the QGIS nightly dev build (which the QGIS dev site said had debug enabled), but dbgview.exe could not read any debug output from the dump file. How can I debug this further?

The issue in question is here: https://github.com/tomchadwin/qgis2web/issues/32


texts: Loading Australian Electoral Commission data (MIF/MID) into QGIS for Ubuntu?

texts: I have installed QGIS for Ubuntu (linux) and downloaded some GIS electoral data from the AEC.  When I try to open the files (they have MID and MIF extensions), the application does not recognise the files.  What can I do to open these files in QGIS?


texts: QGIS - Problem with Print Layout - blacked out sections on map (QQIS 3.2.1-Bonn)

texts: I am experiencing a problem with  Print Layout. When importing my map, random blacked out sections are appearing.


texts: Count the number of polygon features within a polygon

texts:
I am using ArcMap 10.2.2 (Advanced license, no toolbar add ons/no spatial analyst) on a Windows 7 computer (i7 processor, 8GB RAM).  I have used ArcMap for other basic analyses, but am fairly new to the program and have only very basic coding skills.

I am trying to calculate the number of mammal species (TerrestMamm_diss on the image below/light red outlines) that are found within the geographic range of a subset of those species (all species in the order Carnivora, or Carns1 on the map/dark green outlines, one of which is selected).  My ideal output would look like the attribute table below plus an additional field indicating the count (i.e. FID, Shape, binomial, mammal_count).

There are ~250 features on the Carns1 layer and ~5,300 on the TerrestMamm_diss layer (whose attribute table looks just like that of Carns1 shown in the image, since Carns1 is just a subset of TerrestMamm_diss).  Both layers have overlapping polygons.

I do not need to keep the identification/binomials of the mammal species found within each Carnivora range, I just need the total count of species.  I've searched on this forum and elsewhere and I've found ways to count point features within a polygon (and I'm familiar with the model and procedure to count overlapping polygons using the "spaghetti and meatballs" method), but I'm not sure how to do this.

My initial thought was to turn all the terrestrial mammal species (TerrestMamm_diss) into points (feature to point tool) and count them that way, but I realized putting a single point within a large species range would mean that point wouldn't necessarily lie within each of the overlapping polygons.


texts: QGIS Get maximum sum of a feature subset in Atlas-View

texts:
I have a data-layer called data wich looks like this:

I also have an atlas-layout with pagenames defined by the column Jahr from the another Layer, representing the unique values from the Data-Layer.
In the atlas-layout I have a text-field, which shall contain the following, build up by an expression: the maximum of all sums from the column Anzahl, grouped by Datum, on each atlas-page.
So the feature-subset for each atlas-page would be all Datum-fields where Jahr= atlas_pagename. From this subset the maximum sum of Anzahlshall be displayed.
For Exmaple on the atlas-page 2014 the result would be "165"( the sum of all values from the column Anzahl from the Datum 2014-10-01, which is the highest result in this example for 2014).
I already tried the expression sum("Anzahl", group_by:="Datum"), which works well in the field calculator, but I cannot manage to access the Data-Layer from the expression-builder in the atlas-view.


texts: Finding geographic median in ArcGIS Desktop?

texts:
Is it possible to split the total population in 501 polygons, each polygon with a different number of inhabitants into two North/South sub-areas with roughly the same number of inhabitants in each sub-area?

The split line, trending east-west should be along administrative boundaries of the polygons (it is called geographic median).i.e. no polygon divided/shared between sub-areas.

The picture shows the region with 501 villages (polygons) with different numbers of inhabitants each.


texts: Prevent QGIS user to change layers without saving?

texts:
I'm looking for a way of preventing/warning users to select another layer without saving the current one they were using.

The only potential solution I found is an autosave plugin but it won't exactly fit to my question.

I'm using QGIS las palmas.


texts: Arcpy field calculation with user selected field

texts:
I could really use some help with an arcpy field calculation.I need to concatenate a float and a string to populate a TEXT field with distances and units to feed to the buffer function. I will need to buffer multiple points at various distances but with a single unit of measure.
I am reading in the field name from the script tool, and I'm allowing the user to select a units string from a value list so I do not know what the actual field name or unit string will be.

Dist_Field = arcpy.GetParameterAsText(0)
    # In my test data this field is called "Dist"
    # and is sourced to an input table in the script tool
Units = arcpy.GetParameterAsText(1)
    # selection from a value list of valid unit strings -
    # METERS, MILES etc...

# What goes in the expression argument?
arcpy.CalculateField_management (CopiedTable, "DistUnits", ? ,"PYTHON")
    # DistUnits is a TEXT field that I created so I know its field name

    # I attempted to use the following as the expression
      # "'str(Dist_Field)+" Kilometers"'"
      # I don't have the actual line of code handy, but variations on this
      # return the field name + Kilometers:  ( Dist Kilometers )
      # I've also tried "'str(Dist_Field) + str (Units)'"
    # a literal field name string or units string is no use to me as I only
      # know the variable name

How do I concatenate a float field with a string variable and write it to the new TEXT field when I only know variable names?

I'm flailing in my attempts figure out how this works.

This is my first attempt at both Python and ArcPy.


texts: Display two-dimensional color gradient for two different quantities

texts:
I have two overlapping raster layers. One contains quantity A from 0 to 100, and the other contains quantity B from 0 to 100.  I'd like to create a map displaying both sets of data at once based on a color box with two different axes.

So, for an area where:

A is 100 and B is 100, the map would be red
A is 100 and B is 0, it would be yellow
A is 0 and B is 100, it would be green
A is 0 and B is 0, it would be blue

For intermediate values the hue would change accordingly, so A=100, B=50 would be orange; A=50, B=0 would be blueish-green.

I have seen maps display this type of data before, and I have a couple of things I'd like to display this way.  What is the easiest way to do this? I can use ArcGIS, QGIS and R.


texts:
I'm trying to render images of a dome lit from multiple angles in GRASS using the m.nviz.image render plugin.

Whilst the DEM displays nicely in the map display, I'm getting weird fractal patterns present on the rendered image (see shot)

They seem to be a function of:
 1. the rounding of the elevation values of the dome
 2. the region resolution at the time of running m.nviz.image
 3. the actual resolution of the dome raster map when produced

I've had varied success but would like to know what is causing them and the optimum combination for the cleanest output. I attach my python script for generating a dome DEM and nviz image below

import math

size = 200
res = 0.1
gridres = 10

grass.run_command('g.region',n=size, s=0, e=size, w=0, res=res)

# make centre point vector
grass.write_command('v.in.ascii', input='-', output='centre', \
                  stdin=str(size/2)+'|'+str(size/2), overwrite=True)

# rasterise point file
grass.run_command('v.to.rast', overwrite=True, input='centre', \
                  output='centre_rast', use='cat')

# max no. of buffers
numdiv = float(100)
numdist = float(100)
# no. of buffers required
iterations = float(numdist / numdiv)

# raster buffer distance list creation
buffdistlist = []
currentdist = 0
for dist in range(1, int(numdiv) + 1):
    currentdist = dist * iterations
    buffdistlist.append(str(currentdist))

# generate buffers
bufferedpoint = 'centre_rast_elev_buff'
grass.run_command('r.buffer', overwrite=True, input='centre_rast', \
              output=bufferedpoint, distances=buffdistlist)
# make dome
formula = 'dome = if(isnull({0}),null(),{1} - ({2} - sqrt(({3}*{4}) - ({5}*{6}))))' \
          .format(bufferedpoint,numdiv,numdiv,numdiv,numdiv,\
                  bufferedpoint,bufferedpoint)
grass.raster.mapcalc(formula, overwrite=True)

grass.run_command('m.nviz.image',elevation_map='dome', color='black',\
                  position='0.5,0.5',height='1200',twist='90', perspective='90', \
                  shininess_value='255',\
                  light_position='0.68,-0.68,0.8',light_brightness='90',\
                  light_ambient='0', output='test', format = 'tif',\
                  size='4840,4840',resolution_fine='1', mode='fine',\
                  style='surface',shading='gouraud',resolution_coarse='1')


texts: Coordinate transformation issue in postgis

texts: I have a cad drawing that is suppose to be in FL83-EF that I imported into opengis. However I think that the coordinates were wrong. Is there some sort of transformation that I can run on my data given I know some points on the cad drawing and can use as points of reference.


texts:
In QGIS I'd like to add scale to my PyQGIS image exporter.
Similar to the Project -> Import/Export -> Export map to image
where you can add the scale to the image as in the picture below.

vLayer = QgsProject.instance().mapLayersByName("a_g")[0]
sExtent=vLayer.extent()
sOutputDpi= 300
sZoomScale = 2500
sQSize=QSize(10000, 10000)

settings = QgsMapSettings()
settings.setOutputSize(sQSize)
settings.setExtent(sExtent)
settings.outputDpi = sOutputDpi
settings.scale = sZoomScale

render = QgsMapRendererParallelJob(settings)

def finished():
    img = render.renderedImage()
    img.save(image_location, "png")

render.finished.connect(finished)

render.start()

In the output I'd like to include the scale of the map as well. At the moment the image exporter exports the image scaleless (I think) and only takes note of the setExtent value.
Does anyone know if it's possible to include scale in this script?


texts: Running PyQGIS algorithm on selected features in layer

texts:
I try to run algorithms on layers but only on selected features with QgsProcessingFeatureSourceDefinition in a loop and then write the result with QgsVectorFileWriter :

layer1 = QgsVectorLayer(path1, 'layer 1', 'ogr')
layer2 = QgsVectorLayer(path2, 'layer 2', 'ogr')

idx = layer2.dataProvider().fieldNameIndex('colL2')
col2values = layer2.uniqueValues(idx)

fields = layer1.dataProvider().fields()
writer = QgsVectorFileWriter(pathOutput, "CP1250", fields, QgsWkbTypes.LineString, layer1.crs(), "ESRI shapefile")

for value in col2values:
    # Layer 1 and Layer 2 have the same column name 'colL2' but Layer 1 values are longer
    selectLayer2 = QgsExpression(" \"colL2\" LIKE '" + str(value) + "' ")
    selectLayer1 = QgsExpression(" \"colL2\" LIKE '" + str(value) + "%' ")
    selectedLayer2 = layer2.getFeatures(QgsFeatureRequest(selectLayer2))
    idsLayer2 = [f.id() for f in selectedLayer2]
    layer2.selectByIds(idsLayer2)

    selectedLayer1 = layer1.getFeatures(QgsFeatureRequest(selectLayer1))
    idsLayer1 = [f.id() for f in selectedLayer1]
    layer1.selectByIds(idsLayer1)

    parameters = {'INPUT': QgsProcessingFeatureSourceDefinition(layer1.id(), True),
                  'OUTPUT': 'memory:'}
    layer1_line = processing.run('qgis:polygonstolines', parameters)

    parameters = {'INPUT': QgsProcessingFeatureSourceDefinition(layer2.id(), True),
                  'OUTPUT': 'memory:'}
    layer2_line = processing.run('qgis:polygonstolines', parameters)

    parameters = {'INPUT': layer1_line['OUTPUT'],
                  'OVERLAY': layer2_line['OUTPUT'],
                  'OUTPUT': 'memory:'}
    differenced = processing.run('qgis:difference', parameters)

    writer.addFeatures(differenced['OUTPUT'].getFeatures())

    layer1.removeSelection()
    layer2.removeSelection()

del writer

But I get following error :
raise QgsProcessingException(msg)
_core.QgsProcessingException: Unable to execute algorithm
Could not load source layer for INPUT: layer_1_XXX not found


texts: Cutting a polygon using the shape of another polygon layer?

texts:
I have a layer with polygons I have to cut smaller. I also have another layer with the shapes of the polygons I have to use for cutting out some pieces of the polygons on the first layer.

I tried to read the help on "cutting holes in existing polygons" and for some reason I can't do it.

If you're already answering, maybe you could tell me also how can I just cut a hole in existing polygon by freehand?


texts: How to Group Parameters in PyQGIS Processing Plugin

texts: The PyQGIS Processing Plugin I'm writing has a lot of parameters, most of which are advanced and won't be necessary for most users. I'm concerned that a user will be confused if they see all these parameters that they don't know what to do with. Is there a way I can group parameters in a collapseable "advanced parameters" box in the plugin's GUI similarly to how the 'Round raster' built-in processing program does?


texts: Cluster points into evenly numbered groups

texts: I need to take x number of points and equally/spatially divide them into groups. Meaning that if I had 40 points like in the image below and wanted to have 4 groups of points. Each group should have 10 points and all the points should be spatially near each other. I have tried using normal K-means cluster but that only spatially correlates data while not taking into account the even number of points in each group. I have tried writing this in JavaScript as well as Python. Does anyone have any thoughts on what would be the best way to go about evenly and spatially splitting a set of points?


texts: Using Space Syntax Analysis plugin for QGIS?

texts:
I'm trying to run the Space Syntax Analysis plugin for QGIS.

I've added the metadata file that was missing, but it does not appear in the Plugin toolbar when I click on it.

Can someone help me to figure out what is going on?


texts: Cutting holes in custom polygon made with Leaflet.Draw

texts:
I would like to be able to create holes in polygons made with Leaflet.Draw plugin. For example, I currently have code like this: https://jsfiddle.net/4gkrfu5b/

If you draw a polygon on top of another polygon, there will be two polygons overlapping each other. I would like to get an effect, similar to this example:
https://jsfiddle.net/bnue5tjg/5/

The second example nicely implements the cutout possibility with polygon drawn on top of another polygon. The only downside it has, it is made with one existing polygon pre-rendered on the map. I prefer to start with a clean workspace and then be able to create polygons and holes in them. I have tried to implement the second example into my code, but for now I have been unsuccessful with this. Any help with cutouts implementation will be welcome. Hope that the stackexchange community can help me with implementation of cutouts / holes on polygons we make on the map.


texts: Script for selective dissolving of multiple features

texts:
I am kinda new to Python and arcpy and have done some basic online courses on this. So I do have an idea how these scripts are structured/written and how they work. Now I have to perform a seemingly complicated analysis for my master thesis which goes something like this:
I have an excel file in matrix form of around 400 polygon features (municipalities). Each municipality has a relationship with others around it (no. of people commuting daily between the municipalities). I have to somehow make a script/model which dissolves each municipality with only that municipality (out of 400) which has the highest no. of commuters travelling towards it (not from it). This might make one municipality dissolve with multiple municipalities also which is fine.
Anything like this has been done before? I'm thinking about the script and I can't think how I make this selective dissolve?

This is how the table looks like its not really a matrix but I guess it can be converted into one like in mentioned by Richard. The first two columns are origin and destination IDs and the last column is the no. of commuters.


texts:
I have a shapefile with 2 points.

I am using the XY to Line tool in the ShapeTools plugin to create a Geodesic, Great Circle, and Simple Line. However, in each case, the tool is running successfully, but producing empty outputs.

The following screenshots shows the parameters:

The straight line connecting the 2 points is getting successfully created if I use the Points to Path tool.

The aim of this exercise is to be able to see the difference in geodetic distance and Euclidean distance. I also want to show how the projection of the geodetic line connecting the two points look like on a flat map.

If I do not check the optional check boxes, Use the point geometry for the line starting point and Use the point geometry for the line ending point, the tool show the error Please select valid starting and ending point columns.

I am using QGIS 3.4.1


texts: Python code error when using QGIS plugin Openlayers

texts:
I am getting this error message whenever I try to use the Openlayers plugin. I can not read or write Python particularly well at the moment, though I am trying to learn, so I am struggling to identify the problem. I have tried un-installing and reinstalling openlayers plus using it on previous versions of QGIS with no success. Basic attempts but really haven't got the time to delve deeper into more complex procedures with out really knowing what I'm doing.

"An error has occurred while executing Python code:

Traceback (most recent call last):
  File "C:/Users/*****/.qgis2/python/plugins\openlayers_plugin\weblayers\weblayer.py", line 97, in addLayer
    self._addLayerCallback(self)
  File "C:/Users/*****/.qgis2/python/plugins\openlayers_plugin\openlayers_plugin.py", line 130, in addLayer
    layer = OpenlayersLayer(self.iface, self._olLayerTypeRegistry)
  File "C:/Users/*****/.qgis2/python/plugins\openlayers_plugin\openlayers_layer.py", line 330, in __init__
    self.olWebPage = OLWebPage(self)
  File "C:/Users/*****/.qgis2/python/plugins\openlayers_plugin\openlayers_layer.py", line 48, in __init__
    proxy = getProxy()
  File "C:/Users/*****/.qgis2/python/plugins\openlayers_plugin\tools_network.py", line 42, in getProxy
    proxy.setPort(settings.value("/proxyPort", 10, type=int))
TypeError: unable to convert a QVariant of type 10 to a QMetaType of type 2

Python version:
2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]
QGIS version:
2.10.1-Pisa Pisa, d20c5b7
Python path: ['C:/PROGRA~1/QGISPI~1/apps/qgis/./python/plugins\\processing', 'C:/PROGRA~1/QGISPI~1/apps/qgis/./python', u'C:/Users/*****/.qgis2/python', u'C:/Users/*****/.qgis2/python/plugins', 'C:/PROGRA~1/QGISPI~1/apps/qgis/./python/plugins', 'C:\\PROGRA~1\\QGISPI~1\\apps\\Python27\\lib\\site-packages\\matplotlib-1.3.1-py2.7-win-amd64.egg', 'C:\\PROGRA~1\\QGISPI~1\\apps\\Python27\\lib\\site-packages\\nose-1.3.3-py2.7.egg', 'C:\\PROGRA~1\\QGISPI~1\\apps\\Python27\\lib\\site-packages\\tornado-4.0.1-py2.7-win-amd64.egg', 'C:\\PROGRA~1\\QGISPI~1\\apps\\Python27\\lib\\site-packages\\backports.ssl_match_hostname-3.4.0.2-py2.7.egg', 'C:\\PROGRA~1\\QGISPI~1\\apps\\Python27\\lib\\site-packages\\certifi-14.05.14-py2.7.egg', 'C:\\PROGRA~1\\QGISPI~1\\bin\\python27.zip', 'C:\\PROGRA~1\\QGISPI~1\\apps\\Python27\\DLLs', 'C:\\PROGRA~1\\QGISPI~1\\apps\\Python27\\lib', 'C:\\PROGRA~1\\QGISPI~1\\apps\\Python27\\lib\\plat-win', 'C:\\PROGRA~1\\QGISPI~1\\apps\\Python27\\lib\\lib-tk', 'C:\\PROGRA~1\\QGISPI~1\\bin', 'C:\\PROGRA~1\\QGISPI~1\\apps\\Python27', 'C:\\PROGRA~1\\QGISPI~1\\apps\\Python27\\lib\\site-packages', 'C:\\PROGRA~1\\QGISPI~1\\apps\\Python27\\lib\\site-packages\\PIL', 'C:\\PROGRA~1\\QGISPI~1\\apps\\Python27\\lib\\site-packages\\jinja2-2.7.2-py2.7.egg', 'C:\\PROGRA~1\\QGISPI~1\\apps\\Python27\\lib\\site-packages\\markupsafe-0.23-py2.7-win-amd64.egg', 'C:\\PROGRA~1\\QGISPI~1\\apps\\Python27\\lib\\site-packages\\pytz-2012j-py2.7.egg', 'C:\\PROGRA~1\\QGISPI~1\\apps\\Python27\\lib\\site-packages\\win32', 'C:\\PROGRA~1\\QGISPI~1\\apps\\Python27\\lib\\site-packages\\win32\\lib', 'C:\\PROGRA~1\\QGISPI~1\\apps\\Python27\\lib\\site-packages\\Pythonwin', 'C:\\PROGRA~1\\QGISPI~1\\apps\\Python27\\lib\\site-packages\\Shapely-1.2.18-py2.7-win-amd64.egg', 'C:\\PROGRA~1\\QGISPI~1\\apps\\Python27\\lib\\site-packages\\six-1.3.0-py2.7.egg', 'C:\\PROGRA~1\\QGISPI~1\\apps\\Python27\\lib\\site-packages\\wx-2.8-msw-unicode', 'C:\\PROGRA~1\\QGISPI~1\\apps\\Python27\\lib\\site-packages\\xlrd-0.9.2-py2.7.egg', 'C:\\PROGRA~1\\QGISPI~1\\apps\\Python27\\lib\\site-packages\\xlwt-0.7.5-py2.7.egg', u'C:/Users/*****/.qgis2//python', 'C:\\PROGRA~1\\QGISPI~1\\apps\\qgis\\python\\plugins\\fTools\\tools', 'C:/Users/*****/Documents/*****/*****']


texts: Difference between CRS (Coordinate Reference System) and SRS (Spatial Reference System)

texts:
I could see that CRS and SRS are referred as different components in this question
Is WGS84 itself a Coordinate Reference System?

Could anyone please explain whether CRS and SRS are same or different?


texts: Building attribute table for raster building footprints?

texts:
I am doing a project at work looking at rooftop solar resources. I have a set of building footprints for a few small cities, which "become" a set of raster building footprints when I clip our solar isolation raster to the building footprints.

I have completed the whole process, but now my supervisor would like me to go back and find the top 10 buildings (rooftops) in terms of their potential solar resources for each city. The problem is that the raster footprints, obviously, do not have attribute tables and (despite being clipped to the shape of the footprints) are just a collection of individual data points.

I need to somehow aggregate/match the 350,000 unique raster data points (yes, I converted the raster to integer and created an attribute table) back to their respective building footprints, so I can then find the top 10 buildings in terms of their solar resources.

Anyone have an idea how to go about this?


texts: Select n cells from a raster maximizing their variance

texts:
I work in R, especially using 'raster' package. It is a computational problem too hard for my syntax capacities.

STEP 1.
I have a one band raster, and what I want to select is a set of n cells (in my case, n = 15) as different as possible different regarding their values. In mathematical terms, in my selection I want to maximize the variance of values. Obviously, there are a huge amounts of combinations, but I would select the best one solution.

STEP 2.
One further methodological complication is that at the same time I would like to maximize the distance between cells. I mean, once reached a set of possibilities from step1 (variance optimization), I would like to select the best solution in terms distance between cells (for example, the sum of distance between cells).

Is it possible? (at least the step1, but preferentially also step2)


texts: How to convert QGIS Field Calculation to ArcGIS?

texts:
I have a field calculator expression I use in QGIS to generate unique IDs that are sequential, but start with a text string.  I now need to perform the same thing in ArcGIS, but I am unsure of how to start.  Here is my QGIS expression:

 'cos' || ( tostring(  $rownum  + 1000000))

I was thinking that in ArcGIS it would be something like this:

 "cos"+InStr( [FID] + 1000000 )

But it returned an error.

Any help would be appreciated.
Thanks...


texts: Dynamic SLD styling based on min / max value in a PostGIS column

texts:
I have a PostGIS table that contains the GDP of countries together with their polygon boundaries. I am using GeoServer to display this data via OpenLayers.

When styling the data, I want to avoid having to hardcode numeric range / colour pairs in my SLD file.

Is it possible to create an SLD style file that dynamically determines the minimum and maximum value of GDP that are present in a PostGIS column, for example through CQL, and then uses that information to colour the map?


texts: Reading backscatter value of each pixel of Sentinel-1A images in Java?

texts: I am new to working with Sentinel-1A images but have better knowledge of Java. I want to get backscatter value of each pixel of Sentinel-1A image, store it and then get the frequency of each value. Can you please tell how to get it and which data structure or any database should I use to store backscatter values? Please tell the whole process. Please also tell which software I have to use for it?


texts: Determing, whether WMS service supports CQL filtering

texts: Is it possible to query, whether WMS service supports CQL filtering (like Geoserver does) or not? I didn't find any hints from GetCapabilities query result.


texts: how to create graphs while running python with ArcGIS 10

texts:
As Matplotlib seem not to work (what a drag it's such a nice tool) I'am looking for alternatives to create graphs while running python with ArcGIS 10. Besides regular pie and column charts I need in particular the possibility to create star plots/ radar charts / amoeba diagram - most likely something like:

http://matplotlib.sourceforge.net/examples/pylab_examples/polar_bar.html

Most likely with a tool with as less as possible dependencies.
Any recommendations are welcome

Thanks!


texts: Updating point geometry based on geometry of another feature class using arcgisscripting?

texts:
I'm looking to update the point geometry of lets say "oldFC" by obtaining the shape geometry of "newFC" and applying that to the "oldFC".  I also want it to grab the value in the ID field in "newFC"and use that to search for the same ID in "oldFC" and then set the shape value.

I believe I have the loop set up and it appears to be grabbing the ID just fine but it doesn't update the shape geometry.

I'm very much a beginner with using Python so my code is rough.

Here's the code snippet...  running 9.3 (no access to 10):

rows = gp.SearchCursor(newFC):

row = rows.next()

gp.MakeFeatureLayer_management(oldFC, oldLayer, "")
gp.MakeFeatureLayer_management(newFC, newLayer, "")

while row:

newgeo = row.getValue("Shape")
newID = row.getValue("ID")

print newgeo
print newID

where =  "ID" + '= '   + str(newID) #my sql expression

print where

new_rows = gp.UpdateCursor(oldFC, where)
new_row = rows.next()
while new_row:
    new_row.setValue("Shape", newgeo)
    new_rows.UpdateRow(new_row)
    new_row = new_rows.next()
row = rows.next()


texts: Use Display Expression for KML Popup, SHP to KML ArcGIS

texts:
In Arcmap I've created a complex label expression in the Display expression for my layer.  I'd like to have that expression displayed when I convert the layer to .kmz file for display in Google Earth.  I'm using the ESRI Layer TO KML conversion tool.

What happens in the popup right now is the display expression jammed, without the correct formatting, into the feature name box and also the feature description box, followed by the complete list of fields.

I've been looking at the HTML popup tab, but don't see any way to use the Display Expression I already have in that tab.  Is it possible for this to happen?  It looks like otherwise I'll have to write a XSL template.

EDIT:  Thanks to @KHibma for the comment.  What I've ended up doing is

removing all of the fields from display in the shapefile properties in ArcMap before the export so they are not included in the kml popup
Labeling the features in arcmap with one field so that field value is the "title" of the feature in the popup.
adding HTML formatting tags to the Display Expression (which I was writing in Python).  this expression appeared, formatted, in the blue box below the feature title.

I can't add a screenshot of the result right now, but will be able add it later.

This worked fine on a small scale, but when trying to export an entire parcel fabric for a county I would get a C++ runtime error.

One interesting thing I found was that, even before the export, if there was an error in the FindLabel function and it was not able to create a label, it would come back as empty.  So I added a try/except statement within the function like:

def FindLabel([field1],[field2]):
    try:
        #code for function
    except:
        return "there was a problem creating this display"

which worked great.  However, when I tried to import sys or traceback in the function and set the error handling message as the return in the except statement, it still returned blank (i.e. didn't return the text of the error handling message).  In the end, I fixed all the problems with the expression, but still got the runtime error when making the big export.


texts: Leaflet->how to draw a polygon with reference to another geographical position

texts:
I have a basic leaflet map in which I draw a circle representing a city.
Next, I would like to locate a rectangle with reference to that circle. (more specifically drawing a rectangle not with reference to its centroid or any corner x,y; but with reference to the circles x,y. So, I need to define the distance with these two objects, and dimensions of rectangle.)

Please see the image below as a graphical expression of this issue:


texts: Creating line-strings from points in MySQL

texts: I am extracting the lat/long information from Geo-tagged images and storing them in MySQL. I want to get a line-string from the points, I looked into the GROUP_CONCAT() function but I am not able to get it to work with the linestring. Is there a way to store the concatenated geometries of the individual points as a linestring?


texts: Joined table fields don't work in print composer?

texts:
I am trying to use the attributes of the current atlas feature, which is from a joined table. In the expression builder it shows the correct value, but the result is blank when I actually try to insert it. Even using the attribute function, it doesn't let me pull back the data of those joined fields.

Am I doing something wrong? Do I need to flatten the data to the shapefile? I am using version 2.14


texts: Is it possible to have Feature Dataset in another Feature Dataset?

texts:
Understand from ArcGIS website that we can use feature dataset to group together related feature classes.

Is it possible to have another feature dataset in another feature dataset? Understand that there is a performance issue if we have too many things under one feature dataset, but at this point, it is very unlikely to have a lot of accessing the dataset once it is added.

The scenario that I tried to achieve:
1. Have a Report and archive feature dataset to group reports and archives.

2. In each of this feature dataset, it will be further grouped by different category.

For example:

Report Dataset
  Report Type A Dataset
    Report October 2014 Feature Class
    Report November 2014 Feature Class

Archive Dataset
  Report Type A Dataset
    Report October 2013 Feature Class
    Report November 2013 Feature Class

If it is not possible, any other suggestion to achieve this? For example: Create a folder instead of Feature dataset. (Heard that there is a way to do this using ArcGIS command line tools?)


texts: Use Values from Raster Statistics (e.g. sum of all cells) in Raster Calculator

texts:
I need to do calculations on raster layers which require the sum of all cells of the raster. I've used the zonal statistics tool for that but it gets a little annoying to do those extra steps (run the plugin, open the Attribute table, copy the sum...).

Is there a way I can calculate and use the values of raster statistics directly inside a raster calculator?

Even better would be to calculate those values from a raster created in the same step:

  A*B*C / Sum (A*B*C)


texts: Check whether table has overlapping polygons, in PostGIS?

texts:
Using PostGIS, what is an efficient way to check whether a table with a geometry column has geometries that overlap each other?

This is my table:

  Column  |            Type             | Modifiers
---------+-----------------------------+-----------
 sts     | character(1)                |
 geom    | geometry(MultiPolygon,5070) |

I guess I could see whether:

SUM(ST_Area(geom))

is the same as:

SUM(ST_Area(ST_Union(geom)))

But that doesn't feel like it's likely to be the most efficient way to do it, especially since my table has about 40k features.

Any ideas?

One comment on a related question suggests one could "search for overlaps by their bounding boxes" - is that a possibility?


texts: Unable to view local geoserver in that dropdown menu of Open Geosuite?

texts:
I'm using GeoExplorer(Web Map Client) in OpenGeo Suite 4.5 version. I use to upload the layers from my LocalGeoserver by clicking on Add Layers Button. But now i'm unable to view my local geoserver in that dropdown menu.

How can I troubleshoot this problem?


texts: How to handle plugins of two different QGIS installations on the same computer

texts:
I have a computer with 2 versions of QGIS. Both versions are running on Win 7 64 bits with 16 Mb of ram. One is Dufour 32 bits and the other one is version 2.1 OSgeo4W 64 bits.

I am experimenting using the processing menu on both versions.  Nevertheless I noticed that the paths for Grass in my 2.01 32 bits version were the ones of the OSgeo4W install (Processing/Options and configuration/Providers). I changed to the Dufour ones.

Then I tried to replicate a procedure in my 64 bit install but the paths were the Dufour ones!

My question is, then, how I can handle the settings for each install in order that they are independent from each other and I do not have to change my provider paths each time that I want to use the processing menu.


texts: Using ModelBuilder to Iterate Row Selection with ArcMap?

texts:
I have a model that does a lot of linear calculations, then eventually runs into an ITERATE ROW SELECTION tool.  I only want the iterator to cycle through the tools downstream of itself, but instead it runs the entire model for every iteration.

Since the other calculations take several hours to perform and will not change between iterations, is there any way to have the model only loop through the processes downstream of the iterate tool?

Basically I want to initialize a for loop, run only part of the model in the loop, then close the loop, but I can't figure out how to do that.


texts:
I can not find the CARTO dataset to connect the arcgis service:
http://ags.cuzk.cz/arcgis/rest/services/ZABAGED_TOPO/MapServer/

Can you help me?


texts:
Let's say I have two point layers that contain locations of wells on a property.

layer A has 50 points.
layer B has 35 points.

Some of layer B's points already exist on layer A, and some of them do not.

What do you guys think would be the best way to go about removing any redundancy and creating a newly compiled layer?

Is there a pre-built tool for this I am unaware of?


texts: Tool for linear interpolation with a mask?

texts:
I'm asking if anyone knows a tool to do a linear interpolation that supports masks / barriers.
I need to do interpolation on a archipelago and no points on opposite shore should affect on a other shore. So the water area should behave as a mask polygon.

The known points are access times to the location and I want values between points to be linear. Are there other linear interpolation methods than TIN interpolation that use barycentric interpolation?

Here is an example image of my situation

gdal_grid handles the interpolation but is missing the mask option. Are there some other tools?

I would be glad for any help and advise.


texts:
I wanted to know what drawing methodology is used by qgis renderer like gdi,gdi++ or opengl. I have searched its documentation http://qgis.org/api/classQgsMapRendererJob.html
http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/composer.html
but didnt found any methodology used by its renderer for drawing.


texts: How to add layers to MXD programmatically on linux or windows without ArcGIS for Desktop?

texts: Is it possible for me to programmatically create in linux centos 6.4 or a windows with only arcGIS Engine runtime (single use license without sdk)?


texts: Pasting Point goemetry with Z and M values to Simple Point Geometry

texts:
This should be kind of easy but I'm looking for a best practice scenario. I just tried adding 14 points from a point shapefile into another layer in our ArcSDE database. I assumed the layer in SDE was also a point file since it displays as such. I got an error saying there were no files with the appropriate geometry to paste into. I looked and the geometry type of the SDE layer is listed as "Point" but the feature type is "Simple". I inherited this database so I don't know how the layer was created. Is this "simple" geometry something like an event theme (e.g. was it likely just created from a spreadsheet?)?

More importantly, how should I proceed to get the new points from my shapefile into this "Simple geometry" layer?

I also learned just now that the input shapefile has z and m values so I'm exploring stripping those off using "Environment Settings" and the "Copy Features" too.


texts: How to convert a complex project of shapefiles to a SpatiaLite database?

texts:
We have a QGIS project made of ~10 shapefile layers. The layers are saved on a shared server, and the project is used on different computers.

Recently polygons and polylines have been been disappearing when the layer is saved (still present in the attribute table but not in the mapping window). We have been advised that we would be better to move to a SpatiaLite database.

What is the best and quickest way to convert this number of shapefiles (with complex styling) to a SpatiaLite database?

Is it possible to put all of these shapefiles onto one SpatiaLite database, or will they have to be different databases?


texts:
Im using QGIS 2.18.2 and have tried to export layers as dxf using the project tab and by clicking the layer and using save as. Both have given me error message:

Export to vector file failed.
Error: creation of data source failed
OGR error:Failed to open 'dxf-parcels-export.dxf' for writing