texts: Python script crashes at regular intervals when summing time series raster data

texts:
arcpy/python beginner here.

I am trying to write a script that uses the Cell Statistics function to derive annual rainfall data for an entire continent  from 1983 to 2016. The original data is made up of slightly over 12,000 individual tiff files (dimensions 751 x 801 pixels, 1 degree/11km  spatial res.) measuring rainfall for every single day for the entire period. The filenames include a consistent time stamp in the form of *name.yyyymmdd. The annual aggregates I seek are the rainfall totals from June to October (5 months) of each year. I set up my script with the monthly loop nested within the annual loop as follows:

    import arcpy
    arcpy.CheckOutExtension("Spatial")
    arcpy.env.workspace = r'C:/data/rainfall'
    arcpy.env.overwriteOutput = True
    years = range(1983,2017)
    months = ["06","07","08","09","10"]
    yearlist = []
    for year in years:
        for month in months:
            monthlist = arcpy.ListRasters("*"+str(year)+month+"*","TIF")
            yearlist.append(monthlist)
        #convert list of monthly lists to simple yearly lists and sum
        ylist = [j for i in yearlist for j in i]
        annualsum = arcpy.sa.CellStatistics(ylist, statistics_type = "SUM")
        annualsum.save(r'C:/data/rainfall/aggr/'+str(year)+'.tif')
        print str(y)+" annual rainfall total done!"
        del annualsum
     print "ALL DONE!"

The code runs fine initially but crashes while saving the output raster ("annualsum") at the '1988' loop point. No specific error message,  just a pop up that says "pythonw.exe has stopped working." When I alter the begin date to restart after '1988' it would crash again at the '1994', rinse repeat with 2000, 2006 and 2012 points - exactly 6 year intervals. I was convinced it was a memory issue since I am using 32 bit python/ArcGIS but I thought deleting the raster object after saving would solve the problem. While I finally got it all done by restarting after each crash, and doing the 'crash' years separately.

I use 32 bit because I only have the 1 year student/evaluation license for ArcGIS. My pc is however an i5 with 16GB of RAM.


texts: bringToFront(), bringToBack() on L.geojson problem

texts:
I'm trying to make a simple interactive map of Germany. I'm using mapbox.js (2.4.0) and have it set up like this:

var map;
var currLayer, citiesLayer, riversLayer;

$(document).on("ready", function() {
L.mapbox.accessToken = '';
map = L.mapbox.map('map', 'vantuch.d55cb2a1', {zoomControl: true, scrollWheelZoom: true}).
setView([51.0621,10.9863], 6);

    // CALL AJAX REQUEST FOR TESTING
    findCities();
    findRivers();
    findCityPlaces();
});

I am calling AJAX functions to set up a L.geojson layer after success like this:

/* functions for setting up MARKERS (find_river_places) */
function showMarkers(map, geojson, selectType) {

    geojsonLayer = L.geoJson(geojson, {
        // setting up markers
    });

currLayer = geojsonLayer;
map.addLayer(geojsonLayer);
}

function showCitiesMarkers(map, geojson) {
    geojsonLayer = L.geoJson(geojson, {
    pointToLayer: function(feature, latlng) {
    var marker;

        //return new L.CircleMarker(latlng, {radius: 5, fillOpacity: 1.0}).bindPopup(popup);
        marker = new L.Marker(latlng, {
        icon: new L.DivIcon({
            className: 'circle-border-city',
            html: '
' +
          '' + feature.properties.name  +
          '' + '
'
        })
    });
    return marker;
    }
});

citiesLayer = geojsonLayer;
map.addLayer(geojsonLayer);

}

What I want to do is to have layers for rivers and cities always on the top of other layers (user should be able to turn them on or off), no matter which function is called the first. But they're always rendered at the bottom.

I thought that citiesLayer.bringToFront() after adding other layers should work. But after testing it in the console it does nothing. So what am I doing wrong, or is there a better way to do this ?


texts:
Within arcpad (edit form) I want to get the current location of the vbs script running.  This simple vbs script outside of arcpad works fine.

wscript.echo " " + Wscript.ScriptFullName

However, when I call this statement within a button click event in arcpad is gives and error and says object wscript required.

I also tried this:

Dim sCurPath
sCurPath = CreateObject("Scripting.FileSystemObject").GetAbsolutePathName(".")
Application.MessageBox(sCurPath)

However, I'm getting

  C:\Users\username\Desktop

directory (not location of vbs script).

Any ideas of how to find current directory of vbs script running in arcpad edit form?


texts: Is Django's LayerMapping capable of putting a static value into all features?

texts:
I'm using GeoDjango's LayerMapping utility to try and import a set of features. The catch is that aside from the attributes in the source layer, I want all the destination features to have the same value of another attribute. Specifically:

> sourceLayer = theDataSource[5]
> sourceLayer.fields
['DataSources_ID', 'Source', 'Notes', 'Shape']
> from django.contrib.gis.utils import mapping
> fieldMapping = mapping(theDataSource, layer_key=5)
{'datasources_id': 'DataSources_ID', 'source': 'Source', 'notes': 'Notes', 'shape': 'MULTIPOLYGON' }

The destination model has another field, 'owningmap', which is an integer field. I would like to load data from the OGR source via...

> lm = LayerMapping(DestinationModel, theDataSource, fieldMapping, layer=5)
> lm.save()

... but I would like every feature to have the same integer value put in to the 'owningmap' field. As though I could somehow add a "static" mapping to fieldMapping...

> fieldMapping['owningmap'] = 2

Is this possible? Otherwise my last resort is to load via LayerMapping, then calculate the value for all the just-loaded features in a separate database transaction...

Thanks!


texts: How does Personal Geodatabase work from Microsoft Access perspective?

texts:
I am trying to get a better general understanding of how ArcGIS works with MS Access when using a personal geodatabase. If I create a new pgdb and then import a very simple shapefile and finally look at the result through Windows explorer (rather than ArcCat), the only visible item from the pgdb creation process is one mdb file. Opening the mdb file in Access shows there are 11 tables (in the case of a shp file that had one point).

Generally, I can find very little documentation aboutlooking at a pgdb through Access, and if there is any sense in doing so (even if it is just for the benefit of trying to find out how it all works). Most resources about pgdb design don't go say too much compared the the other two recommended gdb alternatives. However, Access is still extremely popular, so I think there is some benefit to looking at the topic in more detail.

So is it possible to know how the mdb tables that ArcGIS creates are related, or are the relations defined by files that neither Access nor Windows can see?

Is it possible to know how ArcGIS uses these tables to provide the "spatial" component to Access?

Is it ever beneficial to use MS Access to work with a shp file stored as a pgdb or is this all futile because pgdb technology is proprietary and no detailed publicly available knowledge exists?


texts: Create QML Widget with a picture in Attributes Form

texts:
I would like a Qmlwidget in the Attributes Form to show a picture.
So far I have created the basic qml for a blue rectangle:

import QtQuick 2.0

Rectangle {
    width: 100
    height: 100
    color: "steelblue"
    Text{ text: "A rectangle" }

}

I'd like to replace the Rectangle dictionnary by a Picture dictionnary with the picture path generated from project folder path and an attribute.

maybe something looking like that?:

Image {
    width: 100
    height: 100
    source: "expression.evaluate("concat(@project_folder,
                                      '/Photos/',
                                      attribute( $currentfeature, 'name'),
                                      '_1.jpg')"
}


texts: How to estimate distance from points to line feature with v.distance?

texts:
I'm trying to use the v.distance module in GRASS 6.4.3 to estimate the distance from a series of point features (animal tracking locations) to the nearest point along a line feature (a coastline). I load in the two shapefiles (both UTM/WGS84),create a new column in the attribute table of the point features to take the estimated distances and start v.distance. I set the 'from' vector map to be the point layer and the 'to' vector map to be the coastline, tick "distance to the nearest point on the 'to' feature along that linear feature", and set the column to upload to be my new distance column.

The script runs ok but the estimated values in the attribute table are complete nonsense. At least half the values are identical, even though the animal was progressively moving away from the coastline and therefore distance should have been increasing. When the animal was 2000km from the coast the distance estimate is 166m, but right nect to the coast it is often >3000m! Funnily though, if I select the option to create an output vector map containing lines connecting the nearest elements, the lines perfectly connect each point to the nearest part of the coastline and are all different lengths.

What is going wrong here?! The frustrating thing is that I ran it once when testing and it worked perfectly, but all subsequent attempts have resulted in the same problems. I can't remember if I used any different settings on the first go. I'm sure it's something simple but it is driving me nuts trying to work out what it is! I know it may be hard to comment without seeing the data, but i'm really hoping someone has run into this before or has some idea where my mistake might be.

All help greatly appreciated! Many thanks,


texts: Script tool: Batch downgrade MXD version

texts:
I am trying to take the guts of this code (Batch convert ArcGIS 9x files to 10) and build it into a tool that we can more flexibly use and share around our organization.  I think I'm almost there but when I run the tool as a script in Arc, nothing happens.  Thoughts?  Thanks for your help.  I'm fairly new to Python but am pretty good on where the guts of it all is.

mxdinloc is set to a folder

mxdoutvers is set to a string - value list where the list choices are (10.1, 10.0, 9.3, 9.2) - no quotes

EDIT: I ran through, got it working and wanted to share the final bit that I brought into my Script Toolbox:

#import modules
import arcpy
import os
import glob

#set parameters for inputs
mxdinloc = arcpy.GetParameterAsText(0)
mxdoutvers = arcpy.GetParameterAsText(1)

#set the values for input versions of ArcMap
if mxdoutvers == "10.1":
    version = "10.1"
    suffix = "_(10.1)"
elif mxdoutvers == "10.0":
    version = "10.0"
    suffix = "_(10.0)"
elif mxdoutvers == "9.3":
    version = "9.3"
    suffix = "_(9.3)"
elif mxdoutvers == "9.2":
    version = "9.2"
    suffix = "_(9.2)"

for pathname, directories, filenames in os.walk(mxdinloc):
    for filename in filenames:
        if filename.lower().endswith(".mxd"):
            print "\nUpdating: " + filename
            mxd = arcpy.mapping.MapDocument(mxdinloc + "\\" + filename)
            filename, ext = os.path.splitext(mxdinloc + "\\" + filename)
            print filename
            if mxd.dateSaved:
                print mxd.dateSaved
                mxd.saveACopy(filename + suffix + ext, mxdoutvers)
                currentMxd = filename + suffix + ext
                mapDoc = arcpy.mapping.MapDocument(currentMxd)
                print "Updated: " + filename + "\n"
        else:
            print "\n" + filename + " is not a MXD and was not updated. \n"


texts:
What defines a Class in the Heat map plugin? As in:
"Black holes like that are usually caused by raster values exceeding the range of the color ramp. Try adding another class with a higher associated value. That should fix it."
This quote comes from an  earlier anwered question:
Why is there a black hole in the Heatmap?

