texts: Color points by value and size them by another one using Arcade in ArcGIS Pro

texts:
I am totally new to the Arcade language and just want to get the position a value has in the list of all values from that column.
More precise: I have a Point Layer in ArcGIS Pro with the columns "Species" and "Diameter". I already used the "Diameter" variable to define the symbology (size of the points in m at ground scale). Now I want to color the points by "Species". Therefore, it seems I have to transform the "Species" values to numeric values (I can only select numeric fields in the color menue or insert an expression). I decided to insert an expression that should give the number that the current "Species" value has in a list of all "Species" values that occur in the layer.
Something like
IndexOf(Distinct($layer, 'Species'), $feature.Species)

Which doesn't work. It fails with the message "Object not found $layer".
What went wrong here?
I also seek other ideas on how to color the points by "Species" (while keeping the diameter variable as the feature size in m at ground scale(!) I am well aware that one can simply set Unique Values in Primary Symbology in order to color by "Species", but I didn't see an option to set color and size in the Primary Symbology menue).


texts: Merging shapefiles with same attributes using QGIS?

texts:
I have two shapefiles, a.shp and b.shp, containing point features with the following attribute tables:

a
-
id    att1   att2
-----------------
1     1      1
2     2      2

b
-
id    att1   att2
-----------------
3     1      1
4     2      1

If I use Vector->Geoprocessing Tools->Union to create a new shapefile ab.shp I end up having an attribute table like this:

id     att1    att2    id_2    att1_2   att2_2
----------------------------------------------
1      1       1       null    null     null
2      2       2       null    null     null
null   null    null    3       1        1
null   null    null    4       2        1

But what I would like to have (and how it is in ArcGIS Desktop) is the following:

id   att1  att2
---------------
1    1     1
2    2     2
3    1     1
4    2     1

The fields are all of the same type (QString). So what would be the right way to produce the desired output?


texts: How to edit attribute with my result from the count

texts:
I need this to loop through all the OBJECTID for the Fitting layer. Then it get the count for every line feature touching it. Now I'm at the point where I need the attribute table to be edited with that count. I have calculate field and I know arcpy.CalculateField_management is wrong and need to change I just don't know what it should change to.

Fitting = r'P:\YU_Restricted\Yuma_Admin\Projects\FY18\WO18_100_199\18-190_Gas_Utilities_Update_(Plicka_Joffre)\Gas_Utilites.gdb\UTILITIES_GAS\GasUtilNode_gControlFitting'
Pipe = r'P:\YU_Restricted\Yuma_Admin\Projects\FY18\WO18_100_199\18-190_Gas_Utilities_Update_(Plicka_Joffre)\Service\Service_Pipe.gdb\GasUtilNode_Segment'

         # Process Selection Loop
 with arcpy.da.SearchCursor(Fitting, ["OID@"]) as cursor:
     for row in cursor:
         # Process: Make Feature Layer
         arcpy.MakeFeatureLayer_management(Fitting, "FittingLyr", "OBJECTID = " + str(row[0]))

         # Process: Select Layer By Location
         arcpy.SelectLayerByLocation_management('GasUtilNode_Segment', 'BOUNDARY_TOUCHES', "FittingLyr", '#', 'NEW_SELECTION', 'NOT_INVERT')

         # Process: Copy Features
         arcpy.CopyFeatures_management('GasUtilNode_Segment', r'P:\\YU_Restricted\\Yuma_Admin\\Projects\\FY18\\WO18_100_199\\18-190_Gas_Utilities_Update_(Plicka_Joffre)\\results', '#', "0", "0", "0")

         # Process: Get Count
         arcpy.GetCount_management("results")

         # Local variables:
         lyrfile = "results"
         result = arcpy.GetCount_management(lyrfile)
         count = int(result.getOutput(0))
         calcuFld = count

         # Process: Delete
         arcpy.CalculateField_management(Fitting,"mediaID",calcuFld)
         arcpy.DeleteFeatures_management(r'P:\\YU_Restricted\\Yuma_Admin\\Projects\\FY18\\WO18_100_199\\18-190_Gas_Utilities_Update_(Plicka_Joffre)\\results')
         arcpy.DeleteFeatures_management("FittingLyr")


texts: send a GetFeatureInfo to Geoserver over a time range

texts: I want query geoserver about FeatureInfo from a set of images temporally spaced. How I should build this route?


texts: Why is a 3D marker symbol blacked out in the legend when saved to a .pdf file?

texts: I must have been needing coffee at the time I posted that, but I am using ArcMap 10.8.1 and I was referring to 10.6 for the earlier episode.  Now, does anybody have any ideas for a solution?


texts: How to update new column with v.distance from GRASS GIS?

texts:
I would like to calculate the distance between two lines and update this value to the attribute table of one of the line layers.So the mean distance between the blue lines and each of the red line segments in the following image:

I have tried solving this with v.distance from GRASS. I have created the two following columns in the attribute table of my red lines:

and run the following command:

This creates a new vector layer containing my red lines without any quantitative information in the attribute table.

Can anyone tell me what I am doing wrong?


texts:
I am trying to export my data to CSV files on Google Drive using this action:

Export.table.toDrive({
  collection: ndviFeatures,
  description: 'ndvi',
  selectors: ['ndvi'],
  folder: 'TRYING'
})

This action creates a task which I need to actively activate.
Is there a way for this thing to happen automatically?
I wish to run a for loop on several dates and change some parameters and have it all downloaded and separated into different files.

Any ideas?


texts: Reading KML Simple Data to a popup in Open Layers

texts:
I have created a pop up for markers read from a KML file, however The KML file data is not being read and I assume due to it being in a <SimpleData> tag and in a <name> tag
How would I be able to read this to get the name of the listed buildings from the file as well as the Grade and Broad Class from this file to the PopContent Variable
Pop up Script
Omap.on('click', function (evt){
    var feature = Omap.forEachFeatureAtPixel(evt.pixel,
        function (feature,layer) {
            return feature;
            });
        if (feature) {
            var Popgeom = feature.getGeometry();
            var Popcoord = Popgeom.getCoordinates();

            var PopContent = '
' + feature.get('name') + '
';
            PopContent += '
' + feature.get('Grade') + '
' + '
The listed property is a ' + feature.get('BroadClass') + ' building
';

            content.innerHTML = PopContent;
            overlay.setPosition(Popcoord);

            console.info(feature.getProperties());
        }
    });

An extract from the KML file:
<Placemark>
    <name>Kelsterton Hall</name>
    <ExtendedData><SchemaData schemaUrl="#LSB">
        <SimpleData name="id">Cadw_ListedBuildings.fid--3eda8af5_1772a727fac_-6786</SimpleData>
        <SimpleData name="RecordNumber">1</SimpleData>
        <SimpleData name="DesignationDate">1949/08/08 00:00:00</SimpleData>
        <SimpleData name="Grade">II</SimpleData>
        <SimpleData name="Location">Set back slightly from the road at the junction of Kelsterton Road and Kelsterton Lane.</SimpleData>
        <SimpleData name="BroadClass">Recreational</SimpleData>
        <SimpleData name="BroadClass_cy">Hamdden</SimpleData>
        <SimpleData name="Report">http://cadwpublic-api.azurewebsites.net/reports/listedbuilding/FullReport?lang=en&amp;id=1</SimpleData>
    </SchemaData></ExtendedData>
      <Point><coordinates>-3.08246357229187,53.2285595916199</coordinates></Point>
  </Placemark>


texts: Publishing PostGIS layer with EPSG:5514 not working in GeoServer of version 2.16.2

texts:
I can not using EPSG: 5514 for publishing layers. Data was stored in PostGIS.

GeoServer returns this error:

  org.apache.wicket.WicketRuntimeException: Method onRequest of
  interface org.apache.wicket.behavior.IBehaviorListener targeted at
  org.apache.wicket.ajax.markup.html.form.AjaxSubmitLink$1@6237d421 on
  component [AjaxSubmitLink [Component id = computeNative]] threw an
  exception

  Caused by: java.lang.reflect.InvocationTargetException    at
  sun.reflect.GeneratedMethodAccessor1324.invoke(Unknown Source)    at
  sun.reflect.DelegatingMethodAccessorImpl.invoke(DelegatingMethodAccessorImpl.java:43)
    at java.lang.reflect.Method.invoke(Method.java:498)     at
  org.apache.wicket.RequestListenerInterface.internalInvoke(RequestListenerInterface.java:258)

  Caused by: java.lang.NoClassDefFoundError:
  org/geotools/filter/function/FilterFunction_arrayAnyMatch     at
  org.geotools.data.postgis.FilterToSqlHelper.createFilterCapabilities(FilterToSqlHelper.java:186)
    at
  org.geotools.data.postgis.PostgisFilterToSQL.createFilterCapabilities(PostgisFilterToSQL.java:100)
    at
  org.geotools.data.jdbc.FilterToSQL.getCapabilities(FilterToSQL.java:431)
    at
  org.geotools.jdbc.JDBCDataStore.getFilterCapabilities(JDBCDataStore.java:541)
    at org.geotools.jdbc.SQLDialect.splitFilter(SQLDialect.java:1281)   at
  org.geotools.jdbc.JDBCFeatureSource.splitFilter(JDBCFeatureSource.java:378)
    at
  org.geotools.jdbc.JDBCFeatureSource.splitFilter(JDBCFeatureSource.java:365)
    at
  org.geotools.jdbc.JDBCFeatureSource.getBoundsInternal(JDBCFeatureSource.java:464)
    at
  org.geotools.jdbc.JDBCFeatureStore.getBoundsInternal(JDBCFeatureStore.java:181)
    at
  org.geotools.data.store.ContentFeatureSource.getBounds(ContentFeatureSource.java:413)

How can I resolve this problem?

Geoserver logs contains next error:

  Illegal access: this web application instance has been stopped
  already. Could not load

Table description:

I tried other table with EPSG: 3857 and problem is same.


texts:
Can't get advanced python label expression to work.  I am trying to get my if else expression to either concatenate or not based on if two fields are blank or not.

If SUB_UNIT_TYPE and SUB_UNIT_ID are not null then concatenate them together on a new line with a space between the values from the two fields.  I can't get python to find the null value.

I tried all the suggestions from this post.  Using if/elif for passing null values to Python of label expression in ArcGIS for Desktop?

def FindLabel ( [Add_Num], [PROP_UNIT_TYPE], [PROP_UNIT_ID],[Sub_Unit_Type],[Sub_Unit_ID]):
   if (not([SUB_UNIT_TYPE] is None) and not([SUB_UNIT_ID] is None)):
       L = ([Add_Num] +  '\n' + [Sub_Unit_Type] + [Sub_Unit_ID])
   else:
      L = ([Add_Num])
   return L

Tried this.  Still need to add a space between the four subaddress fields and write an elif statement for if the PROP_UNIT_TYPE and PROP_UNIT_ID is null but the SUB_UNIT_TYPE and SUB_UNIT_ID fields are not null.  Is this really the best way to go?  Its easier in C# and SQL Server but python is pretty much the only language accepted in labeling expressions.

def FindLabel ( [Add_Num], [PROP_UNIT_TYPE], [PROP_UNIT_ID], [Sub_Unit_Type], [Sub_Unit_ID]):
   if ([PROP_UNIT_TYPE] is None and [PROP_UNIT_ID] is None and [SUB_UNIT_TYPE] is None and [SUB_UNIT_ID] is None):
      L = ([Add_Num])
   else:
      L = '{0}\n{1}{2}\n{3}{4}'.format([Add_Num],[PROP_UNIT_TYPE],[PROP_UNIT_ID],[Sub_Unit_Type],[Sub_Unit_ID])
   return L


texts: Overpass Query to find ways of a route geometry

texts:
Problem
I have a route geometry in the format:

OSM node-1 coordinates
OSM node-2 coordinates
OSM node-3 coordinates
...

And I need to know which OSM ways were "used" by the route.
Question
Is there an Overpass query which can return me the ways which are required to "form" the route geometry with one request?

fyi: The geometries' route is always "connected" by shared way nodes.

e.g.:

input: geometry = [node102, node103, node104, node303, node304]
output: wayIds = [way1, way2, way3]
where way1 is e.g. [node101, node102, node103, node104, node105]
where way2 is e.g. [node104, node303]
where way3 is e.g. [node303, node304, node305]

Context
I'm using a routing API which only returns the Geometry of the Route (OSM node coordinates) but no ids.
The problem is that I need to know which OSM ways were used to construct this route.
What I've tried so far

I can get all the node id from the node coordinates using node(around:0,lat,lon)

Thanks to this I can get the ways which contain a set of nodes using

node(around:0,latNode1,lonNode1,latNode2,lonNode2)->.allnodes;

way(bn.allnodes)->.ways;

foreach .ways -> .singleway (
  node.allnodes(w.singleway);
  way.singleway(bn)(if:count(nodes) == allnodes.count(nodes));
  out geom;
);

Unfortunately I don't know how many nodes of a way were used by the route. Right now I only can think of:
a. find all ways for node1, node2
b. if more than one way found --> add another node to (a) and try again
c. when I found the way, search the next way with the geometries next node pair which is not part of the way I just found (i.e. continue with (a))
But this approach would require me to make many request to the Overpass API. I can't figure out how to write this logic as Overpass Query. I'm running my own overpass server.


texts: Problem with stretched value and classified value in ARCGIS desktop? They are different than each other

texts:
I want to classify the value of an image into 3 groups. when I check properties of a raster image, stretched value is different than classified value. Here I have uploaded the images as follow.

I have no idea where this problem comes. Why they are different? How to fix them?


texts: removeAllFeatures + destroyFeatures too much?

texts:
I'm doing the following in my code and am wondering if it's too much or is there any benefit?

        myLayer.removeAllFeatures();
        myLayer.destroyFeatures();

Reference:

http://dev.openlayers.org/releases/OpenLayers-2.10/doc/apidocs/files/OpenLayers/Layer/Vector-js.html#OpenLayers.Layer.Vector.removeAllFeatures


texts: Why does Save as "Rendered Image" change the raster styling?

texts: Rasters exported from QGIS (via Save As > Rendered Image) do not retain original styling.  See attached.  Any thought on how to prevent this error?  The top image is exported from QGIS as a .bmp, the bottom is a screen grab from tilemill.  Thanks.


texts: Is there a way to create raster by interpolating point values using GeoTools 8.2?

texts:
im new in using GeoTools (8.2). I was able to create shp from asc file which contains my data. Now my question is is there a simple way to create (interpolate) a simple shapefile (only points) to raster image?

Thx in advance for all your answers.

Best regards


texts: Clipping OpenStreetMap basemap using QGIS

texts:
I'm using QGIS, and I am trying to clip my metro region data to only show the OSM basemap in the the metro region, using this thread:
Clipping OSM Standard Basemap using QGIS
I created an inverted polgygon shapefile to mask the metro area, but when it does so it also masks the OSM layer beneath it.  The screenshots are from toggling this layer on and off.
How would I get the metro region to appear, along with the base map beneath it?
Do I need to create a separate shapefile of the OSM basemap I need in order to be able to work with it?


texts: Include ID in an spatial weight matrix R

texts:
I have a data frame like this one:
dataSp <- read.table( text = '
  ID   LATITUDE    LONGITUDE
  A     -85           134
  B      34             2
  C      42             3
  D      45             5
  E      -2            80
  F      -5            79',
  header = TRUE )

My main objective is to generate a spatial weight matrix
This is my code so far:
data_sf <- st_as_sf(dataSp, coords = c("LONGITUDE","LATITUDE"), crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
st_is_longlat(data_sf)
coords <- st_coordinates(data_sf)
col.rel.nb <- graph2nb(relativeneigh(coords), sym=TRUE)
listaw <- spdep::nb2listw(col.rel.nb, style="W")

The problem is that listaw doesn't include the information in ID. How can I identify each neighbour with the ID (i.e: A,B,C,D,E,F)?
PS: This is a repost from a general stackoverflow question (not solve, so far).


texts: Popup for each feature when using Leaflet marker.cluster plugin

texts:
I have leaflet map to display locations of points. Firstly, I use leaflet.ajax plugin to load the points data in GeoJSON format, then I want to enable cluster view using the Leaflet.markercluster plugin.
The cluster view works fine, but the popup only shows "siteid" of the last site in my dataset.
Below is my code, can anyone advise on it?
function map_viewer(map, options){

        var my_data = new L.GeoJSON.AJAX("http://127.0.0.1:8000/my_data/",{
            onEachFeature: function(feature,layer){
            clusters.on('click', function (e) {
            this.bindPopup(feature.properties.siteid);
            });
            }

        });

        var clusters = L.markerClusterGroup();
        my_data.on('data:loaded', function()
        {
        clusters.addLayer(my_data);
        map.addLayer(clusters);
        });

        var groupedOverlays = {
          "Layers": {
            "cluster view":  clusters
          }
        };

        L.control.groupedLayers(groupedOverlays).addTo(map);
    }

Updated on 2021-02-01:
I have modified my code but the popup window completely disappear.
function map_viewer(map, options){
            var clusters = L.markerClusterGroup();
            var my_data = new L.GeoJSON.AJAX("http://127.0.0.1:8000/my_data/",{
                onEachFeature: function(feature,layer){
                layer.bindPopup(feature.properties.siteid);
                }

            });

            my_data.on('data:loaded', function()
            {
            clusters.addLayer(my_data);
            map.addLayer(clusters);
            });

        }

Below is a small part of the importing geojson dataset.
{"type": "FeatureCollection", "crs":
{"type": "name", "properties": {"name": "EPSG:4326"}}, "features": [
{"type": "Feature", "properties": {"siteid": 1, "latitude": -28.004959, "longitude": 153.428117, "pk": "1"}, "geometry": {"type": "MultiPoint", "coordinates": [[153.428117, -28.004959]]}},
{"type": "Feature", "properties": {"siteid": 2, "latitude": -33.870436, "longitude": 151.225013, "pk": "2"}, "geometry": {"type": "MultiPoint", "coordinates": [[151.225013, -33.870436]]}},
{"type": "Feature", "properties": {"siteid": 3, "latitude": -33.92677, "longitude": 151.21356, "pk": "3"}, "geometry": {"type": "MultiPoint", "coordinates": [[151.21356, -33.92677]]}},
{"type": "Feature", "properties": {"siteid": 4, "latitude": -33.854711, "longitude": 150.987657, "pk": "4"}, "geometry": {"type": "MultiPoint", "coordinates": [[150.987657, -33.854711]]}},

Conclusion:
I have fixed the issue just by converting the geometry type in my geojson dataset from "Multipoint" to "Point". It seems this plugin Leaflet.markercluster can only cluster view for Multipoints, but not able to display the bindPopup of each layer.


texts: Is it possible to serve Complex GML Features with Mapserver?

texts:
I want to serve complex GML-Features through WFS 2.0.0 with MapServer. While support for WFS 2.0 was added to Mapserver, I cannot find any current information whether it is possible to serve complex features in Mapserver. Mapserver-Documentation under http://mapserver.org/ogc/inspire_dl.html only talks about metadata.

I know other map servers can handle Complex Features (e.g. GeoServer through its appschema-extension), but is it currently possible in Mapserver (or will it be in the near future)?


texts:
I have created a land classification using QGIS in Geotiff format and need to convert this to Mapinfo format.

I tried to used the conversion tool file type (Northwood Classified Grid Format .grc /.tab) but the output file is .grc and .aux.

Please advise on how to covert this files to a Mapinfo *.grc *.tab files in QGIS.

Thanks,


texts: I have created an elevation contour (topography) with TIN using a spreadsheet of coordinates and elevations of some points. now I want to find the elevation of a new point with these contours. how can I do that ?


texts:
I have two layers in ArcGIS:

Highway layer (9000 data points (or values))
Safe zone (4 types of values i.e. safe zone, medium zone, bad zone, and worst zone).

Now I want to do analysis to see how many highways data points fall into zones and how many do not fall in any zone.

I tried to use the Intersect tool to get the result but it didn't produce the desired result. It only gave those points that lie in a safe zone, it didn't return the points that do not lie in a safe zone.

What function should I use?


texts: How is ST_Azimuth (in PostGIS) calculated?

texts:
What algorithms are internally calculated in PostGIS?
Can you explain it with simple code?
(I want to know the standard of true north.)


texts: Checking whether layer is vector or not in OpenLayers 2?

texts: How do I dynamically check whether a layer object belongs to the OpenLayers.Layer.Vector class?


texts: Openlayers example not working with coordinates

texts:
Hey guys i was just playing around with openlayers 3, and i created THIS fiddle, take a look at the javascript initilization code below:

var styles = [
    'Road',
    'Aerial',
    'AerialWithLabels',
    'collinsBart',
    'ordnanceSurvey'];
var layers = [];
var i, ii;
for (i = 0, ii = styles.length; i < ii; ++i) {
    layers.push(new ol.layer.Tile({
        visible: false,
        preload: Infinity,
        extent: [77.6013, 12.9044, 77.6885, 12.9536],
        source: new ol.source.BingMaps({
            key: 'Ak-dzM4wZjSqTlzveKz5u0d4IQ4bRzVI309GxmkgSVr1ewS6iPSrOvOKhA-CJlm3',
            imagerySet: styles[i]
            // use maxZoom 19 to see stretched tiles instead of the BingMaps
            // "no photos at this zoom level" tiles
            // maxZoom: 19
        })
    }));
}
var map = new ol.Map({
    layers: layers,
    // Improve user experience by loading tiles while dragging/zooming. Will make
    // zooming choppy on mobile or slow devices.
    loadTilesWhileInteracting: true,
    target: 'map',
    view: new ol.View({
        center: [12.9370699, 77.6244159],
        maxZoom: 20,
        minZoom: 13,
        zoom: 13
    })
});

var select = document.getElementById('layer-select');

function onChange() {
    var style = select.value;
    for (var i = 0, ii = layers.length; i < ii; ++i) {

        console.log(layers[i].get('source').get('') + " " + styles[i]);

        layers[i].setVisible(styles[i] === style);
    }
}

select.addEventListener('change', onChange);
onChange();

I took this code from an online fiddle, all i wanted to do was change the co-ordinates , for the center and extent , so i only changed the below code:

center: [12.9370699, 77.6244159],

and extent

extent: [77.6013, 12.9044,77.6885,12.9536],

Now for the center i got the coordinates from google maps and for the extent i got the co-ordinates from the openlayers maps. Now the problem is the map just does't display , Why ??


texts: System policies prevent QGIS installation - false alarm?

texts:
I get the following message during install Quantum GIS 1.8.0-Lisboa on Windows 7:

I tried doing the install as administrator with the same result. The install finishes after the 'ok' button is clicked and QGIS seems to be functioning after the install.

What is causing the installer popup and how do I fix it? Or should I leave as it is?


texts:
I am using this code to try fill a polygon with a regular grid of points, however it returns an empty point layer that is not visible on the map. i don't think it is interacting with the polygon layer i am using.

New to python and qgis

layer = iface.activeLayer()

feats = [ feat for feat in layer.getFeatures()]

points = []

spacing_y = 50
spacing_x = 50

for feat in feats:
    extent = feat.geometry().boundingBox()
    xmin, ymin, xmax, ymax = extent.toRectF().getCoords()

    rows = int(ymax - ymin)/spacing_y
    cols = int(xmax - xmin)/spacing_x

    x = xmin
    y = ymax

    geom_feat = feat.geometry()

    for i in range(rows+1):
        for j in range(cols+1):
            pt = QgsPoint(x,y)
            tmp_pt = QgsGeometry.fromPoint(pt)
            if tmp_pt.within(geom_feat):
                points.append(tmp_pt.asPoint())
            x += spacing_x
        x = xmin
        y -= spacing_y

epsg = layer.crs().postgisSrid()

#points
uri = "Point?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'point',
                           'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(points)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromPoint(points[i]))

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)


texts: OpenLayers3, WFS ajax request doesn't work

texts:
I'm creating a geoportal, which works only on my localhost. I have a problem with adding to my ol.map a vector layer from Geoserver using WFS. I'm using Geoserver v2.7, OL v3.15.1 and jQuery v2.2.4. I had a cross-domain problem, since my gesoerver is hosted at 8080 port and the app is displayed using 9080 port, so I decided to add this layer using JSONP. I added to geoserver\WEB-INF\web.xml the following code:

 <context-param>
 <param-name>ENABLE_JSONP</param-name>
 <param-value>true</param-value>
 </context-param>

Here's my code for adding a layer through AJAX request:

   var vectorLoader = function() {
   $.ajax('http://localhost:8080/geoserver/ows',{
     type: 'GET',
     data: {
      service: 'WFS',
      version: '2.0.0',
      request: 'GetFeature',
      typename: 'brzesko:wydarzenia',
      srsname: 'EPSG:3857',
      outputFormat: 'application/json',
      bbox: '2293539.9663390242,6439346.050046678,2294549.5570354667,6442264.769222975'
     },
     dataType: 'jsonp',
     jsonpCallback: 'loadFeatures',
     jsonp:'format_options'
     });
   };
  var geoJSON = new ol.format.GeoJSON();
  var loadFeatures = function(response) {
   vectorSource.addFeatures(geoJSON.readFeatures(response));
  };
  var vectorSource = new ol.source.Vector({
    format: geoJSON,
    loader: vectorLoader
  });

  var wydarzenia= new ol.layer.Vector({
    source: vectorSource (...)

(can't use outputFormat:'text/javascript', my WFS doesn't support it...)
Then I get an error in the console.log:

  Uncaught SyntaxError: Unexpected token : http://localhost:8080/geoserver/ows?(...)

However, when I open the link in new tab I get a working JSON with all features and information. Does anybody know how to fix this? I've seen several similar posts but none of the given solutions works for me.


texts: How to calculate the accuracy of georeferenced aerial images in meters if the mean error is in pixels?

texts: I have scanned aerial images taken shortly after WWII and have to georeference them and calculate the mean error. I use the Projective transformation so the error after transformation is shown in pixels. I want to know how big is the error in meters. When I go to the properties of the transformed image there is field showing the size of the pixels (see picture). I want to know how reliable this value is and if for example my error is 20 pixels, can I just multiply the error by the given pixel size?


texts: How do I search all primary, secondary, and tertiary street intersections via Overpass?

texts:
I'm trying to identify the coordinates of all of the intersections where a primary, secondary, or tertiary highway meets another primary, secondary, or tertiary highway. I searched in the API docs here but haven't been able to solve it yet. Here's what I have so far but clearly it returns all nodes associated with a primary secondary or tertiary road (I realize the second_road is redundant).

[bbox:{{bbox}}];
way[highway~"(primary|secondary|tertiary)"]->.first_road;
way[highway~"(primary|secondary|tertiary)"]->.second_road;
node(w.first_road)(w.second_road);
out;

I could take this data and reduce it on my end via javascript to find if a node shows up twice or more but I was wondering if there is a better way to do this via the Overpass query language?


texts:
I've recently pulled some of the TIGER data onto my local machine but, I'm at a loss about what it all means. I'm hoping to find some technical documentation on what all the tables and columns mean.

For example, what does the column tlid in the table addrfeat represent?

Specifically, I'm trying to pull up address information by area, i.e, what are all the possible addresses in downtown Albany, NY?

I have it sort of working but I don't feel confident about the results since I just pulled data without regard to the meaning of the columns or tables.


texts: OpenLayers 4 not connecting to GeoServer WMS

texts:
Just started learning OpenLayers 4 and I'm trying to display a simple WMS layer from GeoServer as a web map. I'm not having any luck. If I use a request for OpenStreetMap layers (ie. ol.source.TileWMS) it works, but using ol.source.TileWMS to connect to my GeoServer layers doesn't seem to work. Teh web page is blank and there's not even a box or zoom function. I'm not sure if its my code or if I need to configure GeoServer somehow. Its on a local host, same machine.

The layer properties are:

workspace: basic, layer:roads, EPGS:28356, GeoServer URI for the workspace (not sure if this matters):http://localhost:8080/geoserver/basic

The code is:

<!DOCTYPE html>

    <title>roads</title>
    <link rel="stylesheet" href="https://openlayers.org/en/v4.0.1/css/ol.css" type="text/css">

Some more information; the layer preview URL shows http://localhost:8080/geoserver/basic/wms?.....: when http://localhost:8080/geoserver/basic/wms?request=GetCapabilities is entered there is a long screen of XML: and when the developer console is used to check for JavaScript errors, the line source: new ol.source.TileWMS({ is highlighted. If I delete the extents line out of the code, that error is gone (but the map still doesn't work).


texts:
I would like to build a cartogram map using Web GIS.

I have a basemap and a set of raw data in excel.

I wondered what are the steps to build up a cartogram map using Web GIS?

Note: Originally this Question solicited Answers for both Web GIS and ArcGIS for Desktop.  A new Question has been asked focussed on ArcGIS for Desktop (Creating contiguous pixel based cartogram using ArcGIS for Desktop?) and so this one has now been focussed on Web GIS (which itself is still quite broad).


texts:
I'm trying to merge large raster files (DEM's) together but I run into memory errors using the 'Raster Merge' tool in QGIS 3.8.

I suspect the issue is mainly the high amount of precision in the float32 data format for each pixel. Therefore if I could translate/convert the format into half/float16 format I should be able to do this.

Unfortunately I don't see an option to export GeoTIFF's into this format in the 'Raster Translate' tool or otherwise.

Edit: I need a few decimal places of precision - ideally 3. The rasters are DEM's that are accurate to about 8-10cm, so any decimals past 3 are nonsense and just eat up memory.


texts:
I have recently modified QFIELD into a fit for purpose mobile mapping solution for rural address assignment initiative. Now I need to change the QFIELD icon to reflect our organizational identity and to host the environment internally and locally for confidentiality purposes.

I had successfully downloaded the code on github but Android Studio cannot compile that code. It is giving some errors messages referring to the code itself. I also tried to import the code as an eclipse project and that has failed too.

How can I get it right?


texts: I am looking for a Python package that would allow me to input a UTM zone code and output different geo_referencing strings (for example, WKT or Proj4). Are there any lightweight, widely used packages out there that perform this function?


texts: CORINE land cover 2000: coordinate reference system

texts:
The PDF in this link states that the CORINE land cover 2000 used ETSR89 as its coordinate reference system: http://www.eea.europa.eu/publications/technical_report_2006_7

I wanted to check I am reading this correctly: did CORINE 2000 use ETSR89 as its coordinate reference system?

In addition, what would be an appropriate coordinate reference system to re-project CORINE 2000 over the UK?


texts: QGIS 2.14. minidump crash when closing program

texts:
Every time I close QGIS, I get an error message like this:

  minidump written to C:\Users\name\AppData\Local\Temp\qgis-20160412-100443-7416-5376-ea85bef.dmp

This happened also in the former versions, and reappears when installing the new one (i.e. delete the complete OsGEO4W64 file and new installation). It still remains when I deactivate all plugins. What can be the reason?


texts: new layer from calculations with fields of other layers in QGIS

texts:
I've two layers, both with the same fields: A, V, where A is the shape's area and V a value given by me. I want a new layer with shape created by the intersection of the other layers' shapes and a value given by the calculation: (A1*V1+A2*V2)/(A1+A2).
Is there a command to make it automatically? Obviously in reality I have many layers with many shapes, otherwise I made it manually.

I used intersection to subdivide all shapes from all layers; in this way for every shapes I have 2 weights (V values). Now I can use vector-field calculator and calculate:($area*("x"+"y"))/(sum(x+y), where x and y are the 2 weights. I hoped there was a command that make it automatically… It's possible at least calculate summation of x and y values in vector-field calculator?


texts: Find in between coordinates for the given latitude and longitude from a CSV file

texts:
I have a CSV file which has the latitude and longitude values for the start and the end coordinates as follows. In these the first lon and lat values are the start point then the second coordinates are the end points for one value. Then the third point is the start and the forth is the end and so on.
The csv is in the following format

lat    lon
1.25   6.45
2.45   5.36
3.67   9.67
3.56   6.25

The code which I tried is:

def main(interval,azimuth,lat1,lng1,lat2,lng2):
    '''returns every coordinate pair inbetween two coordinate
    pairs given the desired interval'''

    d = getPathLength(lat1,lng1,lat2,lng2)
    remainder, dist = math.modf((d / interval))
    counter = float(interval)
    coords = []
    coords.append([lat1,lng1])
    for distance in range(0,int(dist)):
        coord = getDestinationLatLong(lat1,lng1,azimuth,counter)
        counter = counter + float(interval)
        coords.append(coord)
    coords.append([lat2,lng2])
    return coords

if __name__ == "__main__":
#point interval in meters
interval = 10
#direction of line in degrees
#start point

lat1 =  18.750000000
lng1 =     85.416666667
#end point
lat2 =  17.650000000
lng2 =  83.250000000
azimuth = calculateBearing(lat1,lng1,lat2,lng2)
#print(azimuth)
azimuth
coords = main(interval,azimuth,lat1,lng1,lat2,lng2)
#print(coords)
coords

This code takes the coordinates manually but I need it to take values from CSV and generate the output as CSV.


texts:
Converting a .mxd to .dxf is something we need to do, and we do not only need the layers but also the layout of a map on our .mxd to be transfered to .dxf. I've done some research on the internet but it seems like an very difficult feat.

In the ArcGIS Discussion Forums, some people mention the possiblity of going from .mxd to .ai (Adoble Illustrator) followed by a conversion from .ai to .dxf, have any one tried it ?

Others indicate that they are working with .mxd to .pdf then .pdf to .dxf, I gave this a try but the layout is not well converted.

Is there any other way we might proceed to have better results ?


texts: Get the list of map services where table "xyz" is used

texts:
I have an SDE feature class:

And I'm looking to find, to how many map services it is pointing.

There are like duplicate map services pointing to same table, so I was trying to figure out how many duplicate services we have.

Is there any script/method to start with?


texts:
I am trying to import CAD files to ArcGIS Pro. I set my basemap and used the add data button to add polylines, polygons, points and annotation to ArcGIS Pro from a CAD drawing.  It imported all the data, but it puts it in the middle of the Atlantic ocean instead of NJ where it belongs.

Any ideas how I can get the drawing to project where it is supposed to be?


texts:
Windows 10, PostgreSQL 11

I successfully managed to export a .tif raster into PostgreSQL using raster2pgsql:

raster2pgsql -s 27700 -I -C -M -Y "D:\Data_Location\country.tif" -F -t 50x50 public.country_rast | psql -U username -d dbname -h hostaddress -p 5432

Now I'd like to query the raster using a function which takes an XY coordinate and returns the value of the raster cell it intersects.

I use this:

SELECT ST_VALUE (rast,1,(ST_GeomFromText('POINT(319006 212230)',27700)))
FROM public.country_rast

However the query returns ALL 1 million rows of the raster, wheras I only want the value that the given XY intersects. How can I do this?


texts: Add drop-down list in parameter input using script editor QGIS 2.18

texts:
In the script editor in QGIS, when choosing string parameter type, a default value must be added according to https://docs.qgis.org/2.8/en/docs/user_manual/processing/scripts.html

What happen is that the parameter line in the GUI will be populated with the default value. 'average' in my case.

What I want to do is to add several different values so that the user can choose from. In my case, I want to add other statistics, so the user will have a drop-down list containing several options.

The solution in Letting user choose from 2 options for QGIS Processing script input parameter?
does not work for me since I'm using QGIS 2.18.2, also the solution for check-box is not relevant.

Is it possible to create a drop-down list for the user in QGIS 2.18.2 las palmas?

In case it is needed, the full code is here:

##Cell Statistics=name
##Select_directory=Folder
##statistic_name=string average
##Output_raster=Output raster
import glob, os
from PyQt4.QtCore import QFileInfo
from qgis.core import QgsRasterLayer, QgsRectangle

os.chdir(Select_directory)
raster_list = []
extent = QgsRectangle()
extent.setMinimal()

statistic_dic = {'average':0, 'count':1, 'median':2,'mode':3, 'minimum':4, 'min_raster':5,'maximum':6, 'max_raster':7, 'std':8, 'range':9}

for raster in glob.glob("*.tif"):
    fileInfo = QFileInfo(raster)
    baseName = fileInfo.baseName()
    rlayer = QgsRasterLayer(raster, baseName)
    # Combine raster layers to list
    raster_list .append(rlayer)
    # Combine raster extents
    extent.combineExtentWith(rlayer.extent())

# Get extent
xmin = extent.xMinimum()
xmax = extent.xMaximum()
ymin = extent.yMinimum()
ymax = extent.yMaximum()
# Run algorithm and set relevant parameters
processing.runalg("grass7:r.series",
                {"input":raster_list ,
                "-n":False,
                "method":statistic_dic[statistic_name],
                "range":'-10000000000,10000000000',
                "GRASS_REGION_PARAMETER":"%f,%f,%f,%f"% (xmin, xmax, ymin, ymax),
                "GRASS_REGION_CELLSIZE_PARAMETER":0,
                "output":Output_raster})


texts: pgr_dijkstra ERROR: invalid memory alloc request size 1080000000

texts:
I am using pgr_version 2.6.3

I tried googling, most of the posts suggest the data is corrupted. So I tried pg_dump and select row by row to try finding the corrupted row, but did not encounter any problem. Some posts also suggest the output row size is too large and suggest to reduce the query size, but I don't see how I can achieve this. The error is identical with every try even changing the source and target id.

I setup the database using osm2po with planet-latest.osm. I have loaded the .sql to my table pl_2po_4pgr, but I cannot get routing to work.

My query is

select * from pgr_dijkstra('select id, source, target, cost from pl_2po_4pgr', 110491754, 110490812, false);

and the error shows:

ERROR:  invalid memory alloc request size 1080000000
CONTEXT:  PL/pgSQL function pgr_dijkstra(text,integer,integer,boolean,boolean) line 17 at RETURN QUERY

I also tried query as https://live.osgeo.org/en/quickstart/pgrouting_quickstart.html suggests, but the error still exist

SELECT seq, node, edge, cost
        FROM pgr_dijkstra('
                SELECT id, source, target,
                        length as cost FROM pl_2po_4pgr',
                100, 600, false
        );

I also went on http://localhost:8888/Osm2poService  and was also be route prefectly within 1 second. The source and target id is right next to each other within 100 meters.

How do I fix this problem? or is there a workaround? or maybe there is a way to output by batch?

I also saw this pgr_dijkstra invalid memory request size but doesn't seem like there is a solid solution

Edit

I tried limiting by using st_buffer and st_dwithin, but neither of them works, something like:

SELECT *
        FROM pgr_dijkstra('
                SELECT osm_id as id, osm_source_id as source, osm_target_id as target,length as cost FROM pl_2po_4pgr
                WHERE geom && st_buffer(st_makeline(
st_setsrid(st_makepoint(-122.39212 46.95), 4326),
st_setsrid(st_makepoint(-122.360618 47.186731), 4326)
                ), 100
             )',
                110491754, 110490812, false
        );

SELECT *
        FROM pgr_dijkstra('
                SELECT osm_id as id, osm_source_id as source, osm_target_id as target,length as cost FROM pl_2po_4pgr
                WHERE st_dwithin(geom,
st_setsrid(st_makepoint(-122.39212 46.95), 4326)
                ), 1000
             )',
             110491754, 110490812, false
        );


texts: Adding new vector layer. Should .shp file only be added or .shp plus corresponding .shx, .dbf, etc.?

texts:
I am a very definite QGIS newbie.  I have been creating a map using files that had been used in ArcGis.  Each time I wanted to add a new vector layer I opened an .shp file and specified CRS as WGS 84.  The resulting layers look fine but when I try to use the "measuring tool" to measure line length I get an error message as follows:

This map is defined with a geographic coordinate system (latitude/longitude) but the map extents suggests that it is actually a projected coordinate system (e.g., Mercator). If so, the results from line or area measurements will be incorrect.
To fix this, explicitly set an appropriate map coordinate system using the Settings:Project Properties menu.

any advice on how I might determine what might be an "appropriate map coordinate system" ?
I am also wondering if I should have also opened the corresponding .shx, .sbx, .dbf, .sbn and .qix files also ?

Thanks


texts: Why are color settings not saved when exporting to a new file?

texts: I extracted with QGIS a municipality from a state map, then divided it with a grid. I don't want the borders of the grid features to be seen so I changed the color of the map and the color of the borders to be the same, then I saved it as a new file (not a new project). Does anyone know why when I open it again it does not save the changes? (the divisions of the grid are visible and the color settings I saved are gone).


texts: ArcMap keeps loading old CSV file after I've updated it

texts:
I have a CSV I'd like to join to a shapefile. I brought the CSV into ArcMap, realized it had an error, fixed it and saved the file with the same name. I removed the CSV from my map document, refreshed the CSV's enclosing folder, and brought the CSV into the map document again. It is still the old CSV, with error. I verified that I had saved it correctly. I had. I closed Excel, and reopened the CSV: still correct. I 'reloaded' the CSV in ArcMap's table of contents. Still the old version.

I know I can resave the CSV with a new name to get around this issue but I want to figure it out. What am I doing wrong? How can I get ArcMap to recognize the updated CSV file?


texts:
I have a project where I need to use the OpenStates API to retrieve information. The code is written and works, I just need to know what the rate limiting is before I turn it on.

I have read through the documentation and did not find it, does anyone know what the OpenState API rate limits are?


texts:
I want to use a shapefile of the Corine Land Cover. This covers all of Europe, but I am only interested in the area of Scandinavia.

How can I cut this file using free software (preferably a command line tool)? I tried ogr2ogr

C:\Program Files\GDAL>ogr2ogr.exe -f "ESRI Shapefile" c:\Users\Christoph\geo\corine\clc06_c312_coniferous\clc06_c312_dalarna.shp c:\Users\Christoph\geo\corine\clc06_c312_coniferous\clc06_c312.shp -clipdst 0 60 20 70

but without any usable result (just almost empty files).


texts: Counting and rasterizing polygon overlaps in ArcGIS Desktop?

texts:
How do you take a vector layer with overlapping polygons and turn it into a raster where each cell counts the number of polygons containing that cell?

We are looking at polygons of fire burn areas between 1976 to 2000 and want to see how many fires would be contained in each raster cell.


texts:
I want to change tileSize from default 256 to 512. I am using WMS version 1.1.0 so i cannot add WIDTH and HEIGHT parameters directly to source.

There are some examples where custom tileGrid is created where we can specify tileSize. In my case it is working only if i know the resolutions

var tileGrid = new ol.tilegrid.TileGrid({
    origin:ol.extent.getTopLeft(projectionExtent),
    resolutions:resolutions,
    tileSize :512
});

How to handle case when i use default resolutions (i do not specify them in view). Is there a way to retrieve the resolutions for projection before creating map? Or is there a better way to handle this kind of situation?


texts: Avoid text/shields from being too close using Mapnik xml

texts:
I'm having the problem of text/shields from different layers appear too close to each other on the map generated from Mapnik xml.

The screenshot below shows the problem

I have three layers in this map.

A layer for shields, appeared as "rectangle 972" in the image
A layer for major roads, which has "SW 24th St" appear on the map
A layer for secondary roads

The problem is the "rectangle 972" and "SW 24th St" are too close to each other, and the street names for the secondary roads are also concatenated with each other (I'm less concerned with this though).

I couldn't find the right configuration option in Mapnik for this to look better. I have tried setting the minimum-distance in ShieldSymbolizer and TextSymbolizer, but I don't think that's the right way to do, and it didn't work of course.


texts: How to perform viewshed analysis for each pixel in raster?

texts:
I am having a problem making a viewshed analysis for every pixel on a raster.
To perform the viewshed analysis from every single pixel in a raster I will need:

A viewshed plugin (which I am willing to use either the one from GRASS or QGIS plugin for visibility analysis

This is what've done
 - STEP 1: Get the points layer from the raster

This is my raster:

What I've done was used Raster Pixels to Points tool:

and then I got this layer:

which I think it's was I need for step 2.

STEP 2: Perform Viewshed Analysis for each Point

As I said I would use one of those plugins/tools, but since the one use in GRASS requires to select one pixel at a time, I was trying to used the other one (QGIS plugin for visibility analysis).
So here is the configuration of the tool that I used:

--> From this point on I have problem achieving what I wanted <--

After I run it I have an error associated with the fact that I am not passing some fields ( which I think are also fundamental, observ_hgt and radius) and I don't know either where to put them...

I even tried to add the fields observ_hgt with the value 1.5 and radius 2000 to each pixel (by clicking Toggle Editing > Open Field Calculator) but afterwards I got another error (as can be seen below):


texts: Deleting pixels of a certain class that are connected to less than 2 other cells in R?

texts:
I want to process my data (classified single-band raster) according to a specific idea. Pixels that do have less than 2 neighbors of the same class (eg. with the value 1)  shall be removed and classified as NA.
It seems a lot trickier than I expected this part to be. I imagined it might be similar to the GameOfLife (Cellular Automata). But I could not implement it with a similar methodology of GameOfLife (I tried it like in the examples of the R package 'focal'). I also tried to use focal in general, which seems not to work good as well, as there are restrictions to the function that can be used (several number as input, one number as output) and I also did not manage to think through how to fix the problem with the outer cells of the moving window.
 Now I tried to use the function adjacend in which I can get the cell-index of 4 or 8 neighbor-cells. This makes it very inefficient though and can not stop thinking that there must be a better alternative.

Has anyone a good idea how I could tackle the problem to change the value of a cell of a specific class in my raster to NA, if it has less then 2 cells with the same classification that it is connected to?

I also tried to work with the mean function in the focal-package, but that approach is way too coarse for my needs.

Original classification
(white = NA, green = 0, red = 1, blue =2):

I want to get rid of the single pixels of only one class, lets say class one which is shown in red in my plot. I don't want the function to delete any values of other classes though, like class 0 (green) or class 2 (blue).

After trying several methods (game of life, adjacent of the raster package, focal of the raster package, clump of the raster package) I figured out that the clump command was the easiest to implement. I came quite far with a function that could be called in focal (from @Spacedman), but it did not work without problems. Clump seems to be an easy way and the calculation is fast. If you want to perform the command only on special classes, like me, you have to define them as a band first and you can later use raster calculations to get the result you want to achieve.

Classification after using the clump command with clumps < 3 of the red class deleted:

The difference might look small, but it's important for my computation.


texts: Produce line from components, length and angle

texts:
I have a CSV file, or point data file which includes information about line vectors which I have previously calculated. This file is outlined as such:

x, y, dx, dy, length, angle(east_from_north)
1,1,1,1,1,1
2,2,2,2,2,2
...
n,n,n,n,n,n

Where x, y are the co-ordinates of each point; dX, dY are the Cartesian components of the line; length is the length/magnitude, and angle(east_from_north) is the bearing of the line (east 0degrees, north 90degrees).

I also have the the original set of x, y data from which I calculated the dX, dY components, length and angle using some python trig, if that would help this process. For example:

x1, y1, x2, y2
1,1,1,1
2,2,2,2
...
n,n,n,n

I would like to produce lines from each point using the x, y values as the originating point with the length and angle creating the length and direction of the line.

I have used a QGIS plugin, Vector Field Renderer to produce lines with arrowheads depicting the 'magnitude' and direction using this information which has been useful for one part of my analysis; however this does not transform the point data into polylines, only produces a line symbol -- for example:

Effectively I would like to reproduce this output, but as a polyline file for which I can do further analysis on. The line does not need to have arrow heads.

Edit: Solutions involving QGIS, GRASS GIS, GDAL, ArcGIS or Python/ ArcPy would be greatly received. Though, as I am currently learning (and a bear of little brain) a well-commented/ walkthrough-esque solution would be amazing!


texts: Online POSTGRESQL database error "undefined symbol: GEOSFrechetDistanceDensify"

texts:
Using QGIS 3 on mac book air and set up a POSTGRESQL database with Acugis 6 days ago. Everything's been working fine and no issues. The yesterday everything went off and was unable to see the layers on POSTGIS window. Lets me connect and login to the schema but I've been getting the error below.

ERROR: could not load library "/usr/pgsql-11/lib/postgis-2.5.so": /usr/pgsql-11/lib/postgis-2.5.so: undefined symbol: GEOSFrechetDistanceDensify
SQL state: XX000

I've researched it but am unable to find a resolution on how this could be fixed.


texts: How to specify relationships (foreign keys) in a spatially enabled database for census geography?

texts:
Lets use the Canadian census geography design as an example:

In my experience, a spatially enabled database for census geography can be created in two ways:

Each geography contains a foreign key of the first order geography above it (FK-Province in EconomicRegion table must exist in Province table; FK-CensusDivision in CensusConsolidatedSubdivision must exist in CensusDivision table)
Each geography contains a foreign key of every geography above it (FK-Province will be in every geography below it... e.g. FK-Provinces in DisseminationBlock must exist in Province table)

This is likely a battle between storage space and computation time:

Storing all related/nested geography FKs in every table increases database size
Storing all related/nested geography FKs allows us to extract data without having to create many joins

One could argue that PK/FK relationships are nonsensical and one should rely on spatial joins instead of creating mundane long-winded join statements, or redundantly storing FKs all over the place.

Can anyone shine any light as to how geography of this nature would be stored in a spatially enabled database?


texts: Mapnik OGCServer doesn't work if theres no active internet connection

texts: Didn't find anything similar so: I have a mapnik OGCServer on an ubuntu 12.04 virtual machine (vmware player) which works fine. The problem occurs  when i disconnect all network connections, when ogc server kind of hangs. On some occasion there was a fault saying something about a broken pipe. I actually noticed this on the airport where i had no connections. Any ideas what might be going on?


texts: Does software exist to create Flow Maps (for migration of people) on Mac?

texts:
Does anybody know a good, and easy to use, software to create flow maps?

I want to create flow maps for showing the flow of people between places (migration).

The software has to be compatible with Mac computers.

I have tried the plugin for QGIS, but it was designed for Windows.


texts: Optimizing a large number of small Geotiffs

texts:
I have a large number - more than 60 thousand - of relatively small Geotiffs, usually from 2 to 12 Mb. I am currently serving them using Geoserver, but things are quite slow even inside our own network. I have a hunch that both Geoserver and these Geotiffs are not configured in the best way possible.

Here is the output from gdalinfo for a typical one:

Driver: GTiff/GeoTIFF
Files: se502c_R9_C1_3_T_2531_sec.rsl.tif
Size is 3396, 2271
Coordinate System is:
PROJCS["ETRS89 / LAEA Europe",
    GEOGCS["ETRS89",
        DATUM["European_Terrestrial_Reference_System_1989",
            SPHEROID["GRS 1980",6378137,298.2572221010002,
                AUTHORITY["EPSG","7019"]],
        TOWGS84[0,0,0,0,0,0,0],
        AUTHORITY["EPSG","6258"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4258"]],
    PROJECTION["Lambert_Azimuthal_Equal_Area"],
    PARAMETER["latitude_of_center",52],
    PARAMETER["longitude_of_center",10],
    PARAMETER["false_easting",4321000],
    PARAMETER["false_northing",3210000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","3035"]]
Origin = (4645325.000000000000000,3992425.000000000000000)
Pixel Size = (25.000000000000000,-25.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 4645325.000, 3992425.000) ( 15d37'30.88"E, 58d55' 2.56"N)
Lower Left  ( 4645325.000, 3935650.000) ( 15d32'42.77"E, 58d24'30.92"N)
Upper Right ( 4730225.000, 3992425.000) ( 17d 5'21.12"E, 58d50'57.39"N)
Lower Right ( 4730225.000, 3935650.000) ( 16d59'19.12"E, 58d20'29.67"N)
Center      ( 4687775.000, 3964037.500) ( 16d18'44.44"E, 58d37'52.11"N)
Band 1 Block=512x512 Type=Float32, ColorInterp=Gray
  NoData Value=-9999

As you can see they are EPSG:3035. They should be internally tiled. In the Coordinate Reference Systems for Geoserver, it says that both Native SRSand Declared SRS are EPSG:3035, and that the handling should be to reproject native to declared.

My web client overlays these Geotiffs on a standard OpenStreetMap layer in Web Mercator. All the geotiffs we have tried appear correctly, so I guess reprojection is still happening at some stage.

Would I gain any speed changing the declared SRS to Web Mercator? Or reprojecting the original Geotiffs to Web Mercator? What other properties of the Geotiff could I tinker with?

Bonus question - if a Geotiff is changed, reprojecting or retiling or what have you, does the Geoserver layer associated with it also change automatically?


texts:
I'm trying to extract multiple shapefiles from a folder, all these shapefiles are zipped and some of them have subfolders under those zips (most of them have shapefiles right under the zip structure)

I'm trying to automate unzipping these shapefiles (and accompanying files into one main folder without keeping the structure.

my_dir = "/Users/myusername/Documents/My_Dataset/new"

for subdir, dirs, files in os.walk("/Users/myusername/Documents/My_Dataset/"):
    for file in files:
        #print os.path.join(subdir, file)
        filepath = subdir + os.sep + file

        if filepath.endswith(".zip"):
            print (filepath)

This code allows me to print every single zip file name regardless of they are in the folder or in a subfolder. I'm trying to unzip all these shapefiles into a folder but some of these zipfiles has subfolders and I couldn't find a way to go into those subfolders and extract shapefiles from there. All I was able to achieve is unzipping these zipfiles while keeping the subfolder structure. How can I get rid of that?


texts: The resulting values are not what they indicate?

texts:
I am using ArcMap 10.3.

When I create a watershed, using the Watershed tool in the Hydrology tool set, I am receiving values that should really be representing the number of water splitting lines, correct ?

However, I did not receive those results, what does the value I'm receiving then indicate?

Here is the resulting graph:

here is the tool set I'm working with:


texts:
I try to find topology error in this polygon layer using QGIS 2.18 with topology rule of "must not overlap":

but i don't understand why with the Topology checker, only several areas marked in red, whereas i expect that also other areas (i marked it with blue circle) will be defined as topology error:

this is the attribute table of the layer:


texts: Map symbols fade after exporting to PDF, JPG, BMP from ArcGIS for DEsktop?

texts: Any solutions how to make the legend symbols not blurry after exporting as an image?


texts:
I'm new to QGIS, and have what seems like basic requirements; however I'm struggling to find a solution. I have GPS data in excel showing start and end points of a boat survey. In a row, it has start x, start y, end x, end y.

I would like to plot these points on a map, and draw a line between each start and end point, but without drawing lines between points in different rows.

I've seen a few similar questions in the forums, but none of them seem to work when I try. Points to path seem to plot lines between rows, and only used one point per line, rather than the two required. I've also seen WKT mentioned, but I don't know how this works.


texts: Simple iterate SetNull in ArcPy not working

texts:
For a large number of rasters, I need to define a pixel value less than zero as being Null and output the rasters. However, my simple SetNull loop (below) only processes the first raster in my list as Null, before stopping and not looping through the rest.

Triple commented ### is a dummy ExtractByMask function which I used to test the loop. When commented in, the ExtractByMask does loop successfully through all my rasters, indicating the loop itself seems okay. Can someone please advise, is there something specific I am not doing with SetNull please?

import arcpy, os, sys
from arcpy import env
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = True

arcpy.env.workspace = "C:\\GIS\\in" # set up a folder containing rasters to clip
outws = "C:\\GIS\\out" # set up a folder to output clipped rasters to
#mask = "C:\\GIS\\dump\\pan.shp" # shapefile that defines clip

rasterlist = arcpy.ListRasters("*", "tif")
print(len(rasterlist))

for i in rasterlist:
    print i
    temp_nulled = SetNull(i,i,"VALUE < 0")
    ####temp_extracted = ExtractByMask(i, mask)
    rastername = str(i) # converts loop(i) raster name into a string
    rastername_crop = rastername[:7] # takes first n characters of raster name
    outname = os.path.join(outws, rastername_crop + '_NULL.tif') # Create the full out path
    ####outname = os.path.join(outws, rastername_crop + '_EXTRACT.tif')
    print outname
    temp_nulled.save(outname)
    ###temp_extracted.save(outname)


texts: Creating multipoints from points and setting incrementing M value in spatialite

texts:
I have a layer of points, which I want to group into multipoints grouping by some column. This part is easy:

SELECT name_pre, CastToMulti(GUnion(geometry)) AS geometry
FROM sp_survey
GROUP BY name_pre;

But the trick is that I'd like to assign to each part of multipoint incrementing m value (0, 1, 2 etc.) based on the "depth" column in sp_survey table, ordered from lowest (0) to highest (n) for each multipoint feature.

Is it possible in spatialite?


texts: Trying to Create a Global Choropleth Population Density Map then add Biodiversity Level

texts:
I am trying, quite unsuccessfully, to create a global choropleth population density map for a specific 120 countries, using ggplot2.  Once I have created this choropleth map, I want to add number of threatened species (biodiversity) on top. I have gathered the population of all 120 countries and have the area of each country also. I also have the total number of threatened species per country.  I am trying to create something like this:

Here is my code so far:

(Import Data set called "data")

library(ggplot2)
library(splines)
library(mgcv)
library(plyr)
library(maps)

Creating data sets

humpopden2010 <- (data$X2010.Census)/(data$Area) # 2010 Census Population Divided by Country area

specthreat2010 <- (data$ThreatMammals2010 + data$ThreatBirds2010)/(data$Area/10^6)

specrich <- (data$TotalMammals2004 + data$TotalBirds2004)/(data$Area/10^6)

ATTEMPT 1 (Just trying to make pop density map)

world_map <- map_data("world")
names(humpopden2010) <- tolower(names(humpopden2010))
humpopden2010$Country <- tolower(rownames(data$X2010.Census))
choro <- merge(world_map, population, by = "Country")
choro <- choro[order(choro$order), ]
qplot(long, lat, data = choro, group = group, fill = humpopden2010,
  geom="polygon")

ATTEMPT 2

world_map <- map_data("world")
names(world_map)[5] <- "country"

p_country <- ddply(data, "country", summarise, total = length(Country))

p_map <- merge(p_country, world_map, by = "country", all = T)
p_map <- p_map[order(p_map$order), ]

ggplot(p_map, aes(long, lat)) +
geom_polygon(aes(group = group, fill = log(total)), colour = "grey60", size = .3) + ylim(-55, 85)

I can't even get the pop density map to create... need to do that before I can try to put on the number of threatened species.

Made a map with nothing on it here...

world <- map_data("world")
worldmap <- ggplot(world, aes(x=long, y=lat, group=group)) +
geom_path() +
scale_y_continuous(breaks=(-2:2) * 30) +
scale_x_continuous(breaks=(-4:4) * 45)

print(worldmap)


texts: How to output coordinates from a shapefile to a GMT or CSV file in QGIS?

texts:
I want to get the coordinates for a line shapefile as csv or gmt files.

I click create "new layer" - "new shapefile layer" - I check the "line" box and then "OK". Then I make my shapefile. When I move the cursor I can see the coordinates at the lower bar but when I create a new shapefile and look in the attribute table there is no coordinates.

How to get the all the coordinates for the shapefile?


texts: ModelBuilder: How to force a specific transformation in the "Project" tool?

texts:
At 10.7.1, here's the first part of a model I am creating.  Please note that I have parameterized the input and output:

My problem is within the Project tool: I want it to always use the same datum transformation, but I am unable to get that transformation to appear.  When editing the model, I can drop-down the tool's transformation list, but I cannot make the selected transformation "stick" in the transformation box; the box remains blank. Weird, because I can select the desired Output Coordinate System just above.
I think that the problem has to do with the parameterized input, because if I hard-code an input shapefile, the transformation sticks.

How do I "force" the desired transformation to appear in the model while still allowing parameterized input?
EDIT:
For this tool, the input data will always be NAD83 (specifically:
NAD_1983_StatePlane_Washington_North_FIPS_4601_Feet
and the output data will always be NAD83 HARN (specifically:
NAD_1983_HARN_StatePlane_Washington_South_FIPS_4602_Feet / VCS:NAVD_1988)
When I Project these two manually, Arc offers a variety of transformations.  I always select:
NAD_1983_To_HARN_WA_OR.
This is the transformation that I want to force my tool to use.


texts: Is kriging an exact interpolation method?

texts:
Can kriging be an exact* method? For example, when optimized with the Geostatistical Wizard? Or in any other way? (as far as I know, "normal" kriging isn't an exact method).

*exact: predictions at known observed values yield exactly the same values again.


texts: Layer doesn't turn back on after clicking in OpenLayers 3

texts:
I have the following code. I have setup a default function of turning the showing all my layers and on click the layers disappear. But they don't turn back on when I click on them.
So what have I done wrong here?

 <input type="checkbox" onClick="vector_summer.setVisible(false);"  id="ch_summer" name="summer" checked>Summer
 <input type="checkbox" onClick="vector_winter.setVisible(false);"  id="ch_winter" name="winter" checked>Winter
 <input type="checkbox" onClick="vector_relay_route.setVisible(false);"  id="ch_SLC2002" name="SLC2002" checked>Salt Lake City 2002
 <input type="checkbox" onClick="slc_2002_torch.setVisible(false);"  id="ch_none" name="none" checked>Torch

Then this is how I am calling the layers later on:

var map = new ol.Map({
  layers: [rasterLayer, us_states_vec, slc_2002_torch, vector_summer,
           vector_winter,vector_relay_route],
  overlays: [overlay],
  target: 'map',
  view: new ol.View({
    center: ol.proj.transform([34.59, 31.59], 'EPSG:4326', 'EPSG:3857'),
    zoom: 2
  })
});
vector_summer.setVisible(true);
vector_winter.setVisible(true);
vector_relay_route.setVisible(true);
us_states_vec.setVisible(true);
slc_2002_torch.setVisible(true);

What do I do so that the layers get turned back after I uncheck them?


texts: Are there any exisiting .NET implementations for WFS that are extendable?

texts:
I've been told its a bad idea to implement WFS from scratch:
Implementing WFS within .NET

Does anybody know of any exisiting .net implementations that I can extend?


texts: Change the default output width and height in QGIS "Save map as image" feature

texts:
I need to output many snapshots of the project I'm working on in QGIS (v3.10) into reports I create. I'm currently using the 'save map as image' function which is great as it allows copy to clipboard and I can then paste maps into my report, this is great as it's really quick and I've got reports that have had up to 50 maps in previously and foresee creating even more voluminous work in the future.

Trouble is I'm having to enter the output width and output height each time, I know the exact size to perfectly fit the map into my tables in my report template (1610 x 1040px), can I set these to replace the default 977 x 920px values?

Only other workaround I can think of is setting my project canvas to the size of my report paste area (1610 x 1040px) and using the 'map canvas extent' setting in the 'save map as image' feature, can you set a custom project canvas size (I could save it as a template for future work going forward)


texts: What projection should be used along Gulf Coast?

texts: What projection should I use along the Gulf coast - spanning the west coast of Texas to the Florida panhandle. I need to perform spatial analysis for a data-set along the coast. I see a similar question has been asked before but my query is different as I am specifically searching for a specific projection in ArcGIS which allows for spatial analysis in the region described above The answer to the similar question is to adopt an Orthographic projection near the center of the Gulf Coast Region. I am new to ArcGIS and wanted to confirm a projection which can be used for spatial analysis and will not generate wide error margins.


texts: Raster data like heatmap from polygons w.r.t their own values

texts:
I have polygons as below, and each polygon is assigned to some value.
What I'd like to do are ...

Smoothen boundaries (Values should change more smoothly around boundaries)
Fill N/A part properly
Generate raster data which met the above

How can I get such a raster data with QGIS ?


texts: create individual buffers from vector points in QGIS

texts:
I am using QGIS and the python console to achieve the following: my problem is simple in that all I want to do is to create individual buffer shapefiles for each vector point.

The code I am using produces one shapefile with all the buffers in the that one file. What I'd like it to do is: for 5 points in the vector file I want 5 individual buffers.

from qgis.utils import iface
from qgis.analysis import QgsGeometryAnalyzer
import processing

mc = iface.mapCanvas()
layer = mc.currentLayer()

QgsGeometryAnalyzer().buffer(layer, "path\to\output\output.shp", 50000, False, False, -1)

I understand that this might be achievable using a for loop so that would require a small change like so

for feature in layer:
    QgsGeometryAnalyzer().buffer(layer, "path\to\output\output.shp", 50000, False, False, -1)

but this produces the following error

Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "c:/users/lazarusa/appdata/local/temp/tmprs5egx.py", line 9, in <module>
    for feature in layer:
TypeError: 'QgsVectorLayer' object is not iterable


texts: How you get values within intersect layers and autofill attributes to another layer form when creating a new object

texts:
When im creating a point that intersect on an polygon, how you can get values from attribute polygon and autofill the attribute in the point layer.

There will be some code?

I´m using qgis 3.12

Thanks everyone, im so sorry for my bad english


texts:
Actually, when I'm trying to do a landuse supervised classification with Google Earth Engine, a strange message appears :

  Collection.geometry: Unable to perform this geometry operation.
  Please specify a non-zero error margin.
      in , line 11
      in , line 24

What is wrong with my code?

var image = ee.Image(ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
    .filterBounds(roi)
    .filterDate('2018-05-01', '2018-06-30')
    .sort('CLOUD_COVER')
    .first());
Map.addLayer(image, {bands: ['B4', 'B3', 'B2'],min:0, max: 3000});
var classNames = turban.merge(twater).merge(tforest).merge(tagriculture).merge(tsoil);
print(classNames);
var bands = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7'];
var training = image.select(bands).sampleRegions({
  collection: classNames,
  properties: ['landcover'],
  scale: 30

});

print(training);
var classifier = ee.Classifier.cart().train({
  features: training,
  classProperty: 'landcover',
  inputProperties: bands,
  });
var classified = image.select(bands).classify(classifier);
Map.centerObject(classNames, 11);
Map.addLayer(classified,
{min: 0, max: 4, palette: ['red', 'blue', 'green','yellow','orange']},
'classification');


texts: Do any geographers on here have an estimate of the mean and median city radius sizes? I'm trying to optimize some Google Places API queries [e.g.] and this information would be helpful for that purpose.


texts: How to approximate point to segment distance on sphere?

texts:
I would like to approximate the distance of a line segment from a reference point on the map (coordinates are gives as lat/lon).

In the Cartesian coordinate system I would simply find the closest point of the line segment to the reference point using the dot product and calculate the distance.

I am new to working with geo-coordinates, and from what I have learnt, segments are called great circles on the Earth's surface. This site gives a clear formula for calculating what I would like (which is I think called the cross-track distance). However it seems an overkill for what I want to do. Firstly, I am interested in shorter distances only, maximum a few hundred kilometers (say 700 km). Second, I can do with an approximate value (say maximum 2-3% error). Third, I want the calculation as quick as possible, and the exact formula contains around 20 trig calculations.

I would need something similar to equirectangular approximation:

x = (lon2-lon1) * cos((lat1+lat2)/2)
y = lat2-lat1
d = sqrt(x*x + y*y) * R

Is there a way to generalize the above formula for point-segment distances? Or, alternatively, is there a handy approximation to find out the distance of a point and a great circle for small-to-mid distances with reasonable accuracy?


texts: Counting features into intersections of Shapely Polygons

texts:
I have a geopandas GeoDataFrame containing hundreds of shapely Polygon and MultiPolygon geometries. The polygons overlap in lots of places. I'd like to make a new geometry that contains a count of how many of them overlap. Something like this:

Does anyone have any ideas for how to approach this? I can't even see a way in.

Eventually I'd especially like to be able to weight the polygons, so that some polygons might be worth 2 on their own. Doing this with shapely's Z field might be nice.

Aside: I'm not particularly tied to any of these libraries, it's just where I've ended up. The coordinates in these geometries are actually pixel coordinates — I'm stumbling towards making a raster to overlay on another image. I'd prefer to keep my footprint as small as possible, since I'd like to be able to deploy this stuff to cloud servers etc, where I might not be able to install random stuff.


texts: Professional Indemnity Insurance for freelance GIS?

texts:
I am about to go freelance and I know there are a few freelancers on the forum so I'd welcome recommendations/experiences relating to Professional Indemnity Insurance.  Are there companies (brokers or providers) that specifically cater for GIS and would you recommend them?   I am based in the UK - so, while all replies are welcome, I am particularly looking for responses from UK freelancers.

Following Mapperz comments, to be clear: This question is about insurance for GIS professionals and NOT GIS for insurance analysis.  I am looking for information on companies who PROVIDE professional indemnity insurance to freelancers rather than work opportunities. I am also interested in knowing the sort of level of cover people have and whether the level of cover has ever been an issue when tendering for work.


texts: How to convert/typecast from Polygon to Multipolygon in PostGIS?

texts:
Is there any way to typecast Polygons to Multipolygon type?

For example, I have created a circle using ST_Buffer() on a point. The return type of that function is of Polygon type. The problem is that on the geometry column there is a restriction limiting it to Multipolygon type. So, I cannot insert a simple polygon.

I have used ST_MPolyFromText() function, but as usual it is returning null for that and ST_Dump() is not giving the proper output.

The constraint is:

CONSTRAINT enforce_geotype_geom CHECK (geometrytype(geom) = 'MULTIPOLYGON'::text OR geom IS NULL),


texts:
I've got a few on demand scripts via the web that run ArcGIS in the background via Arcpy. Now, every one of these scripts takes a solid 15-30 seconds to run, which I find okay. After doing some trial and debugging on it, I find that 90% of the time spend running the script is taken importing Arcpy.

I've broken down a script I wrote here

Import standard python modules like: os, cgi, time (less than 1 second)
Grab web query and disseminate into a dictionary (less than 1 second)
Import ArcPy (15 seconds)
Import utm (less than 1 second)
Convert lat long to utm (less than 1 second)
Arcpy insert cursor, insert a new point into a feature class (less than 1 second)
Arcpy update cursor, add attributes to a the new point (less than 1 second)
Print web page based on outcome of query and Arcpy cursors (less than 1 second)

This whole script take about 16-20 seconds to run and most of that is spent importing Arcpy. The PC I am using is good, so what is the issue?

Is Arcpy so big that it actually takes this much time to import? Is there any other reason why the import takes so long? Also is there any work around to speed this up?


texts: rectifying KML file to match GoogleEarth imagery

texts:
Ubuntu 14.04, QGIS 2.12:

In QGIS I have a polygon (33 segments) which I have successfully exported to GoogleEarth as a KML file.
Due to rectification of GE's image, the polygon doesn't quite line up, of course.

The one certain point I can line up in both mappings requires a shift of KML data 5 metres westward.
So I have carefully edited the KML file to reduce the Easting of every point by 0.00004° .
And it works, but it is very laborious and error-prone, so I'm thinking there must be a better way.

Really I want to do a rectification on the polygon's coordinates, to stretch/squeeze/rotate to match what Google has done to their image.

I can certainly write a PHP-CLI script to manipulate the file, but is there a better way, or has someone this before?


texts:
I have bug with QGIS  version 3.10 (also happenes in 3.8) that I don't know how to fix.
Since few days ago, the desktop screen looks like this:

It might look fine on first sight, but the qgis window covers the start menu and also I cant close the screen or minimize it because I "lost" the acess to the '-' or 'x' buttons on the top.

I have been told that it might relate to plugins but I couldn't find any solution or solve it.

My end goal: to have the screen on the right size again


texts: Should OGR2Layers plug-in respect layer transparency?

texts: When I export a graduated style layer from QGIS Dufour to OpenLayers using the OGR2Layers plug-in and QGIS rendering the colors and categories come across correctly. However, the opacity in the HTML is set to 1, where it should be .3 in my example . Am I doing something wrong, or does OGR2Layers not handle transparency? It's easy enough to edit the transparency in the HTML, but if there's a direct way, that would be nice.


texts: ArcGIS Online - Programmatically Upload Tabular Data Source

texts:
I have some ArcGIS Online maps and want to reference a simple data source like a CSV or XML file.  I know that I can upload the file, but is there a way I could programmatically upload the data file on a scheduled basis?  Or in other words, is there an API that I can use in a program or script that allows me to send data files to my ArcGIS Online storage location?

The first thought I had was to write a program that connects as me and sends the file by POSTing to the web service, but there might be a better way.


texts: Generate a grid through an extensive area without killing my machine

texts:
I need to generate a grid of 10x10m on a certain extent so I can extract only the squares that have some area inside of a certain shape.

My problem is that the shape I need to extract those squares are just 60m wide, while being 400km long. Since it's the shape of an actual existent road and its vicinity, if we check it's extent, it's a rectangle of about 150x400km, and running a fishnet on it through ArcGIS comes back negative because the file size would be too big, and tessellation is taking so far 50 minutes and I'm afraid my CPU is overheating because of this.

Is there any way I can do this without having to make multiple instances of fishnet or tessellation, on ArcGIS 10.x, to merge them together later?


texts: Precision loss when trying to select the latitude and logitude coordinates from an sde.ST_GEOMETRY using SDE.STX() and SDE.STY() functions

texts:
In an oracle 11g enterprise geodatabase I am trying the following:

INSERT INTO MySTATION (STATION_ID, DB_VERSION_NU, CREATED_DT, MOD_DT, SHAPE)
VALUES (1, 0, sysdate, sysdate, SDE.ST_GEOMETRY('multipoint z(1.2345 1.2345 10)', 4326));

INSERT INTO MySTATION (STATION_ID, DB_VERSION_NU, CREATED_DT, MOD_DT, SHAPE)
VALUES (2, 0, sysdate, sysdate, SDE.ST_GEOMETRY('point(1.2345 1.2345)', 4326));

SELECT STATION_ID, SDE.ST_X(SHAPE), SDE.ST_Y(SHAPE) FROM MySTATION ;

For both of the above inserts the x, y coordinates returned with the select are 1,23450000000003 instead of 1.2345. Is there a way to read/write data to an ST_GEOMETRY column without loss of precision?

Also, can anybody explain me why is the precision lost?


texts: Clip vector on Geoprocesing QGIS 2.18.13 "Algorithm Clip starting... global name 'features' is not defined"

texts:
I'm using QGIS version 2.18.13 and I'm getting the error message while clipping a layer line within a polygon.

Error message:

  "Algorithm Clip starting... global name 'features' is not defined"


texts: Use a python expression in Field Calculator to sum each row across several columns

texts:
this seems like it should be easy but I am new to python and am having trouble... I have four columns. I would like perform a field calculation that sums the values by row for three of the columns and inserts that sum into the associated row of the 4th column. In vb script is it very straight forward. I would like to know how to do it in python.
ie.

Col1    Col2    Col3    Col4(SummedValues)
  1       3     1              5
  2       1     1              4

Here is what I have:

I have tried adding a for statement (for rows in !resuly!:)
but it returns a parsing syntax error (Error 000989 Python Syntax Error: Parsing error Syntax Error: invalid syntax (line 2)).

I am in ArcMap V.10.4

I feel like it is an iteration issue- maybe I need to tell it to go through each row i just don't know how to do that correctly in the field calculator...

Per Phillip and JR comments: Code updated to and it worked perfectly.


texts: Find locations of raster maximum within polygons: Basin Delineation

texts:
I'm trying to identify the Location (lat/lon) and upstream area of every raster point that flows into the Arctic Ocean (there are a lot so it has to be fully automated).  I'm starting with the following rasters:

Flow Direction
Accumlated Flow Area

My approach so far has been to use the above rasters to create a delineated basin shapefile.  From there I was hoping I could use the shapefile to look at the Flow Accumulation Raster and identify the location of the maximum value in each basin.  In terms of output, I just need a table that shows [Basin#, Longitude, Latitude, UpstreamArea] for every basin.

I've already looked at How to find the locations of the highest values in a raster? but I don't know how to restrict the raster to only the area bounded by the polygons in the shapefile]

**Using: ArcGIS10, Python, etc.*


texts: Difference between ST_Simplify function in PostGIS 3.0.0 and PostGIS 2.5.3

texts:
Faced such a problem (or maybe it's a bug).
I get different results when using the ST_Simplify function in PostGIS 3.0.0 and PostGIS 2.5.3.

There are source data that are simplified by using a query:

SELECT ST_Simplify(geometry, 20) FROM source_table

For verification, I created two databases on different OS. PostgreSQL11-PostGIS-3 and PostgreSQL11-PostgiGIS-2.5 with the same data mapping and SRID.

Result:

Postgis 2.5.3

Postgis 3.0.0

PostGIS 3 lost a lot of small hole geometries, although the source data and the request are the same. Why is this happening and how to fix it in PostGIS-3?

Update

Yes, SRID 3857 in both cases.

I try ST_Simplify(geometry, 20, false) and ST_Simplify(geometry, 20, true) in Postgis 3 it did not help.

Postgis_full_version();
POSTGIS="3.0.0 r17983" [EXTENSION] PGSQL="110" GEOS="3.7.2-CAPI-1.11.2 b55d2125" PROJ="Rel. 5.2.0, September 15th, 2018" LIBXML="2.9.4" LIBJSON="0.13.1" LIBPROTOBUF="1.3.1" WAGYU="0.4.3 (Internal)"


texts: Installing Whitebox in QGIS without broken repo link

texts:
I need to install Whitebox tools in QGIS. I have tried to follow the instructions provided in here but the repo link is broken (https://plugins.bruy.me/plugins/plugins.xml and https://plugins.bruy.me/processing-whitebox.html). The same problem is also mentioned here in this post
I have tried to install the plugin from zip but it did not work. The zip file was download from Github repo.
Is there a working link to the repo or another method of adding the plugin offline in QGIS?


texts: Does ArcGIS Pro support personal geodatabases (mdb)?

texts:
Does ArcGIS Pro support mdb or not?

In ArcGIS Pro,I can't find anywhere can load data from mdb;
I have not find the help page.


texts: Using Postgresql/PostGIS/Pgrouting in PYQGIS Processing Scripts

texts:
I have started exploring PYQGIS scripts and learning writing PYQGIS user scripts in python.

I want to use pgsql/postgis/postgresql in the script. Please suggest
Also Can anyone suggest That the script starts with a double Python comment symbol (##) and have the following structure

  [parameter_name]=[parameter_type] [optional_values]

Can the parameter_type be a database connection  So once the script is executed the pop window ask for the connection with the database etc.

How can one use pgrouting algorithms inside scripts.


texts:
How do you configure TauDEM 5.1.2 to run in QGIS 2.10.1?

TauDEM plugin activated, and visible in the Processing Toolbox (Advanced Interface).

However, TauDEM tools do not run.

I get the message '"mpiexec"' is not recognized as an internal or external command, operable program or batch file.

NOTE: TauDEM DOES run from the command line interface (mpiexec is recognized).

System paths have been updated for MSMPI and TauDEM. OS is Windows 8.1' TauDEM 5.1.2 and QGIS 2.10.1 are each 64-bit installations (confirmed).


texts: Opening gkpg in geopandas throws: UnsupportedGeometryTypeError: 12

texts:
I wanted to open a gkpg file in geopandas 0.8.2, on windows, using conda mith channel conda-forge
obj = r'path.gpkg'
layer="layername"
gdf = gpd.read_file(obj, layer=layer)

Throws the errors:
fiona/ogrext.pyx in fiona.ogrext.Session.get_schema()
fiona/_geometry.pyx in fiona._geometry.normalize_geometry_type_code()
UnsupportedGeometryTypeError: 12

other programms using fiona can open the data, so I tried it in fiona:
features = []
with fiona.open(obj, 'r', layer=layer) as src:
     for feature in src:
            features.append(feature)

features

Which made me wonder why, I cannot open it with geopandas using fiona but with fiona itself?
What even works is using the features to build a geoDataFrame
gdf = gpd.GeoDataFrame.from_features([feature for feature in features], crs=25833)


texts: Automated way to create polygons from coloured map

texts:
I have a scanned map that has areas of different colours.  I want to create polygons for each area and would like to avoid doing it manually (tracing over the shapes).

Is there a way to do this?  I do not have ArcGIS - am using QGIS.

I have tried various feature extraction options to no success (probably choosing the wrong options).

I have tried creating a layer of just the e.g. blue sections with everything else transparent, importing that, and then polygonising this - but I just get a polygon the size of the whole area.

I have also tried to follow How to digitize a raster automatically (Processing toolbox no longer exists, options I tried in Orfeo unsuccessful); and https://docs.qgis.org/2.2/en/docs/training_manual/forestry/stands_digitazing.html (polygons whole area).

I have something like this

except the colours aren't quite so clear - it's a shaded map overlaying satellite imagery, so there's lots of tiny "gaps", but I can fill these without too much hassle if needed.


texts:
In the layer tree, I want to move a few layers at root level to a group also at root level.
From what I have gathered so far, I have to clone the layer first, add it to the group and then delete the original layer.

https://docs.qgis.org/3.4/en/docs/pyqgis_developer_cookbook/cheat_sheet.html#advanced-toc

So I am trying to fill the variables with the group and layers I want to work with.
The for-loop finds the group, but for some reason, I am getting None for the value of LyrOne.

from qgis.core import QgsLayerTreeGroup, QgsLayerTreeLayer
root = QgsProject.instance().layerTreeRoot()

for child in root.children():
    if isinstance(child, QgsLayerTreeGroup):
        print ("- group: " + child.name())
    elif isinstance(child, QgsLayerTreeLayer):
        print ("- layer: " + child.name())

grpQS = root.findGroup("QS")
LyrOne = root.findLayer("One")
print(QS_group)
print(LyrOne)

The results of the for-loop and the print statements are:

- layer: One
- group: QS
<qgis._core.QgsLayerTreeGroup object at 0x0000000008032A68>
None

Why can I get a reference to the group, but not for the layer?

EDIT:
I managed to create a dictionary with the names of the layers as keys

dictLayer = {}
for child in root.children():
    if isinstance(child, QgsLayerTreeLayer):
        lyrName = child.name()
        lyrID = child.layerId()
        dictLayer[lyrName] = lyrID

LyrOne = dictLayer["One"]

Unfortunately, this solution
Add layer to a QGIS group using Python

toc = self.iface.legendInterface()
groups = toc.groups()
groupIndex = groups.index(u'myGroup')
toc.moveLayer(newLayer, groupIndex)

doesn't work in QGIS3.
How can I move a layer to a group in QGIS 3?


texts: Direct connect to PostgreSQL without ArcSDE?

texts:
I have postgreSQL installed on my machine.My machine already has ArcGIS 10. No Sde installed. Can we create geodatabase in PostgreSQL with out installing ArcSDE?

Jay Cummins is right. I am intended to use create esri geodatabase. Can we do this with direct connect by simply installing postgreSQL on same/remote machine? Can we connect to Postgres with ArcCatalog and create featureclasses. Can we load raster in it? I have seen in many articles mentioning post-installer to setup the ArcSDE objects in the database even for direct connect. Should we need ArcSDE?


texts: Earth Engine: combining reducers does not give same output as running reducers separately?

texts:
I am trying to combine reducers, which according to my understanding, is running parallel yet independent reducers when using sharedInputs: true (see doc). But actually, I notice that with combine(), the output is different than compared to running each reducer separately. Is this a bug, or I did not understand properly the combine() operation?

For the example, I simply re-run the EE example on weighted reducers, which compares the output of weighted sums versus non-weighted sums. Now I just do the same thing, but combining those. We get:

Separate reducers:

weighted: -9079.394626741121
unweighted: -9086.503929115412

Combined reducers:

weighted: -9438.194106051067
unweighted: -9079.394626741121

What is going on? Why does combining reducers not give the same output as running them separately?

Code

Link: https://code.earthengine.google.com/8a4fffe6b446f908d6f1db5b33304048

// First lines are exactly from: https://developers.google.com/earth-engine/reducers_weighting

var image = ee.Image('LANDSAT/LC08/C01/T1/LC08_044034_20140318');

// Creat an arbitrary region.
var geometry = ee.Geometry.Rectangle(-122.496, 37.532, -121.554, 37.538);

// Make an NDWI image.  It will have one band named 'nd'.
var ndwi = image.normalizedDifference(['B3', 'B5']);

// Compute the weighted mean of the NDWI image clipped to the region.
var weighted = ndwi.clip(geometry)
  .reduceRegion({
    reducer: ee.Reducer.sum(),
    geometry: geometry,
    scale: 30})
  .get('nd');

// Compute the UN-weighted mean of the NDWI image clipped to the region.
var unweighted = ndwi.clip(geometry)
  .reduceRegion({
    reducer: ee.Reducer.sum().unweighted(),
    geometry: geometry,
    scale: 30})
  .get('nd');

/// New code combining reducers
var weighted_unw = ndwi.clip(geometry)
  .reduceRegion({
    reducer: ee.Reducer.sum().combine({
      reducer2: ee.Reducer.sum().unweighted(),
      outputPrefix: "weighted",
      sharedInputs: true
    }),
    geometry: geometry,
    scale: 30});

var unwweighted_w = ndwi.clip(geometry)
  .reduceRegion({
    reducer: ee.Reducer.sum().unweighted().combine({
      reducer2: ee.Reducer.sum(),
      outputPrefix: "weighted",
      sharedInputs: true
    }),
    geometry: geometry,
    scale: 30});

// Observe the difference between weighted and unweighted reductions.
print('weighted:', weighted);
print('unweighted', unweighted);

print("weighted_unw", weighted_unw)
print("unwweighted_w", unwweighted_w)


texts: Modifying properties of Attributes Table in QGIS?

texts:
In my QGIS project I'm using vector layers that require a lot of data inputting/editing.

The problem is that some of the data I'm inputting varies greatly in matter of length and structure (various number of lines of text for example). When it comes to edit those lines in many objects or to compare the data in various objects it is painful to do it by selecting each object(line or point) separately. Instead of doing all this work by selecting each object and edit its properties I'm using Attributes Table.

I would like to know is there a way so Attributes Table remembered height and width of each column?

Is it possible to modify those properties in python?

And last but not least, do you think there is reasonable effort in writing plugin that extends Attributes Table for functionalities mentioned above.


texts:
I have an sf point dataframe, with a single incorrect point. I know what it's coordinates should be, but I'm not sure how best to directly edit them. I'm aiming for a conditional edit like the following

new_point <- st_point(c(145.38626, -16.43350)) %>%
               st_sfc(crs = 4283) %>%
               st_transform(., 3577)

allpoints_clean <- allpoints %>%
  mutate(geometry = ifelse(s_id == '1', st_geometry(new_point), geometry)

but this produces the error 'no simple features geometry column present', which is definitely not the case. What should I be doing instead?


texts: Field does not get updated with Geoprocessing Service

texts:
Can anybody tell why the field gets updated when I run this script in ArcMap but not when I run a geoprocessing service?

#Import modules
import os, sys, arcpy, traceback, arcgisscripting
gp = arcgisscripting.create(10.1)
Stewardship = gp.GetParameterAsText(0)
Officevar = "\"College Station\""
# This is to send the data input to the database
with arcpy.da.UpdateCursor(Stewardship, ("City")) as rows:
     # row comes back as a tuple in the order specified here, so Office is row[0], Forester is row[1]
        for row in rows:
            row[0] = Officevar
            rows.updateRow(row)


texts: Multiple warnings when launching QGIS 2.18 Las Palmas on Ubuntu 16.04

texts:
I just installed QGIS 2.18.4 Las Palmas on Ubuntu 16.04 with the UbuntuGIS dependencies. While it seems to be working (for the moment), I get a barrage of warnings when I launch it from the terminal.

I added the key:

wget -O - http://qgis.org/downloads/qgis-2016.gpg.key | gpg --import
gpg --fingerprint 073D307A618E5811
gpg --export --armor 073D307A618E5811 | sudo apt-key add -

I added the following to my /etc/apt/sources.list:

deb     http://qgis.org/ubuntugis xenial main
deb     http://ppa.launchpad.net/ubuntugis/ubuntugis-unstable/ubuntu xenial main
deb-src http://qgis.org/ubuntugis xenial main

I ran:

sudo apt-get update
sudo apt-get install qgis python-qgis qgis-plugin-grass

The following warnings appear after launching QGIS from my terminal:

myuser@myuser:~$ qgis
Warning: loading of qgis translation failed [/usr/share/qgis/i18n//qgis_en_US]
Warning: loading of qt translation failed [/usr/share/qt4/translations/qt_en_US]
Warning: QCss::Parser - Failed to load file  "/style.qss"
QInotifyFileSystemWatcherEngine::addPaths: inotify_add_watch failed: No such file or directory
Warning: QFileSystemWatcher: failed to add paths: /home/jmnoor/.qgis2//project_templates
Warning: QLayout: Attempting to add QLayout "" to QgsPanelWidgetStack "mWidgetStack", which already has a layout
Warning: QVariantMap DBusMenuExporterDBus::getProperties(int, const QStringList&) const: Condition failed: action
Warning: QVariantMap DBusMenuExporterDBus::getProperties(int, const QStringList&) const: Condition failed: action
Warning: void DBusMenuExporterPrivate::addAction(QAction*, int): Already tracking action "&Toolbox" under id 258

Anyone experience this?


texts: Can Google Earth App for the iPad specify a KML using a URL without dealing with Google MyMaps?

texts: Can Google Earth App for the iPad specify a KML using a URL without having to go through  with one's Google My Maps Account?