texts: ArcPy copy structure of group layers then copy feature layers to this structure

texts:
I have ArcGIS project (.mxd) with group structure like shown below with some feature layers inside (which are stored in file geodatabase).
What I need is to copy group structure from "old_layers" to "new_layers" and clone feature layers preserving group structure they were before.

Structure before transformation:

Structure after transformation (layers which are in "new_layers" group should be "new", cloned based on layers from "old_layers" group:

What I achieved so far is cloning layers giving them different suffix (arcpy.Copy_management method)
Then I created new .lyr file which is a copy of "old_feature" layers so the structure of groups is preserved, then I can add new (cloned) feature layers to the .lyr file with arcpy.mapping.AddLayerToGroup method, some loops and if conditions but the last step (removing old feature layers from .lyr - arcpy.mapping.RemoveLayer) - doesnt work, so I have old and new feature layers in the same structure.

I hope you get my point.


texts:
I have a set of data held in a backend database that I wish to render into a heatmap using esri. The data is dynamic based on users search and filters.

I'm creating the feature layer and applying the heatmap renderer but getting no luck with this, one thing I notice is that when I log the layer to the console the graphics geometry's are null whether this is the issue I'm not sure.

I've created a jsfiddle of where I am currently.

https://jsfiddle.net/sgdm9an7/


texts: Cannot create extension postgis_raster on CentOS 7.6

texts:
I am using CentOS 7.6 at work. One one server I have PostgreSQL 12 and PostGIS 3 installed from the pgdg12 repo. This is working fine for vector data, but now I would like to use it for raster data as well, but when I run

-bash-4.2[Kartdata]$ psql <database> -c "create extension postgis_raster;"

I get

  ERROR:  could not load library "/usr/pgsql-12/lib/postgis_raster-3.so": /usr/gdal30/lib/libgdal.so.26: undefined symbol: GEOSMakeValid_r

I have been checking that PostgreSQL, PostGIS, GEOS and GDAL all are the newest versions, but except that PostgreSQL and PostGIS both are from the pgdg12 repo but GEOS and GDAL are from pgdg10 (I have tried to install it from pgdg12, but CentOS is still using pgdg10) everything seems perfect.

Any other ideas where I could look?

bash-4.2[<>]# ls -l /usr/pgsql-12/lib/postgis*
lrwxrwxrwx 1 root root      30 Nov  5 08:28 /usr/pgsql-12/lib/postgis-2.5.so -> /usr/pgsql-12/lib/postgis-3.so
-rwxr-xr-x 1 root root 1067904 Sep 30 12:14 /usr/pgsql-12/lib/postgis-3.so
lrwxrwxrwx 1 root root      37 Nov  5 08:28 /usr/pgsql-12/lib/postgis_raster-2.5.so -> /usr/pgsql-12/lib/postgis_raster-3.so
-rwxr-xr-x 1 root root  862368 Sep 30 12:14 /usr/pgsql-12/lib/postgis_raster-3.so
lrwxrwxrwx 1 root root      39 Nov  5 08:28 /usr/pgsql-12/lib/postgis_topology-2.5.so -> /usr/pgsql-12/lib/postgis_topology-3.so
-rwxr-xr-x 1 root root  552168 Sep 30 12:14 /usr/pgsql-12/lib/postgis_topology-3.so

bash-4.2[<>]# ls -l /usr/gdal30/lib/
total 20664
drwxr-xr-x 2 root root        6 Nov  4 19:34 gdalplugins
lrwxrwxrwx 1 root root       17 Nov  5 08:28 libgdal.so.26 -> libgdal.so.26.0.2
-rwxr-xr-x 1 root root 21156328 Nov  4 19:35 libgdal.so.26.0.2


texts: Fill with @row_number based on group field in QGIS 3

texts:
Is there a way to fill the table with a row number so that count will restart as new group (identified by a separate column) starts? As a result, "row number" column will restart the count when the new attribute value for a group is found, and will produce row number (1,2,3 etc.) record within this group.
This post on stack overflow  seems to answer my question, but SQL in QGIS does not allow such statements (DENSE_RANK) (or it does but I don't know how?).


texts: Geopackage raster layer name is blank in QGIS Layers panel?

texts:
I added two rasters to QGIS: "site01.tif" and "site02.tif" then added them to a geopackage. When I drag the geopackage into QGIS, the geopackage layers display as expected but the layer names are blank in QGIS' Layer panel. I got the idea that they should just display automatically from watching this video: https://www.youtube.com/watch?v=_KuTfTd9RKc

I used the following to create the geopackage and add the two rasters:

C:\>gdal_translate -co APPEND_SUBDATASET=YES -co RASTER_TABLE=site01 -a_srs EPSG:2258 -of GPKG C:\Users\richard\Desktop\geopackage-temp\site01.tif C:/Users/richard/Desktop/geopackage01.gpkg
Warning 1: TIFFFetchNormalTag:Incompatible type for "RichTIFFIPTC"; tag ignored
Input file size is 26294, 40994
0...10...20...30...40...50...60...70...80...90...100 - done.

C:\>gdal_translate -co APPEND_SUBDATASET=YES -co RASTER_TABLE=site02 -a_srs EPSG:2258 -of GPKG C:\Users\richard\Desktop\geopackage-temp\site02.tif C:/Users/richard/Desktop/geopackage01.gpkg
Warning 1: TIFFFetchNormalTag:Incompatible type for "RichTIFFIPTC"; tag ignored
Input file size is 46713, 44301
0...10...20...30...40...50...60...70...80...90...100 - done.

C:\>gdalinfo C:\Users\richard\Desktop\geopackage01.gpkg
Driver: GPKG/GeoPackage
Files: C:\Users\richard\Desktop\geopackage01.gpkg
Size is 512, 512
Coordinate System is `'
Subdatasets:
  SUBDATASET_1_NAME=GPKG:C:\Users\richard\Desktop\geopackage01.gpkg:site01
  SUBDATASET_1_DESC=site01 - site01
  SUBDATASET_2_NAME=GPKG:C:\Users\richard\Desktop\geopackage01.gpkg:site02
  SUBDATASET_2_DESC=site02 - site02
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)

C:\>

I drag the geopackage into QGIS with no issues however, the layer names are blank in the side panel. In the image below, "site01" and "site02" are the original rasters. The two no-name rasters are the two geopackage layers.:

The layer names are also blank when I drag the geopackage into the project:

I'm expecting to see "site01" and "site02" in the Layers panel for both the original rasters AND the geopackage layers. Instead, the original rasters display their names but the geopackage layers are blank.I'm assuming that the layer name would parse from the RASTER_TABLE.

QGIS 3.6.3

GDAL/OGR 2.4.1

Windows 10 Pro


texts: ImportError dll load failed the specified module could not be found (rasterio)

texts:
I want to use the rasterio package to make some test with the Planet Explorer Api (https://www.planet.com/docs/guides/quickstart-ndvi/). I've seen on many forums that I needed to use a wheel file to install rasterio (and gdal), it worked for gdal but for rasterio I got this issue :

>> import rasterio
Traceback (most recent call last):
 File "<stdin>", line 1, in <module>
 File "C:\Users\Hortense\Miniconda3\envs\envpythonfs\Lib\site-packages\rasterio\__init__.py", line 23, in <module>
  from rasterio.__base import (
ImportError: DLL load failed: the module specified can't be found.

I create a virtual environnement with conda, with the version python 2.7, I used those files from this site http://www.lfd.uci.edu/%7Egohlke/pythonlibs/

GDAL‑2.1.3‑cp27‑cp27m‑win_amd64.whl

rasterio‑1.0a9‑cp27‑cp27m‑win_amd64.whl

I also download the Microsoft Visual C++ 2008, but it doesn't work either.
I don't know how to solve this, can someone help me ?


texts: Points not displaying in feature class in file geodatabase

texts:
I was working with shape files when I recently decided to put all my shape files into a single file geodatabase.  I just tried to add some more points to the feature class by inputting the coordinates in the columns I have created for the x y coordinates.  However, the point is not displaying.  I am thinking that the feature class is just not recognizing the xy columns as xy attributes since this is the first time I have tried to add points to the table.  Originally to add the points to the map I just imported an excel sheet with the xyz and attribute data as a table then exported the table as a shapefile.

Is there a way to specify columns as coordinates so when adding points in the attribute table they will display without having to construct a point?


texts:
I want to aggregate points of interest (POI) using FME. The aggregation should make it possible to group all POI within an area based on a buffer distance from the centerpoint of n POI.

This would mean an iterative approach where all POI are subject to a buffer-select analysis counting the number of POI within the bufffer distance. The most "clusteder" POIs will then be used to represent a combination related to the POI types.

If possible the attributes from all POI should be merged so that all atttributes from all aggregated POI are kept as separate attributes in the representative POI.

The procedure is straightforward, but can this be done using FME?

This question is related to my question Group and align icons in QGIS atlas.


texts:
I have NED elevation data in 3 rasters. (That's how it downloaded from USGS NED page.) The files are named 61471006, 61471006_2, and 61471006_3. I succeeded in using the GRASS GUI and menu [File-> Import raster data -> Common import formats] to import the first raster, by browsing into the file folder named 61471006 and selecting the file "w001001.adf". However, I need to also bring in the .adf files from 61471006_2 and 61471006_3, but within their respective folders, the .adf files are named the same as the .adf file in the first folder. Therefore GRASS does not let me bring in the same-name files.
I tried re-naming the .adf files in 61471006_2 and 61471006_3, but then GRASS put up the error message that the file type was not recognized.
How can I get these ostensibly same-name rasters to all three be imported into my GRASS map? (Each raster covers a different geographic area within the overall area.)


texts:
I have a line layer in PostGIS that I would like to add to a Leaflet map. I have the PHP setup, but I can't seem to get the Javascript correct. I know my issue is converting the line to geoJSON.

The following is the JavaScript code I have setup:

function refreshAttractions(){
    $.ajax({
        url:'load_attractions.php',
        type:'POST',
        data:{filter:$("#filter").val()},
        success:function(response){

            lyrAttractions=L.geoJSON(JSON.parse(response), {
                LineString: function(latlngs) {
                    var strPopup = "
"+feature.properties.bikeway_ty+"

";
                    strPopup += "
Sponsor: "+feature.properties.sponsor+"
";
                    strPopup += "
Location: "+feature.properties.location+"
";
                    strPopup += "
From: "+feature.properties.from_limit+"
";
                    strPopup += "
To: "+feature.properties.to_limit+"
";
                    strPopup += "
Name: "+feature.properties.local_name+"
";
                    strPopup += "
Length: "+feature.properties.length_mil+"
";
                    //strPopup += "";
                    return L.polyline(latlngs).bindPopup(strPopup);
                }
            });
            lyrAttractions.addTo(mymap);
    }});
}

And the PHP code that fetches data from PostGIS and returns JSON is:

<?php
    $strQry="SELECT sponsor, location, from_limit, to_limit, local_name, length_mil, bikeway_ty, ST_AsGeoJSON(geom, 5) as geom FROM bikeways_existing";

    $db=new PDO('pgsql:host=*****.com;port=5432;dbname=geozeoco_demo;','geozeoco','******');
    $sql = $db->query($strQry);
    $features=[];
    while ($row = $sql->fetch(PDO::FETCH_ASSOC)) {
        $feature=['type'=>'Feature'];
        $feature['geometry']=json_decode($row['geom']);
        unset($row['geom']);
        $feature['properties']=$row;
        array_push($features, $feature);
    }
    $featureCollection=['type'=>'FeatureCollection', 'features'=>$features];
    echo json_encode($featureCollection);
?>


texts: How do I perform one-step transformation in python?

texts:
I would like to perform transformation for this example data set.
I have here four pairs of adjustment points and want to transform given point coordinates from primary system to secondary system according to adjustment.

primary_system1 = (3531820.440, 1174966.736, 5162268.086)
primary_system2 = (3531746.800, 1175275.159, 5162241.325)
primary_system3 = (3532510.182, 1174373.785, 5161954.920)
primary_system4 = (3532495.968, 1175507.195, 5161685.049)

secondary_system1 = (6089665.610, 3591595.470, 148.810)
secondary_system2 = (6089633.900, 3591912.090, 143.120)
secondary_system3 = (6089088.170, 3590826.470, 166.350)
secondary_system4 = (6088672.490, 3591914.630, 147.440)

#transform this point
x = 3532412.323
y = 1175511.432
z = 5161677.111<br>

eventually, how do I compute helmert transformation parameters, shifts, rotations and scale factor?
EDIT
at the moment I try to average translation for x, y and z axis using each of the four pairs of points like:

#x axis
xt1 =  secondary_system1[0] - primary_system1[0]
xt2 =  secondary_system2[0] - primary_system2[0]
xt3 =  secondary_system3[0] - primary_system3[0]
xt4 =  secondary_system4[0] - primary_system4[0]

xt = (xt1+xt2+xt3+xt4)/4    #averaging

...and so on for y and z axis

#y axis
yt1 =  secondary_system1[1] - primary_system1[1]
yt2 =  secondary_system2[1] - primary_system2[1]
yt3 =  secondary_system3[1] - primary_system3[1]
yt4 =  secondary_system4[1] - primary_system4[1]

yt = (yt1+yt2+yt3+yt4)/4    #averaging

#z axis
zt1 =  secondary_system1[2] - primary_system1[2]
zt2 =  secondary_system2[2] - primary_system2[2]
zt3 =  secondary_system3[2] - primary_system3[2]
zt4 =  secondary_system4[2] - primary_system4[2]

zt = (zt1+zt2+zt3+zt4)/4    #averaging

Is it right way?


texts: Why does an insert cursor only insert rows when it is followed by del statement even though the del statement is outside the loop?

texts:
My insert cursor failed in insert any rows.  After spending time trying to fix something within the main loops, I added the del inCursor after the loop and now the cursor inserts rows.  How can something outside and after a loop impact something above?

#------------THIS DIDN'T WORK:
    inCursor = arcpy.da.InsertCursor(fcEmpty,["Shape@","NAME","CO_WI", "CO_NRI","A_NUMBER","LSE_NO","ORIG_LENGTH"])
    with arcpy.da.SearchCursor(fcPolyline,["Shape@","NAME","ORIG_LENGTH"]) as firstCursor:
        for firstRow in firstCursor: # For each feature in input1
            inputFeat = firstRow[0]
            with arcpy.da.SearchCursor(testingINLAND, ["Shape@","CO_WI", "CO_NRI","A_NUMBER","LSE_NO"]) as intersectCursor:
                for intersectRow in intersectCursor: # for each feature from input2
                    # intersect it with the feature from input1
                    intersectFeat = intersectRow[0]
                    print "check"
                    # dimension ==> polygon geometry (where 1 = point, 2 = polyline)
                    clippedFeature = intersectFeat.intersect(inputFeat, 2)
                    print clippedFeature
                    # setup the insert cursor (for the input you previously created)
                    inCursor.insertRow([clippedFeature,firstRow[1],intersectRow[1],intersectRow[2],intersectRow[3],intersectRow[4],firstRow[2]]) # Insert the resulting geometry from INTERSECT geometry method
    del firstRow, firstCursor
    del intersectRow, intersectCursor

#------------THIS DOES WORK:

inCursor = arcpy.da.InsertCursor(fcEmpty,["Shape@","NAME","CO_WI", "CO_NRI","A_NUMBER","LSE_NO","ORIG_LENGTH"])
with arcpy.da.SearchCursor(fcPolyline,["Shape@","NAME","ORIG_LENGTH"]) as firstCursor:
    for firstRow in firstCursor: # For each feature in input1
        inputFeat = firstRow[0]
        with arcpy.da.SearchCursor(testingINLAND, ["Shape@","CO_WI", "CO_NRI","A_NUMBER","LSE_NO"]) as intersectCursor:
            for intersectRow in intersectCursor: # for each feature from input2
                # intersect it with the feature from input1
                intersectFeat = intersectRow[0]
                print "check"
                # dimension ==> polygon geometry (where 1 = point, 2 = polyline)
                clippedFeature = intersectFeat.intersect(inputFeat, 2)
                print clippedFeature
                # setup the insert cursor (for the input you previously created)
                inCursor.insertRow([clippedFeature,firstRow[1],intersectRow[1],intersectRow[2],intersectRow[3],intersectRow[4],firstRow[2]]) # Insert the resulting geometry from INTERSECT geometry method
del inCursor
del firstRow, firstCursor
del intersectRow, intersectCursor


texts: How to calculate area in sq km of an ellipsoid shapefile

texts:
In QGIS, I often find myself setting up projects with an ellipsoid coordinate system (e.g. WGS84) so that the resulting maps have lat lon scales, etc.

Often in these projects I also need to do some calculations in other units (e.g. calculate the sq km of a feature).

At present I do this my maintaining two versions of a shapefile, one projected for calculations and one ellipsoid for layouts.

This becomes complicated in large projects!

Is there an easier way? Am I overlooking a technique for doing these calculations on-the-fly from spheroid shapefiles?


texts: How to sum 10 raster bands instead of 9 from code?

texts:
I want to sum 10 raster bands but the problem is that the following code will only work for summing 9 bands, which means the output image consists of 9 bands.

from qgis.core import *
from qgis.gui import *
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
from math import *
from PyQt4.QtCore import QFileInfo
import os
from osgeo import gdal
from array import *
import ctypes

entries = []
fileName = "E:/WRA/pet/mothmean81j82m.tif"
fileInfo = QFileInfo(fileName)
baseName = fileInfo.baseName()
rlayer =QgsRasterLayer(fileName,baseName)

boh1 = QgsRasterCalculatorEntry()
boh1.ref = 'boh@1'
boh1.raster = rlayer
boh1.bandNumber = 1
entries.append( boh1 )

boh2 = QgsRasterCalculatorEntry()
boh2.ref = 'boh@2'
boh2.raster = rlayer
boh2.bandNumber = 2
entries.append( boh2 )

boh3 = QgsRasterCalculatorEntry()
boh3.ref = 'boh@3'
boh3.raster = rlayer
boh3.bandNumber = 3
entries.append( boh3 )

boh4 = QgsRasterCalculatorEntry()
boh4.ref = 'boh@4'
boh4.raster = rlayer
boh4.bandNumber = 4
entries.append( boh4)

boh5 = QgsRasterCalculatorEntry()
boh5.ref = 'boh@5'
boh5.raster = rlayer
boh5.bandNumber = 5
entries.append( boh5 )

boh6 = QgsRasterCalculatorEntry()
boh6.ref = 'boh@6'
boh6.raster = rlayer
boh6.bandNumber = 6
entries.append( boh6 )

boh7 = QgsRasterCalculatorEntry()
boh7.ref = 'boh@7'
boh7.raster = rlayer
boh7.bandNumber = 7
entries.append( boh7 )

boh8 = QgsRasterCalculatorEntry()
boh8.ref = 'boh@8'
boh8.raster = rlayer
boh8.bandNumber = 8
entries.append( boh8 )

boh9 = QgsRasterCalculatorEntry()
boh9.ref = 'boh@9'
boh9.raster = rlayer
boh9.bandNumber = 9
entries.append( boh9 )

boh10 = QgsRasterCalculatorEntry()
boh10.ref = 'boh@10'
boh10.raster = rlayer
boh10.bandNumber = 10
entries.append( boh10 )

addop='boh@1+ boh@2+ boh@3+ boh@4+ boh@5+ boh@6+ boh@7+ boh@8+ boh@9+ boh@10'
print addop
calc = QgsRasterCalculator(addop, 'E:/WRA/pet/addoop.tif', 'GTiff', rlayer.extent(), rlayer.width(), rlayer.height(), entries )
calc.processCalculation()


texts: CartoDB not showing title or description (or other options)

texts:
I'm using CartoDB and I can only get one option every to work cartodb_logo. The rest seem to be ignored. Here's the relevant part of my code. What is wrong with this to cause title and description to be ignored (not on my map)?

var layerSource = {
        user_name: 'myname',
        type: 'cartodb',

        sublayers: [{
            sql: "SELECT * FROM mydata",

            cartocss: "#mydata{marker-width: 10; marker-line-color: #FFF; marker-line-width: 0.5; marker-line-opacity: 1; marker-fill-opacity: 1.0; }"
        }]

    };
    var options = {

        https: true,
        infowindow: true,
        cartodb_logo: false, //this works fine
        title: "testing title", //this isn't rendered on map
        description: "testing description", //this isn't rendered on map
        shareable: true, //no sharing buttons shown
        legends: true //no legends either
    };
    var map_object = new L.Map('test_heat', {
        center: [39.8282, -98.5795],
        zoom: 5,
    });
    L.tileLayer('https://cartodb-basemaps-{s}.global.ssl.fastly.net/light_all/{z}/{x}/{y}@2x.png', {

        attribution: '&copy; OpenStreetMap contributors, &copy; CartoDB'
    }).addTo(map_object);

    cartodb.createLayer(map_object, layerSource, options)
        .addTo(map_object)
        .on('done', function(layer) {
            cdb.vis.Vis.addInfowindow(map_object, layer.getSubLayer(0), options, ['cartodb_id', 'number_of_cases'], {
                infowindowTemplate: $('#infowindow_template').html()
            });
            var sublayer = layer.getSubLayer(0);
            sublayer.setInteraction(true);
            sublayer.on('featureClick', function(e, latlng, pos, data) {
                //alert("Hey! You clicked " + data.cartodb_id);
            });
        })
        .done(function(data) {

        });


texts: I am trying to put title on leaflet markers. This explains the way to give title to marker which would appear on hover. Is there any way we can set time on that i.e. make the title appear after 50ms of the mouse hover?


texts: How to add feature to auxiliary layer and commit changes to auxilliary field in QGIS 3.x?

texts:
I can't figure how to do that, here is my code :
auxKeyField = myLayer.fields().field(myLayer.fields().indexOf('id'))
auxStorage = QgsProject.instance().auxiliaryStorage()
myAuxLayer = auxStorage.createAuxiliaryLayer(auxKeyField, myLayer)
myLayer.setAuxiliaryLayer(myAuxLayer)

propAux = QgsPropertyDefinition()
propAux.setName('auxfield')
propAux.setDataType(QgsPropertyDefinition.DataTypeNumeric)
res = myAuxLayer.addAuxiliaryField(propAux)
auxLayerFields = myAuxLayer.fields()

auxFeature = QgsVectorLayerUtils.createFeature(myAuxLayer)
auxFeature.setAttribute('ASPK', 1)
auxFeature.setAttribute('_auxfield', 123)
res = myAuxLayer.addFeature(auxFeature)

myAuxLayer.save()

Auxiliary field is created, addFeature returns true, the auxiliary field column appears in attribute table, but it contains only NULL values, even for id=1.
EDIT 1 : Adding myLayer.startEditing() at the beginning of the code, adding myLayer.commitChanges() at the end, and removing myAuxLayer.save() seemed to solve the problem. Unfortunately, as soon as I save the project, every value in the auxiliary field disappears...
EDIT 2 : Looking at unit tests I ended up with this code :
auxKeyField = myLayer.fields().field(myLayer.fields().indexOf('id'))
auxStorage = QgsProject.instance().auxiliaryStorage()
myAuxLayer = auxStorage.createAuxiliaryLayer(auxKeyField, myLayer)
myLayer.setAuxiliaryLayer(myAuxLayer)

propAux = QgsPropertyDefinition(
    'offset',  # name
    'string',  # description
    QgsPropertyDefinition.Offset,  # type
    'symbol',  # origin
    'decalop'  # comment
)
res = myAuxLayer.addAuxiliaryField(propAux)

req = QgsFeatureRequest().setFilterExpression('id=1')
myFeature = QgsFeature()
myLayer.getFeatures(req).nextFeature(myFeature )
auxFieldName = QgsAuxiliaryLayer.nameFromProperty(propAux, True)
index = myLayer.fields().indexOf(auxFieldName)
myLayer.changeAttributeValue(myFeature.id(), index, 123)

Sadly, every value in the auxiliary field still disappear as soon as I save the project...


texts: Is there a way to use a variable in a rule-based style?

texts:
I have several symbol layers:

current: "from" <= 1832 AND "to" > 1832
former: "from" <= 1832 AND "to" <= 1832
future: "from" > 1832

What I'd like to do is be able to change 1832 to another year without separately editing each rule. Ideally I'd be able to define a variable somewhere and refer to that variable in the rules.

I've found a workable solution by creating an extra column in the attribute table and filling it with 1832, then referring to that in the rules. It's easy to change with the field calculator, but is rather inelegant. I tried using a joined table and editing that instead, but it doesn't update the rendering properly when I save the table with a new year value.

Is there a better way to do this, or have I determined the optimal way?


texts:
In ArcGIS engine I have a polyline with From and To points. I also have a Point. How can I determine if the point is located on the right or left of the line on the plane considering the direction of the line?

I thought of using MeasurementTool to get a point from the line and the point to be related and find the angle between them in relation with the azimuth of the line, but there lies the problem of which point from the line I should choose.


texts: Getting correct WKT result to 4 decimal places using PostGIS?

texts:
I'm using PostGIS version 2.3. I want to get the SW and NE corners of a bounding box, corrected to 4 decimal places.

Right now I have this to retrieve the SW and NE corners from my geometry column:

select
st_astext(st_transform(st_makepoint(st_xmin(st_extent(geom)), st_ymin(st_extent(geom))), '+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs ',4326)) as sw,
st_astext(st_transform(st_makepoint(st_xmax(st_extent(geom)), st_ymax(st_extent(geom))), '+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs ',4326)) as ne
FROM HOLDINGS
WHERE custname='XXX';

Which returns results like this:

                     sw                     |                     ne
--------------------------------------------+--------------------------------------------
 POINT(-0.685239554498251 51.4940418030595) | POINT(-0.680132974825013 51.4962802022361)

But how can I amend this to get the SW and NE corners corrected to 4 decimal places?

The results can be four numbers rather than two WKT POINTs (in fact that would be preferable).


texts: Opening map (*.mxd) using ArcGIS Online?

texts:
I created a map yesterday using ArcMap on my laptop.

I do not have access to the programme right now.

Is it possible to open this map (mxd) file on ArcGIS online?


texts: Local or global Moran's I minimum neighbor connectivity rules of thumb and justification

texts:
I am trying to run a univariate local Moran's I (in GeoDa) and I'm questioning my spatial weights, trying to do different neighbor counts / distances in increments and evaluate my results.

I was told that all of my locations need a minimum of eight (8) neighbors. Using the range of spatial weights I think is most appropriate (based on familiarity with the dataset and the expectation of what I'm looking for) I have between twelve (12) and twenty (20) records with <8 neighbors out of over fifteen-thousand (15,391) total records. These now-neighbor records are larger polygons near the edges of my map.

I tried understanding 'why' but I can't find much online. Only that ESRI lists it as a requirement and several other places reference it. I'm assuming it has something to do with estimating the local covariance, but I'm not really sure.

I want to know how my analysis will be impacted if it violates this rule. What is the concern with having a small or large number of neighbor counts below eight (8) but above zero (0).

What happens if I continue with these records included and my spatial weights as they are?


texts: I need to define the boundary of a number of towns. I have found maps of some of the towns in google maps but have no idea how I can digitize or export them for use in ArcMap. So far i have not seen any tools in Google Maps that could be used.


texts:
I have a problem with my newly created layers. I have drawn a lot of features and most of them disappeared during the weekend while the pc was turned off. They are not visible in the drawing, in the attribute table and in the ArcCatalog. The problem is caused not by wrong saving because it happened twice. At the first time, after problem occured, while drawing the features once again, the old ones started to appear on the drawing and attribute table. Now, I tried to review the layer with another pc, I tried to export it as a copy and I still can't see my features.

Does anyone know what could be the problem of this and what is the solution?


texts: GDAL/OGR Returning wrong reprojected shapefile

texts:
I use the script in this link to reproject a shapefile and save it. The target projection is derived from a georeferenced raster. The problem is while the target coordinate system is UTM 33N (epsg: 32633) the output shapefile has some other coordinate system (I expected to get Central Europe and Africa and I get Australia and New Zealand). I did notice that the output shapefile didn't have a .prj file, just the .shp, .shx and .dbf files. What could be the problem?
Here's the script I use:

def reproject_shapefile(coast, epsg, prj):
import os

"""epsg: target srs
prj: target projection"""
driver = ogr.GetDriverByName('ESRI Shapefile')

# get the input layer
inDataSet = driver.Open(coast)
inLayer = inDataSet.GetLayer()

InSpatialReference = inLayer.GetSpatialRef()
epsg_src = InSpatialReference.GetAttrValue('AUTHORITY', 1)
if not epsg == epsg_src:
    # inSpatialRef = osr.SpatialReference()

    # output SpatialReference
    outSpatialRef = osr.SpatialReference(wkt=prj)

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

    # create the output layer
    outputShapefile = coast.replace('.shp', '_%s.shp' % epsg)
    if os.path.exists(outputShapefile):
        driver.DeleteDataSource(outputShapefile)
    outDataSet = driver.CreateDataSource(outputShapefile)
    outLayer = outDataSet.CreateLayer('coast_%s' % epsg, geom_type=ogr.wkbMultiPolygon)

    # add fields
    inLayerDefn = inLayer.GetLayerDefn()
    for i in range(0, inLayerDefn.GetFieldCount()):
        fieldDefn = inLayerDefn.GetFieldDefn(i)
        outLayer.CreateField(fieldDefn)

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

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

    # Save and close the shapefiles
    inDataSet = None
    outDataSet = None

    coast = outputShapefile

return coast


texts:
I use geoserver and OL3 to display one layer. I have configured my sld in order to show labels after a specific zoom level. It works quite well but there is one specific polygon that doesn't show the label at a specific zoom level (see attachment - upper left). All the others seem to work fine.

It doesn't seem to be a problem of overlapping or tiled layers. Any ideas?

This is how my SLD rule looks like:

<sld:Rule>
      <sld:Name>rule1</sld:Name>
      <sld:Title>Grid</sld:Title>
      <sld:MaxScaleDenominator>1701630</sld:MaxScaleDenominator>
      <sld:PolygonSymbolizer>
        <sld:Fill>
          <sld:CssParameter name="fill">#1B1818</sld:CssParameter>
          <sld:CssParameter name="fill-opacity">0.4</sld:CssParameter>
        </sld:Fill>
        <sld:Stroke>
          <sld:CssParameter name="stroke">#0000FF</sld:CssParameter>
          <sld:CssParameter name="stroke-opacity">0.9</sld:CssParameter>
          <sld:CssParameter name="stroke-width">1.2</sld:CssParameter>
        </sld:Stroke>
      </sld:PolygonSymbolizer>
      <sld:TextSymbolizer>
        <Geometry>
          <ogc:Function name="centroid">
            <ogc:PropertyName>shape</ogc:PropertyName>
          </ogc:Function>
        </Geometry>
        <sld:Label>
          <ogc:PropertyName>map_code</ogc:PropertyName>
        </sld:Label>
        <sld:Font>
          <sld:CssParameter name="font-family">Arial</sld:CssParameter>
          <sld:CssParameter name="font-size">12</sld:CssParameter>
          <sld:CssParameter name="font-style">normal</sld:CssParameter>
          <sld:CssParameter name="font-weight">normal</sld:CssParameter>
          </sld:Font>
        <sld:LabelPlacement>
          <sld:PointPlacement>
            <sld:AnchorPoint>
              <sld:AnchorPointX>0.5</sld:AnchorPointX>
              <sld:AnchorPointY>0.5</sld:AnchorPointY>
            </sld:AnchorPoint>
          </sld:PointPlacement>
        </sld:LabelPlacement>
        <sld:Halo>
          <sld:Radius>0.5</sld:Radius>
            <sld:Fill>
              <sld:CssParameter name="fill">#FFFFFF</sld:CssParameter>
              <sld:CssParameter name="fill-opacity">0.78</sld:CssParameter>
            </sld:Fill>
        </sld:Halo>
        <sld:Fill>
          <sld:CssParameter name="fill">#0616F6</sld:CssParameter>
        </sld:Fill>
        <sld:VendorOption name="goodnessOfFit">0.3</sld:VendorOption>
        <sld:VendorOption name="spaceAround">-10</sld:VendorOption>
        <sld:VendorOption name="conflictResolution">false</sld:VendorOption>
        <sld:VendorOption name="partials">true</sld:VendorOption>
        <sld:VendorOption name="charSpacing">23</sld:VendorOption>
      </sld:TextSymbolizer>
    </sld:Rule>


texts: Get error when I try to display a GeoJson file with Openlayers 3

texts:
I get some errors when I try to add a Geojson file with Opelayers 3.

I did follow the first answer here : Add GeoJSON layer to OpenLayers 3

But I got the following error :

XMLHttpRequest cannot load file:///C:/wamp64/www/Web%20Mapping/fond.geojson. Cross origin requests are only supported for protocol schemes: http, data, chrome, chrome-extension, https.

Uncaught SyntaxError: Unexpected end of JSON input
    at parse (<anonymous>)
    at Fr (ol.js:494)
    at Nr.l.ja (ol.js:494)
    at up.<anonymous> (ol.js:331)
    at ad (ol.js:43)
    at io.l.dispatchEvent (ol.js:41)
    at so (ol.js:322)
    at io.l.en (ol.js:321)
    at io.l.Wg (ol.js:321)

this is my code :

    <link rel="stylesheet" href="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.4/css/bootstrap.min.css">

    <link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/ol3/3.6.0/ol.css" type="text/css">

I tried changing the url to :

http://localhost/Web%20Mapping/fond.geojson

And I get another error :

XMLHttpRequest cannot load http://localhost/Web%20Mapping/fond.geojson. No 'Access-Control-Allow-Origin' header is present on the requested resource. Origin 'null' is therefore not allowed access.


texts:
I have a vectortile layer of points I added to PostGIS using ogr2ogr, which runs on a remote server. The points show nicely on my client in the browser. Now I would like to add a point from my clients' browser and store it in PostGIS. What tools do I need?

I find queries like:

UPDATE sysinst.bio SET geom = ST_Transform(ST_SetSRID(ST_MakePoint(longitude, latitude), 4326),26913);

which I understand. But how do I exactly use this query (from my browser)? My stack is node.js (express.js) based.


texts: MapServer 7 and File Geodatabases and Expressions

texts:
I am trying to define within my map definition document test1.map a layer that comes from a ESRI file geodatabase (from ArcMap 10.2), and display different styles/fill colors for each item based upon a field within the layer, in this case ws_condition.  I have defined it like this:

LAYER
    NAME "analysis_sp__test1__a1__demo2__ws_condition"
    TYPE POLYGON
    STATUS ON
    PROJECTION
        "init=epsg:32049"
    END
    METADATA
        "wms_title" "ws_conditionLayer"
        "wfs_title" "ws_conditionLayer"
        "wfs_srs"    "EPSG:32049"
        "gml_include_items" "all"
        "gml_featureid"  "objectid"
        "wfs_enable_request" "*"
    END
    CONNECTIONTYPE OGR
    CONNECTION "C:\OSGeo4W\apps\demo1\data\demo1.gdb"
    DATA "demo2__a1__test1__merged"
    CLASS
        EXPRESSION('[ws_condition]' < -.6)
        STYLE
            COLOR 255 0 0
            OUTLINECOLOR 0 0 0
        END

    END
    CLASS
        EXPRESSION('[ws_condition]' > -.6 AND '[ws_condition]' < -.2)
        STYLE
            COLOR 240 128 128
            OUTLINECOLOR 0 0 0
        END

    END
    CLASS
        EXPRESSION('[ws_condition]' > -.2 AND '[ws_condition]' < .2)
        STYLE
            COLOR 229 204 255
            OUTLINECOLOR 0 0 0
        END

    END
    CLASS
        EXPRESSION('[ws_condition]' > .2 AND '[ws_condition]' < .6)
        STYLE
            COLOR 142 229 238
            OUTLINECOLOR 0 0 0
        END

    END
    CLASS
        EXPRESSION('[ws_condition]' > .6)
        STYLE
            COLOR 0 0 139
            OUTLINECOLOR 0 0 0
        END

    END

END

Now I receive nothing on the map, and worse, when I try to open it inside of the MapServer Manager, I receive an error saying there is a parsing error.   Is there any way to make this work within the map definition file, or possibly how to create a simple WFS server layer in OpenLayers (keeping it simple, the WFS server will be on the same server and domain as the OpenLayers server).   I have tried many different things, including trying to use CLASSITEM and specifying the field name, but it seems to only be doing a string comparison and not a numeric comparison.


texts: Creating a geoserver layer with curl (but with a bounding box)

texts:
I'm attemping to load a layer into Geoserver from a PostGIS database using cURL.
It works, but I also need to set the bounding box to -180,-90,180,90.  Does anyone know how to add the bounding box to the below command?

curl -v -u xxx:xxx -X POST -H "Content-type: text/xml" -d "<featureType><name>Airports</name></featureType>" http://xxx.xx.com:8080/geoserver/rest/workspaces/PostGIS/datastores/PostGIS/featuretypes

Thanks for your help!


texts: Extracting coordinates of pixels inside polygon in Google Earth Engine as single feature


texts: arcpy Tool works in debug mode but not in normal model

texts:
When I try to run the following simple script from arcmap it fails but when run using right--click debug and then F9 (debug) in python scripter it works fine. I don't understand why???

# Modified by: Georgec Corea (ATGIS) coreagc@gmail.com

# Author:  ESRI
# Date:    July 5, 2010
# Version: ArcGIS 10.0
# Purpose: This script will print one or more map documents to a local printer.
#          The script is intended to run within a script tool.  There are two
#          parameters:
#               1) Select Map Documents to Print,
#               2) Select Output Printer (auto populated using a validation script)
#
#Notes: The print order of the MXDs is based on how they are entered.  The MXD
#       at the top of the list is first followed by those below it.

import arcpy, string, datetime, shutil, os
import arcpy.mapping as MAP

#Read input parameters from script tool
MXDList = string.split(arcpy.GetParameterAsText(0), ";")
#printer = arcpy.GetParameterAsText(1)

count=0
#Loop through each MXD and print
for MXDPath in MXDList:
    count=count+1
    arcpy.AddMessage(str(count)+'. Working on '+str(MXDPath)+' @ '+str(datetime.datetime.now())+"\n")
##    try:
    mxd=arcpy.mapping.MapDocument(MXDPath)
    outPDFpath = mxd.filePath[:mxd.filePath.rfind('.')]+'.pdf'
    ddp = mxd.dataDrivenPages
    #MXD = MAP.MapDocument(mxd)
        #arcpy.AddMessage(ddp, MXD)
    if str(os.path.isfile(outPDFpath))=='False':
            outputPDF = arcpy.mapping.PDFDocumentCreate(outPDFpath)
            mxd.dataDrivenPages.exportToPDF(outPDFpath)
        #arcpy.mapping.ExportToPDF(mxd, outPDFpath)
            name='c:\\temp\\'+outPDFpath[outPDFpath.rfind('\\'):]
            shutil.copy2(outPDFpath,name)
            arcpy.AddMessage('Printed '+str(outPDFpath)+' @ '+str(datetime.datetime.now())+"\n")
    else:
            arcpy.AddMessage('File Exists '+str(outPDFpath)+' skipped pdf creation @ '+str(datetime.datetime.now())+"\n")
##    except:
##        errorm=arcpy.GetMessages()
##        arcpy.AddMessage(' Error creating file...' +str(outPDFpath)+'\n'+str(errorm)+'\n'+' continuing...')

    #MAP.PrintMap(MXD, printer)
    #Remove variable reference to file
del mxd

When run in debug


texts:
I have an OpenLayers cluster map. I want to create a customized popup as follows:

The problem I'm facing is auto-pan doesn't work
var popup = new ol.Overlay(//** @type {olx.OverlayOptions} *//*
({
    element : container,
    autoPan : true,
    autoPanAnimation : {
        duration : 250
    }
}));

i.e if I click on a cluster on the border of the visible view the screen doesn't adjust to fit the popup rather the popup gets cut.

Second issue I'm facing is the popup width, I want to give the popup a responsive width
i.e. it should have a width in relation to its content but it gets scrunched up. I've tried giving it in % but only hard coded px width works.

I want to give it responsive width, with a min width beyond which it should not go

Third issue I'm facing is if my popup has a lot of data I'm showing a scroller by using:
.popover-content {overflow-y: scroll;}

How do I customize the scroll bar??

cluster map js file: http://pastebin.com/XamgaYHC
HTML file: http://pastebin.com/M0ncX3hV
css file: http://pastebin.com/QnjnUv1N


texts:
I have a project where I am making a lot of rule based labels from a polygon layer.  The layer is saved as a shapefile (like the other label-layers (point and linestrings). All other layers in the project is saved in a database (PostgreSQL).

I have made the layer with no symbology and a lot of rule based labels, where only a small part is active.

When I turn the layer on, it overrides all other layers with a black color. The more I zoom out the more of the data frame it covers in black. It looks like this:

And when the layer is of it looks fine:

Does anyone know why I have this problem on a polygon layer, but not on the line or point layers, with the same settings? And does anyone have a solution?


texts:
I am probably making a simple mistake but can't seem to figure it out, only a beginner.The Select By Location Script is working and it exports to excel when it is run from Python GUI, it also works in ArcGIS when i use the Pyton Shell but when i try to make a tool out of it in ArcGIS it does not work. The error was ERROR 000714, which the help says there is an error in the script. If it works in Python GUI and in ArcGIS Python Shell (picture below) , why does this error appear when i try to run the custom tool? any ideas?

import arcpy

arcpy.env.workspace = "Z:\GIS TEST\Select_by_Location"
BOL='Breakout_Location.shp'
arcpy.MakeFeatureLayer_management('Bld_Locations.shp', 'Bld_Locations_lyr')
arcpy.SelectLayerByLocation_management('Bld_Locations_lyr', "WITHIN_A_DISTANCE", 'Breakout_Location.shp', "2000 Meters", "NEW_SELECTION")
arcpy.SetParameterAsText(0, BOL)

arcpy.MakeFeatureLayer_management('Bld_Locations.shp', 'Bld_Locations2_lyr')
arcpy.SelectLayerByLocation_management('Bld_Locations_lyr', "WITHIN_A_DISTANCE", 'RiversITM.shp', "100 Meters", "ADD_TO_SELECTION")

CSVFile = 'Z:\\Excel.csv'

fields = [f.name for f in arcpy.ListFields('Bld_Locations_lyr')]

for i,f in enumerate(fields):
    if f == 'Shape' or f == 'Shape_Length' or f == 'OBJECTID':
    del fields[i]

with open(CSVFile, 'w') as f:
    f.write(','.join(fields)+'\n') #csv headers
    with arcpy.da.SearchCursor('Bld_Locations_lyr', fields) as cursor:
        for row in cursor:
            f.write(','.join([str(r) for r in row])+'\n')

There will only be one input into the Custom Tool, that is why i only used one parameter.


texts: How to combine labelling with limited boundingbox when plotting a geodataframe

texts:
I am trying to plot a limited boundingbox of my geodataframe setting up some x and y limits to the ax. Everything works but the labelling. This is my piece of code and current result:

fig, ax = plt.subplots(figsize=(15, 15))
labels = []

xlim = [-24.213867, 16.171875]
ylim = [31.952162, 47.100045]
ax.set_xlim(xlim)
ax.set_ylim(ylim)

basemap = world.plot(ax=ax, color='grey', edgecolor='white')
grid.plot(ax=ax)
nearest_cities.plot(ax=ax, marker='o', color='pink', markersize=40)
nearest_cities.plot(ax=ax, marker='o', color='red', markersize=5)

for x, y, label in zip(nearest_cities.geometry.x, nearest_cities.geometry.y, nearest_cities["nearest_name"]):
    labels.append(plt.text(x, y, label, fontsize = 8))

If I comment the last two lines, the bbox works. So my expected result would be like the following image but with labels on top of the markers. How can I make it work?


texts: Creating lines to visualize differences between polygon features in PostGIS

texts:
I've a PostGIS table polygon_b with some polygon features. There is also a table polygon_a which contains the same polygons as polygon_b but with minor changes. Now I want to create lines to visualize the differences between the polygon features.

I suppose that ST_ExteriorRing and ST_Difference will do the job but the WHERE clause seems to be quite tricky.
CREATE VIEW line_difference AS SELECT
row_number() over() AS gid,
g.geom::geometry(LineString, yourSRID) AS geom
FROM
    (SELECT
    (ST_Dump(COALESCE(ST_Difference(ST_ExteriorRing(polygon_a.geom), ST_ExteriorRing(polygon_b.geom))))).geom AS geom
    FROM polygon_a, polygon_b
    WHERE
    -- ?
    ) AS g;

Can anyone help me?

As posted by 'tilt' I've tried ST_Overlaps(polygon_a.geom, polygon_b.geom) AND NOT ST_Touches(polygon_a.geom, polygon_b.geom) but the result is not as expected.
CREATE VIEW line_difference AS SELECT
row_number() over() AS gid,
g.geom::geometry(LineString, your_SRID) AS geom
FROM
    (SELECT
    (ST_Dump(COALESCE(ST_Difference(ST_ExteriorRing(polygon_a.geom), ST_ExteriorRing(polygon_b.geom))))).geom AS geom
    FROM polygon_a, polygon_b
    WHERE
    ST_Overlaps(polygon_a.geom, polygon_b.geom) AND NOT ST_Touches(polygon_a.geom, polygon_b.geom))
     AS g;

workupload.com/file/J0WBvRBb (example dataset)

I've tried to turn the polygons into multilines before using ST_Difference, but the results are still strange.
CREATE VIEW multiline_a AS SELECT
row_number() over() as gid,
ST_Union(ST_ExteriorRIng(polygon_a.geom))::geometry(multilinestring, 4326) AS geom
FROM
polygon_a;

CREATE VIEW multiline_b AS SELECT
row_number() over() as gid,
ST_Union(ST_ExteriorRIng(polygon_b.geom))::geometry(multilinestring, 4326) AS geom
FROM
polygon_b;

CREATE VIEW line_difference AS SELECT
row_number() over() as gid,
g.geom
FROM
    (SELECT
    (ST_Dump(COALESCE(ST_Difference(multiline_a.geom, multiline_b.geom)))).geom::geometry(linestring, 4326) AS geom
    FROM
    multiline_a, multiline_b)
As g;


texts: Is there a way to perform a relative XPath expression filter in Geoserver?

texts:
I have a feature collection that follows the generic structure displayed below:

<Activities>
  <Run>
    <id>1234-5678</id>
    <GeneralInfo>
      <StartTime delay="0">2017-10-18T04:10:00Z</StartTime>
    </GeneralInfo>
  </Run>
  <Walk>
    <id>5678-9012</id>
    <GeneralInfo>
      <StartTime delay="0">2017-10-18T02:10:00Z</StartTime>
    </GeneralInfo>
  </Walk>
</Activities>

I would like use an XPath filter that will select features whose StartTime was at a specific time, regardless of the type of activity.  I've come up with a simple XPath that will target any element named StartTime

XPath=//StartTime

Which lead me to try the following fes Filter:

<fes:Filter xmlns:gml="http://www.opengis.net/gml/3.2" xmlns:gco="http://www.isotc211.org/2005/gco" xmlns:gmd="http://www.isotc211.org/2005/gmd" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:om="http://www.opengis.net/om/2.0"  xmlns:fes="http://www.opengis.net/fes/2.0" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <fes:PropertyIsEqualTo>
    <fes:ValueReference>//StartTime</fes:ValueReference>
    <fes:Literal>2017-10-18T02:10:00Z</fes:Literal>
  </fes:PropertyIsEqualTo>
</fes:Filter>

Which results in an exception being thrown from Geoserver.

<ows:ExceptionReport xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:ows="http://www.opengis.net/ows/1.1" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="2.0.0" xsi:schemaLocation="http://www.opengis.net/ows/1.1 http://internal.geoserver.local/geoserver/schemas/ows/1.1.0/owsAll.xsd">
    <ows:Exception exceptionCode="NoApplicableCode">
        <ows:ExceptionText>java.lang.StringIndexOutOfBoundsException: String index out of range: 0
String index out of range: 0</ows:ExceptionText>
    </ows:Exception>
</ows:ExceptionReport>

Is there a way to perform these relative XPath expression filters in Geoserver?

I have tried not including the leading '//', but I then get an empty feature collection back from geoserver.  Essentially it is unable to find the attribute I am attempting to filter on.


texts: What is the correct way to apply a SpatialFilter on a Layer using GDAL & Python?

texts:
I'm looking at Selecting those Polygons in One Shapefile which are contained inside the Polygon which is stored in another Shapefile.

I have managed to make some progress, and I have modified the sample given here: https://pcjericks.github.io/py-gdalogr-cookbook/layers.html#spatial-filter

I'm trying to figure out exactly what this function does, but the documentation says that this might give incorrect results.

What the correct way to apply a SpatialFilter on a Layer using GDAL & Python?


texts: Add multiple pages within ArcGIS Pro layout?

texts:
I'm creating a 6 page Map-Series enabled PDF document. The first page has a map frame and dynamic text boxes all linked to a map series layer/attributes. The remaining 5 pages contain static text which will never change.

How can I add multiple pages to a layout? I can't see any options.


texts:
When attempting to display the national chargepoint registry from the following URL:

http://chargepoints.dft.gov.uk/api/retrieve/registry/format/json

The following error returns:

http://chargepoints.dft.gov.uk/api/retrieve/registry/format/json is not a valid or recognized data source

I've tried following the instructions in another similar post without success: How to display geoJSON from an URL in a layer. Any help much appreciated.


texts:
I'm new to this field of study.

my teacher ask me to do a hybrid MLC and ISODATA classification for my bachelor thesis. so I got one journal that combine both method, but I think it too complex for bachelor degree because it's also using a decision tree, which I don't understand at all.

so I was thinking, what if I made one for myself, a simple one but useful, accurate, and easy to apply. here how it goes:

I classify image with ISODATA, and then get the result, if accuracy assessment value is low, we redo the classification, if high, we use that result as training area for MLC. Then after we classify the image with MLC, we compare the result also from MLC but with training area from ground truth.

The conclusion will be, if both result are similiar, that mean we don't have to do ground truth anymore to get training area for MLC, we can just use ISODATA to get the training area and then use the result as training area for MLC.

I'm planning on doing all of this with ERDAS IMAGINE 2014

Can you give me your opinion for this method?


texts: Using a custom UI for a QGIS processing algorithm script

texts:
In QGIS 2.x I have a PyQGIS script that utilises a custom form created in Qt Designer. I could simply call the form and connect a button to a function with the following:

class MyWidget(QtGui.QWidget):
    myForm = uic.loadUi('ui_file.ui')

    def doSomething(self):
        print 'Hello World'

myWidget = MyWidget()
searchButton = myWidget.myForm.findChildren(QtGui.QWidget, 'searchButton')[0]
searchButton.pressed.connect(myWidget.doSomething)
myWidget.show()

I'm now trying to update this script for use in QGIS 3.4, which now requires processing scripts to be written by extending the QgsProcessingAlgorithm class. I've used the template script as described in the documentation, and used createCustomParametersWidget() to show my custom ui in place of the default form.

However, I can't figure out how to get the signal/slot code to work correctly. I've tried defining the signal in the createCustomParametersWidget function and in the initAlgorithm function. When the script is run, the custom form shows, but the button is not connected.

class MyTool(QgsProcessingAlgorithm):

    INPUT = 'INPUT'
    OUTPUT = 'OUTPUT'

   class MyWidget(QWidget):
        myForm = uic.loadUi('ui_file.ui')

        def exec_ (self):
            exec_ = self.myForm.show()

        def doSomething(self):
            print('hello World!')

    def createCustomParametersWidget(self, QWidget):
        myWidget = self.MyWidget()
        searchButton = myWidget.myForm.findChildren(QToolButton, 'searchButton')[0] #also tried this in initAlgorithm()
        searchButton.pressed.connect(myWidget.doSomething)

It seems the initAlgorithm takes a configuration parameter (but I can't find any documentation on how to set this up, and this may not be the solution anyway).

Could anybody offer a code snippet on how to make this work in QGIS 3.4, or some links to some further information on how to the this up correctly? If possible, I want to avoid making this into a fully fledged plugin.


texts: Expressions Plus plugin and aggregate functions QGIS 2:18

texts:
I used the function IsSelected('nome_layer') in group_by of the aggregate function Sum:

sum ($ area, group_by: = IsSelected ('layer'))

I notice a strange behavior of the labels, sometimes labeled the selected features (fig. 1) sometimes complementary feature (fig.2);

what am I doing wrong?

fig.1

fig.2


texts: QGIS field-calculator time formating produces NULL

texts:
I'm struggling to format a QString attribute field to a DateTime or Date field.
I've been trying the built in functions to_date to_time to_datetime of QGIS 3.10.14 and none of them seem to do the trick. The layer unit is meter not degrees. I even created a spatial index, not knowing whether it would help. I saved the changes on the layer, of the project.
The result is always the same - the created field is NULL.  I seem to overlook something fundamental. The weirdest part that the Field Calculator shows a valid result, as displayed in the image below.

Any suggestions?


texts: What is the intended Behavior of -skipfailures in ogr2ogr?

texts:
I am trying to run ogr2ogr to convert some geojson to shapefile. All is good until any of the geojson features have an invalid geometry (for example, the coordinate array of point geometry is empty).

I included the -skipfailures flag, but ogr2ogr quits as soon as it hits the first invalid feature. Am I not doing something wrong here?

Example ogr2ogr call:

ogr2ogr -skipfailures -f "ESRI Shapefile" myshape.shp mygeojson.geojson


texts: "Error in Parcel Selection" in SaeX in ArcGIS 9.3

texts:
I am using SaeX extension in ArcGIS 9.3 and whenever i try to select a parcel it gives:

  Error occured in parcel selection.

Other tools are working fine.


texts: Error when saving vector geometry with crs definition in PyQGIS

texts:
I came up with the following issue in PyQGIS (2.8.9-Wien on Windows) when trying to save a created point to a vector file: the point was saved as (inf,inf) in the vector file if the coordinates in Lambert 72 were used with the EPSG code 31370 to create the crs. The point was properly saved when using WGS84 and the respective EPSG code 4326. See code example below.

# create a point and add to project
v_path = r'C:\tmp\test.shp'
crs_epsg = 31370
pt = (161000, 224000)
nlayer = QgsVectorLayer("Point", "poi", "memory")
pr = nlayer.dataProvider()
# add fields and update vector layer
pr.addAttributes([QgsField("name", QVariant.String)])
nlayer.updateFields()
# add a feature
fet = QgsFeature()
fet.setGeometry(QgsGeometry.fromPoint(QgsPoint(pt[0], pt[1])))
fet.setAttributes(["poi"])
pr.addFeatures([fet])
# Commit changes
nlayer.updateExtents()

# create the crs
crs = QgsCoordinateReferenceSystem()
crs.createFromId(crs_epsg, QgsCoordinateReferenceSystem.EpsgCrsId)
print crs.isValid()

# save to vector file
_writer = QgsVectorFileWriter.writeAsVectorFormat(nlayer, v_path, 'system', crs, "ESRI Shapefile")
del _writer
del nlayer

# load created vector file and check feature
nlayer = QgsVectorLayer(v_path, 'poi', "ogr")
f = QgsFeature()
features = nlayer.getFeatures()
for f in features:
    print "F:", f.id(), f.attributes(), f.geometry().asPoint()

This returns the following output:

True
F: 0 [u'poi'] (inf,inf)

But with pt = (4.52656, 51.3257) and crs_epsg = 4326 it works fine?

True
F: 0 [u'poi'] (4.52656,51.3257)

The solution for me could be easy (converting the points to WGS84 before saving).

Can anyone could point me to the error in my code?


texts: Now that TileMill (specifically TileMill-v0.10.1-291-g31027ed) supports FileGDB https://www.mapbox.com/blog/tilemill-adds-support-for-esri-filegdb-format/ ...Is it possible to add label annotations that were created in ArcGIS to TileMill? I added a gdb containing some highway shields as annotations and they imported into TileMill but show up as square polygons (probably from the default styling option). I just wasn't sure if I'm supposed to be styling them based on a marker, shield, etc. or if styling them is a lost cause because TileMill fails to support annotated labels.


texts: How to retrieve highlighted dynamic map image in ArcGIS rest API for given polygon geometries?

texts:
I am using ArcGIS rest api export map operation in dynamic map service to retrieve png for given extent. I want to highlight specific polygon (already have geometry of polygon) with any colour using rest api and retrieve png from server. Is it possible to highlight (add graphics) specific geometry using rest api only?

Note that client is in Java and I don't want to introduce any new technology (ArcObjects etc).


texts: Import shapefiles on my PostGIS table using R

texts:
Using R, I would like to import shapefiles that are in a directory on a geometric table I have on PostgreSQL/PostGIS.

What is the instruction that allows me to do it?

EDIT: here is my code to generate contour lines from the PostgreSQL data, the result is a shapefiles I would save .shp on my table PostGIS  already create before:

library(ggplot2)
library(gstat)
library(sp)
library(maptools)
library(rgdal)
library(zoo)
library(xts)
library(RPostgreSQL)
library(spacetime)
library(raster)
library(foreign)
dbname = "postgis"
user = "postgres"
password = "***"
#password = ""

drv <- dbDriver("PostgreSQL")
con <- dbConnect(drv, dbname="BDDMeteo", user=user, password="***",host='localhost', port='5432')
estonia_air_temperature_2=dbGetQuery(con,"select lat,lon,temperature from observation, station, date where id_s=nom_ville AND id_d=id_date")

estonia_air_temperature_2_test <- estonia_air_temperature_2  # duplicate air temp. data file
estonia_air_temperature_2_test$x <- estonia_air_temperature_2_test$lon  # define x & y as longitude and latitude
estonia_air_temperature_2_test$y <- estonia_air_temperature_2_test$lat
coordinates(estonia_air_temperature_2_test) = ~x + y

plot(estonia_air_temperature_2_test)
x.range <- as.numeric(c(-8.5, 9.5))  # min/max longitude of the interpolation area
y.range <- as.numeric(c(26.97, 37))  # min/max latitude of the interpolation area
grd <- expand.grid(x = seq(from = x.range[1], to = x.range[2], by = 0.1), y = seq(from = y.range[1],
                                                                                to = y.range[2], by = 0.1))  # expand points to grid
coordinates(grd) <- ~x + y
gridded(grd) <- TRUE
plot(grd, cex = 1.5, col = "grey")
points(estonia_air_temperature_2_test, pch = 1, col = "red", cex = 0.4)

idw <- idw(formula = temperature ~ 1, locations = estonia_air_temperature_2_test,
           newdata = grd)  # apply idw model for the data
residual_grid = raster(idw, "var1.pred")
plot(residual_grid)
x <- rasterToContour(residual_grid)
class(x)
plot(residual_grid)
plot(x, add=TRUE)
writeSpatialShape(x, "/interpolate/contour/contourline")

Edit:

> ogrDrivers()
             name write
1      AeronavFAA FALSE
2          ARCGEN FALSE
3          AVCBin FALSE
4          AVCE00 FALSE
5             BNA  TRUE
6             CSV  TRUE
7             DGN  TRUE
8             DXF  TRUE
9          EDIGEO FALSE
10 ESRI Shapefile  TRUE
11     Geoconcept  TRUE
12        GeoJSON  TRUE
13       Geomedia FALSE
14         GeoRSS  TRUE
15            GML  TRUE
16            GMT  TRUE
17       GPSBabel  TRUE
18  GPSTrackMaker  TRUE
19            GPX  TRUE
20            HTF FALSE
21         Idrisi FALSE
22            KML  TRUE
23   MapInfo File  TRUE
24         Memory  TRUE
25   MSSQLSpatial  TRUE
26           ODBC  TRUE
27            ODS  TRUE
28        OpenAir FALSE
29    OpenFileGDB FALSE
30         PCIDSK  TRUE
31            PDF  TRUE
32            PDS FALSE
33         PGDump  TRUE
34           PGeo FALSE
35            REC FALSE
36            S57  TRUE
37           SDTS FALSE
38       SEGUKOOA FALSE
39           SEGY FALSE
40            SUA FALSE
41            SVG FALSE
42            SXF FALSE
43          TIGER  TRUE
44        UK .NTF FALSE
45            VRT FALSE
46           Walk FALSE
47           WAsP  TRUE
48           XLSX  TRUE
49         XPlane FALSE


texts: extracting statistics from a raster within the extent of a polygon

texts: I have a polygon representing the extent of a study area. This polygon is overlaid to a raster representing the quality of land, which is classified into 5 classes (from bad [1] to good [5]). I would like to have arcGIS (10.1, advanced license) calculate how much land within the polygon is bad, how much is mediocre, and so on. How can I accomplish this?


texts: Dissolving adjacent parcels by same owner but include those parcels located across street with same owner?

texts:
In trying to determine the total contiguous stream length per landowner, I dissolved parcels by parcel owner name attribute and intersected these dissolved parcels with streams.  However, dissolving by parcel owner doesn't include parcels located across the street (right-of-way) when 'Create multipart features' is unchecked since they don't share common edges; and when I do check 'Create multipart features', it creates the desired result across the street but also includes piecemealed parcels per owner that are spread out in space which are not considered contiguous.  I was even thinking of modeling this somehow in Spatial Analyst.  I'm using ArcGIS 10.1 with a basic license.

Any ideas of how to achieve the desired result?


texts: Change annotation font size on a IGeoFeatureLayer

texts:
I would like to change the font size for the annotations on an IGeoFeatureLayer.  I tried to set it in the lpLabelEngine.Symbol.Font.Size below but the set doesn't persist.

     IGeoFeatureLayer geoLayer = thisLayer as IGeoFeatureLayer;
            if (geoLayer != null)
            {
             geoLayer.DisplayAnnotation = annotationsOn;
            IAnnotateLayerPropertiesCollection propertiesColl = geoLayer.AnnotationProperties;
            IAnnotateLayerProperties labelEngineProperties = new LabelEngineLayerProperties() as IAnnotateLayerProperties;
            IElementCollection placedElements = new ElementCollectionClass();
            IElementCollection unplacedElements = new ElementCollectionClass();
            propertiesColl.QueryItem(0, out labelEngineProperties, out placedElements, out unplacedElements);
            ILabelEngineLayerProperties lpLabelEngine = labelEngineProperties as ILabelEngineLayerProperties;
            if(geocode != null) lpLabelEngine.Expression = geocode;
            if(annotationLabelColor != null) lpLabelEngine.Symbol.Color = annotationLabelColor;
            labelEngineProperties.AnnotationMinimumScale = minScale;
            labelEngineProperties.AnnotationMaximumScale = maxScale;
                if (fontsize != null) lpLabelEngine.Symbol.Font.Size = Convert.ToDecimal(fontsize);

Any ideas?

ArcEngine 10, C#, VS2010


texts: When did the word "geospatial" first come into usage?

texts: For historical reasons, I'm trying to figure out when the term "geospatial" came into use? The earliest usages I can find are from around 1994, but I would be really surprised if it wasn't used earlier. Please share references, if you have them.


texts: ESRI IdentityManager.registerToken unable to register Multiple servers

texts:
I have multiple secured service that needs to be added to my map for which I need to register multiple servers using ESRi IdentityManager.registerToken method

      IdentityManager.registerToken({
            'server': "https://testurl1",
            'token': "abcasdasasd"
           });

           IdentityManager.registerToken({
            'server': "https://testurl2",
            'token': "MV5VF5SeK9pigO5094rC9MMEyY39SnmB"
           });

The First IdentityManager.registerToken for the testurl1 is able to register the token successfully but the  IdentityManager.registerToken for the testurl2 throws the following error

Cannot read property '3' of null
    at new k (init.js:457)
    at b.hasSameOrigin (init.js:731)
    at init.js:731
    at Object.some (init.js:70)
    at Object.b.canUseXhr (init.js:731)
    at Object. (init.js:1999)
    at Object.forEach (init.js:71)
    at Object.registerServers (init.js:1998)
    at Object.registerToken (init.js:1999)
    at Object.onReceiveData (Widget.js?wab_dv=2.10:48) "TypeError: Cannot read property '3' of null
    at new k (https://js.arcgis.com/3.31/init.js:457:489)

Please provide some thoughts on how to fix this issue..

Thanks
Muthu


texts: Generate 3D elevation points grid from 3D contour lines!

texts:
I have a shapefile with 3D contour lines.

I want that based on this information to generate for all this area a 3D grid of points with certain distance from each-other like every 50cm.

I use Arcgis Desktop Advanced 10.2.2


texts: How to getfeatureinfo results in a column than a row from geoserver postgis openlayers

texts:
Below is a snippet of a code I modified from openlayers. Basically, I click on a wms map and data from postgis is queried by geoserver and served on a single row table. There are 11 variables in the postgis layer called var1 through var11 if which one is dynamically chosen by the user from 20 available options. I would like to arrange the variables in a column instead of the default single row. In addition, I’d like to place the resulting table in a fix location (e.g., 900 pixels from the left edge and 1pixel from the top) within the screen instead of showing up at where one clicks or offsetted with some constant from where the click was done. Lastly, var1 to var 4 columns contain an html link to an image and I’d like to have a word or two that can be clicked to open the linked image. I suspect it got something to do with the style and popover sections of the code but have very minimal coding background so I’m stumped.

<style>
  html, body, #map {
    height: 100%;
  }
  #map {
    padding-top: 50px;
  }
  .legend {
    position: absolute;
    z-index: 1;
    left: 10px;
    bottom: 10px;
    opacity: 0.6;
  }
  .popover {
    max-width: 440px;
  }
  .popover-title {
    float: right;
    background: none;
    border: 0;
    cursor: pointer;
  }
  .popover-content iframe {
    width: 400px;
    height: 120px;
    border: 0;
  }
</style>

and the javascript:

var popup = new ol.Overlay({
  element: $('#popup')
});
olMap.addOverlay(popup);

// Handle map clicks to send a GetFeatureInfo request and open the popup

olMap.on('singleclick', function(evt) {

  var view = olMap.getView();

  var url = wmsLayer.getSource().getGetFeatureInfoUrl(evt.coordinate,
      view.getResolution(), view.getProjection(), {'INFO_FORMAT': text/html'});

  popup.setPosition(evt.coordinate);

  $('#popup-content iframe').attr('src', url);

  $('#popup')

    .popover({content: function() { return $('#popup-content').html(); }})

    .popover('show');

  // Close popup when user clicks on the 'x'

  $('.popover-title').click(function() {

    $('#popup').popover('hide');

  });

});


texts: QGIS python plugin: PyQGIS for spatial join

texts:
For a QGIS 2.0 and 2.2 Python plugin, I am trying to update one layer's field attribute with another layer's field attribute based on geometry function QgsGeometry.within(). My first layer is a point layer, the second layer is a buffer of a vector polyline layer containing azimuth measurements. I would like to update the point layer to include azimuth information of the buffer polygon it is within (essentially a spatial join). It is to automatize the process described here: How to horizontally align irregularly shaped point symbol to vector layer in QGIS?. Currently, I'm not getting an error, but nothing is happening either.

rotateBUFF = my buffer polygon layer
pointLayer = my point layer to obtain azimuth data

rotate_IDX = rotateBUFF.fieldNameIndex('bearing')
point_IDX = pointLayer.fieldNameIndex('bearing')
rotate_pr = rotateBUFF.dataProvider()
point_pr = pointLayer.dataProvider()
rotate_caps = rotate_pr.capabilities()
point_caps = point_pr.capabilities()
pointFeatures = pointLayer.getFeatures()
rotateFeatures = rotateBUFF.getFeatures()

for rotatefeat in rotateFeatures:
    for pointfeat in pointFeatures:
        if pointfeat.geometry().within(rotatefeat.geometry()):
            pointID = pointfeat.id()
            if point_caps & QgsVectorDataProvider.ChangeAttributeValues:
                bearing = rotatefeat.attributes()[rotate_IDX]
                attrs = {point_IDX : bearing}
                point_pr.changeAttributesValues({pointID : attrs})


texts: How to Geocode a list of address into a map in qgis?

texts: I have a list of country and city capitol  on a Excel, and I want to upload them into a map in QGIS?


texts:
I'd like to open a popup in on the midpoint of a polyline/linestring in Mapbox GL, like it's the default case in Leaflet.

How can I do this via turf.js? If I use centroid/center/centerOfMass it will get a point usually not on the line, but somewhere next to it.

My only idea is to:

Measure the lenght of the line in some units
Use along() with the half of the measured distance.

Is this really the best option?


texts: drawing concentric circles in arcmap 10.1

texts: I am using ArcMap 10.1 and I would like to know if there is a way to draw circles at a specific distance from a georeferenced coordinate. I have a property and I would like to know how far you can walk from it within different amounts of time, and I would like to portray this with multiple circles.


texts:
I have a raster with 500 meter resolution that was in sinusoidal projection, so I first converted it to WGS84using the raster package in R using the following code:
library (raster)
library(maptools)
#Load the tif file
tif <- raster("~/Sinu.tif")
# Extract the SRS of the original tif
CRS <- "+proj=sinu +lan_0=0 +x_0=0 +y_0=0 +R=631007.181 +units=m +no_defs"
#Transform the projection to WGS84
llcrs <- "+proj=longlat +ellps=WGS84 +datum=WGS84"
#Create a data frame of the original sinusoidal tif
DF <- as.data.frame(Sinu, xy=TRUE)
data(DF)
#Create the Transform
rll <- projectRaster(Sinu,crs = llcrs)
#ras_sinu <- raster(nrow=391, ncol=695, ext=extent(-9031818, 4773511, -870816, 4954666), crs=incrs)
#rs <- rasterImage(w, ras_sinu)
plot(rll)
Sin2WGS84 <- writeRaster(rll,filename = "~/Sinu2WGS84.tif")

The conversion was done to align the pixels of this finer resolution raster with another raster which is in WGS84and has a spatial resolution of 1/16 (approximate to 6 kilometers). The coarser resolution raster has a representative gridded polygon file. I then imported the transformed tif file to ArcGIS and converted it to a Polygon in ArcGIS. This was done to see whether 12x12 (12 rows by 12 columns of pixels) 500M meter pixels fall in a coarser pixel. But as you can see in the following image I am getting either 10x10 or 11x11 pixels. What can I do in R or ArcGIS to solve this alignment problem? If anyone wants to access the data, please click here.


texts: How to hide certain views from ArcCatalog?

texts:
I have users that connect to an enterprise geodatabase via OS authentication. They use a couple of methods:

An .SDE connection in ArcMap
A ODBC connection for MS Office

The users need to be able to connect to all of the views in the GDB using the ODBC connection.

However, the users only need to connect to some of the views in ArcMap. The rest of the views are not needed in ArcMap (they just create clutter--there are hundreds of them).

Is there a way to hide certain views from ArcCatalog?


texts:
I have a map service which is using data from a local repository, viewed using the ESRI Silverlight API.  I am dynamically adding a tiled map service for imagery from ESRI.

The map service can be viewed here:
http://server.arcgisonline.com/ArcGIS/rest/services/ESRI_Imagery_World_2D/MapServer

When the tiled map service is added, I can not zoom in closer than 0.002 Decimal Degrees, regardless of whether the map service is visible or not.  This is despite the fact that according to the ArcGIS Services Directory, the imagery layer has levels of detail at much smaller resolutions.

Is there a reason why I would not be able to zoom in closer than 0.002 Decimal Degrees?


texts: Reading, modifying and writing a geotiff with GDAL in python

texts:
I'm trying the learn the ropes of Remote Sensing image processing using Python GDAL bindings and numpy. As a first attempt, I'm reading a Landsat8 geotiff file, do a simple manipulation and write the result to a new file. The code below appears to work fine, except that the original raster is dumped in the output file, rather than the manipulated raster.

Any comments or suggestions are welcome, but particularly notes on why the manipulated raster does not show in the result.

import os
import gdal

gdal.AllRegister()

file = "c:\~\LC81980242015071LGN00.tiff"
(fileRoot, fileExt) = os.path.splitext(file)
outFileName = fileRoot + "_mod" + fileExt

ds = gdal.Open(file)
band = ds.GetRasterBand(1)
arr = band.ReadAsArray()

[cols, rows] = arr.shape
arr_min = arr.Min()
arr_max = arr.Max()
arr_mean = int(arr.mean())

arr_out = numpy.where((arr < arr_mean), 10000, arr)

driver = gdal.GetDriverByName("GTiff")
outdata = driver.Create(outFileName, rows, cols, 1, gdal.GDT_UInt16)
outband = outdata.GetRasterBand(1)
outband.WriteArray(arr_out)
outdata = None

print arr_min
> 0
print arr_max
> 65535
print arr_mean
> 4856

I use Python 2.7.1 on a Windows 7 32 bit machine.


texts: How to merge polygons but keeping all the inner polygons?

texts:
I am merging polygons using ST_Union, but what I get is the "exteriors" only of the polygons. How can I merge them, but keeping all the eventual inner polygons?

For example given two polygons A and B with relative inner rings, how can "stack" them when they overlap (ST_Intersects):


texts:
An isoscope shows me how far I can get within a specific amount of time. Is there a similar visualization that shows me everything within X distance by roads of location?

If so what would be the easiest way to create my own visualizations of this type? Is there an R package which provides this, and allows me to fetch map data?


texts: Verifying that account used to set up ArcGIS server can connect to server?

texts:
I have ArcGIS Server and IIS installed on //GISSERVER01 and on //GISSERVER02 I have my raster data with ArcGIS Desktop installed from which I create my services. Until this morning everything was working fine and I was able to publish services from //GISSERVER02. There is no internal firewall.
In a web browser I am able to connect to arcgisserver manager as Administrator and have full access to all functions. I can stop services but I can´t start any of the services. The following errors are shown when trying to start a service. I can create the service file but I cannot publish any ArcGIS Server 10.3.1 Services and get the following Error:

  001359: Failed to connect to the server.

  Description

  The attempt by the server-side publishing tool (systems tools) to
  make an admin connection to the server has failed.

  This server system tool is invoked by the publishing client (Manager
  or ArcGIS for Desktop) when you create a
  service.

  Solution

  Verify that the account you used to set up the server can connect to
  the server.

What steps can I take to verify that the account I used to set up the server can connect to the server?


texts: Statistics from one image to another in Google Earth Engine

texts:
I now have two images of the same place, image1 stores image values and each pixel on image 2 stores time (jday).I want to count the sum of image values of image1 at the same time point on image2.For example, there are 10 time points on January 1, 2000 on image2, and I want to get the sum of these 10 positions image1.The final result is presented in the form of a table or chart.
Here is my example code https://code.earthengine.google.com/0db5219ac3e642d9ac28c73e20c6e7cb
Hope you can answer this question


texts: Dividing a geographical area (bbox) into smaller equal sized square miles

texts: I'm working on an automated system for processing geographical information by bbox using PostGIS with an Openlayers 3 UI.  The user will use dragbox to 'draw' the bbox they want processed. The catch is the area can be no longer than 2 sq miles.  Is there a method for dividing up a bbox into equally sized smaller bbox's that are no larger than 2 sq miles?


texts: Prevent newly created feature class being automatically added by ArcPy (Python AddIn) to Map Document?

texts:
I have a python addin, which as part of its function checks to see if a template point feature class exists on the user's machine.  If it doesn't, it creates this file.  This template file is used elsewhere in the addin to create other files.  I do not want this template file to be added to the map document.

When the addin is run, if the template file doesn't exist it gets created and automatically added into the currently active data frame.  Is there something I can add to the scrip to prevent this from happening?  Or is the neatest way to let the feature class be added to the data frame and then add a line to the script to remove it from the map document?

My code is currently as follows:

coordsys = arcpy.SpatialReference(27700)

#Define the schema file for the output
schemafile = defgdb + "/Template"

#Check to see if the schema file exists, if it doesn't, create it.
if arcpy.Exists(schemafile):
    pass
else:
    arcpy.CreateFeatureclass_management(defgdb, "Template", "POINT", "", "", "", coordsys)

My current method for removing the layer from the map document (rather than preventing the layer from being added) is:

for layer in arcpy.mapping.ListLayers(mxd, "Template"):
    arcpy.mapping.RemoveLayer(df,layer)


texts: How to export from ModelBuilder as script and then run using Python. 3.4?

texts:
I'm currently studying for a Msc in GIS. In my assignment I have been asked to produce a model, export that model as a script file and then run that script as a stand alone in IDLE.
So I have created the model in model builder and exported it as a script and have even had a go at cleaning the script. But whenever I try to run the script I get an error. I have attached a jpeg of my model and have also attached the script file. Does anyone know where I'm going wrong? Im using Arcmap 10.1 and python 3.4

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# floodrisk.py
# Created on: 2014-04-21 11:34:26.00000
#   (generated by ArcGIS/ModelBuilder)
# Usage: floodrisk <River_thames> <POI> <SOA> <POI_clipped> <SOA_clipped__6_>
# Description:
# ---------------------------------------------------------------------------

# Import arcpy module
import arcpy
arcpy.env.workspace = "C:\Users\Imran\Documents\MASTER GIS\ASSIGNMENT4\Floodrisk.gdb"
arcpy.env.overwriteOutput = True

# Script arguments
River_thames = arcpy.GetParameterAsText(0)
if River_thames == '#' or not River_thames:
    River_thames = "River_thames.shp" # provide a default value if unspecified

POI = arcpy.GetParameterAsText(1)
if POI == '#' or not POI:
    POI = "POI.shp" # provide a default value if unspecified

SOA = arcpy.GetParameterAsText(2)
if SOA == '#' or not SOA:
    SOA = "SOA.shp" # provide a default value if unspecified

POI_clipped = arcpy.GetParameterAsText(3)
if POI_clipped == '#' or not POI_clipped:
    POI_clipped = "POI_clipped.shp" # provide a default value if unspecified

SOA_clipped__6_ = arcpy.GetParameterAsText(4)
if SOA_clipped__6_ == '#' or not SOA_clipped__6_:
    SOA_clipped__6_ = "SOA_clipped.shp" # provide a default value if unspecified

# Local variables:
River_thames_buffer1 = "River_thames.shp"
SOA_clipped = "River_thames_buffer1.shp"
SOA_clipped__2_ = "SOA_clipped.shp"
SOA_clipped__3_ = "SOA_clipped__2_.shp"
SOA_clipped__4_ = "SOA_clipped__3_.shp"
SOA_clipped__5_ = "SOA_clipped__4_.shp"

# Process: Buffer
arcpy.Buffer_analysis(River_thames, River_thames_buffer1, "300 Meters", "FULL", "ROUND", "NONE", "")

# Process: Clip
arcpy.Clip_analysis(POI, River_thames_buffer1, POI_clipped, "")

# Process: Clip (2)
arcpy.Clip_analysis(SOA, River_thames_buffer1, SOA_clipped, "")

# Process: Add Field
arcpy.AddField_management(SOA_clipped, "AreaB", "DOUBLE", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")

# Process: Calculate Field
arcpy.CalculateField_management(SOA_clipped__2_, "AreaB", "!shape.area! ", "PYTHON", "")

# Process: Calculate Field (2)
arcpy.CalculateField_management(SOA_clipped__3_, "AreaB", "[AreaB]/1000000", "VB", "")

# Process: Add Field (2)
arcpy.AddField_management(SOA_clipped__4_, "areaP", "DOUBLE", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")

# Process: Calculate Field (3)
arcpy.CalculateField_management(SOA_clipped__5_, "areaP", "[AreaB]/ [AreaSQKM]", "VB", "")


texts:
I'm new to QGIS.

I'm trying to add a map just of the State of NJ. I figured out to load Google street map, but that's as far as I got.

My goal is to get a Map of NJ with its counties. I found a few maps online that are PDFs. Not sure how to bring it into QGIS and still be able to  geocode with that.


texts: Create a projection to emphasize a given direction in ArcMap?

texts:
Suppose one had a dataset (centered on a given coordinate) of lines that were very close together in terms of longitude but also had important turns/bends far separated by latitude.  In order to map this, I want to over-emphasize the x direction and de-emphasize the y-direction.

Can one do this by creating a projection in ArcMap? If so, how?


texts: Is it possible to create dynamic print templates in Composer? (QGIS 2.8)

texts:
I'm trying to create a universal template for different map sizes.

Typically, our maps have to fold up to A4 size, which means that the page size is a multiple of A4 in lots of different combinations. The template I want to create should automatically adjust the position and size of map elements (map frame, north arrow, scale bar, ...) in relation to the chosen page size.

Using the Expression string builder, is it possible to refer to the values of Width an Height in the Paper and quality tab?

I know that something similar can be be done using Atlas and the Expression string builder as shown in this example: Multiple format map series using QGIS 2.6 – Part 2


texts:
I see quite often in Python GDAL code that people close datasets at the end of their script. Why does it makes sense to close a dataset in Python GDAL? Are there any consequences if I don't do it?

import gdal

# open dataset
ds = gdal.Open('test.tif')

# close dataset
ds = None


texts: QGIS Modeler - Select by expression error

texts:
QGIS Modeler works fine with my data if I store the layers as temporary layers, but soon as I will save the output in the selected folder (as shapefile or geopackage) I get a wrong result. The select by expresion is not working fine anymore.

The expresssion is: "AREA" != maximum("AREA")

The select by expresion is not working fine anymore


texts: Weighted flow accumulation: using deforestation as weight?

texts:
Right off the bat: if a small river has 10% of its upstream catchment deforested, and a big river has 5%, I want that small river to have higher output values. But now it won't because of how much more flow accumulation the big river has.

I need to figure out which rivers in a certain area have the most potential influence from upstream deforestation.
So, I used ArcGIS Pro to calculate flow accumulation with my deforestation raster with pixel values ranging from 0 - 1.
The result, of course, still is biased towards bigger rivers with more flow accumulation.
Any idea on how I can separate flow accumulation from deforestation values?


texts:
I have a Point

RW     HW     EPSG
3011.65;341146.38;31256

which ST_Transform transforms into projection
EPSG=4326, Lat;Lon
48.2086988651298;16.3686458354975

However, when using http://cs2cs.mygeodata.eu/ and checking with Google Maps and openstreetmap, the actual result should be
48.2087130847;16.3726546685;

This is much to inaccurate.

Am I doing sthg wrong when using Transform?

EDIT1: That's the way I use it:

select ST_Y(a.b), ST_X(a.b) from
  (select ST_Transform(ST_SetSRID(ST_MakePoint(3011.65, 341146.38), 31256), 4326) as b) as a


texts:
I have a requirement to add a 'heat map' style layer to a Javascript/jQuery based application that I am working on.  So far, I have used heatmap.js with no problems at all.

My only concern with this approach is with how well it will handle a large number of points. I have asked for some type of upper bound, and the answer I got was 'a-lot', most likely going to be in the thousands.

So I am thinking I may try to create a WMS style web service - where a bounding box, image height & width are the request parameters, the server queries the points and does the rendering, and the resultant image sent back to client.

Does anyone know any such open source apis/sdks that contain this functionality?  I know such tools exist in qgis and arcgis desktop - but i am looking to do this in a service.


texts: Why are features added to layers via Python script not being displayed?

texts:
I am using PyQGIS with QGIS version 2.18.11.

Can layer.addFeature() be used on a layer which is a clone of a layer which was imported from a .xml file?

This code throws no errors:

def updateDestination(self, importLayer, masterName):
    masterLayer = Utils.getLayerByName(masterName)
    self.iface.setActiveLayer(masterLayer)
    self.iface.actionDuplicateLayer().trigger()

    tempLayer = Utils.getLayerByName(masterName + " copy")
    tempLayer.startEditing()
    tempLayer.setLayerName(masterName + "_temp")

    for feature in importLayer.getFeatures():
        print "Adding feature", feature["ID"]
        tempLayer.addFeature(feature)

    tempLayer.commitChanges()
    tempLayer.updateFields()
    tempLayer.updateExtents()

But, when I go to check the attributes table, it still only has the original data. If the action sequence I require is

Clone layer
Add features to the end of that layer
UI button to trigger continuing
Delete clone layer
Add features to the end of the original layer

Am I limited to using something that isn't an imported vector file? Data was originally entered into an Access Database. xlm was used to transfer some test data across.


texts: Accessing processing algorithm results from Whitebox Tools in PyQGIS

texts:
When I run a GDAL processing algorithm from the Python Console in QGIS that loads a result raster:
resultClip = processing.runAndLoadResults("gdal:cliprasterbymasklayer", {
                                          'ALPHA_BAND' : False,
                                          'CROP_TO_CUTLINE' : True,
                                          'DATA_TYPE' : 0,
                                          'EXTRA' : '',
                                          'INPUT' : 'E:/geodata/nsw100k.vrt',
                                          'KEEP_RESOLUTION' : False,
                                          'MASK' : 'C:/Users/brennant/Downloads/data.geojson',
                                          'MULTITHREADING' : False,
                                          'NODATA' : None, 'OPTIONS' : '',
                                          'OUTPUT' : 'TEMPORARY_OUTPUT',
                                          'SET_RESOLUTION' : False,
                                          'SOURCE_CRS' : None,
                                          'TARGET_CRS' : None,
                                          'X_RESOLUTION' : None,
                                          'Y_RESOLUTION' : None })

I can access the temporary raster output file via
resultClip['OUTPUT']
However, when I do the same thing with one of the Whitebox Tools algorithms (eg BreachDepressions):
resultBreach = processing.runAndLoadResults("wbt:BreachDepressions", { 'dem' : resultClip['OUTPUT'], 'max_depth' : '20.0', 'flat_increment' : 0.00001, 'fill_pits' : 1, 'output' : 'TEMPORARY_OUTPUT' })

resultBreach['output'] returns a QgsProcessingOutputLayerDefinition variable:
<QgsProcessingOutputLayerDefinition {'sink':TEMPORARY_OUTPUT, 'createOptions': {}}>}

The temporary layer is created and loaded in the layers window, but from a PyQGIS processing point of view, I can't see how to access it for use in future steps.
And if I look at both algorithms using processing.algorithmHelp
eg
processing.algorithmHelp("gdal:cliprasterbymasklayer")

the output type for both is the same (QgsProcessingOutputRasterLayer)
Suggestions welcomed.


texts: QGIS 3 plugin not able to detect Python packages

texts:
I have installed the custom plugin which was created through plugin builder 3. I just added an import pymango statement and have deployed the plugin to QGIS. This gave below error with stacktrace.

  ModuleNotFoundError: No module named 'pymongo'

Traceback (most recent call last):
  File "/usr/lib/python3/dist-packages/qgis/utils.py", line 335, in startPlugin
    plugins[packageName] = package.classFactory(iface)
  File "/home/ubuntu/.local/share/QGIS/QGIS3/profiles/default/python/plugins/sample/__init__.py", line 35, in classFactory
    from .sample import custom_plugin
  File "/usr/lib/python3/dist-packages/qgis/utils.py", line 672, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/home/ubuntu/.local/share/QGIS/QGIS3/profiles/default/python/plugins/sample/custom_plugin.py", line 33, in
    import pymongo
  File "/usr/lib/python3/dist-packages/qgis/utils.py", line 672, in _import
        mod = _builtin_import(name, globals, locals, fromlist, level)
    ModuleNotFoundError: No module named 'pymongo'

Python Path:
/usr/share/qgis/python
/home/ubuntu/.local/share/QGIS/QGIS3/profiles/default/python
/home/ubuntu/.local/share/QGIS/QGIS3/profiles/default/python/plugins
/usr/share/qgis/python/plugins
/usr/lib/python36.zip
/usr/lib/python3.6
/usr/lib/python3.6/lib-dynload
/home/ubuntu/.local/lib/python3.6/site-packages
/usr/local/lib/python3.6/dist-packages
/usr/lib/python3/dist-packages
/home/ubuntu/.local/share/QGIS/QGIS3/profiles/default/python

The pymango python library is installed in /home/ubuntu/.local/lib/python3.6/site-packages folder. You can see this folder in python path, still QGIS3.6 desktop application is not able to locate the package. It is not detecting any external Python package that I installed. Any solution?

System details:

Python version: 3.6.7 (default, Oct 22 2018, 11:32:17) [GCC 8.2.0]

QGIS version: 3.6.0-Noosa Noosa, 5873452

Ubuntu: 18.04


texts: Learning resources for beginning differential topology for a programmer?

texts: With this question I'm looking for suggestions about how and where to start learning differential topology starting from a good base of linear algebra and geometry, also consider that this know-how is supposed to be used in a programming context and for reading and writing algorithms.


texts: Accessing ORACLE table from QGIS gives TNS-12545 error?

texts:
I have to access one ORACLE tab from QGIS but it does not work.

From Instantclient_12_1 the connection and parameters are ok, while from QGIS we always get TNS-12545 error.

ORACLE_HOME and TNS_ADMIN var are setted. But from trace log we noticed that the QGIS client seems to be 11_2.
Maybe there is a different client in QGIS dir? And this client does not know our env var?

I don't know how to start to search solutions.


texts: Export from scene layer in ArcGIS Hub to file 3D

texts:
I want make 3D model using 3D buildings and trees from this layer: https://hub.arcgis.com/datasets/esrinl-content::3d-bomen-rd-1
There is any way to export a part from that scene layer to a 3D file?
I tried with ArcGIS Pro making a clip but I only got errors. My knowledge of ArcGIS Pro is very limited.
Edit: I created a new layer as my workspace (rectangle) and made a clip between the new layer and the 3D layer. I got a "layer not found" error. I don't know if there is a problem with the 3D format (maybe other tools are needed to deal with 3D data) or a problem with the data source. My goal is to get a kml with the trees and buildings and use it in SketchUp.


texts:
How would I write a proper geoJson, readable by a GIS, from the results of the rasterstats package?

stats = zonal_stats("somePolygons.shp",
                     "someTif.tif",
                    geojson_out=True)

Something like

import geojson
with open("newfile.geojson", "w") as myFile:
     myfile.write(geojson.dump(stats))

does not work.
The object I am getting from rasterstats is a geoJson. There are 200 Polygons in my geoJson object.
Do I really have to write my geoJson "manually" like suggested in an answere here, or is there a way to do that automatically?


texts: Create an arc of lines from a line and a value

texts:
I'm trying to recreate an Origin-Destination plot like this:

I have managed to munge the data into a MSOA to LAD table and can draw a map like this for one of the origin MSOA.

Which once you allow for the (now ridiculous) distances people in the Peak District commute to work is close.

But I quite like the effect that the author has achieved by "splaying" out the lines. Obviously, with flows of 522 and 371, I can't go for a single line per commuter but it would be nice to produce a proportional arc of lines to show the number of people making the trip.

I thought that I would be able to use the Geometry Generator but without a loop construct, I can't seem to make headway.


texts:
I have a shapefile of concentration values, they range from 0 7000. I want to either contour or display these in log scale (0.1, 1, 10, 100 and 1000) what is the best way of doing this?

I have about 200 files to batch process.

We have ArcGIS Desktop 10.3.1 Basic license with 3D Analyst and Spatial Analyst.


texts: unable to change the work space using arcpy

texts:
I try to create a folder and Geodatabase at the beginning of my code, so all results will be saved in it.
My code is:

class btn_WF(object):
        """Implementation for Addins_addin.button (Button)"""
        # this button allows the user to determine where he/she wants to store all output files
        # The model generates a folder called (), and a File GeoDataBase called ()
        # The model makes the generated GeoDataBase is the current workspace
        def __init__(self):
            self.enabled = True
            self.checked = False
        def onClick(self):
            import pythonaddins, arcpy
            # Let the user select a folder
            Work_Location = pythonaddins.OpenDialog("Please select a folder to store the project and files in", False, "C:/ModelBuilder/AAA","Select a Folder")
            # create a new folder to save the project in
            Project_Folder = arcpy.CreateFolder_management(Work_Location, "Flash Flood Risk Assessment Projcet")
            # create a GeoDataBase to store files and outputs
            Output_File_GDB=arcpy.CreateFileGDB_management(Project_Folder, "FFRA_WRRI.gdb")

            # the problem is the next line
            # the "arcpy.env.workspace" doesnt accept to equal the variable "Output_File_GDB"
            arcpy.env.workspace = Output_File_GDB

            print (arcpy.env.workspace)

it gives me the follwoing error

Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\geoprocessing\_base.py", line 537, in set_
    self[env] = val
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\geoprocessing\_base.py", line 597, in __setitem__
    ret_ = setattr(self._gp, item, value)
RuntimeError: Object: Error in accessing environment <workspace>

Could you help me


texts: CAD To GeoJSON after publishing it shows the job failed

texts:
I have written the code for CAD to GeoJSON. When I ran ArcMap or ArcCatalog it's working fine, after publishing the code it's giving the error. The code is too long. Can you provide the code

ArcGIS Manager Error description:

Tool Properties and Result in Catalog


texts: Get the width and height of a rectangular polygon in PostGIS

texts:
I have a table with point features saved in the my PostgreSQL database, and I can run the following query to get the extent of the point features and I turned the extent back into a geometry:

SELECT ST_AsText(ST_SetSRID(ST_Extent(pts.geom), 4026)) AS Extent_Geom

FROM Table_w_Point_Features pts

Here is how the output geometry of a rectangle polygon that represents the extent of the point features looks like:

POLYGON((-117.536012078 33.7130048820001,
         -117.536012078 34.3860049530001,
         -116.825011941 34.3860049530001,
         -116.825011941 33.7130048820001,
         -117.536012078 33.7130048820001))

Is there any way in PostGIS that I can get the Width and Height of the extent that I created above?


texts: Android WebMap Popup Viewing not showing my layer info

texts:
I used the code for WebMap Popup Viewing from ArcGIS for Developers (https://developers.arcgis.com/android/sample-code/webmap-popup-viewing/) and I just changed the map URL with a map I created and contains 4 layers.

However, every time I tap on a point where I'm supposed to get a Pop Up Dialog containing my marker's details, I see the querying message for a moment and it disappears, I don't see any info.

The application doesn't report any errors, I guess it has something to do with the way I fill in information in the layer? What is the right way to do it?