I'm getting all sorts of clunky results and I'm not sure whether I need to find the right settings, alter my data, or both.  I'm working on identifying all important parameters: obviously Radius, and cell size, their relation, extents (I'm foggy on), projection (I'm still not positive how it affects the heat map as long as it sit on top of the data.  And it does.), projection units,...

I'd like to know better specifically what affects the curve of the histogram. Also is the color ramp synonymous with the histogram curve?

In addition to lots of other issues, I too have black holes.  But I have to say it is awesome fun ferreting all this out.  And Thank you, Underdark for pointing the way.  I'm starting to find some light with the raster histogram.


texts: Which hardware for running FME 2015 on linux

texts:
I am running FME on Xubuntu. It works very fine save for some minor glitches with the inspector. I am planning on setting up a couple of rather heavy vector calculations (areaonarea,buffering and dissolves) as well as raster calculations. For the latter I will probably be using Orfeo Tool Box from FME.

I asked IT for new hardware (I am running some pretty ancient hardware on my linux desktop) and they told me to specify a new computer.

I am therefore wondering:

Can FME 2015 on Linux can make use of more than 16 Gb of memory?
Safe software suggests Intel core i7. Which processor speed would be
a good option?
If SSD disks are necessary?
If the software can also make use of graphics cards?


texts: Geoserver strange inverse polygon rendering

texts:
I found strange behaviour of rendering one particular polygon on Geoserver. Red polygon in the link. If you compare SHP data from ZIP file in some desktop software like QGIS and rendering from Geoserver, you can see that Geoserver renders that polygon differently.

I am trying to find a reason why this one polygon is wrongly rendered on Geoserver. Does anybody know why this is happening?

I have tested it with Geoserver 2.4 and Geoserver 2.5.1 and data is stored in PostGIS not in SHP.


texts: How to add ECW support to GDAL bundled with PostGIS Windows

texts:
I'm currently trying the new raster functionalities of PostGIS 2.0 and I successfully imported some GeoTIFF file into PostGIS (although it only find one band so it's all black and white for some reason). I used raster2pgsql.exe bundled with PostGIS 2 to store it.

Most of the raster data i have is in ECW format and I would like to import them directly with raster2pgsql. I know i can convert my ECW to GeoTIFF but i would like to do it without conversion.

When i list the supported input format of raster2pgsql (raster2pgsql -G), it does not list ECW as it can not be distributed with the install. So i downloaded the ERDAS SDK following these instructions:
http://trac.osgeo.org/osgeo4w/wiki/pkg-gdal-ecw

And copied the DLL in my SysWOW64 system folder (i'm using Windows 7 64-bits). But it still won't list ECW as an input format. When looking up online, i found that PostGIS uses a bundled GDAL found in

C:\Program Files (x86)\PostgreSQL\8.4\bin\libgdal-1.dll

I don't know where to copy the DLLs so that the GDAL bundled with PostGIS will grab them and accept ECWs rasters. Or if anyone know how to tell PostGIS to use the GDAL install that comes with FWTools (C:\Program Files (x86)\FWTools2.4.7\bin) which is directly compatible with ECWs (path variables?)

Thanks a lot for your help!

Fabien


texts: Creating new empty memory layer with fields scheme from other layer in QGIS?

texts:
Being new with python so far, I got inspired by this answer: https://gis.stackexchange.com/a/155988/8202
to dive into it a bit.

When creating features in a new target layer from an existing source layer, all "solutions" I found so far to create that new layer, were to either

prepare the new target layer manually in QGIS beforehand, adding all fields (with field types, lengths etc).
Hard-code the field settings into the python script when creating that new target layer.
How to create a new empty vector layer programmatically?

It seems a bit odd to me that there should be no programmatic way to clone the attribute scheme of a layer, and create a new (memory) layer with this scheme, so you wont have to change your script for each new layer, or fumble things together manually.

I'm focused on the creation of memory-layers, cause I want to have the opportunity to decide by myself if the outcome is satisfactory and save it manually when it fits my needs.


texts: Combining wildcard command along with %s in where clause statement for ArcPy?

texts:
I am trying to write a script where the user can search for a five digit zipcode using an input statement and a search cursor. The zipcodes in the file though are nine characters long (i.e. xxxxx-xxxx) so I'm trying to use the wildcard character after five digits in the where clause  but I can't figure out either where to put it or the syntax I'm supposed to use. Currently I can get the use input zipcode to print out but when I try putting the % wildcard in I get an error. Can anyone help me?

Here is my code so far:

import arcpy

#set workspace
work=raw_input("work= ")
arcpy.env.workspace=work
zip5=input("zip5= ")
strZip= str(zip5)

def searchHospitals(work,zip5):
    fc= "Hospital.shp"
    whereClause= '"ZIPCODE" LIKE'+"'%s'"% strZip
    print whereClause


texts: How to automatically draw lines with a given angle from one single point in QGIS?

texts:
I need to draw a bundle of lines all originating from the same point and with a 10° angle between each line (e.g. lines bearing at 0,10,20 ..90..110..220..350).
I can't decide a priori the length of each line since I need to intercept the contour of another layer (i.e. coastlines) therefore I guess I will use a really long length.

Is there a way to automatically draw these lines in a shapefile? Do I have to write a specific script or are plugins available? I do not know python but I can work in R if needed.

I've seen similar question have already been asked but always with different softwares or in slightly different ways (e.g. lines from one point to another) so that none of them was useful to me.


texts: Make part of legend label italics either for layer header or items in QGIS Print Composer

texts:
I am using QGIS 2.18.  I want to make PART of a legend label/heading italic within print composer when creating a map layout (for example make part of my label italic if its a botanical name and keep the non botanical part normal text). I know it is possible to format entire text for Titles, Subgroups Groups and Items but this apply the style to the entire label.  Is label expressions like ArcGIS implemented for legend text and headings to bold/italic/superscript part of a legend text?

I tried applying HTML formatting to the legend labels, without any success in the style properties and the print composer legend properties because I know it can be used for text boxes in Print Composer.

I can go back to ArcGIS but this will defeats my purpose.

Note: A similar question was asked in 2015. Is it possible to print one item legend item in italic in print composer? but it does not cover the issue that I need only a part of the label in italics.


texts: Creating Left and Right Image pair to view with 3D goggles in stereo from DEM or Hillshade?

texts:
How does one create the stereo pair for viewing in a program such as StereoPhoto Maker?

I have processed lidar that I'd like to turn into hillshades and then view in this kind of software, with 3D scope goggles, two images on my screen, etc.

I see a lot of questions/answers on how to create a DEM from stereo pairs, but am having trouble finding techniques to actually make the left and right images from terrain data in or outside of a GIS software...unless I'm missing something obvious, any ideas?


texts: I have several layers on my map. I need to set minimum and maximum scale level to prevent zooming. What's the best way of implementing it?


texts: Selecting Feature Class from drop down list of all Feature Classes in Geodatabase In ModelBuilder?

texts:
Is there a feature of ModelBuilder that allows a user to select a feature class based on a list of all feature classes in a geodatabase, ie upon opening the tool, the user will select the geodatabase and then will see a list of all feature classes available to run the tool on.

Iterate doesn't seem to work, because I can't find a way to input the results into the a list in the parameters. Multivalue doesn't work either if I have to manually populate the list.


texts:
So I have a list like this:

list = [[[-3, -2, 2], [3, -2, 2], [3, 2, 2], [-3, 2, 2]], [[3, -2, -2], [-3, -2, -2], [-3, 2, -2], [3, 2, -2]], [[-3, -2, 2], [-3, 2, 2], [-3, 2, -2], [-3, -2, -2]], [[-3, 2, 2], [3, 2, 2], [3, 2, -2], [-3, 2, -2]], [[-3, -2, -2], [3, -2, -2], [3, -2, 2], [-3, -2, 2]]]

Now I want to create a 3D Shapefile with Z-Values out of it. I tried several approaches. The one that works best is this:

import arcpy
arcpy.env.overwriteOutput = True
arcpy.env.outputZFlag = "Enabled"

features =[]

for feature in list:
    # Create a Polygon object based on the array of points
    # Append to the list of Polygon objects
    features.append(
        arcpy.Polygon(arcpy.Array([arcpy.Point(coordPair[0], coordPair[1], coordPair[2], coordPair[2]) for coordPair in feature])))

# Persist a copy of the Polyline objects using CopyFeatures
arcpy.CopyFeatures_management(features, "C:/Users/test_shape.shp")

What I get are flat Polygons which don't show Z order M Values in ArcScene. But at least they show up in 2D.

I already tried this one:

arcpy.Polygon(arcpy.Array([arcpy.Point(coordPair[0], coordPair[1], coordPair[2], coordPair[2]) for coordPair in feature]), None, True, True))

But then no polygon shows up.


texts: Enable load extension in Python 3.6 mac osx

texts:
As explained here : https://gis.stackexchange.com/a/277071/20802 the standard Python framework on Mac OS X is not compile with the load extension function. Therefore it is impossible to use sql queries with this package and within QGIS (2.18 or 3.0) on Mac OSx.

The kyngchaos version of QGIS is compile to work with the standard version of Python (http://www.kyngchaos.com/software/qgis) but as this problem arise could someone tell me how to point out the pip version of Python to QGIS rather than the standard one? Or any other solution

https://gis.stackexchange.com/users/55203/steven-kay gave a link to an older topic which might was working with older python version. but wasn't successful with Python 3.6


texts: MapBox GL - filtering GeoJSON by feature index

texts:
I'm using MapBOX GL JS with GeoJSON data source. The content of the GeoJSON is changing dynamically. It is a FeatureCollection which contains polygons. I would like to update my MapBox data source very frequently (60 ms). Because the size of the GeoJSON is growing I want to filter the data to maintain the performance of my application. Thus, I would like to pre-filter the data and use only the last 100 polygon from my GeoJSON.
This is original GeoJSON with only one element:
var geojson = {'type':'FeatureCollection','features':[{'type':'Feature','geometry':{'type':'Polygon','coordinates':[]},'properties':{'c':'0'}}]};

I tried to pre-filter the data before updating the source. I can do the filtering based on the category of the features:
 var newGeoJSON = geojson;
newGeoJSON2.features = geojson.features.filter(feature => feature.properties.c ==='1');

And set the map source:
map.getSource('geojson').setData(newGeoJSON);

It works, however, I would like to filter the GeoJSON by the index of the features. For example filter out the last 100 elements of the array.
Is it possible? Do you have eny suggestions how to do it?


texts:
Is there a way to use a text editor to edit the various layers of an ArcGIS MapServer?

This is the service I'm working with.

Here's my scenario: I'm using ArcGIS Online to create a custom WebMap that will be viewed in ArcGIS Explorer for iOS. The MapServer layers already have symbol and label info included, but because the iOS device has a Retina display, the symbol widths and text sizes are unreadably small. My solution was to globally double these values, but I can't do so in ArcGIS Online without losing the existing symbology (things like highway shields that aren't editable in ArcGIS Online).

As a non-professional Mac user, I don't have any other ArcGIS applications at my disposal, but can I download the layers somehow, edit them in a text application, and re-upload to the WebMap from a file?


texts: Copy specific cell value from one feature class and append to another feature class

texts:
I want to copy specific cell value from one feature class and append to another feature class. I have the following code, however, when I run this, it does not produce the results, nor shows any errors.

import arcpy
import os
import math
from arcpy import *
from arcpy import env

env.workspace = "C:\RG\Kings County"
outworkspace = "C:\RG\Kings County\Test.gdb"

arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = 1 # this code overwrites existing output

InFC = "C:\RG\Kings County\LuChange.gdb\NVF_050413"
OutFC = "C:\RG\Kings County\Test.gdb\Summary"

try:
    rows = arcpy.da.SearchCursor(InFC, ("Field_ID"), '"Final" = 1')
    for row in rows:
        CopyFeatures_management(InFC, row[0], OutFC)
    del rows, row

except Exception as e:
    print e
    #print arcpy.GetMessage()


texts: "Other" category in combo box replacement by line edit on QGIS Custom Form?

texts:
I am trying to create a custom form for QGIS 2.0 with comboboxes containing an "other" category. I wonder how to provide a line edit to enter info replacing the "other" choice if chosen. Not sure if I can sync a combo box and a line edit at the same time under field properties.

In pseudo code:

myfield = dialog.findChild(QComboBox, "Myfieldname")
otherfield = dialog.findChild(QLineEdit, "Myfieldname")

if myfield.currentIndex() == -1:
    set myfield to otherfield.text()

This thread's idea seems very similar, but I'm having trouble adapting it: How to automatically populate fields instantly?

my custom form code so far:

Naturefield = None
Nat_autrefield = None

def formOpen(dialog, layerid, featureid):
    global myDialog
    myDialog = dialog
    global Naturefield, Nat_autrefield
    Naturefield = dialog.findChild(QComboBox, "Nature")
    Nat_autrefield = dialog.findChild(QLineEdit, "Nat_autre")
    Naturefield.currentIndexChanged.connect(autreChanged)

def autreChanged(value):
    #intending to take on form of line edit only if 'Autre' chosen
    if Naturefield.currentText() == 'Autre':
        Naturefield.setLineEdit(Nat_autrefield)

setLineEdit() currently transforms my combo box to make the list editable for all choices if 'autre' is chosen and the combo box stays a line edit once 'autre' is chosen--not exactly what i'm aiming for. Instead, I'd like the combo box to receive the input from a line edit only if 'autre' is chosen, and revert to normal combo box if other choices are chosen. This is all in the interest of data quality control.

I'm new to QtDesigner and Python.


texts: QGIS changes CRS for new layers using GRASS toolbox functions. Is that normal?

texts:
I noticed that QGIS changes the CRS of newly generated layers, that is, I have a layer in my desired CRS, then I use a toolbox function with that layer as input, and the resulting layer is often in a different CRS. In my concrete example, I downloaded National Map elevation data and added it as a first layer, so the project CRS was EPSG:4269. I then reprojected that layer to EPSG:3719. I open a new QGIS project and add the reprojected layer, because I want EPSG:3719 as my project CRS. Then I use r.fill.dir on the EPSG:3719 layer and the resulting layer CRS is EPSG:2956. Then I use r.watershed with the depressionless DEM and the resulting flow direction and accumulation layers are in EPSG:2956. I have the setting checked that new layers should get the CRS of the project, but apparently it's not doing that.

Please note that this is not the case when I use SAGA functions, i.e., the SAGA fill algorithm comes up with the correct CRS (EPSG:3719).

Am I missing something obvious? I could of course warp/reproject all new layers, but that is not the point...


texts:
I have two PostGIS tables with geometries. One with number attributes and one with letter attributes. Each table has two features (rows).

I would like to apply a operation on both stables so that I create the following result:

The resulting table has 7 rows and contains all attributes from the input tables. Some features have both the letter attribute and number attributes.

I struggle to obtain this result using PostGIS / SQL. So far I've only been able to create the features that have both letter and number attributes.

SELECT
    table1.letter,
    table2.number,
    ST_Intersection(table1.geom,table2.geom)
FROM
    table1,
    table2
WHERE
    ST_Intersects(table1.geom,table2.geom)

Is there a way to include the remaining features in the same query? I could also create multiple tables and perform a join but this is a non-parsimonious solution. I've also tried ST_Union to no avail.

In ArcMap and Geopandas this is often referred to as "Union" without dissolving however as far as I can tell ST_Union in PostGIS does something different.
http://geopandas.org/set_operations.html
https://desktop.arcgis.com/en/arcmap/10.3/tools/analysis-toolbox/how-union-analysis-works.htm


texts: Can LiDAR differentiate between rock and soil?

texts:
Can LIDAR by itself, tell the difference between rock (large and small) and soil?  I'm under the assumption that the answer is no. However, I wanted to check here with the experts as I was told by someone else (not an expert) that it is possible.

If it is not possible and you know of what sensor data for me to look into, that would be very helpful.


texts:
I come from the Java world and go straight to C# and ArcGis 10 (I have to make a simple task for the university). So this means I don't really know much of C# but I can live with it and don't know ArcGis API at all. All I have is an example from the classes. It is about how to filter the cities by coordinates X and Y and a Buffer. So we have layers City(Point) and States(Polygon). My task is to make other functionality, and by functionality I mean method, that takes the name of 1 state and get all the cities from state's area.
In the example code that we have some ISpatialFilter that we use with the buffer like this:

public List<City> SearchCities(...)
{
    IFeatureClass citiesClass = MyHelper.GetLayer(Layers.CITIES);

    ISpatialFilter spatialFilter = getSpatialFilter(centerX, centerY, distance);

    IFeatureCursor citiesCursor = citiesClass.Search((QueryFilter)spatialFilter, true);
    .....
    .....
}

private static ISpatialFilter getSpatialFilter(double centerX, double centerY, double distance)
{
    IPoint center = getCenter(centerX, centerY);
    IGeometry buffer = getBuffer(distance, center);

    ISpatialFilter spatialFilter = new SpatialFilterClass();
    spatialFilter.Geometry = buffer;
    spatialFilter.SpatialRel = esriSpatialRelEnum.esriSpatialRelIntersects;

    return spatialFilter;
}

private static IGeometry getBuffer(double distance, IPoint center)
{
    return (center as ITopologicalOperator).Buffer(distance);
}

private static IPoint getCenter(double centerX, double centerY)
{
    IPoint center = new Point();
    center.PutCoords(centerX, centerY);
    return center;
}

So I decided to reuse some of the code and only change the way the data is filtered. So, I came up with something like this:

private IQueryFilter GetQueryFilter(string stateAbbr)
{
    IGeometry stateGeometry = getGeometry(stateAbbr);

    ISpatialFilter spatialFilter = new SpatialFilterClass();
    spatialFilter.Geometry = stateGeometry;
    spatialFilter.SpatialRel = esriSpatialRelEnum.esriSpatialRelIntersects;

    return (QueryFilter)spatialFilter;
}

// TODO fix it
private IGeometry getGeometry(string stateAbbr)
{
    IFeature state = GetStateByAbbr(stateAbbr);

    return (state as ITopologicalOperator).Boundary;
}

private IFeature GetStateByAbbr(string stateAbbr)
{
    IFeatureClass statesClass = MyHelper.GetLayer(Layers.STATES);
    IQueryFilter queryFilter = GetStateQueryFilter(stateAbbr);
    IFeatureCursor stateCursor = statesClass.Search(queryFilter, false);

    return stateCursor.NextFeature();
}

private IQueryFilter GetStateQueryFilter(string stateAbbr)
{
    ISpatialFilter spatialFilter = new SpatialFilterClass();
    spatialFilter.WhereClause = StatesWhereBuilder(stateAbbr);

    return (QueryFilter)spatialFilter;
}

private string StatesWhereBuilder(string stateAbbr)
{
    return (MyHelper.isEmpty(stateAbbr)) ? "" : MyHelper.CreateEqualStatement(StateConst.STATE_ABBR, stateAbbr, true);
}

I know it is probably very amusing for some of you, but this is all I came up with :]  The problem is that I get NullReferenceException after calling (state as ITopologicalOperator).Boundary but the debug shows that state is not null. So after the as something happens. I might be wrong with the ITopologicalOperator. If so, what should I use instead? Thanks in advance.

EDIT: Please tell me what do you think about this question (how it is exposed, etc.)

I changed the GetGeometry() method to:

private IGeometry getGeometry(string stateAbbr)
{
    IFeature state = GetStateByAbbr(stateAbbr);
    IPolygon statePolygon = state.Shape as IPolygon;

    return (statePolygon as ITopologicalOperator).Boundary;
}

and now I don't have NullPointerException, but I don't get any result.


texts: ArcGIS Server WMS content URL for reverse proxy

texts:
I configured my ArcGIS Server behind the IIS reverse proxy. My ArcGIS Server machine named as gisserver . And I changed my ArcGIS Server manager system property:

WebContextURL : https://myservices.domain.com/arcgis

So, when I get the ArcGIS WMS capabilities XML document from the internet using reverse proxy, the WMS inner URLs are following:

<OnlineResource xmlns:xlink="http://www.w3.org/1999/xlink" xlink:type="simple" xlink:href="http://gisserver:6080/arcgis/services/Countries/MapServer/WmsServer?"/>

How can I change urls as my domain as following?

<OnlineResource xmlns:xlink="http://www.w3.org/1999/xlink" xlink:type="simple" xlink:href="http://myservices.domain.com/arcgis/services/Countries/MapServer/WmsServer?"/>

ArcGIS Server version is 10.5.


texts:
I would like to use the google maps downloader to save satellite images from google maps (http://www.allallsoft.com/gmd/index.html), but is it legal to use?

In the terms of service you can read that: "Unless you have received prior written authorization from Google (or, as applicable, from the provider of particular Content), you must not: (a) copy, translate, modify, or make derivative works of the Content or any part thereof; (b) redistribute, sublicense, rent, publish, sell, assign, lease, market, transfer, or otherwise make the Products or Content available to third parties...(d) use the Products in a manner that gives you or any other person access to mass downloads or bulk feeds of any Content, including but not limited to numerical latitude or longitude coordinates, imagery, and visible map data..."

This program downloads the tiles and match them. The internet explorer does the same. When I view google maps in internet explorer, it downloads the tiles to "C:\Users\username\AppData\Local\Microsoft\Windows\Temporary Internet Files" and put them together. So it copies imagery to folder and mass downloads them.


texts: Is there a Raster Equivalent of the Split tool?

texts:
I have the DEM of a State in ERDAS Imagine (.img) format. I have a shapefile containing the counties of that state.

I want to split the whole raster into DEMs of the individual counties. The counties are irregular polygons. So this isn't as easy as splitting the Raster into tiles.

How do I split a raster with a vector polygon shapefile?

I am looking for what is basically the Raster equivalent of the Split tool.

I tried using Clip tool, and the Extract by mask tool, but both of them are unhelpful since they work only on single polygons.

I have looked at the answer given here: Need to clip raster based on field name in ArcGIS 9.3
, but Hawth Tools aren't available for ArcGIS 10.0

I have ArcGIS 10.0+ Spatial Analyst Extension, and QGIS 1.8. A solution in either will be acceptable.


texts:
I created a little plugin to geocode addresses using the HERE API and QGIS 3.

Especially in batch mode, I can see that the plugin is iterating over my features in a table and my messageBar progress is working as expected.

After a few seconds, QGIS gets unresponsive and the progressbar stops. This happens in a non-reproducible way.

The relevant code section is:

    progressMessageBar = iface.messageBar().createMessage("Looping through " + str(layer.featureCount()) +" records ...")
    progress = QProgressBar()
    progress.setMaximum(layer.featureCount())
    progress.setAlignment(Qt.AlignLeft|Qt.AlignVCenter)
    progressMessageBar.layout().addWidget(progress)
    iface.messageBar().pushWidget(progressMessageBar, level=1)
    for id in range(0, layer.featureCount()-1):
        urlPart=""
        oldAddress=""
        for key in addressLists.keys():
            if key != "oldIds":
                urlPart+="&" + key +  "=" + addressLists[key][id]
                oldAddress += addressLists[key][id] + ","
        url = "https://geocoder.api.here.com/6.2/geocode.json?app_id=" + appId + "&app_code=" + appCode + urlPart
        r = requests.get(url)
        if r.status_code == 200:
            if len(json.loads(r.text)["Response"]["View"])>0:
            #as the response may hold more than one result we only use the best one:
                responseAddress = json.loads(r.text)["Response"]["View"][0]["Result"][0]
                geocodeResponse = self.convertGeocodeResponse(responseAddress)
                lat = responseAddress["Location"]["DisplayPosition"]["Latitude"]
                lng = responseAddress["Location"]["DisplayPosition"]["Longitude"]
                ResultFet = QgsFeature()
                ResultFet.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(lng,lat)))
                ResultFet.setAttributes([
                    addressLists["oldIds"][id],
                    etc...

                ])
                ResultFeatureList.append(ResultFet)
        progress.setValue(id)
        iface.mainWindow().repaint()
    pr.addFeatures(ResultFeatureList)
    iface.messageBar().clearWidgets()

Here is the behaviour in a gif:

The main functionality below is embedded in a function which is connected with the button. The connection is defined in the main run function of the plugin:

    def run(self):
        self.dlg.show()
        self.dlg.LayerSelect_2.currentIndexChanged.connect(self.searchFieldsPopulate)
        self.dlg.batchGeocodeFieldsButton.clicked.connect(self.batchGeocodeFields)
        result = self.dlg.exec_()

        if result:
            pass

UPDATE
When adding a short sleep prior and after progressbar update the progress bar seems to update better and the system is responding and not showing any unresponsive messages in the plugin/QGIS application.

    time.sleep(0.2)
    progress.setValue(id)
    time.sleep(0.2)


texts: ArcGIS 10.5 - Filter raster datasets by it's ISO 19115 metadata

texts: I'd like to use search engine of ArcGIS 10.5 to filter raster datasets by it's ISO attributes (for example resolution). Search engine by default can filter only by tags, summary, description and credits of datasets. Do you have any idea how to solve that problem?


texts: OpenLayers Hide/Show Layers from external checkboxes

texts:
I'm using OpenLayers with mapserver to display a lot of raster images and queryable wms layers. Everything works fine, the layers can be turned on/off in the OpenLayers layer switcher...But now I want to turn on/off the layers by external checkboxes and I don't know how to do this.

I'm little bit new to Javascript-the only thing I know is that I have to pass the checkbox layer value to an OpenLayers function. Till now I tried to call a function with an onClick event

<input type="checkbox" name="ckeckbox_layer1" onClick="layervisibility(this.checked);">

But I don't know how the function 'layervisibility' has to look like so that my layers that are actually placed in the openlayers init() function turn on/off correctly when I check/uncheck the checkboxes.


texts:
Is there a way to do a "left join" when using the "merge" command on a geopandas df to merge by attribute?

The problem I am running into is that I am joining a geopandas df with CA counites to a pandas df that does not contain all the counites. However, I want the final geopandas df to contain all of the counties, but by default the merge command does an inner join. There is no information on the geopandas merge documnetation page: http://geopandas.org/mergingdata.html


texts: Removing property from Feature (or FeatureCollection) using Google Earth Engine?

texts:
How do I remove a property from a Feature (or FeatureCollection), by only specifying the name of the properties I want to remove? I know I can "remove" by selecting only properties I want, but in my case it is quite cumbersome...
For the minimal example in code below (see also link), how to remove property color?
var features = [
  ee.Feature(ee.Geometry.Point(-73.96, 40.781), {name: 'Thiessen', color: "a"}),
  ee.Feature(ee.Geometry.Point(6.4806, 50.8012), {name: 'Dirichlet', color: "a"})
];

// Create a FeatureCollection from the list and print it.
var FC = ee.FeatureCollection(features);
print(FC);

// Now I want to remove property color in each feature?


texts: How to mix two raster in R (same extent)?

texts:
I want to produce one PNG image for each country in the world. The pixel colors come from a raster that describes the vegetation, and the boundaries of the countries come from a shapefile.

I managed to do almost all of it using the raster and rgdal libraries in R:

library(raster)
r <- raster('gm_lc_v3.tif')
library(rgdal)
p <- readOGR('ne_10m_admin_0_countries','ne_10m_admin_0_countries')
i <- 66 # Dominican Republic, a small country for fast example
country <- crop(r,p[i,]) # produces the rectangular crop below
plot(country)

m1 <- mask(country,p[i,]) # crops on the country boundary
plot(m1)

m2 <- mask(country,p[i,],inverse=T) # crops outside the boundary

m2[m2 < 20] <- NA # turns all land into NA

What I need now is to combine the two rasters (m1 and m2) into a single raster, with all data from the country, and only land/water data from outside the country. Like this (made in Gimp):

I have tried all the following so far, but without success:

m3 <- cover(m1,m2)
m3 <- cover(m1,m2,identity=T)
m3 <- merge(m1,m2)
m3 <- overlay(m1,m2,fun=function(x,y){return(x+y)})
m3 <- mask(m1,m2)
m3 <- stack(m2,m1)

They return all black, or different colors from the original. What am I missing?


texts: Removing part of string in ArcMap using Python?

texts:
I have a Field in this table which reads, and this is a sample because there are over 82,000 tracks in this county, "40 acres".

I do have this Script which I know took the "Acres" portion out of the field and the number portion is the only part of the field I need, and this looks like this.

This is the way I have entered this Script:

But I keep getting this error message:

What am I doing wrong?


texts: I am new to geographical statistics so please pardon this question if it is a duplicate. I would like to test if points that are similar are more or less likely to be near each other. I am latitude and longitude data for several alligator nests. I need to know if a singly sired nest is more likely to be near another singly sired nest or not. What would be the most appropriate way to formally test this?


texts: Multi Labelling with missing values in QGIS

texts:
The code for multi labels eg.:
Field1 || '\n' || Field2

works fine, unless you have a missing value (NULL) in one of your fields, in which case QGIS doesn't show any labels at all (for the specific object).
I use QGIS 3.10.14-A Coruna.
Any suggestions?


texts: "SDE error 6" when publishing GeoServer layer

texts:
I have successfully created an ArcSDE datastore in GeoServer 2.2.4. When I try to publish a feature class in the SDE database, the "Compute from data" link returns nothing, and the "Compute from native bounds" link produces the following error.

Error computing the geographic bounds:[SDE error 6][Error desc=UNKNOWN ERROR CODE. ][Extended desc= ]

The Native SRS says UNKNOWN in the text box but has the correct system to the right "GCS_WGS_1984...".

Because I do not have either of the bounding boxes filled in, obviously I can not save the layer.

I have tried feature classes with other coordinate systems, deleting and re-creating the datastore and creating whole new feature classes. Extensive searching on the internet has not turned up any kind of solution.

Any help would be great! Thanks!


texts:
Is there a way to preserve the mask of a numpy array when reprojecting it with rasterio?
See the reproducible example with real data below:

import os
import tempfile

import rasterio
from rasterio.crs import CRS
from rasterio.warp import reproject, calculate_default_transform
import matplotlib.pyplot as plt
import numpy as np
import requests

print(rasterio.__version__) #1.0a10

# Download function, only downloads the file once
def download(url):
    filename = os.path.join(tempfile.gettempdir(), os.path.basename(url))
    if not os.path.isfile(filename):
        r = requests.get(url, stream=True)
        with open(filename, 'wb') as dst:
            for chunk in r.iter_content(chunk_size=1024):
                dst.write(chunk)
    return filename

# Download test data (14MB)
scene_id = 'LC08_L1TP_026047_20170616_20170629_01_T1'
url = 'https://earthexplorer.usgs.gov/browse/gisready/landsat_8/%s.zip' % scene_id
archive = download(url)

# Read data as numpy array and apply mask
rgb = '/vsizip/%s/%s.tif' % (archive, scene_id)
mask = '/vsizip/%s/%s_QB.tif' % (archive, scene_id)
with rasterio.open(mask) as src:
    m = src.read(1)
with rasterio.open(rgb) as src:
    src_crs = src.crs
    src_transform = src.transform
    src_bounds = src.bounds
    src_width = src.width
    src_height = src.height
    red_array = np.ma.array(src.read(1), mask=m!=0)

# VIsualize masked array (clouds and cloud shadows are properly masked)
# plt.imshow(red_array)
# plt.show()

# Prepare destination array for warping
dst_crs = CRS({u'proj': u'laea', u'lon_0': -106, u'lat_0': 24}) # Equal area crs more or less centered on mexico center
dst_transform, dst_width, dst_height = calculate_default_transform(src_crs,
                                                                dst_crs,
                                                                src_width,
                                                                src_height,
                                                                *src_bounds,
                                                                resolution=(30, 30))
red_laea = np.ma.zeros((dst_height, dst_width), dtype=np.uint8)

# reproject array
reproject(source=red_array,
          destination=red_laea,
          src_transform=src_transform,
          src_crs=src_crs,
          dst_transform=dst_transform,
          dst_crs=dst_crs,
          dst_nodata=0)

# Visualize (clouds are no longer masked)
plt.imshow(red_laea)
plt.show()

What I would like/expect is that the mask (cloud and shadows in the example above) is transferred to the destination array. I was sort of expecting that behaviour from the description of the src_nodata argument (rasterio.warp.reproject function), but I may have misinterpreted it.


texts: How to test whether two countries are adjacent using ArcGIS Desktop?

texts:
I have an interesting question:
Is there a dataset by which I can use to test whether two countries are adjacent? For example, if inputs are Mexico and U.S., they are adjacent; if inputs are England and Russia, they are not. Is it easy to produce an adjacency matrix for all the countries?

Edited: This question might be extended as, if we get a political boundary map, is there any ways to produce an adjacency matrix for all the political entities (e.g. Given a map of Florida, can a counties matrix be built to indicate the adjacency of all the counties)


texts: Getting findAndReplaceWorkspacePaths to work?

texts:
I am copying and pasting a gdb and an mxd into new folders using arcpy.  I then want to change the data sources of the layers in the mxd to point to the new version of the gdb instead of the old one.

My code is:

mxd = arcpy.mapping.MapDocument("pathtomynewmxd")
mxd.findAndReplaceWorkspacePaths("pathtomyoldgdb", "pathtomynewgdb")

this runs successfully, but nothing actually happens.
Additonal info: arcGIS 10.3, the mxd does store relative path names already.


texts:
I have some DEM files that use a coordinate system that doesn't align well with WGS84. I can align them through the data frame properties -> Coordinate Systems Transformation, where i type Only the x,y,z correction in meters.

I want to make this transformation permanent (create a new DEM file that aligns well with WGS84) but dont seem to find a tool to do this like i do in the data frame properties (by using the 2 coordinate systems and the x,y,z correction in meters).

Any ideas?


texts: gdal.VectorTranslate returns an empty file

texts:
When testing the following code:

from osgeo import gdal
import json

def read_file(filename):
    f = gdal.VSIFOpenL(filename, "rb")
    if f is None:
        return None
    content = gdal.VSIFReadL(1, 10000, f).decode('UTF-8')
    gdal.VSIFCloseL(f)
    return content

gdal.VectorTranslate('out.json', """{
  "type": "FeatureCollection",
  "features": [
      { "type": "Feature", "geometry": { "type": "Polygon", "coordinates": "coordinates":[[[-188.56933593749997,61.58549218152362],[-173.9794921875,61.58549218152362],[-173.9794921875,66.00015035652659],[-188.56933593749997,66.00015035652659],[-188.56933593749997,61.58549218152362]]] } },
  ]
}""", format = 'GeoJSON', layerCreationOptions = ['RFC7946=YES', 'WRITE_BBOX=YES'])

got = read_file('/vsimem/out.json')
gdal.Unlink('/vsimem/out.json')

print json.loads(got)

from this answer, the file out.json is empty. The folder is with all permissions for writing allowed. When testing the method gdal.OpenEx, is possible to load an object, but still the method gdal.VectorTranslate returns an empty geojson when trying to pass this new loaded object as a srcDS object. Any hint on solving this issue?


texts: How to reset the primary key (sequence) for a PostGIS table in PostgreSQL

texts:
I am uploading 64 features from 64 shapefiles into the postgis server (I use INSERT TO to append them one by one), for some reasons I made some mistakes at the first attempt so I deleted them (except where id = 1), then I re-upload them. However, the id (primary key generated automatically when uploading) does not start from 2 but from 65.

Why is postgis doing this? How do I change it back to start from 2? Also is there any way that I can arrange the table by the id in an ascending order (so that when the user drag it into QGIS and open the attribute table, they can see that the features are arranged ascendingly)?


texts: Transforming Shapely Polygon and MultiPolygon objects

texts:
Is there an easy way of transforming Shapely objects (namely, Polygons and MultiPolygons) from one projection to another without having to dig around and extract coordinates by hand?

In fact, I don't even care if they are Shapely objects at this point, I just want to pass features and a projection, and get a reprojected set of features back.

Does this sort of functionality exist, or must it be hand coded?


texts: Offline SLIPPY Map (Tiles) Database for Leaflet

texts: How can create an offline cache of a large number of SLIPPY tile (PNG, 256x256) inside a web browser and then display them from that offline cache. For example, can one use HTML5 IndexedDB as an storage location for offline map tiles?


texts: Assigning UNC Path to arcpy.env.workspace

texts:
For some reason my python script crashes ever time it reaches the following code which attempts to assign a unc path as my workspace.  Can any of you please help me solve this problem?  Also, I'm not getting any error code when PythonWin crashes, it just does and yes I'm running my code in PythonWin.

arcpy.env.workspace = r'\\ARSENAL.LAMB.CRAYON.NET\SHARE\LEY\CTPCW3\DATA\GIS_TOTA\Projects\TOTA\2014\Data_for_Submission\2014_MASTER_TOTA.gdb'

Currently using Arcmap 10.0
PythonWin 2.6.5 (r265:79096, Mar 19 2010, 21:48:26)


texts: How to create selective output from a QGIS layer, using QGIS2leaf?

texts:
I want to create a web map on which only a selection of data from a layer is displayed. To the best of my knowledge the tool seems to output the full layer.

In Mapinfo I'd run a query to create a new table which would display only the selection of data I want. I've not done anything similar with QGIS - and I don't know if it would work if I did.

I realise I could pursue a non-elegant solution in the short term - extracting data from the layer and saving it as a new file... but there must be a more flexible/elegant way (because this layer would need to be updated each time the parent layer changed).


texts: Getting value set in binding input widget in QGIS expression?

texts:
I would like to set dynamic line width and it should be scaled by the value "Stroke width".

When I set the value to 0.15, it will be less than 0.15 if $length is too small. I need a function to get the value set in "Stroke width" so that I can modify it directly in "Stroke width" instead of in Expression String builder.

CASE
  WHEN ($length < 1) THEN ($length * 0.15 / 1)
  WHEN ($length >= 1) THEN (0.15)
END


texts:
I hesitated a bit before asking this, as it isn't necessarily a GIS specific question, though I'm using it to update geodatabases. I have found plenty of hits on a Google search, but frankly, many of the websites out there for .bat files are hacky and not real trustworthy.

Do you have a resource that you use for database manipulation through batch files?

Any GIS/Geodatabase specific resources for .bat's?

A little background, I'm just starting to dabble in batch files. I'm aware that I can automate many things that are done manually, but am terrified of doing anything I shouldn't be doing/cannot undo. I have found introductory sites that don't give much more information than printing to the console. I'm looking for something that goes into detail about databases.


texts: Processing vector to raster faster with R

texts:
I am converting vector to raster in R. However the process was too long. Is there possibility to put the script into multithread or GPU processing in order to do it more faster?

My script to rasterized vector.

r.raster = raster()
extent(r.raster) = extent(setor) #definindo o extent do raster
res(r.raster) = 10 #definindo o tamanho do pixel
setor.r = rasterize(setor, r.raster, 'dens_imov')

r.raster

class       : RasterLayer
dimensions  : 9636, 11476, 110582736  (nrow, ncol, ncell)
resolution  : 10, 10  (x, y)
extent      : 505755, 620515, 8555432, 8651792  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0

setor

class       : SpatialPolygonsDataFrame
features    : 5419
extent      : 505755, 620515.4, 8555429, 8651792  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=24 +south +ellps=GRS80 +units=m +no_defs
variables   : 6
names       :    ID,CD_GEOCODI,   TIPO,  dens_imov,area_m,domicilios1
min values  : 35464, 290110605000001,  RURAL, 0.00000003,100004,1.0000
max values  : 58468, 293320820000042, URBANO, 0.54581673,99996, 99.0000

Print of setor


texts: ArcGIS Online: User manually add object/feature to see effect on surroundings

texts:
I would like to make an web-app, I know the basics, but is it possible for users to  (1) add a temporary feature (in this case: physical building) and (2) directly solve/update the map to see the effects on the surroundings.

Ive done this in excel. Its a bit conditional. So whenever you add a feature, values around the grid will be affected and will directly update/solve in the web-page.

Is this possible, if so: are there examples or ways to do this?


texts: Make WMS layer background transparent in Leaflet

texts:
I have inserted the following layer from a GeoServer WMS in the following Leaflet map

:

var nexrad = new L.TileLayer.WMS("URL", {
    layers: 'Wrecks:WrecksGreaterNorthSea',
    format: 'image/png',
    transparent: true
});

but the background is not transparent. I have tried to set it by modifying the tiff image with QGIS creating a transparency band but it does not work. I have a second WMS layer in the map but I do not have this problem with this one.

Could you help me solve this issue, please?


texts: How to start WFS with layers that are stored in a particular workspace?

texts: I would like to start WFS with layers that are stored in particular workspace. Now my WFS address (ex. domain.com/geoserver/wfs) gives access to layers saveed in all workspaces. Wher to change it?


texts: In ArcGIS raster calculator can you access the individual cells surrounding a central cell like you can in GRASS r.mapcalc

texts:
Is there a way to do the following in Raster Calculator in ArcGIS?
In GRASS r.mapcalc you can access the cells surrounding the current cell by writing  inputRaster[a,b] where a and b represent the number of cells away from the current cell. So you can write complex map algebra expressions using the surrounding and central cell, not just what is directly above or below it in the raster stack. For example if I wanted to estimate the 3d surface area of a DEM, I could add up the the area of each of the 8 triangles formed between the central cell and each pair of its 8 neighbors (each triangle of course would only go half way out to the surrounding cells).
Is there a way to do this in Raster Calculator in ArcGIS? Or some other tool, or in python maybe?


texts:
following this tutorial I set up a OpenTripPlanner (version 0.11.0) stand alone server.
It works very well at

 http://localhost:8080/index.html

and for example at

http://localhost:8080/otp/routers/default/index/routes

returns a routes json.

Now I'd make rest api request (I'm interesting to Isochrone api).
I try return isochrone geometry as a set of GeoJSON polygons at the url specified at this link but doesn't work.

Can you suggest the right url to get Isochrone GeoJSON?


texts: Using a field to define label font size in MapInfo?

texts:
I have a tab file of ordnance survey MasterMap data, this was obtained from a client so I have no way of re-downloading it.

Usually when I have these files the text is automatically displayed according to a number of fields (e.g size, angle etc.). However this file I have everything is displayed as lines.

I have used style override to remove the lines and have turned labelling on.

However each label has the same text size and makes the map quite unreadable when zoomed out to the level I require.

Is there any way in MapInfo to set the font size/height to that defined in a field in the table for the data?


texts: Filtering duplicate features on mapbox vector tile source

texts:
I am using mapbox and have a vector tiles source.
How can i filter out duplicate  features in order to display only unique feature?
I checked the expressions in the filter and i cannot see how can this be done
Thanks


texts: Using Windows folders as ArcPy Parameters?

texts:
I am having problems trying to perform simple folder management operations in an ArcGIS python toolbox.  First here's some code snippets (I have left out all of the stuff that does not pertain to this problem):

import os

def getParameterInfo(self):
    #define the parameters for the toolbox

    param = arcpy.Parameter(displayName="Destination Folder",
        name="in_destination",
        datatype="DEFolder",
        direction="Input")

    params = [param]
    return params

def execute(self, parameters, messages):

    #get the destination root folder parameter
    #I have also tried using .valueAsText but it does not solve the problem
    dest_root = parameters[0].value

    #if the folder does not exist create it
    if not os.path.exists(dest_root):
        os.makedirs(dest_root)

    #join path to store shapefiles to the root folder
    curr_path = os.path.join(dest_root,"/Shapefiles")

    #create folder to store shapefiles
    if not os.path.exists(curr_path):
        os.makedirs(curr_path)

    #start doing a bunch of stuff to create shapefiles

The problem I am having is that the folder parameter that ArcGIS creates is an unescaped string, so if I wanted to use the C:\Temp folder as my destination root the dest_root variable in the execute function would have a value of 'c:\temp'.  This is causing problems for all of the calls to the os methods because (I think) they are escaping the characters that come after the backslash.

I have tried to replace the backslash character with a forward slash by doing this:

dest_root = dest_root.replace("\\","/")

but this does not solve the problem because the organization I'm working at uses employee ID numbers for the Windows profile user names.  So if I want to use the Documents folder for a destination the escapes are being treated like hex codes:

dest_root = "C:\Users\678910\Documents"
dest_root = dest_root.replace("\\","/")
#dest_root is now equal to "C:/Users78910/Documents"

Have I missed something basic here?  This has been driving me crazy all morning.


texts:
We are using OpenGeo suite 4.6 (GeoServer ver. 2.8) using PostreSQL ver 9.6 on Windows Server 2012 R 2.

We are publishing layers on GeoServer, while publishing layers found some strange layers (attached as screenshot)

Does any know what it is exactly?
is there any kind of layers or views?


texts: Creating layer in ArcMap from SQL Server database table?

texts: Is it possible to create a layer from an SQL Server table instead of a shapefile?


texts:
I have a TIFF file that is 2880x1440 pixels. It is a global map in simple cylindrical projection (i.e., each pixel is a 1/8th degree square):

  If the uppermost leftmost pixel has coordinates x = 1, y = 1, then

  north latitude (deg) = (720.5 - x) / 8

  east longitude (deg) = (y + 0.5) / 8

I'd like to use the projection information to make a GeoTIFF that I can load into ArcMap.

Is that possible or do I need to create a world file instead?

Note also that it is a map of Mars, so I will use the Mars2000 coordinate system (available in ArcMap).


texts: How can I resample statistics from one set of districts to another?

texts:
I have statistical data for a certain region based on a subdivision of that region into districts (e.g. number of schools per district). Now I want to resample that data to a different subdivision of the same region. The two subdivisions don't have anything in common aside from covering the same overall region.

My specific use case is to resample statistics based on city districts to voting districts. I have the outlines of the districts in two GeoJSON files (one for the city districts, one for the voting districts). The former also contains the statistics I want to resample (as properties of the GeoJSON features). Resampling based on the area shared between the two types of districts would be a good start, although I'd love to take the population count of the voting districts into account, too (that information is available in the second GeoJSON file).


texts:
I am currently trying to get my update cursor to populate every row in each attribute table for 175 different feature class fields. I have written the script and it is populating the first row but all the next rows are coming out "NULL" (blank).

This is my script

import arcpy, csv

geodatabase = "C:\Users\kd16342\TSB_Values_Project\Tool\RasterTool.gdb"
csvfile = "C:\Users\kd16342\TSB_Values_Project\Documents\Scripts\ServicesforMosaic.csv" #change "Ecosystem Services" to "Ecosystem"

fields = ["Ecosystem","Cur_Con","Scale_Efft","Eff_Con","Recover","Typ_of_ser","Services","Eff_on_ser","Dur_s_aff","Reason","Eff_Durat","Unique_ID","Serv_ad_af","Serv_pe_af"]

arcpy.TableToTable_conversion(csvfile, geodatabase, "ServicesforMosaic")

for x in range (9,176):
    with arcpy.da.SearchCursor(geodatabase+"//ServicesforMosaic",fields) as scursor:
        for rw in scursor:
            if rw[11] == str(x):
                with arcpy.da.UpdateCursor("C:\Users\kd16342\TSB_Values_Project\Tool\RasterTool.gdb\Unique_ID_"+x,fields) as ucursor:
                    for row in ucursor:
                        row[0] = rw[0]
                        row[1] = rw[1]
                        row[2] = rw[2]
                        row[3] = rw[3]
                        row[4] = rw[4]
                        row[5] = rw[5]
                        row[6] = rw[6]
                        row[7] = rw[7]
                        row[8] = rw[8]
                        row[9] = rw[9]
                        row[10] = rw[10]
                        row[11] = rw[11]
                        row[12] = rw[12]
                        row[13] = rw[13]
                        ucursor.updateRow(row)
                    del cursor, row

All helped welcomed!!!


texts: Why is my ArcPy script deleting its workspace?

texts:
I wrote an ArcMap script to go through a feature layer, identify unique values in a given field, and output a feature class from each unique value to its own feature class or shapefile. As an exercise, I decided to experiment with allowing the user to set a limit on the number of outputs, so the script would only take the first 10 or 20 unique values and create a feature class from each in a designated workspace.

However, if there are more unique values than the limit, the script deletes whatever geodatabase or folder is designated as its output location, along with any files in that location. Occasionally, the script will throw a 000210 error when creating its first file for no reason that I can determine (the workspace is not open in ArcCatalog or ArcMap, the file name is valid, and nobody else in my agency is using it) and if this happens, the output location is also deleted.

Below is my script in its entirety, except for a header that is all comments. Has anyone experienced similar behavior from a script? Any ideas why this is happening? I've written a number of scripts, both stand-alone and added to ArcMap, and I'm terrified to write any more until I know why this one is acting so oddly.

import arcpy, time, os
from arcpy import env

#set parameters
# SplitFeat - input features
SplitFeat = arcpy.GetParameterAsText(0)
# SplitField - input field
SplitField = arcpy.GetParameterAsText(1)
# OutLocation - output folder
OutLoc = arcpy.GetParameterAsText(2)
#MaxSplit = Maximum number of files to create
MaxSplit = arcpy.GetParameterAsText(3)
# Overwrite output

arcpy.env.workspace = OutLoc
arcpy.env.overwriteOutput = True

IndNum = 0

arcpy.MakeFeatureLayer_management(SplitFeat, "Split_lyr")
CatList = []
CodeList = []
FldList = (SplitField)
with arcpy.da.SearchCursor("Split_lyr", SplitField) as cursor:
    for row in cursor:
        if row[0] not in CatList and IndNum < MaxSplit:
            AttString = str(row[0])
            print AttString
            CatList.append(AttString)
            CodeString = ''.join(e for e in AttString if e.isalnum())
            CodeList.append(CodeString)
            IndNum = IndNum + 1
        elif IndNum == MaxSplit:
            arcpy.AddMessage("{0} classes identified. Proceeding with script...")
    else:
        "All classes identified. Proceeding with script..."

#splitting features
#new iteration number
IterNum = 0
for Cat in CatList:
    CodeName = CodeList[IterNum]
    ClassName = CodeName[0:30]
    ClassFile = os.path.join(OutLoc, ClassName)
    arcpy.AddMessage("Creating file {0}...".format(ClassFile))
    #select by attributes
    SelClause = "{0!s} = '{1!s}'".format(SplitField, Cat)
    arcpy.SelectLayerByAttribute_management("Split_lyr", "NEW_SELECTION", SelClause)
    #Copy features, using CodeList index as key
    arcpy.CopyFeatures_management("Split_lyr", ClassFile)
    arcpy.AddMessage("{0} has been created.".format(ClassFile))
    IterNum = IterNum + 1


texts:
I am using MapInfo to display some Bing aerial photos. It seems as if it is using the wrong scale of imagery - they look blurred, as if the image has been enlarged.

Is there a way to force it to use a larger scale of imagery and therefore get more detail?


texts: Correcting coordinates of CAD drawing imported into ArcGIS?

texts:
I am having some issues with some CAD drawings that I have imported to ArcGIS. I have a number of CAD drawings where by the features within the drawing have an error in their coordinates. For example, the CAD file has a coordinate like 308092578,52639792281 where as it should read 308092.578, 52639792.281.

The result of this error is obviously causing the features to be placed far outside its actual location.

Is there any way to correct this error outside of manually moving the features?


texts: How to make both hover and click popup work?

texts:
I am pretty new to all of this.

I want to have a hover/mouseover AND click popup work on my geojson layer.
Here is the code so far, but I don't know where to include the hover without messing up the click popup. The hover should show the content of a text field and highlight the circle marker.

$.getJSON('http://soundgoods.cartodb.com/api/v2/sql?format=GeoJSON&q=SELECT * FROM soundgoods_mixtape_map_1', function(data) {
var geojsonMarkerOptions = {
    radius: 8,
    fillColor: "#FFCC00",
    color: "#323232",
    weight: 2,
    opacity: 0.5,
    fillOpacity: 0.4
    };

var geojson = L.geoJson(data, {

    pointToLayer: function (feature, latlng) {
        var popupOptions = {maxWidth: 500};
        var popupContent = '' +
                feature.properties.soundcloud +
                '
' + "Post & tracklist" + '
' + '';
        return L.circleMarker(latlng, geojsonMarkerOptions).bindPopup(popupContent,popupOptions,{

        });

    }

});
markers.addLayer(geojson);

// CONSTRUCT THE MAP
var map = L.map('map').setView([0,0],3);
baseLayer.addTo(map);
markers.addTo(map);


texts: How to obtain the value of a feature programatically, based on a coordinate

texts:
I have a point layer, and a polygon layer, and I wish to calculate the value of an attribute from the polygon layer at each point's location. Eg, find the area of the polygon which is found at each point, and write it to the point's attribute table:

One simple approach is to perform a spatial join - but this creates a new output point layer. I would prefer to update the existing point layer.

I could iterate through the points layer and:

obtain each point's OBJECTID
create a layer from the selected point, using the OID as the Where clause
perform a Select By Location, using the above layer as the select_features
iterate through the selected polygon(s) and write the Area value to the point

Steps 2 and 3 seem a bit extraneous - is there a programmatic way to obtain the value of a layer based solely on an XY coordinate?

This would be akin to using the Identify tool in ArcMap, whereby clicking on a location reveals the attributes of the layers at that location.


texts: Transformation accuracy or offset in ArcMap?

texts:
I'm trying to transform a map from UTM ED 50 zone 33N to WGS84 and vice versa in ArcMap 10.2.2.

I'm doing this by using the Project tool in the Arc toolbox and choosing the ED_1950_TO_WGS_1984_12 transformation (since the area in question is Malta).  However when I compare the completed transformation with another map of the same data in the same projection, there is an offset of 22m.

Can anyone please indicate why this is happening, and how I can improve the accuracy of my transformation?


texts: Feature Service from plain PostgreSQL/PostGIS in ArcGIS Server 10.2

texts:
Using 10.2 ArcMap/ArcGIS Server on Windows I've successfully published some map services from my PostgreSQL/PostGIS data and now I want to enable editing the features through a JavaScript web app. This is NOT through SDE tables (i.e., creating a geodatabase). Yes, I am using PostgreSQL 9.3.1 and PostGIS 2.1... but everything else seems to work fine but this PK error. How do I get ArcGIS Server to trust me that the incrementing PK I've created is OK?

  "00166: Layer's ObjectID is not maintained by the database and the feature service is being published with Create capability enabled"

So, following the JS API examples I am enabling both the geometry service and the Feature Access service with all the Operations and Capabilities enabled.

I began using the serial data type for the OBJECTID primary key which appears to be a shortcut for creating an integer, auto-incrementing key as follows:

CREATE TABLE test_schema.volmonwgs84e
(
  objectid serial NOT NULL,
  volmon24_i integer,
  site_id integer,
  waterbody character varying(30),
  geom geometry(MultiPoint,4326),
  CONSTRAINT volmonwgs84e_pkey PRIMARY KEY (objectid)
)
WITH (
      OIDS=FALSE
);
ALTER TABLE test_schema.volmonwgs84e
  OWNER TO postgres;

CREATE INDEX volmonwgs84e_geom_gist
  ON test_schema.volmonwgs84e
  USING gist
  (geom);

However, when I go to publish the service (from ArcMap) the Analyzer shows the error above.

Right-clicking the error kindly offers to "Add an auto-incrementing ID" for me, to which it later responds with the following. AFAIK, I ALREADY have an auto-incrementing field anyhow.

"Executing: AddIncrementingIDField testing.test_schema.volmonwgs84e BOO
Start Time: Tue Dec 10 15:55:06 2013
ERROR 000152: You cannot add an incrementing ID field to a table registered with the geodatabase.
Failed to execute (AddIncrementingIDField).
Failed at Tue Dec 10 15:55:07 2013 (Elapsed Time: 0.70 seconds)[/CODE]

So, what data type and rules need to be on this auto-incrementing field, and what is AddIncrementingIDField trying to do for me?

UPDATE:

Thinking that my issue was the "new" PostGIS/PostgreSQL approaches to things I've also tried creating a table with these older approach as follows:

manually setting the sequence on integer PK column instead of using the SERIAL type,
using "SELECT AddGeometryColumn" to create the geometry column, and therein using the old constraints approach with the final argument = false (http://postgis.net/docs/AddGeometryColumn.html)

CREATE SEQUENCE test_schema.mypoint_const_b_seq
INCREMENT 1
MINVALUE 1
MAXVALUE 9223372036854775807
 START 1
CACHE 1;

 ALTER TABLE test_schema.mypoint_const_b_seq
 OWNER TO postgres;

CREATE TABLE test_schema.mypoint_const_b
(
  objectid integer NOT NULL DEFAULT nextval('test_schema.mypoint_const_b_seq'::regclass),
  volmon24_i integer,
  site_id integer,
  waterbody character varying(30),
  CONSTRAINT mypoint_const_b_pkey PRIMARY KEY (objectid)
)
WITH (
  OIDS=FALSE
);

SELECT AddGeometryColumn ('test_schema','mypoint_const_b','geom',4326,'MULTIPOINT',2,     false);

ALTER TABLE test_schema.mypoint_const_b
  OWNER TO postgres;

CREATE INDEX mypoint_const_b_geom_gist
  ON test_schema.mypoint_const_b
  USING gist
  (geom);

But I still get a fatal error when trying to create a service with Feature Access "Create" capability enabled with the same message: [CODE]"00166: Layer's ObjectID is not maintained by the database and the feature service is being published with Create capability enabled"
[/CODE]

I know that the database IS maintaining the PK because I can add features using QGIS just fine. AND, ArcMap's properties for the indexes on the layer shows objectid is a unique index.

PS, buried in this excellent link (http://proceedings.esri.com/library/userconf/proc13/tech-workshops/tw_192.pdf)  I see that Arc doesn't like fields with quotes (which is how you force CAPS on OBJECTID. Lowercase objectid appears to be just fine to ArcMap, but perhaps not ArcServer.


texts: Vector to grid - Computing statistics by grid cells

texts:
I have a vector point data from drifting buoys (speed, direction, sst, slp) that I want to analyse in a grid. What are the right tools for this?

Example of my points data below. Each point has the variables lat, long, date, time, SST, SLP, direction and velocity.

First I need to create a grid of the desired size (im thinking 1/4 degree)and then do basic statistics within each cell to create an overall picture. Basic stats would be a number of points in a cell, mean, max and min of the variables.

I tried to rasterize the vectors but can't figure out what the results are showing. I want to compute the means, max and mins for different variable in a cell. Any suggestions?


texts:
I have been using QGIS 2.8.1 to map an archaeological project and I have now installed Midvatten 1.4.6 with the aim of visualizing the stratigraphy from various test pits using qgis2threejs.

After creating the Midvatten database, I have made up some dummy test data for a single pit, loaded the pit name as "obsid" and its position into table "obs points" and the test stratigraphy data into "stratigraphy".  By selecting the test pit, I can obtain a stratigraphy plot from Midvatten.

My problem arises when I use "Midvatten> Utilities> Prepare 3D-data for Qgis2threejs plugin".  I see that the data layers are generated but the stratigraphy does not appear in the qgis2threejs visualization and on inspection I find that the 3D-data tables are completely empty of any attributes.


texts:
I have TanDEM-x data and according to the documentation:
"The elevation values represent the ellipsoidal heights relative to the WGS84 ellipsoid in the WGS84-G1150 datum."
I would like to convert this to EGM2008 but I cannot find the right starting EPSG code for the WGS84-G1150 datum. I have tried the following:
gdalwarp  -s_srs EPSG:7661 -t_srs EPSG:6893 TDM1_DEM__04_N04E006_DEM.tif TDM1_DEM__04_N04E006_DEM_EGM08.tif

While this runs and produces an output the elevation pixel values have not changed and so it has not done the datum conversion.
I have also tried:
gdalwarp -srcnodata -32767 -dstnodata -32767 -of GTiff -co COMPRESS=LZW -s_srs "+proj=longlat +datum=WGS84 +no_def" -t_srs "+proj=longlat +datum=WGS84 +no_defs +geoidgrids=egm2008.gtx" TDM1_DEM__04_N04E006_DEM.tif TDM1_DEM__04_N04E006_DEM_EGM08.tif

but this also does not do the datum conversion.
Does anyone know how to do this conversion in gdal or python?


texts: Requested operation is invalid on closed state using ArcSDE for Oracle?

texts:
I have a version on an Oracle SDE that I am working on using ArcMap / ArcPy 10.2.1.

I started an edit session and inserted a row in to a feature. Saved the edit and stopped editing.

I then tried inserting data through an arcpy script.

val = [uuid.uuid4(),'bryan', dt.strftime('%d-%m-%Y'),
       dt.strftime('%d-%m-%Y'),dt.strftime("%Y-%m-%d %H:%M"),
       'ins','qa','0','','done','xyz','test','','206901','5997594',
       'id1','id2','terr','co1','flag','issue','bid','official',
       'area','kv','sub','0','99',uuid.uuid4(), arcpy.Point(206901, 5997594)]
edit = arcpy.da.Editor(myconnect)
edit.startEditing(False, True)
with arcpy.da.InsertCursor("V.DATA/V.TABLE",fieldnames) as iCur:
    iCur.insertRow(val)
edit.stopOperation()
edit.stopEditing(True)

When the section of code executes below it fails at the iCur.insertRow(val) line.

RuntimeError('The requested operation is invalid on a closed state [V.TABLE]',), <traceback object at 0x0DB535D0>)

What am I doing wrong?


texts: Confidence intervals on raster produced by kriging spatial interpolation?

texts:
I have used spatial interpolation to generate a raster that estimates percentage of a variable of interest (e.g. expected percentage at unsampled locations given frequencies at nearby sites). As I am ultimately interested in identifying regions where this percentage is high, I have reclassified the raster so that values below a cutoff (e.g. 90%) are set to 0.

Initially I used inverse distance weighting spatial interpolation in ArcGIS but I'm wondering whether kriging would allow me to estimate some sort of confidence interval/estimates of error of the location of the boundaries where estimates of my variable drops to less than 90%? Does ArcGIS provide output that can be used for such purposes?


texts:
I am using QGIS version 2.14.1 (Essen) and trying to update data in my table using field calculator.

I am inputing the following statement:

update "table_1"
set  "table_1" ='y'
where  "table_2" ='test'

But I get the following error message:

syntax error, unexpected COLUMN_REF, expecting $end

Anyone know where i'm going wrong?


texts:
I have converted a shapefile into a geodatabase which I use for geocoding via the ILocator interface.  It has the extension (.mdb) and contains address fields with a point which is stored as a long binary object (I think).  I have queried it directory using libraries like OleDB. I'm wondering if there's a way to do a SELECT * FROM mytable using the ArcEngine interfaces.

I'll provide some context as well.  I'm trying to get around the inflexibility of the SoundEx implementation and some limitations with the locator styles. The locator style interace works like this for address search results:

IField addressField = candidateFields.get_Field(fieldIndex);

But the fields returned only match the fields in the locator, so I get "Zone" and not city, state, zip, unit number, owner name, parcel price, etc.


texts: Displaying HTML widget in the map generated by QGIS2web plugin

texts:
I would like to render on my map the HTML widget prepared in the QGIS. The widget looks like below.

Unfortunately, the issue is not simple, because in the popup generated by QGIS 2 web all the elements come from separate files: qgis2web.css qgis2web.js, separate layer file and geojson file.

The standard popup provided by QGIS2web plugin looks very simple. I would like to have the header in different color and table like in my QGIS HTML widget. Is it possible?
UPDATE:
In the qgis2web.js file I tried to embed this whole HTML code inside the following section:
           if (doPopup) {
                popupText += '

';
       /*.....
           popupText += '
';
            }

but it didn't work


texts:
I need to get the outer polygon of the LSOAs within 10 miles of each postcode centroid. There are 180k LSOA polygons in table1 and 2.6 million postcode centroid points in table2.

I have created 2 geography columns (polygon.geog, pois.geog) and one geometry column polygon.coordinates. All of them are indexed using GIST.

I tried running for 10 using the code below which took 1 minute.

SELECT pois.objectid, ST_Multi(ST_Union(polygon.coordinates)) as coords
              into table3
              FROM table1 as polygon
              INNER JOIN table2 as pois
              ON ST_DWithin(polygon.geog,pois.geog, 10* 1609.34,false)
              group by pois.objectid
            limit 10;

Query plan,

[
  {
    "Plan": {
      "Startup Cost": 0.84,
      "Plans": [
        {
          "Partial Mode": "Simple",
          "Startup Cost": 0.84,
          "Plans": [
            {
              "Startup Cost": 0.84,
              "Plans": [
                {
                  "Startup Cost": 0.43,
                  "Scan Direction": "Forward",
                  "Plan Width": 39,
                  "Node Type": "Index Scan",
                  "Plan Rows": 3011500,
                  "Relation Name": "table2",
                  "Alias": "pois",
                  "Parallel Aware": false,
                  "Parent Relationship": "Outer",
                  "Total Cost": 11352214.75,
                  "Index Name": "postcode_pkey"
                },
                {
                  "Filter": "(((pois.geog)::geography && _st_expand(geog, '16093.4'::double precision)) AND _st_dwithin(geog, (pois.geog)::geography, '16093.4'::double precision, false))",
                  "Startup Cost": 0.41,
                  "Scan Direction": "NoMovement",
                  "Plan Width": 11560,
                  "Node Type": "Index Scan",
                  "Index Cond": "(geog && _st_expand((pois.geog)::geography, '16093.4'::double precision))",
                  "Plan Rows": 1,
                  "Relation Name": "table1",
                  "Alias": "polygon",
                  "Parallel Aware": false,
                  "Parent Relationship": "Inner",
                  "Total Cost": 78.31,
                  "Index Name": "gidgeo"
                }
              ],
              "Node Type": "Nested Loop",
              "Plan Rows": 181100,
              "Join Type": "Inner",
              "Parallel Aware": false,
              "Parent Relationship": "Outer",
              "Plan Width": 5783,
              "Total Cost": 247219664.75
            }
          ],
          "Node Type": "Aggregate",
          "Strategy": "Sorted",
          "Plan Rows": 181100,
          "Parallel Aware": false,
          "Group Key": [
            "pois.objectid"
          ],
          "Parent Relationship": "Outer",
          "Plan Width": 39,
          "Total Cost": 247223286.75
        }
      ],
      "Node Type": "Limit",
      "Plan Rows": 1,
      "Parallel Aware": false,
      "Plan Width": 39,
      "Total Cost": 1365.96
    }
  }
]

Is it possible to improve the speed of this code?


texts: GDAL_Translate corrupts large bitmap image

texts:
I am trying to translate a large bitmap (60000 x 25000) with GDAL_Translate but it appears to corrupt the image. This works fine for smaller images but has strange results with larger ones.

I can view the image in WildBit viewer prior to translating and it has no corruption. After the translation the top of the image appears to be a repeat of the bottom of the image. Is there a way to prevent this?

I have tried the -co BIGTIFF=YES option and that has no effect.

My arguements:

gdal_translate.exe -of GTiff -a_nodata 0 --config GDAL_DATA "*data path*" -a_ullr tlX tlY brX brY -a_srs EPSG:3857 InputPath OutputPath

The end result

GDALInfo of input image

Size is 63264, 24960
Coordinate System is `'
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,24960.0)
Upper Right (63264.0,    0.0)
Lower Right (63264.0,24960.0)
Center      (31632.0,12480.0)
Band 1 Block=63264x1 Type=Byte, ColorInterp=Red
Band 2 Block=63264x1 Type=Byte, ColorInterp=Green
Band 3 Block=63264x1 Type=Byte, ColorInterp=Blue

GDALInfo of output image

Size is 63264, 24960
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],
    AUTHORITY["EPSG","3857"]]
Origin = (-8500416.405268350600000,4435398.925636629600000)
Pixel Size = (0.149300716298692,-0.149018012050865)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-8500416.405, 4435398.926) ( 76d21'37.94"W, 36d58'24.22"N)
Lower Left  (-8500416.405, 4431679.436) ( 76d21'37.94"W, 36d56'48.10"N)
Upper Right (-8490971.045, 4435398.926) ( 76d16'32.49"W, 36d58'24.22"N)
Lower Right (-8490971.045, 4431679.436) ( 76d16'32.49"W, 36d56'48.10"N)
Center      (-8495693.725, 4433539.181) ( 76d19' 5.21"W, 36d57'36.17"N)
Band 1 Block=63264x1 Type=Byte, ColorInterp=Red
  NoData Value=0
Band 2 Block=63264x1 Type=Byte, ColorInterp=Green
  NoData Value=0
Band 3 Block=63264x1 Type=Byte, ColorInterp=Blue
  NoData Value=0


texts: Storing a large amount of vector layers on PostGIS

texts:
I'm developing a web application where PostGIS will be used as a vector tile server backend to serve mvt tiles. For each user of the application at least one vector layer will be created. So, potentially there are going to be thousands of layers. Since these layers will be used to serve tiles, read-queries performance is rather critical. My question is, whether I should store all geometries in a feature table with a foreign key referencing layer table, or should I have a separate table for each layer, eg. layer_1, layer_2, etc.

The second approach feels more natural and common for relational databases. However, I have arranged a little experiment, importing 5000 layers to a single table and saw that there is a significant fall in performance.

What is the general approach to this problem? Searched the web but couldn't find the answer.


texts: How to apply Log Transform to Flow Accumulation?

texts:
I have run a flow accumulation for a DEM that I'm working with. I'm trying to apply a log transform because I've heard it could give a better visual representation.

How would I go about doing this process.... Can this be done without Python?


texts: Get geometry of Line Callout using ArcObjects?

texts:
How can I extract the geometry of a line callout and the border of an annotation featureclass in ArcGIS desktop. I want to transfert those to a separate featureclass.

I'm able to get the Anchor point of the callout, but it's not enough for my needs.

I have try to use following code, but it doesn't return nothing. In the code, I try to get a polyline, but I've also try to get a polygon with the same results. I'm using ArcMap version 9.3 but I also tried with 10 and I have same results. Quite strange also that the method QueryBoundary on IQueryGeometry says at run-time that it's not implemented!

Dim pActiveView As IActiveView = pMap
Dim pScreenDisplay As IScreenDisplay = pActiveView.ScreenDisplay
Dim pDisplayTransform As IDisplayTransformation = pScreenDisplay.DisplayTransformation

Dim pFDOGraphicsLayer As IFDOGraphicsLayer = pMap.Layer(0)   'The Annotation featureclass
Dim pFDOGraphicsLayerRead As IFDOGraphicsLayerRead = pFDOGraphicsLayer

pFDOGraphicsLayerRead.StartGeneratingGraphics(Nothing, pScreenDisplay, True, True, False)

Dim pTextElement As ITextElement = pFDOGraphicsLayerRead.NextGraphic

Do Until pTextElement Is Nothing
      Dim pFormattedTextSymbol As IFormattedTextSymbol = pTextElement.Symbol
      Dim pTextBackground As ITextBackground = pFormattedTextSymbol.Background

      If TypeOf pTextBackground Is ILineCallout And pTextBackground IsNot Nothing Then
            Dim pLineCallout As ILineCallout = pTextBackground
            Dim pQueryGeometry As ESRI.ArcGIS.Display.IQueryGeometry = pLineCallout
            Dim pElement As IElement = pTextElement
            Dim pPointGeometry As IGeometry = pElement.Geometry
            pScreenDisplay.StartDrawing(pScreenDisplay.WindowDC, ScreenDisplay.ActiveCache)
            Dim pPolyline As IPolyline = New Polyline

            pPolyline = pQueryGeometry.GetGeometry(pScreenDisplay.WindowDC, DisplayTransform, pPointGeometry)

            If pPolyline IsNot Nothing Then
               'At this point the polyline is always nothing***
               MsgBox(pPolyline.Length)
            End If

            pScreenDisplay.FinishDrawing()
        End If
    pTextElement = pFDOGraphicsLayerRead.NextGraphic
Loop

Source: http://forums.esri.com/thread.asp?c=93&f=993&t=275411


texts: ST_Centroid in view produces geometry without further definition from geometry(polygon, ssid)

texts:
I'm a Postgresql novice.

In a postgresql 9.2 + PostGIS instance I have one table with a geometry(polygon, ssid) column. From this table I am trying to create a view using ST_Centroid  like so:

CREATE OR REPLACE VIEW schema.polys_points AS
   SELECT poly_table.name,
          ST_Centroid(geom) as geom
   FROM poly_table;

This produces a view with a geometry() column where I would like it defined as geometry(point, ssid) column. Is this possible?


texts:
I am evaluating Hurricane data that has a radius which 64kt wind speeds were observed for a specific event. The properties of a Hurricane dictate that the wind speeds change in each part of said hurricane. This data has specified where the maximum radius from the eye of the storm the 64kt winds are found in each quadrant: NE, NW, SE, & SW.

I have created four buffers for showing the different sizes (i.e. one buffer with the NE quadrant's data but it is buffering all sides. However, I am trying to figure out if there is a way to cut the buffers into the needed cardinal quadrants?

The goal is to then dissolve the four buffers for a point into one buffer that shows how the radius changes for each quadrant. I know you can buffer right left sides of lines but this is a little more complicated.

(I am using 3.6 and have access to 3.4)


texts:
What is the most direct and straightforward, hopefully interactive, method to conduct a For each item Do... loop, where the input list is a selection of polygons?

A word about "interactive". I'm not opposed to a script solution, but it's not the ideal. The basic idea is that the whole process be as dynamic and immediate as possible. Point and click at some stuff, do something with it, point and click at the results, do something else, and so on for an hour or so. Then save the final results and forget about everything in the middle, it was all adhoc and doesn't need to be repeated exactly.

Software at hand is ArcGIS 10 (with 3D and Spatial Analyst), Arcinfo Workstation, and QGIS 1.7 (though don't feel like a solution from a different product is off-topic). In pseudo code the idea is:

for each selected poly in "index_layer" do:
   clip "raster_layer" by poly geometry to "d:\out\raster_poly[$ID].tif"

In the following image the selected set is outlined in heavy black and with halo text.


texts: How to use Arcpy to set Extent to "Union of Inputs"

texts:
When I try to append two polygon feature classes, I get a geometry domain error. When using just ArcMap, I know that setting the processing extent to union of inputs fixes the issue.

However, how do I do that using arcpy, I know you can set a specific extent, but am not sure how to turn on the Union of Inputs option found under processing extent in ArcMap itself using only python.

Using Python 2.7.5 and ArcGIS 10.2


texts: Identifying which singlepart non-overlapping polygons are islands using ArcGIS for Desktop?

texts: Is there a simple method to select singlepart non-overlapping polygons in a single feature-class which are not contiguous with other polygons ie they are island polygons?


texts: gdal.jar and .dll placement with gdal.Allregister();

texts:
I built gdal and java bindings for gdal.
My environment variables I added (Windows 7 64 bit):

CLASSPATH = C:\Program Files\Java\JDK1.8.0_40\lib; C:\gdal-1.11.1\swig\java\gdal.jar

the dlls (gdal111.dll (I run version 1.11.1), gdaljni.dll, ogrjni.dll and gdalconsthni.dll) are in the same folder as gdal.jar

GDAL_DATA = C:\gdal-1.11.1\data

That is the file where compdcs.csv etc. are.

JAVA_HOME = C:\Program Files\Java\jdk1.8.0_40;

Path = C:\gdal-1.11.1\bin; C:\Program Files\Java\jdk1.8.0_40\bin;C:\ProgramData\Oracle\Java\javapath;

In C:\gdal-1.11.1\bin is gdal_contour.exe, gdal111.dll etc.

Anything suspicious so far?
When I use gdal.Allregister(); I get the

java.lang.UnsatisfiedLinkError: no gdaljni in java.library.path
Exception in thread "main" java.lang.UnsatisfiedLinkError: org.gdal.gdal.gdalJNI.AllRegister()V
      ...

and now comes the interesting part: When I copy gdal.jar, gdal111.dll (I run version 1.11.1), gdaljni.dll, ogrjni.dll and gdalconsthni.dll in to the same folder as MyClass.class and run java MyClass it works.

And in Eclipse I have to put them into the bin file.. What is going on?

Oh and I should mention, that I have Windows 7 64 bit


texts: What Geotransformation (NAD27 to WGS84) to use for Gulf of Mexico?

texts:
I'm working with datasets that are spread out across the northern part of the Gulf of Mexico, stored in NAD 27.  I've been told the polygons shown above were obtained from BOEMRE (formerly known as MMS).

I'm projecting it into Web Mercator (WGS 84), and would like to choose the most suitable geotransformation.  Arcmap presents a long list of choices, but I'm not sure which one would be the overall best fit.  Is there a map somewhere showing how the error varies for various geotransformations?

Better yet, is there an industry standard geotransformation for GoM whose parameters I could plug into a custom geotransformation?


texts: How do I bypass class extensions when reading from ArcSDE with FME 2014?

texts:
I'm trying to read Feature Class from our Enterprise Geodatabase and because we have custom class extensions registered on all our our data.

I am not able to set the Reader to bypass the class extensions.

I get this warning in the log

"Unable to read from the table/feature class 'YK_SPATIAL.Buildings'. If 'YK_SPATIAL.Buildings' contains a class extension, make sure it is registered on this computer"

Thanks


texts: How to blend multiple semi-transparent GeoTIFFs

texts: How to blend multiple overlapping semi-transparent GeoTIFFs to single semi-transparent GeoTIFF on Linux (GDAL...) ?


texts: How to control QGIS from within Anaconda's Spyder (Continuum Analytics Python distribution) in Windows

texts:
I would like to control QGIS from within Anaconda's Spyder, in particular use the QGIS core modules.

My config:

Windows 7
Anaconda 64x for Python 2.7
Spyder 2.3
Postgresql 9.4
PostGIS 2.1
QGIS 2.8 (Wien)

Unfortunately it seems (from several answer I read) that one needs to recompile the distribution prior to be able to do that.

So I came to read this documentation on the github:
https://raw.githubusercontent.com/qgis/QGIS/master/INSTALL

I don't quite see where it is that one can modify the python distribution enclosed into the install. Could anyone give me some comments or a starting point to build an Anaconda compatible distribution of QGIS?


texts: AttributeError: 'QgsVectorLayer' object has no attribute 'geometry'. How to fix?

texts:
Ok... So I am trying to add a column (Length) in attribut table that shows the length of each line of a line feature. A part of my code is this:

            geom = layer.geometry()
            len = geom.length()
            processing.runalg('qgis:fieldcalculator', layer_output2, 'Length', 0, 10, 2, True, len, layer_output2)

And the error says that QgsVectorLayer has no attribute 'geometry'...What am I missing?

QgsVectorLayer is here:

    for layer in layers:
        if layer.type() == QgsMapLayer.VectorLayer:
            self.dlg.layerCombo.addItem( layer.name(), layer )

Maybe my data has a problem...but everything else I try to do, worked.
So, maybe my code is wrong...(probably, I am not good at coding).
Any tips?


texts:
I've some problems concerning Point Distance in ArcGIS 10.5.
 Whenever I try it to open, this tool appeared "require licensed tools".