texts: How to cartographically generalise buildings/settlements in QGIS?

texts:
I have a set of buildings in a polygon layer which I'm going to use in a 50k map. Since the buildings are tiny I cannot show the actual form of the building, instead I'm taking the centroid of the polygon and then make a point symbol (a square) which is a lot larger than the buildings original form.

When I do this in rural areas this is fine, but when there are a lot of tiny buildings near each other I have to do some kind of generalisation.

I have managed to remove smaller buildings that are in a certain distance from larger buildings (using buffer), but when I have a lot of buildings with the same size I'm not sure how to solve this.

Here are some images:

This image shows the original size of the building polygons

This image shows building symbols created from centroid of each building. As you can see this doesn't look good

This is somewhat what I'd like to achieve. I have removed some of the buildings manually.

Any ideas?

I've seen that ArcGIS has a function called Simplify Building which is similar to whad I'd like to achieve. Unfortunately Open Source applications on OSX is the only possibility right now.


texts: Use a map service exposed via ArcGIS REST API as a data source in Data Integration Extension ETL tool

texts:
I am trying to use the SpatialETL tool in ESRIs Data Integration Extension to cut out rasters from a map service exposed via a REST URL on the following format

http://utility.arcgis.com/usrsvcs/servers/REDACTED/rest/services/REDACTED/MapServer

How can I use this data as a data source in the SpatialETL tool?

It works fine as a layer in ArcGIS Pro but I am failing to get it to load from the ETL tool.


texts:
I am trying to creat a webmapping application to make the client to change the style of a specific Layer like this :

var ab =new ol.layer.Vector({
  name:'ab',
  source: vectorSource,
  style: new ol.style.Style({
    stroke: new ol.style.Stroke({color: '#41f46e', width: 1}),
     fill : new ol.style.Fill({ color:'#4286f4'})
  })
});
var map = new ol.Map({
  layers: [ ab],
  target: document.getElementById('map'),
  view: new ol.View({
    center: [-1095791.453557, 3422374.879112],
    maxZoom: 19,
    zoom: 5
  })
});

I did found a code from https://github.com/weskamm/ol3vectortiles to change color style after selecting a feature in the map. but what i want is to change style directly to the layer ab:

    var styleFeaturesBox = document.getElementById('stylefeatures'),
        styler = document.getElementById('stylecontainer');

    styleFeaturesBox.onchange = function(evt) {
        if (styleFeaturesBox.checked) {
            if ($('.minicolors-input').length === 0) {
                $('#colorpicker').minicolors({
                    control: 'wheel',
                    inline: true,
                    change: function(hex) {
                        handleFeatureSelection(hex);
                    }
                });
            }
            styler.style.display = 'block';
        } else {
            styler.style.display = 'none';
        }
    };

var handleFeatureSelection = function(hex) {
    var selectCtrl = map.getInteractions().getArray().
        filter(function(i){return i instanceof ol.interaction.Select;})[0],
        selectedFeature = selectCtrl.getFeatures().getArray()[0];

    if (selectedFeature) {
        var type = selectedFeature.get('type'),
            lArr = map.getLayers().getArray()[0].getLayers(),
            layer;

        // determine the layer
        lArr.forEach(function(l) {
            if (l instanceof ol.layer.Image && l.getVisible()) {
                var feats = l.getSource().getSource().getFeatures();
                for (var i = 0; i < feats.length; i++) {
                    if (feats[i] === selectedFeature) {
                        layer = l;
                        break;
                    }
                }
            }
        });

        if (layer === 'ab') {
            layer.getSource().setStyle([new ol.style.Style({
                fill: new ol.style.Fill({
                  color: hex
                }),
                stroke: new ol.style.Stroke({
                  color: '#B3B3B3'
                })
            })]);

       }
        layer.getSource().changed();
    }
};


texts: Difference between Merge Raster Layers and Build Virtual Raster in QGis?

texts: Could you explain the difference between those two functionalities. I have three rasters I want to merge into one. What is the best option ?


texts: How to get extent of rasters in folder ArcGIS

texts:
I have 100 raster image in a folder and want to get extent of all raster in ArcGIS. How to get extents of all raster in one go or write in text file?

For single raster I used this is working fine:

import arcpy
from arcpy.sa import *
elevRaster = Raster('C:/arcgis/ArcTutor/3DAnalyst/Exercise2/workspace2/congrd')
myExtent = elevRaster.extent
print myExtent


texts: Showing Select By Location dialog window when it is open but behind other ArcMap windows?

texts:
My "Select By Location" dialog window keeps falling behind other windows(such as my table window) that I have on my second monitor.

So far, the only way to fix it is to move the table window to the first monitor which allows me to get the the "Select By Location" dialog window, and then window juggle.

The "Select By Location" dialog window doesn't create an item in the taskbar which would be useful (I could right-click and move).  Or even just having the window raised to the top when I select the "Select By Location" menu-item again would help.

Are there any clever workarounds or maybe a shortcut key to deal with this issue?


texts:
I'm trying to add the Draw Control availble in with the Leaflet.Draw Plugin

I've downloaded the code and stored in htdocs and point my page towards the js and css libraries needed. The page loads ok but minus the draw control panel as seen here

Here's the code I'm using for the page.

    Leaflet Draw

<link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-0.5.1/leaflet.css" />
<link rel="stylesheet" href="http://localhost/leaflet/leaflet/Leaflet.draw/dist/leaflet.draw.css" />


texts: Connecting waypoint pairs and calculating length?

texts:
I want to make a path out of points and calculate the distances of each single line.
I have a point layer with waypoints of a research vessel. The way points have in the attribute table a date and a time and are labeled either as on effort or as off effort. Now I want to make a line layer with a path from point to point according to the recorded time. This line layer should have the information about the associated waypoints in the attribute table. I that line layer I want to calculate the distances from one waypoint to the next in the attribute table.

I have tried the extensions "points to path" and "points2one", but in both cases I only got one polyline. That is not suitable for me, because for further calculations I have to neglect the distances between off effort waypoints.

In ArcView it works with the extension "Path with distances and bearings". Is there a way I can solve this in QGIS?

This is how the point layer attribute table looks like:

This is how the attribute layer of the line layer should look like:

I am using QGIS 1.8.0 on a Mac OSX 10.6.8


texts: Road number boxes dependent on string length?

texts:
I would like to have a coloured box behind each road number on my map, matching the colour of the road it relates to and scaled to the length of the road number.

I understand this could be done with custom SVG symbols but it seems like it would be a time consuming process.

The road numbers in question are labels attached to points. Could I then work a solution where a box is created around that point as some of kind of data representation, basing the width on the string length of the road number?

EDIT: I think I may have found a way to do this with Ellipse Marker symbols. If I work it out I will post an answer.


texts:
parcelsFC = os.path.join('temp.gdb',JsonDirectoryLocation)

arcpy.JSONToFeatures_conversion(parcelsFC)

Created gdb with above code.

How can I add attributes to gdb with another JSON file?


texts: Raster to Vector reduces Segmentation accuracy too much (Earth Engine)

texts:
When I try to convert a segmented Image to a Feature Collection I loose accuracy, the new polygons are much bigger than after the segmentation.

import ee
import json
import eeconvert
import geopandas as gpd
ee.Initialize()
from IPython.display import Image

neuburg_area = ee.Geometry.Polygon([[11.135076548789584,48.761643884690855],
[11.137136485313022,48.68078409548941],
[11.25901606294974,48.6821441576126],
[11.254896189902865,48.76571740675962],
[11.135076548789584,48.761643884690855]])
neuburg = neuburg_area.toGeoJSON()

img_col = ee.ImageCollection("COPERNICUS/S2")
geometry = neuburg
img_col = ee.ImageCollection("COPERNICUS/S2")
image = ee.Image(img_col.filterBounds(geometry).filterDate('2018-05-01', '2018-07-01').sort('CLOUD_COVERAGE_ASSESSMENT').first())
image = image.clip(geometry)

seeds = ee.Algorithms.Image.Segmentation.seedGrid(10) #36
snic = ee.Algorithms.Image.Segmentation.SNIC(image= image,
  size= 32,  #32
  compactness= 5,    #5
  connectivity= 8,    #8
  neighborhoodSize=2,   #2
  seeds= seeds
).select(["B2_mean","B3_mean","B4_mean","B8_mean", 'clusters'], ['B2', 'B3', 'B4', 'B8', 'clusters']) #32 size
clusters = snic.select('clusters')
clustervis = clusters.randomVisualizer()
Image(url=clustervis.getThumbUrl({"region":geometry,"bands":["viz-red","viz-green","viz-blue"]}))

That's how the later polygons should look like (Image after Segmentation):

# Image to Feature Collection with single polygons
col = clusters.reduceToVectors()
Image(url=col.reduceToImage(["label"],"first").randomVisualizer().getThumbUrl({"region":geometry,"bands":["viz-red","viz-green","viz-blue"]}))

And that's how they look like:


texts: How to get all intersections between two simple polygons in O(n+k)?

texts:
Can someone help with a following problem: report all intersections between two simple polygons in O(n+k) where n - is a total number of polygon edges and k - is number of intersections. As a result - algorithm should output pairwise edges that interesect. It is very important that algorithm should be linear.

I've read that this can be solved in given time with linking to U. Finke and K. Hinrichs. Overlaying simply connected planar subdivisions in linear time. But I couldn't find this paper and I'm not even sure that this is my case.


texts: Recovering pixels at intersection of polygon and raster image in PyQGIS?

texts:
I am beginner in PyQGIS.

I would like to recover the pixels that intersect with my shapefile (polygon) using a loop for ... apparently, the loop is incorrect (no result).

Below the code

from PyQt4.QtCore import *
from osgeo import gdal
from osgeo.gdalconst import *

#driver = gdal.GetDriverByName ('GTiff')
Driver.Register # ()

rlayer = "F:\ MASTER GEOM \ MasterGeomatic_data \ PLEIADES \ PLEIADES_20130415_COLNAT.TIF"
rlayer = gdal.Open (rlayer, GA_ReadOnly)
provider = QgsVectorLayer ("F: try1.shp", 'try', 'ogr')
extent = provider.extent ()

collars = rlayer.RasterXSize
rows = rlayer.RasterYSize
#print (passes)
#print (rows)
val = []
for range (cols) in extent:
     print (range)


texts: Changing text element font size using ArcPy?

texts:
I need my script to change the font size of a text element if the string is too long to fit in the title section of the map.

I've looked into it using ArcPy, but haven't had any luck.

Does anyone know how to do this?


texts: Using ArcGIS for Desktop when learning via QGIS?

texts:
I'm currently trying to teach myself GIS.

I use a Mac and have limited access to ArcGIS for Desktop at a library, so I'm thinking about going the QGIS route.

My concern is that most job vacancies that are looking for a proficiency in GIS will be using ArcGIS for Desktop.

I'm not too familiar with either yet, so I'm wondering how different the program interfaces really are.

If I work on becoming proficient using QGIS, will I still be totally lost if I need to use ArcGIS for Desktop?


texts: Creating an OpenLayers map in Node.js View

texts:
I'm new to node.js and am attempting to put a simple OpenLayers map on a page. However, the map is not showing up.
Code:

app.js

var express = require('express');
var app = express();

app.set('view engine', 'jade');
app.use(express.static(__dirname + '/public'));

app.get('/', function(req, res){
    res.render('index', {title: 'Map Viewer'});
});

app.get('/map', function(req, res){
    res.render('view');
});

app.listen(3000, function(){
    console.log('Server listening on port 3000...');
});

layout.jade

doctype html
html
    head
    title= title
    link(rel='stylesheet', href='/stylesheets/style.css')

    body
        block content

view.jade

extends layout
block content
    #map

script(type='text/javascript').

var map = new ol.Map({
    target: 'map',
    layers: [
      new ol.layer.Tile({
        title: 'Global Imagery',
        source: new ol.source.TileWMS({
          url: 'http://demo.opengeo.org/geoserver/wms',
          params: {LAYERS: 'nasa:bluemarble', VERSION: '1.1.1'}
        })
      })
    ],
    view: new ol.View({
      projection: 'EPSG:4326',
      center: [0, 0],
      zoom: 0,
      maxResolution: 0.703125
    })
  });

In package.json I'm using "openlayers": "^3.15.1"
Result:
The above code is producing a blank page with the following html:
<html>
    <head>
    </head>
    <title>
    </title>
    <link rel="stylesheet" href="/stylesheets/style.css">
    <body>
        <div id="map">
        </div>
    </body>
</html>

Any idea what I'm doing wrong?


texts: Understanding relationship between heading, pitch, roll to phi omega kappa?

texts:
I am trying to find some good online reference (formulas, perhaps a script) on the relationship between heading, pitch, roll to phi omega kappa. In my current situation i receive camera orientations (and positions in world frame) in the hybrid notation (1x extrinsic then 2x intrinsic) of heading pitch roll, which refer to the local tangent plane. I have some other data i want to superimpose into the images which are provided in the UTM "world" frame so the data is not directly combinable. For this i would like to have the camera orientations in the phi omega kappa notation which correspond to the world UTM frame. There is a whole mess of different notations and definitions of angles in this topic which makes everything very confusing.

Does someone have a good reference on this?

Particulary conversion of HPR to OPK?


texts: How to get native bounding box (minx,miny,maxx,maxy) using OpenLayers?

texts:
I am trying to develop a GIS application using GeoServer and OpenLayers.  I want to get the bounds from GeoServer Native Bounding Box data of a layer instead of manually hard code:

var map;
var key = 0;
var fromNo = "";
var toNo = "";
var colorId = "";
var count = 0;
function init(){
format = 'image/png';
if(map)
map.destroy();
var bounds = new OpenLayers.Bounds(
    8568700.4305088, 1677704.30921461,
                    9195694.68892519, 2070938.3484273
);
var options = {
    controls: [],
    maxExtent: bounds,
    maxResolution: "auto",
    projection: "EPSG:900913",
    units: 'degrees',
    minZoomLevel: 0, numZoomLevels:15,maxZoomLevel:15
};
var  click_ap ;
var values = new Array();
map = new OpenLayers.Map('map', options);
    $('.selectionTable tr').each(function() {
        fromNo  =  $(this).closest("tr").find(".fromNo").val();
        toNo    = $(this).closest("tr").find(".toNo").val();
        colorId = $(this).closest("tr").find(".colorSelect option:selected").val();
        apZoomLevel1 = new OpenLayers.Layer.WMS(
                   "ap",
                "http://localhost:8090/geoserver/scmapWithSelectColor/wms??request=getCapabilities ",
                    {'layers': 'scMapQuery',
                    viewparams: 'fromNo:'+fromNo+';toNo:'+toNo+';colorId:'+colorId+';',
                    'transparent': true,isbaselayer:false},{
                       singleTile: true,opacity: .5
                    }
              );
              values.push(apZoomLevel1);
});
click_ap = new OpenLayers.Control.WMSGetFeatureInfo({
 url:'http://localhost:8090/geoserver/scmapWithSelectColor/wms',
 title: 'Identify features by clicking',
 layers: [apZoomLevel1],
 queryVisible: true,
 eventListeners: {
 getfeatureinfo: function(event) {

if (ContentFilterAp(event.text,'hos').length>0){
//alert("localitieslevel3");
    map.addPopup(new OpenLayers.Popup.FramedCloud(
        "apZoomLevel1",
        map.getLonLatFromPixel(event.xy),
        null,
        ContentFilterAp(event.text,'hos'),
        null,
        true
            ));}
        }
        }
    });
gmap = new OpenLayers.Layer.Google
(
    "Google Streets"
);
ghyb = new OpenLayers.Layer.Google
(
   "Google Hybrid",
   {type: google.maps.MapTypeId.HYBRID}
);
    map.addControl(new OpenLayers.Control.Navigation());

    map.addControl(new OpenLayers.Control.LayerSwitcher({'div':OpenLayers.Util.getElement('layerswitcher')}));
    map.addLayers([gmap,ghyb]);
    var length = values.length;
    if(length > 0)
    {
        for(i=0;i<length;i++)
        {
            //debugger;
            map.addLayers([values[i]]);
        }
    }
    map.addControl(new OpenLayers.Control.PanZoomBar({
    }));
    map.addControl(new OpenLayers.Control.Navigation());
    map.addControl(new OpenLayers.Control.Scale($('scale')));
    map.addControl(new OpenLayers.Control.MousePosition({element: $('location')}));
    map.zoomToExtent(bounds);
    map.addControl(click_ap);
    click_ap.activate();
    map.events.register('zoomend', this, function (event) {
}); }

Actually above bounds are hard coded but actually apzoomlevel1 have native bounds how to get those one? actually those bounds useful in map.zoomToExtent(bounds); based on that map zoom will change


texts: Construct Sight Lines 999998 : Unexpected Error

texts:
My problem occurs when I want to produce Construct Sight Lines. When the process is completed I get an 999998 : Unexpected Error. I see the feature class in my geodatabase. But I can't see the Sight Lines on the Map when I add the Data in the Table Of Contents.
I am using ArcMap 10.5.1
Is there anyone that knows what is going wrong?


texts: Performing Row (Element or feature)-wise Union using ArcPy with ArcGIS Pro?

texts:
I have a number of features which I would like to individually Union with itself to remove internal gaps. The best way I've thought of so far is to use cursors and iterate through, essentially performing the Union as many times as I have rows in my feature class.

Is there a tool or other approach that might be more efficient?

Note I do not want to Union everything in the feature class together - the individual results should stay separate and should have as many rows output as I originally input.

I'd like to perform this in one of Esri's products/ArcPy (ArcGIS Pro preferred over ArcMap).

Per @Keagan Allan I tried the dissolve tool but wound up with the same polygon - the interior gap was not removed. A polygon example is below. (Note the view is filtered to only show this polygon - see second image.)

The second image shows that the polygon is overlapping several other polygons which need to be kept separate (and if they have gaps, those gaps should be removed as well.)

UPDATE: The "Eliminate Polygon Part" (Advanced license) tool worked. I don't know whether something similar to this tool is available in other license levels or open source GIS, as replicability is a concern. (It may be better to use the repeated union for that reason.)

arcpy.management.EliminatePolygonPart("input_layer", "output_layer_eliminate", "PERCENT", "0 SquareMeters", 99, "CONTAINED_ONLY")

UPDATE2: I tried @FelixIP's method in arcpy (Python 2.7) and see a change in the outer border but it's no longer following the line. In the below screenshot, the greenish line is the resulting polygon from the code below, and the cyan is the original polygon.

I did the "check for gaps" and "remove gaps" separately because I wanted to examine the results of the first operation before proceeding, but theoretically could combine the two.

# Identify whether a polygon has an interior gap
with arcpy.da.UpdateCursor("infc", ["uniqueid", "SHAPE@", "HasInnerRing"]) as cur:
    for row in cur:
        print(row[0])
        for part in row[1]:
            partlist = [pnt for pnt in part]
            #print(partlist)
            if None in partlist:
                print("{} has inner ring".format(row[0]))
                row[2] = 1
            else:
                row[2] = 0
        cur.updateRow(row)

# Remove gaps, if they have been identified in previous step
with arcpy.da.UpdateCursor("infc", ["uniqueid", "SHAPE@", "HasInnerRing"]) as cur:
    for row in cur:
        if row[2] == 1:
            part = row[1].getPart(0)
            pgon = arcpy.Polygon(part)
            row[1] = pgon
        cur.updateRow(row)


texts: Finding the location of raster cell based on selected value in ArcGIS

texts:
How to get the location (X and Y) of the raster's cell based on selected value in ArcGIS.

For example, I have a raster, which has the min value 100 and maximum value 300, If I want to get the location of a pixel which has the value 200, how to reach on a particular cell in ArcGIS.

I have tried it using properties of the raster, but the function is returning the Mean, Max, Min value.


texts:
Mapbox showcases a beautiful "light" style: https://www.mapbox.com/maps/light
However, this style is not available in the style gallery: https://www.mapbox.com/gallery/
Would anyone know where/how to get the "light" style json in order to import it in Mapbox Studio?


texts: Adding a Spatialite 2.3.1 layer in QGIS 3: invalid layer

texts:
I am running QGIS 3.0.3 and I am trying to add a Spatialite layer from this link

http://www.gaia-gis.it/spatialite-2.3.1/test-2.3.zip

The program says that it is an invalid layer and that it can't be opened. I've tried different ways to opening this file.
I got this link from a course held in Udemy.
He is running QGIS 2.0.3, and I am wondering if that could be the reason for the error.


texts: Georeferenced map lines up better with different CRS – why?

texts:
I have an old map that I'm trying to georeference. The map says that it's a Bonne projection, with the origin at 34ºN, 130º,30'E.
I downloaded a DEM (in WGS84), and reprojected it using a custom Bonne projection +proj=bonne +lat_1=34.00 +lon_0=130.5 +ellps=bessel +units=m. I then georeferenced the map in QGIS (Transformation type: Projective, Resampling method: Nearest neighbor, Target SRS: the custom Bonne projection).
The georeferenced map is not horrible, however not perfect. I tried reprojecting the original DEM into EPSG: 30162, a projection that's a good fit for my target area, and it seems to line up better (although not perfectly) with my georeferenced map. The project CRS is set to the custom Bonne projection, so I understand that the EPSG:30162 DEM is being reprojected into the Bonne projection on the fly, but why does it work better than the DEM that is in the Bonne projection natively?
Tutorials I've read about georeferencing have told me to define a project CRS that is the same as the map, import a DEM in that projection, and then georeference in that CRS for best results. Can anyone shed some light on why that process doesn't seem to be delivering the best results in my case?
Did I make a mistake georeferencing? Did I make a mistake defining the Bonne projection? Or is there some other step or option that I failed to take into account?
Here is a section of the the DEM with the Bonne projection:

Here is a section of the DEM in EPSG:30162 (project CRS still set to the Bonne projection):

Here are my GCP points in the georeferencer:

I'm still new to georeferencing.


texts: sr = dsc.spatialReference AttributeError: DescribeData: Method spatialReference does not exist

texts:
I have to assign shapes without a spatial reference a coordinate system. Therefore I build a model which goes through different folders. But I need this Python script for if then else: if the shape has no spatial reference it should get one, else should be done nothing.

I nearly can take the script of the ESRI support
http://desktop.arcgis.com/en/arcmap/10.3/analyze/modelbuilder/using-if-then-else-logic-for-branching.htm

but there is this error: "AttributeError: DescribeData: Method spatialReference does not exist"

I use ArcMap 10.2.2 and Python 2.7.5 shell.
Here is my script:

# Import modules
import arcpy
import sys
import traceback

# Set local variables
prj = ""
indata = "Y:/Koordinatenumstellung_Test/mit_und_Ohne_Koo"
dsc = arcpy.Describe(indata)
sr = dsc.spatialReference
prj = sr.name.lower()

try:

   # check if indata is in StatePlane, has no PRJ, or one other than StatePlane
   if prj.find("_unknown_") > -1:
       # Set the Is Unknown parameter to FALSE, and the Is StatePlane parameter to TRUE
       arcpy.SetParameterAsText(1,"true") #The first parameter refers to the "Is Unknown" variable

       arcpy.AddMessage("Coordinate system unknown")

   else:
       # Set the Is Unknown parameter to FALSE, and the Is StatePlane parameter to FALSE
       arcpy.SetParameterAsText(1,"false")

       arcpy.AddMessage("Coordinate system is known")

except Exception as e:
    AddPrintMessage(e[0], 2)


texts: Looping through two raster folders to perform raster calculation gives ERROR 999998?

texts:
Last week I asked Looping through two raster folders to perform raster calculation? about a script to loop through two geodatabases and perform a raster calculation by matching the two rasters up based on their coordinates. This is my post for clarification:

"I have two geodatabases with rasters containing information about "treecover" and "loss". The file names are slightly different but can be matched based on the last 8 characters (the coordinates of each raster granule).See image

I can perform this simple conditional statement using either CON or RASTER CALCULATOR to retrieve all "Loss" data situated on pixels which have a value for "treecover" greater than or equal to 50, with 9999 being assigned as the output value if the statement is FALSE. The output raster is exactly what I want. I just need to do it for 137 pairs of rasters

RASTER CALC STATEMENT: con("treecover raster" >= 50, "Loss raster", 9999)"

The script so far looks like this:

    import arcpy, os, math
... arcpy.CheckOutExtension("Spatial")
... from arcpy.sa import *
... CVR = r"F:\Work\Aquifers_recalc\Final_recalc\Continent_Data\TEST\Treecover_TEST.gdb"
... LSS = r"F:\Work\Aquifers_recalc\Final_recalc\Continent_Data\TEST\Lossyear_TEST.gdb"
... # get first set of rasters
... arcpy.env.workspace = CVR
... cvr_rasters = arcpy.ListRasters()
... cvr_rasters = [os.path.join(CVR, r) for r in cvr_rasters]
...
... # get second set of rasters
... arcpy.env.workspace = LSS
... lss_rasters = arcpy.ListRasters()
... lss_rasters = [os.path.join(LSS, r) for r in lss_rasters]
...
... arcpy.env.workspace = r"F:\Work\Aquifers_recalc\Final_recalc\Continent_Data\TEST"
... for cvr_ras, lss_ras in zip(cvr_rasters, lss_rasters):
...     r1 = arcpy.sa.Raster(cvr_ras)
...     r2 = arcpy.sa.Raster(lss_ras)
...     result = Con(r1 >=50, r2, 999)
...     out_name = "loss_{}.tif".format(cvr_ras[-8:])
...     result.save(out_name)

This script worked for two geodatabases with 4 rasters in each, but when i tried it with 9 rasters in each it failed on the 6th with this error. I will need to do this for two geodatabases with 137 rasters in each.

    Runtime error
    Traceback (most recent call last):
      File "<string>", line 20, in <module>
      File "c:\program files (x86)\arcgis\desktop10.6\arcpy\arcpy\sa\Functions.py", line 269, in Con
        where_clause)
      File "c:\program files (x86)\arcgis\desktop10.6\arcpy\arcpy\sa\Utils.py", line 53, in swapper
        result = wrapper(*args, **kwargs)
      File "c:\program files (x86)\arcgis\desktop10.6\arcpy\arcpy\sa\Functions.py", line 256, in Wrapper
        ["IfThenElse", in_conditional_raster, in_true_raster_or_constant, in_false_raster_or_constant])
    RuntimeError: ERROR 999998: Unexpected Error.

UPDATE: I deleted the pair of rasters causing the issue and the script still fails. It seems like it always fails after 5 pairs of rastes, without completing the 6th.

UPDATE:By adding the loop I recieved this error, but I got 6 rasters this time:

Tool failed in loop: 6 while executing rasters F:\Work\Aquifers_recalc\Final_recalc\Continent_Data\TEST\Treecover_TEST.gdb\Hansen_GFC2015_treecover2000_00N_070W and F:\Work\Aquifers_recalc\Final_recalc\Continent_Data\TEST\Lossyear_TEST.gdb\Hansen_GFC2015_lossyear_00N_090W
Tool failed in loop: 8 while executing rasters F:\Work\Aquifers_recalc\Final_recalc\Continent_Data\TEST\Treecover_TEST.gdb\Hansen_GFC2015_treecover2000_00N_090W and F:\Work\Aquifers_recalc\Final_recalc\Continent_Data\TEST\Lossyear_TEST.gdb\Hansen_GFC2015_lossyear_00N_070W
>>>


texts: Creating dictionary for area calculation of Resolve Ecoregions (vector)

texts:
I'd like to calculate area by Ecoregion dataset and create a dictionary for each class as well.
I'd like to also select particular metadata within this dataset (e.g: BIOME_NAME)
However, I haven't found how to clip and calculate area using vector files in Earth Engine.
//** Center on the Region of Interest
Map.centerObject (roi,10);

//** Import Ecoregion
var ecoRegions = ee.FeatureCollection("RESOLVE/ECOREGIONS/2017")
                .filterBounds (roi);
print (ecoRegions);

//** Set up color
var colorUpdates = [
                  {ECO_ID: 204, COLOR: '#B3493B'},
                  {ECO_ID: 245, COLOR: '#267400'},
                  {ECO_ID: 259, COLOR: '#004600'},
                  {ECO_ID: 286, COLOR: '#82F178'},
                  {ECO_ID: 316, COLOR: '#E600AA'},
                  {ECO_ID: 453, COLOR: '#5AA500'},
                  {ECO_ID: 317, COLOR: '#FDA87F'},
                  {ECO_ID: 763, COLOR: '#A93800'},
                  ];

//** loop over all other features and create a new style property for styling
//** later on
var ecoRegions = ecoRegions.map(function(f) {
  var color = f.get('COLOR');
  return f.set({style: {color: color, width: 0}});
});

//** make styled features for the regions we need to update colors for,
//** then strip them from the main asset and merge in the new feature
for (var i=0; i < colorUpdates.length; i++) {
  colorUpdates[i].layer = ecoRegions
      .filterMetadata('ECO_ID','equals',colorUpdates[i].ECO_ID)
      .map(function(f) {
        return f.set({style: {color: colorUpdates[i].COLOR, width: 0}});
      });

  ecoRegions = ecoRegions
      .filterMetadata('ECO_ID','not_equals',colorUpdates[i].ECO_ID)
      .merge(colorUpdates[i].layer);
}

//** Use style property to color shapes
var imageRGB = ecoRegions.style({styleProperty: 'style'});

//** Visualize the Ecoregion dataset
Map.addLayer(imageRGB.clip(roi), {}, 'RESOLVE/ECOREGIONS/2017');

//** Create List
var list = ecoRegions.toList(10,0)
print ('--Resolve Ecoregion 2017 by AOI--', list)

``
How do I complete this script?


texts:
I'm trying to convert a polygon GeoJSON of countries into a line layer that contains only the borders. It should not contain the edges of the polygons that border the ocean (in my polygon dataset, nothing).

How can I in QGIS 3.8 do that? Basically, "give me the linestring of this polygon data, but only where there is a polygon on both sides of it".

I've checked similar answers and they don't cover this case. Most importantly, I actually want the data, not just set the display fill style to none.


texts: Using OSM data, pgRouting, and dummy data to develop an origin/destination travel matrix

texts:
Further to my previous question I was wondering if it is possible to add dummy node data without to do routing analysis.

ISSUE

I have a polygon layer that represents zonal attributes
I created a centroid (point layer) for each polygon
I would like to add the new points to the database with "dummy" links to the nearest highway or nearest road of a certain classification
Route from between all the nodes to develop a trip matrix

Is this possible with OSM and pgRouting?

NOTE I will add more details to this question tomorrow


texts: I have a hard time understanding what ASCII files are. I have this DEM in asc. I try to import them into ArcMAP 10.3 and they appear in my naked eyes as images however what I found on internet is ascii .txt.....what are they? Is there other extensions called ascii too? I want to use this DEM for HEC-GeoRAS and I understand that HEC-GeoRAS need DEM in TIN as data input. I have no idea what to do with this ascii file.


texts: Adding popup to geojson layer based on feature property?

texts:
I'm using leaflet and trying to add a popup to a point on the map that will return a property of the point from the geojson file. Going off this post has been helpful problem adding pop-up to GeoJSON layer in Leaflet , but it only returns the same text for each point. Here's part of the code:

function eqcolors(mag) {
          return  mag >= 5 ? 'maroon':
                  mag >= 4 ? 'red':
                  mag >= 3 ? 'orange':
                  mag >= 2 ? 'yellow':
                  mag >= 1 ? 'green':
                            'rgb(255,255,255)';
 }

function eqstyle(feature){
    return {
        radius: 8,
        fillColor: eqcolors(feature.properties.mag),
        color: eqcolors(feature.properties.mag),
        weight: 1,
        opacity: 1,
        fillOpacity: 0.8
    }
}

function onEachPoint(feature, layer){
    switch (feature.properties.mag){
        case 1.0 : return layer.setRadius(5);
        case 1.1 : return layer.setRadius(5);
        case 1.2 : return layer.setRadius(5);
        case 1.3 : return layer.setRadius(5);
        case 1.4 : return layer.setRadius(5);
        case 1.5 : return layer.setRadius(5);
        case 1.6 : return layer.setRadius(5);
        case 1.7 : return layer.setRadius(5);
        case 1.8 : return layer.setRadius(5);
        case 1.9 : return layer.setRadius(5);
        case 2.0 : return layer.setRadius(10);
        case 2.1 : return layer.setRadius(10);
        case 2.2 : return layer.setRadius(10);
        case 2.3 : return layer.setRadius(10);
        case 2.4 : return layer.setRadius(10);
        case 2.5 : return layer.setRadius(10);
        case 2.6 : return layer.setRadius(10);
        case 2.7 : return layer.setRadius(10);
        case 2.8 : return layer.setRadius(10);
        case 2.9 : return layer.setRadius(10);
        case 3.0 : return layer.setRadius(15);
        case 3.1 : return layer.setRadius(15);
        case 3.2 : return layer.setRadius(15);
        case 3.3 : return layer.setRadius(15);
        case 3.4 : return layer.setRadius(15);
        case 3.5 : return layer.setRadius(15);
        case 3.6 : return layer.setRadius(15);
        case 3.7 : return layer.setRadius(15);
        case 3.8 : return layer.setRadius(15);
        case 3.9 : return layer.setRadius(15);
        case 4.0 : return layer.setRadius(20);
        case 4.1 : return layer.setRadius(20);
        case 4.2 : return layer.setRadius(20);
        case 4.3 : return layer.setRadius(20);
        case 4.4 : return layer.setRadius(20);
        case 4.5 : return layer.setRadius(20);
        case 4.6 : return layer.setRadius(20);
        case 4.7 : return layer.setRadius(20);
        case 4.8 : return layer.setRadius(20);
        case 4.9 : return layer.setRadius(20);
        case 5.0 : return layer.setRadius(25);
        default: layer.setRadius(100);
    }
}

function pointToLayer(feature, latlng){
    var popupOptions = {maxWidth: 200};
    var popupContent = "This is some content";
    return L.circleMarker(latlng, eqstyle).bindPopup(popupContent, popupOptions);
     }

var geojson;
//map output
geojson = L.geoJson(eq, {
    style: eqstyle,
    onEachFeature: onEachPoint,
    pointToLayer: pointToLayer
}).addTo(map);

and the geojson file:

var eq =
{
"type": "FeatureCollection",
"features": [
{ "type": "Feature", "properties": { "time": "2014-02-16T20:23:35.420Z","depth": 5.0, "mag": 3.2, "magType": "mb_lg", "nst": 0.0, "gap": 52.0, "dmin": 0.47, "rms": 0.45, "net": "us", "id": "usc000mrls", "updated": "2014-02-17T16:34:23.777Z","type": "earthquake"}, "geometry": { "type": "Point"

For polygons layer.bindPopup(feature.properties.popupContent) works. Is a different method used for points?


texts: Validate a supervised classification in Google Earth Engine using ground truth data

texts:
I ran a supervised classification in Google Earth Engine (GEE) using the following script.

// add an image, bounds, filter data and sort by cloud cover

var limitiNapoli = 'PathToMyasset'

var image = ee.Image(ee.ImageCollection(sentinel)
.filterBounds(roi)
.filterDate('2019-07-01', '2019-07-28')
.sort('CLOUD_COVERAGE_ASSESSMENT')
.first());
//Map.addLayer(image, {bands: ['B4', 'B3', 'B2'], max: 0.3}, 'image');

var imageNapoli = image.clip(limitiNapoli)
// Map.addLayer(imageNapoli, "", "ImageNapoli")

print (image)

// training data and merge. Printing training data
var newfc = urban.merge(water).merge(forest).merge(crops);

// 'B2', 'B3', 'B4'

var bands = [ 'B5', 'B6', 'B7', 'B8', 'B8A', 'B11', 'B12' ];
var training = image.select(bands).sampleRegions({
  collection: newfc,
  properties: ['landcover'],
  scale : 10
});
print(training);

var trueColor = {
  bands: [ "B4", "B3", "B2"],
  min: 0,
  max: 3000 };

// train the classifier
var classifier = ee.Classifier.cart().train({
  features: training,
  classProperty: 'landcover',
  inputProperties: bands
});

// run the classification
var classified = image.select(bands).classify(classifier);

// display classification
Map.centerObject(newfc, 12);

Map.addLayer(imageNapoli, trueColor,  'Sentinel image Napoli');

  Map.addLayer(classified.clip(limitiNapoli),
  { min: 0, max: 3, palette: ['CF521F', '39C3EE', '2F9516', '53FD29']},
  'classification');
// Map.addLayer(newfc);

Map.addLayer(limitiNapoli, "", "Limiti Comune Napoli")

Export.image.toDrive({
  image: classified, //     <----- put the image you want to download
  description: "Land Cover Classification",
  scale: 10, // <------ set the pixel size in meters
  region: exportation
});

As output i have a raster in which the pixels are classified into 4 different classes: water, urban, forest and crops.
I would like to validate the forest pixels by running an accuracy assessment using some ground truth data i have from the field (polygons). How could i do it?


texts: Can QGIS projects be opened in ArcGIS for Desktop?

texts:
I am about to start a project with a number of people for a nonprofit organization. I am wanting to create the files in QGIS, but some of the people that I am working with are city officials and might feel more comfortable using ArcGIS for Desktop as opposed to QGIS.

If I were to create project files in QGIS could others open up them up in ArcGIS for Desktop? Or is there a way for me to export it to a file that could be opened in ArcGIS for Desktop? Or maybe there is file conversion plug-in that could be used?


texts: Qgis 2.18 unknown Error when Extract by location

texts:
I hope somebody can help me.

At first I have merge some vectorlayers (lineshape) without problems. Then I've create a bounding box with "polygon from layer extend" (polygon-shape) without problems too.

And then when I tried to "extract by location", I always get the following error (screenshot).

At first i've thought the reason would be the large nummer of datasets (>3 million), but i get the same error with fewer datasets. Now i THink the problem is because of the merging tool, perhaps. But I can't imagine.

I've attached a screenshot from the error text. Perhaps you can understand it better than me


texts:
I'm trying to plot a series of contours in QGIS based on dxf data exported from some specialised mining software.  The attirbute output is pretty ugly and merges a few fields together.

DXFName=C$5 DXFGroup= DXFFeature= DXFValue=50.0 DXFDescription=Major contour at level 50.000 DXFPrimitive= DXFLine=0 DXFPattern=0

I need to extract just the 'DXFValue=' value to allow categorised symbology, and labelling of the contours.

I've played around with rule based symbology using the 'regexp_substr' function, but I don't have the regex knowledge to produce what I need.

Anyone able to assist?
Thanks!


texts: Adding UTM grid (NAD83) with 1x1km spacing and label the grid lines?

texts:
I'm currently running QGIS 2.18.4 (the latest release).

I'm trying to add a UTM grid to my maps, specifically the NAD83 reference system for zone 15N, that hopefully can show UTM coordinates and with 1x1km spacing. Is there an easy way of doing this? I can't find a plugin for it, and I can't find any wms server with utm grids either.

I've gone to view > decorations > grid and enabled a grid at 1000m spacing, I think its correct, however, I can't get line labels for it.

I've searched for a while and I feel like I have a very specific question as I can't seem to find anything that works.


texts: Running a QC check for gaps between vertices on a polyline

texts: I'm checking polylines to make sure power pylons were correctly extracted from imagery.  Is there a way to find all of the gaps or missed pylons in the database?  Currently, I'm manually making a point file from the vertices on the line and checking inconsistencies at a high level.  I thought of making a buffer and looking the areas that didn't overlap, but the power pylons are gapped at different distances so it would take forever to run several different buffers for every line in the database.  Any other ideas for a more efficient process?


texts:
I have a shapefile where there are polygons perfectly fitting, I created  contour lines from all of the polygons so I have a shape containing only the lines. My problem is as the polygons were fitting, most of the lines cover each other. I cannot merge them because I need to keep the attributes of the lines.

What is the best way to get rid of the unnecessary second lines? ArcGIS or QGIS usage is preferred.


texts:
I am trying to select grids (polygon geometry) that intersect with a building layer. I can do it by using "Select by location" tool and then I generate centroids from the selected grids. There is no problem in these steps. Later I generated a model with the sextante modeler but it doesnt work properly. Select by location tool dosn't work in sextante modeler. It cannot select features but interesting, I can run this tool perfectly directly from sextante toolbox!

So I need help on how I can extract the grids that intersect with another layer.

EDIT : Sextante has been renamed to Processing as of QGIS 2.0


texts: On QGIS, with a processing script, how to refresh the canvas?

texts:
In very simple beginners terms... how do I refresh the QGIS map canvas (or a specific named layer) as part of a processing script?

I have an extremely simple processing script (below) set up just to make it really easy to change a project variable. One of my map layers has a conditional styling rule set up based on the project variable. After running the script I currently have to manually refresh the canvas (by clicking the refresh button or moving the canvas) to update the map display. For purposes of simple elegance and efficiency I'd like the script to do this automatically.

Please answer in the simplest terms - I'm well out of my depth here.

The script is as follows:

##variable_for_styling=number 0

from qgis.core import *
from PyQt4.QtCore import *

QgsExpressionContextUtils.setProjectVariable('NameOfMyProjectVariable',variable_for_styling)


texts: What would be appropriate output pixel size while converting contours into DEM/DTM?

texts:
I have contours having 10 meter (32.8 feet approximately) interval. Now, I am using these contours to Interpolate a hydrologically correct raster surface in a bid to calculate slope and other topographic attributes. The 'Topo to Raster' tool available in ArcGIS will be utilized to do so.

Keeping in view of my 10 meters contour interval, what would be optimal output cell size should I choose to get maximum from these contours.

Is there any rule of thumb to get pixel size based on the interval of contours?


texts:
I am building a tool to calculate navigation charges for aircraft flying within Europe.
This function is administered by an organisation called EUROCONTROL.

EUROCONTROL have an application of their own (I have my reasons for wanting to recreate the functionality). Within that application is a file containing some sort of geo file in .txt format with the boundaries for the various zones.

The UK extract is below as an example, EG being the identifier for the UK.
Assuming I'm looking at what I think I'm looking at, its a closed polygon.

Is anyone able to identify the format, and even better, explain it?

44 3259 -198 1 0 660 0 0 0 0 0 0 0 0 EG
3274 -600
3300 -600
3360 -600
3420 -600
3480 -600
3540 -600
3600 -600
3660 -600
3660 -540
3660 -484
3660 -480
3660 -420
3660 -390
3660 -360
3660 -300
3660 -275
3660 -240
3660 -180
3660 -120
3660 -111
3660 -60
3655 0
3600 0
3420 300
3300 300
3090 120
3067 120
3060 88
3040 88
3000 -15
3000 -120
2930 -480
2940 -480
2970 -480
3000 -480
3060 -480
3140 -330
3235 -330
3265 -490
3320 -415
3325 -440
3320 -495
3285 -540
3274 -600

Edit
I used Pythons svgwrite to draw the number pairs as coordinates and got the below. I'm trying to find something that matches it.


texts: Add QGIS Server WMS layer to GeoExt template

texts:
How to add QGIS Server WMS layer to GeoExt template (OpenGeoSuite SDK)?
I used the following code, but unsuccessful display layer

    qgis: {
        ptype: "gxp_wmssource",
        url: "http://localhost/cgi-bin/qgis_mapserv.fcgi",
        title: "Local QGIS MapServer",
        version:"1.1.1"
    },

    {
        source: "qgis",
        title: "QGIS Server Layer",
        map: "/home/hast/UA/qgis-mapnik.qgs",
        layers: "buildings",
        opacity: 1,
        group: "overlays"
     }

I didn't find any information about it....


texts: PyQGIS composer doesn't change atlas feature zoom

texts:
I have found this code in this site and with the some changes I want to create a automate print PyQGIS composer using atlas.

The code looks like it should work, I have an error.

The code doesn't change my atlas feature layer zoom, any idea why?

Stay in the default zoom in my .qpt but my atlas change features correct

Here is the code:

from qgis.PyQt.QtXml import QDomDocument
from qgis.utils import iface

def quick_export():
    alayer=iface.activeLayer()
    # Add all layers in map canvas to render
    myMapRenderer = iface.mapCanvas().mapRenderer()

    # Load template from file
    myComposition = QgsComposition(myMapRenderer)
    template = 's4.qpt'

    myFile = 'C:/Users/name/Desktop/tt/s4.qpt'
    myTemplateFile = file(myFile, 'rt')
    myTemplateContent = myTemplateFile.read()
    myTemplateFile.close()

    myDocument = QDomDocument()
    myDocument.setContent(myTemplateContent)
    myComposition.loadFromTemplate(myDocument)

    # Get map composition and define scale
    myAtlasMap = myComposition.getComposerMapById(1)
    myAtlasMap.setAtlasDriven(True)
    #myAtlasMap.setNewScale(int(scale))

    # Setup Atlas.
    #myAtlas = QgsAtlasComposition(myComposition)
    myAtlas = myComposition.atlasComposition()
    myComposition.setAtlasMode(True)
    myAtlas.setCoverageLayer(alayer)
    myAtlas.setComposerMap(myAtlasMap)
    #myAtlas.setFixedScale(True)
    #myAtlas.fixedScale()
    myAtlasMap.setAtlasScalingMode( QgsComposerMap.Auto )
    myAtlas.setHideCoverage(False)
    myAtlas.setFilterFeatures(True)
    myAtlas.setFilterFeatures(True)

    # Generate atlas
    myAtlas.beginRender()
    for i in range(0, myAtlas.numFeatures()):
        myAtlas.prepareForFeature( i )
        jobs = r"C:/Users/name/Desktop/test/"
        output_jpeg = jobs + str(i)+ "_BMS_plan.jpg"
        myImage = myComposition.printPageAsRaster(0)
        myImage.save(output_jpeg)
    myAtlas.endRender()

quick_export()


texts: Converting long, lat to tile of NASA GIBS - WGS 84 / Lat-lon / Geographic (EPSG:4326)?

texts:
I'm following GIBS API for Developers and I need to convert long, lat of POI to col and row of tile for given zoom. For example to get right tile in endpoint like:

https://gibs.earthdata.nasa.gov/wmts/epsg4326/best/MODIS_Terra_CorrectedReflectance_TrueColor/default/{date}/250m/{zoom}/{row}/{col}.jpg

I'm using python and already tried few libs:

https://github.com/geometalab/pyGeoTile
https://github.com/jswhit/pyproj
https://gist.github.com/tucotuco/1193577

But didn't get right row, col values. As I see the main problem that all of them are basing on idea that lat map direct to y and you should just divide to 180˚ and mult to 2^zoom. But actually it doesn't work and it is very easy to check.

New Zealand location is about 41° 174°. And for zoom = 3 you should calc row that way:

row = int((90˚ + 41˚) / 180˚ * 2 ^ 3) ≈ int(5,8222) = 5
col = int((180˚ + 174˚) / 180˚ * 2 ^ 3) ≈ int(15.73333) = 15

So here you should see New Zealand:

https://gibs.earthdata.nasa.gov/wmts/epsg4326/best/MODIS_Terra_CorrectedReflectance_TrueColor/default/2016-09-03/250m/3/5/15.jpg

but actually we are out of tiles. Image:

Response:

<ExceptionReport xmlns="http://www.opengis.net/ows/1.1" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://schemas.opengis.net/ows/1.1.0/owsExceptionReport.xsd" version="1.1.0" xml:lang="en">
  <Exception exceptionCode="TileOutOfRange" locator="TILECOL">
    <ExceptionText>TILECOL is out of range, maximum value is 9</ExceptionText>
  </Exception>
  <Exception exceptionCode="TileOutOfRange" locator="TILEROW">
    <ExceptionText>TILEROW is out of range, maximum value is 4</ExceptionText>
  </Exception>
</ExceptionReport>

Why 4 and 9 are here?

As well for zoom = 0 it should be a whole planet on one tile, but not in NASA GIBS:

https://gibs.earthdata.nasa.gov/wmts/epsg4326/best/MODIS_Terra_CorrectedReflectance_TrueColor/default/2016-09-03/250m/0/0/0.jpg :

so you could ever ask row=0,col=1 tile:

Why it is so strange and where am I wrong?


texts: An unexpected error occurred processing the request when requesting the layer definition Url

texts:
When I publish my MXD in ArcGIS 10.0, the service url which returns the layer definitions in json format returns an error.

{"currentVersion":10.05,"error":{"code":500,
"message":"An unexpected error occurred processing the request.","details":[]}}],"tables":[]}

The url looks something like this:

http://myserver/arcgis/rest/services/myservice/MapServer/layers?f=json

The MXD contains 150-200 layers. If I remove a few of the layers near the bottom, it renders correctly. These layers at the bottom work fine in other MXDs.  If I remove the bottom layers to get it working, but try adding completely different layers, it fails.

I found the errors being logged in a file called restlog.txt. So it looks like a bug in the code that renders the layer definitions.

ERROR:Index was outside the bounds of the array.
   at ESRI.ArcGIS.REST.MapServer.MapLayerRequestHandler.get_LayerDefinitionExpression()
   at ESRI.ArcGIS.REST.MapServer.MapLayerJsonFormatter.WriteJsonResponse()
   at ESRI.ArcGIS.REST.JsonFormatter.WriteResponse()

Here is a public Url (not related to me) illustrating the error:
http://epamap5.epa.gov/ArcGIS/rest/services/EMEF/TRINAT_2012/MapServer/layers?f=json

Is there a documented or known bug in the 10.0 Arcgis server code that renders the layer definitions.

Any idea how to get around this error and stay on 10.0?


texts: Opening LIDAR DEM in QGIS with several .adf files

texts:
I have a LIDAR DEM dataset of a bay. What I got is a single folder with these files

I found out that the w001001.adf is the one that contains all the raster data and that the dblbnd, hdr, prj.sta provides metadata, however, there are several other (z001001.adf, z001002.adf, z001003.adf, z001003x) files which are as heavy as the w001001.adf and I dont know what they are for? I'm pretty much a beginner with raster data.


texts: ArcGIS polyline dangling endpoints count

texts: Is there a way for count the number of dangling endpoints of a polyline? I need that because I would like to define a feature complexity parameter. For example, on the image below, the number of dangling endpoints of the selected feature is 6


texts: How to check whether an optional FeatureSet parameter was set by the user in a geoprocessing script tool?

texts:
I have created a custom geoprocessing script tool which contains an optional FeatureSet parameter.

When unpacking the parameters in the script, the FeatureSet is defined as a RecordSet object, even though no input was received. It seems to set a default value if it was not entered.

I need to check if this parameter was specified by the user or not.

For example, when I run this through ArcMap without a FeatureSet parameter the following occurs:

exportAOI = arcpy.GetParameter(3)
arcpy.AddMessage('AOI is: ' + str(exportAOI))
# AOI is: <geoprocessing record set object object at 0x14133578>

I was expecting exportAOI to be None.

How can I determine whether the user actually set this type of parameter?


texts: dynamically added class not retained on leaflet marker

texts:
I am trying to add an active class to the marker that was clicked by the user. For this, I am storing the clicked marker's icon in a variable and doing the following

let activeMarker = e.target._icon;
activeMarker.classList.add('leaflet-marker-active');

Since we are showing more than 500 points, I use the map.flyTo function to flyTo the selected marker's position and then use map.invalidateMapSize() as well.

With all this, once the map is completely loaded, the active class added disappears from the marker. The markers keep on refreshing and redrawn, as a result of which the active marker class is not retained.

How can we add class to leaflet marker, such that it is not removed even when we move the map or drag the map?


texts: How to store and serve data and its progress over time in GeoServer?

texts:
I recently started working with GeoServer and stored some shapefiles. I successfully can retrieve the data as JSON. The data I want to maintain in the database changes over time. I do not want to update the values but add a new revision with a time stamp of the very same data. Further, in the query reply I would like to see all revisions of one feature being returned.

Question: How can I store and serve various revisions of the same data with GeoServer?

Update: I am also interested in other software that fulfils the requirement to serve geospatial data with the versioning feature.


texts: Which default formats does GeoServer REST API datasource PUT endpoint accepts?

texts:
I'm using the following endpoint to create datasources in Geoserver:

PUT /workspaces/{workspaceName}/datastores/{storeName}/{method}.{format}

I have been able to create vector datastores from an already available .shp at Geoserver data directory, using the method as external and  format option as .shp

My issue is that the documentation only gives the .shp format example, and i can't seem to find whats the correct option for geotiff format.

So, does any one know the list of available options for format and the correct way to use it for geotiff data?


texts: Nearest line to polygon based on equal attribute

texts:
I have a polygon shapefile and a line shapefile with several thousand features each. I want to get:

The nearest distance between a polygon and a line that have the same attribute (e.g. ID). Example here (I want to get the green line as an output).

I have played around with various combinations of all distance tools in QGIS, like Distance to nearest hub, Join attributes by nearest, v.distance.

v.distance SHOULD be able to do what I want with 'upload' = 'to_attr', but I am getting no output.


texts: Plot BNG grid data on Matplotlib Basemap

texts:
I would like to plot UK population density data in Matplotlib Basemap.

I got data in ESRI ASCII using BNG OSGB63 from
https://catalogue.ceh.ac.uk/documents/61f10c74-8c2c-4637-a274-5fa9b2e5ce44

and then used the following code from
https://stackoverflow.com/questions/23264404/gridded-data-basemap-mismatch-what-have-i-done-wrong#comment75355592_23264404

which in my case looks like this:
https://pastebin.com/52GH9FnH

Here is an excerpt:

# Define a BNG coordinate grid for the data, based on information in the header:
   # LLcorner eastings = -200000
   # LLcorner northings = -200000
   # grid box size = 5000
   # All units in metres
lle = -10000 + 2500 # because we want to plot the box point at the centre of each box, not its lower left corner.
lln = 0 + 2500
delta = 1000
nn,ne = grid.shape

ns = np.arange(0.,nn,1)*delta + lln
es = np.arange(0.,ne,1)*delta + lle

# convert to mesh of eastings/northings.
eastings,northings = np.meshgrid(es,ns)

# convert mesh to lon/lat (wgs84) using pyproj:
bng = pyproj.Proj(init='epsg:27700')
wgs84 = pyproj.Proj(init='epsg:4326')
lons,lats = pyproj.transform(bng,wgs84,eastings,northings)

lats = lats[::-1] # because otherwise the map is upside down (differences between imshow and pcolormesh)
fig = plt.figure(figsize=(15, 15))
# Get a basemap of the UK
m = Basemap(projection='merc',llcrnrlon=-15 ,llcrnrlat=47.8,urcrnrlon=7.2 ,urcrnrlat=62.5,resolution='h')

# project the lon/lat grids onto map coordinates
x,y = m(lons,lats)

# Draw the map!
m.drawcoastlines()
m.pcolormesh(x,y,grid,cmap=plt.cm.YlOrRd)
plt.show()

This is the result:

Any idea, what I'm doing wrong and how to make the basemap outline and my grid data align?


texts:
I've got a website where I'd like to insert maps with some placemarks (restaurants, points of interest, ...). I'm tired of Google Maps and like to use ESRI basemaps. Is it possible to do this?

I'd like to create these maps directly from my webpages (with JavaScript) and without embeding ArcGIS online maps (I suppose this would require to create the map in 'arcgis.com' before).

Is it possible with any ESRI JavaScript API?


texts: XY Points incorrect on map (ArcGIS Javascript API 4.7)

texts:
Hopefully something silly I'm overlooking. I started building my map with this sample code

But once I started using my own data, it isn't appearing correctly on the map. Example: These coordinates [38.90309, -77.05070999] in Washington DC show up somewhere in Antartica, as do the other US coordinates.

I've read a little about specifying the coordinate system that is being used, and it seems that maybe that is the problem. But I can't find any documentation about how to do that in the JavaScript API.


texts: shapefile and multiple lines of census data QGIS

texts:
Possible Duplicate:
  How to join multiple records to single feature?

I am hoping to be able to link in the attribute table more than one line of data from a census entry per polygon.  for example 8 separate lines of names and demographic info.  I attached a spread sheet file to the attribute table but ended up with this in one continuous row or line; with "1" representing the id number for the polygon.

1 Peltier  Jean 26 farmer 1 Gautier Louis 28 laborer

in one line,

I want two or more lines like this;

1 Peltier Jean 26 farmer
1 Gautier Louis 28 laborer

Is this possible?  I want to keep demographic data tied to each person not within a continuous single line or row.


texts: Extracting raster pixel data in QGIS for project analysing deforestation rates

texts:
very new to QGIS and GIS software in general. I'm using QGIS to obtain spatial data from raster data for a spatial econometrics project I'm working on.

I have 14 .tif images displaying change in forested area from a given year across Australia. The 14 images combine to form a mosaic of the country, and each pixel value displays the year in which deforestation occurred for a given pixel area (0 if no deforestation, 1 for year t+1 etc).

I need to separate the pixel data so that I can analyse deforestation changes in a given area, for each year. How do I go about separating the data? Currently, I'm leaning towards creating shapefiles for each pixel value in each raster and combining all shapefiles for a given year, before using these shapefiles to calculate deforestation area proportions in administrative areas. Am I going about this the right way?


texts: Downloading Google Maps municipal geometry

texts:
When searching for a city on google maps, google is that smart it displays a geometry line on the boundary off the municipal.
See this link: http://g.co/maps/23dz5

Is there a place in the google maps api where I can download this geometry?


texts: Labeling polygon with several parts with lines converging on one label using ArcMap

texts:
I'm trying to label in ArcMap 10.6
a polygon with three parts. With the labeling manager, and I've arrived to label the three parts with three labels (image 1).

It doesn't look horrible but it takes too much place on the map. I want to label the three parts of the polygon with three lines converging on one label (image 2).

Is this possible to do that with the label manager in ArcMap or is there a "rudimentary" solution?


texts: QGIS 2.0 Plug In Repository - How to Set Parameters to qgis=1.8

texts:
A number of very useful QGIS plug ins available in 1.8 are not available in 2.0. Is there a way to set the parameter page in "Repository Details" to 1.8 so I can access the 1.8 plug ins and check functionality?

The option is grayed out as shown in the screenshot below.


texts: Generating "theme gradient" legend graphics from MapGuide Open Source

texts:
In MapGuide Open Source, I've defined 10 layers, along with a theme for each. As I understand it, a theme represents how different data values are represented. For each theme, I picked two colors and created a gradient from the lowest data value to the highest.

On my map, I use WMS to render the layers. They render properly, with theming visible. However, I'd also like to generate a legend graphic for the theme.

I've tried this request: /mapguide/mapagent/mapagent.fcgi?version=1.0.0&operation=getlegendimage&THEMECATEGORY=-1&TYPE=-1&format=image/PNG&SCALE=44076.203877482876&layerdefinition=Library://OneOfMyLayers.LayerDefinition. However, that only returns a 16x16 pixel representation of a single point.

Instead, I want a graphic that displays the theme gradient. This is a standard WMS function, yet I haven't found any way to do it out of MapGuide Open Source.

Is there any way to generate such a theme gradient legend with MapGuide Open Source?


texts:
Is there an alternative source for 8-digit (vector shapefiles) US watersheds?

The USDA link is down and they are not responsive to my query for getting the data:

https://gdg.sc.egov.usda.gov/


texts: Calculate change in canopy height over 10 years?

texts:
I want to know how many "tall" trees (greater than 30 m) fell down between 2005 and 2015.  I have LIDAR imagery from 2005 and 2015 from an area where many tall trees are known to have fallen due to ice storms.

I created 2 tree height rasters, one for 2005 and one for 2015, by subtracting the ground height from the canopy height (I originally had 4 rasters).  What I want to know is the total number of pixels that declined by at least 20 m in the period 2005-2015.

I've tried a variety of approaches with this and haven't quite figured it out.  Does anyone have any good hints for me?

(using ArcGIS 10.2.1)


texts: Can a Geometric network be hosted in a Spatial Data Service?

texts:
I'm wanting to deploy a geometric network to Sql Azure, where clients can trace through it.

The Spatial Data Services Deployment Utility allows me to deploy each featureclass from a geometric network into its own Azure Spatial Data Service.  However, I don't see a way to perform traces.  Since Sql Server Spatial Tools still can't be deployed to Azure I would not even be able to roll my own tracing.

Is there a way to host a Geometric network in Sql Azure with tracing functionality exposed?


texts: Retrieving tiles in EPSG:4326 format. Is there a service for that like CloudMade?

texts:
Is there a service that can provide custom tiles in the EPSG:4326 projection? I've looked at Cloudmade and it would be the ideal solution but their tiles are using Spherical Mercator.

Basically the end result I'm looking for is to create a square grid overlay on a map. Like this : http://openlayers.org/dev/examples/graticule.html

So retrieving the maps via WMS would be a good choice but I need the map to be custom styled like you would be able to do in cloudmade. But maps like cloudmade and google maps seem to distort the world and so a square grid won't be possible on their tiles if i wanted to create the grid by equal increments of latitude and longitude.


texts: Splitting Polygon based on overlaid vertices using ArcGIS Desktop?

texts:
Is there a way to split a polygon using overlaid vertices using ArcGIS Desktop?

I am trying to determine parcels that may be used for future development. I have both building and parcel shapefiles. The parcels I'd like to split based on the building vertices. An example is below:

From this, I'd like to split the parcel based on the building vertex with the maximum distance to the parcel boundary like so:


texts: How to execute a clip after a buffer in a custom script tool

texts:
I'm working on a script which

Buffers roads in a layer, then
clips the buffered roads from an area

Here is my script thus far:

import arcpy
import os

from arcpy import env

env.workspace = "C:/Code/Module6"

arcpy.env.overwriteOutput = True

Roads_to_Buffer = arcpy.GetParameterAsText(0)
Output_Road_Buffer = arcpy.GetParameterAsText(1)
Road_Buffer_Distance = arcpy.GetParameterAsText(2)
Focus_Area = arcpy.GetParameterAsText(3)
Clipped_Buffer_Roads = arcpy.GetParameterAsText(4)

# step 1 create a buffer around the roads
arcpy.Buffer_analysis(Roads_to_Buffer, Output_Road_Buffer, Road_Buffer_Distance)

# Step 2 Clip the buffered roads from the focus area
arcpy.Clip_analysis(Output_Road_Buffer, Focus_Area, Clipped_Buffer_Roads)

here is a screenshot of my parameters:

However, I get the following error:

Traceback (most recent call last):
   File "C:\Code\Module6\Module6test.py", line 19, in <module>
    arcpy.Clip_analysis(Output_Road_Buffer, Focus_Area, Clipped_Buffer_Roads)
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\analysis.py",     line 56, in Clip
raise e
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Input Features: Dataset C:\Code\Module6\Bufferroads does not exist or is not supported
Failed to execute (Clip).

Failed to execute (Module6Test).

Can anyone help me fix my script?


texts: Good projections for Middle East and North Africa in D3?

texts: I'm a data visualizer who's mostly worked with US data, and I'm familiar with AlbersUSA project. However, I'm going to be doing a series of graphics on the Middle East. Any suggestions or interesting bl.ocks you've come across?


texts: Field Equipment: Geneq SXBlue II with iPad and ArcGIS Collector app

texts: My company just recently purchased a Geneq SXBlue II GNSS receiver to be pair with an iPad and the Collector for ArcGIS app.  I have heard great things about the submeter accuracy of the unit; however, the biologists on my team want the ability to see the PDOP per point and want the GPS to only capture the point if it meets a certain PDOP threshold, just like our Trimble GeoXT unit does.  I believe there is an ability to change the PDOP threshold on the GPS receiver, but does anyone have experience with how or if the Collector app lets you know if the threshold has not been met? Does it simply not take the point? Give you a message?


texts: Need a specific field in FME KMZ to shapefile conversion

texts:
I am converting a KMZ to shapefile and need to see a field called "Type Name", but it is not in any of the shapefiles' attribute tables.

I have every name exposed in Format Attributes in the Reader and the reader's User Attributes are blank:

All the fields that show up in the attributes table are listed in the writer's User Attributes tab and I have tried it with Automatic and Dynamic selected:

Also, it's strange that none of the fields have any real data. Every record in most of them is either blank or 0, and the others have meaningless constants.
Am I doing what I need to do to in the conversion to see all the fields in this KMZ?


texts: How to do geometric correction of remote sensing data in R?

texts:
My challenge is to replace most of the remote sensing commercial tools with R code (I'm teaching assistant at the University and I would like to move to open source R code to show the added values to the students). I'm currently stuck as I cannot find a way to perform geometric correction of remote sensing data. The exercise should be:

Select Ground Control Points on both the georeferenced image and the
image to correct
Use polynomial mapping for image correction (students should try to test different  polynomial order and calculate the RMSE for each step)
Warp the image in a defined projection using the interpolation method (either nearest neighbour, bi-linear interpolation, cubic convolution)

The 2 last steps could be processed using GDAL but it is not really interactive with R.

Does someone has any idea or experience with such processing?


texts: Is there any method to generate the minimum number of intersecting grid squares to a polygon?

texts: I know I can generate a group of fixed sized intersecting square grids for a polygon using the bounding box of the polygon and generate gird squares on that and keep those, which are intersecting with the polygon. But this is not the best minimal fit of the grid squares. So is there any method to generate the minimum number of intersecting grid squares to a polygon?


texts: add marker to map from list of coordinates

texts:
Forgive me if this question has an incredibly obvious answer. I'm very new to ESRI and I'm trying to add markers to my map via gps coordinates. I have not found an example that does something this simple. When I try to adapt the more complicated ones, my code breaks.

Would this be an image layer? a SimpleMarkerSymbol? a PictureMarkerSymbol?

What would really help is an extremely streamlined example of a marker on a map so I can see how it works and then adapt my code from there. It seems all the samples I'm finding have all kinds of fancy extras that I'm not familiar enough to play with.

EDIT

I've actually got some markers on my map! yay!
But now, I'm having a tough time with the coordinates showing up where I expect. As you can see from the screenshot, the highlighted marker uses gps coordinates (42.10560 -88.380801) that should put it somewhere in northern Illinois but for some reason, its all the way off the east coast. The second marker shows up in Illinois but I guestimated the coordinates for that one.

I think it has something to do with the SpatialReference class but I've tried changing the WKID attribute to various things that some of the mapping experts in my office have suggested but nothing seems to change anything.

Here is my code for playing around with.

<!DOCTYPE html>

    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">

    <meta name="viewport" content="initial-scale=1, maximum-scale=1,user-scalable=no">
    <title> Topographic Map with graphics</title>
<link rel="stylesheet" href="http://js.arcgis.com/3.7/js/dojo/dijit/themes/claro/claro.css">
<link rel="stylesheet" href="http://js.arcgis.com/3.7/js/esri/css/esri.css">

I must be missing something. Please help!


texts:
When I try to use GRASS commands in QGIS, I get the following error message:

Missing dependency.This algorithm cannot be run :-(  This algorithm
requires GRASS to be run. Unfortunately, it seems that GRASS is not
installed in your system, or it is not correctly configured to be used
from QGIS Click here to know more about how to install and configure
GRASS to be used with QGIS

I'm not sure if this is related, but in processing options and configuration I have the following paths:

Grass Folder: C:/PROGRA~2/QGISCH~1/apps\grass\grass-6.4.4
MSYS Folder: C:/PROGRA~2/QGISCH~1/apps\msys

Should I reroute these folders? Which are the correct folders that Grass folder and MSYS folder should point towards?


texts: Can't use GDAL with Python for opening *.GML file (vector file) from Sentinel-2 (L2A type) imagery *.zip archive

texts:
I'm trying to use GDAL 3.0.4 to open some of the *.gml files directly from the Sentinel-2 *.zip archive with following lines:

gml_path = f'GML:/vsizip//s2zip_path/path_to_GML_file'
gml_ds = gdal.Open(gml_path, gdal.GA_ReadOnly)

but I receive an exception like:
"Exception: GML:vsizip//zip path//path to file: No such file or directory

I doubled checked several times the used path and it is fine. I also made a try without mentioning the driver, like this:

gml_path = f'/vsizip//s2zip_path/path_to_GML_file'
gml_ds = gdal.Open(gml_path, gdal.GA_ReadOnly)

and it seems now that GDAL saw the file but it through another exception like 'not recognized as a supported file format.'

I would like to do fix that for a script.

What can be wrong ?


texts: (Google Map) auto update of info window content like time

texts: Is there anyway to frequently auto update content of info window like live  local time of particular city like New York?


texts: What are QGIS equivalent functions for ArcPy's Update/Delete Row/Field?

texts:
I am trying to reprogram some scripts from ArcPy to QGIS (1.8 or 2.0) and there are some simple functions that I want to be able to redo but unfortunately documentation in QGIS is lacking in certain areas.

Namely the three most important for me are:

Add Field - Add field

arcpy.AddField_management(Feature, "ID", "SHORT")

Calculate Field Management - Update that field

arcpy.CalculateField_management(Feature,"ID","!FID!")

Update/Delete Rows - Update/Delete rows based on condition (without copying the shapefile)

keep = ["Bob","Janet","John","Mike"]
Counter = 0
rows = arcpy.UpdateCursor(Feature)

for row in rows:
    if row.Name in keep:
        row.ID = Counter
        rows.updateRow(row)
    else:
        rows.deleteRow(row)
    Counter += 1

Now I can iterate through each feature in QGIS using SEXTANTE and obtain its geometry which I should be able to rewrite into a new shapefile and thereby update/delete a row or field. Starting with something along the lines of...

layer = st.getobject(Polygon)
features = st.getfeatures(layer)
for f in features:
    f.geometry().asPolygon()

but I can't find a simple solution for those functions mentioned above?


texts: Create callout from label to edge of polygon with Geometry generator with QGIS

texts:
With make_line(centroid ($geometry), make_point("x","y")) you can draw a line between the label and the center of an polygon, like this:

Is there a way to draw a line just between the border/edge of the polygon and the label?

The polygons have a transparent fill, so cannot just be drawn over the line.

Yes, the polygons are separate from each other or located at the edges. Is there any possibility to change the ankerpoint from endpoint of the polygon to the nearest point between the edge/border of the polygon and the label?

So that i can adjust the Label all around the polygon with a connection line between label and polygon edge/border?


texts: Select by attribute using date field values two days before machine date of ArcGIS for Desktop?

texts:
I have a date field within a file geodatabase of ArcGIS 10.2 for Desktop that I am accessing using ArcMap.

An example of a value in this field: 19/11/2014 08:49:56

I want to select all rows that the date is before two days from machine date.

An example : if the machine date is 24/11, I want to select all rows before 22/11.

I need query to select rows before two days depends on machine date.

Is it possible?


texts: Surface distance on DEM between two coordinates

texts: Given two coordinates on a digital elevation model (DEM), how do I compute the actual distance traveled between the two locations, assuming a straight line route? How will changing the resolution of my DEM affect the results?


texts: Python script for filling  values with previous value in one field

texts:
I have a table (inside a File Geodatabase) that looks like this:

I need the  values of Owner attribute to be automatically filled with the previous non-Null value. As Such:

Is there a Python script, usable in Field Calculator, that can achieve this?


texts:
I need to add in the QGIS a vector layer using a getMap request. I have added other layers too which doesn't require any authentication I am able to do easily. But for one layer I need to pass the header in QgsRasterLayer. I have done a lot of searches but I am not able to find it. Can anybody provide me any solution?
Here is my sample code :
    layer = QgsRasterLayer(path, baseName)
    QgsMapLayerRegistry.instance().addMapLayer(layer)
    if layer.isValid() is True:
        print "Layer was loaded successfully!"
    else:
        print "Unable to read basename and file path - Your string is probably invalid"


texts: Drawing new features with QGIS2web plugin

texts:
I would like to draw shapes with QGIS2web plugin
According to the query here:
Add OpenLayers3 controls to qgis2web project
I should be able to develop the qgis2web.js file further behind the ol.Map section.
I found a good drawing shapes test for Openlayers here:
https://cdn.rawgit.com/HamHamFonFon/ol3-drawFeatures/82f29a3f/examples/basic_use.html
from where I picked up the relevant files. I did the following steps:

Linked the essential scripts and csses into my index.html page:
<link rel="stylesheet" href="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.6/css/bootstrap.min.css">

<script type="text/javascript" src="https://cdnjs.cloudflare.com/ajax/libs/jquery/2.2.0/jquery.min.js"></script>

Added, and next linked the following script and css files:
    <link rel="stylesheet" href="./resources/ol3-drawcontrols.css">
    <script src="resources/ol3-controldrawfeatures.js"></script>

Placed the following section:
  var vector_draw = new ol.layer.Vector({
source: new ol.source.Vector(),
style: new ol.style.Style({
    fill: new ol.style.Fill({
        color: 'rgba(255, 255, 255, 0.2)'
    }),
    stroke: new ol.style.Stroke({
        color: '#ffcc33',
        width: 2
    }),
    image: new ol.style.Circle({
        radius: 7,
        fill: new ol.style.Fill({
            color: '#ffcc33'
        })
    })
})
 });

in the layers.js file

Situated another section:
    var options = {
   "style_buttons" : null, /** @var {string} glyphicon|default */
   "draw": {
    "Point": true,
    "LineString": true,
    "Square": true,
    "Circle": true,
    "Polygon": true
}
 };

before the
    var layerSwitcher = new ol.control.LayerSwitcher({tipLabel: "Layers"});

in the qgis2web.js file

At the finish I placed the following line:
 var buttonsDrawControls = new ol.control.ControlDrawFeatures({vector_draw, options});
 map.addControl(buttonsDrawControls);

but the map is gone.
The console says:
Uncaught TypeError: ol.inherits is not a function
at ol3-controldrawfeatures.js:209
at ol3-controldrawfeatures.js:16
at ol3-controldrawfeatures.js:18
which refers to:
at ol3-controldrawfeatures.js 209 ol.inherits(ol.control.ControlDrawFeatures, ol.control.Control);
at ol3-controldrawfeatures.js:16   root.ControlDrawFeatures = factory(root.ol);
at ol3-controldrawfeatures.js:18   }(this, function(ol) {
The similar problem was raised here:
ol.inherits in OpenLayers 6
Uncaught TypeError: Cannot set property 'Ending' of undefined
at new ol.control.ControlDrawFeatures (ol3-controldrawfeatures.js:24)
at qgis2web.js:135
which refers to:
at ol3-controldrawfeatures      options.draw.Ending = true;
at qgis2web.js        var buttonsDrawControls = new ol.control.ControlDrawFeatures({vector_draw, options});
Is there something I have forgotten here?


texts: Variable buffer radius in ArcGIS JavaScript API?

texts:
I'm trying to combine the following Esri samples, to end up with the ability to create buffers based on a variable radius, that will then query data within the buffer zone.

https://developers.arcgis.com/javascript/jssamples/util_buffergraphic.html
https://developers.arcgis.com/javascript/jssamples/query_buffer.html

I'm new to JS so I'm not sure if I'm on the right lines or not, but using the query_buffer as a template in my latest attempt, I'd thought that adding inputs for the distance and units to the HTML, then replacing the radius and radiusUnit values in the JS with this, would work.

 //when the map is clicked create a buffer around the click point of the specified distance.
    map.on("click", function(evt){
      circle = new Circle({
        center: evt.mapPoint,
        geodesic: true,
        radius: dom.byId("distance").vlaue,
        radiusUnit: dom.byId("unit").value
      });
      map.graphics.clear();
      map.infoWindow.hide();
      var graphic = new Graphic(circle, circleSymb);
      map.graphics.add(graphic);

Click on the map to select census block points within buffer.

Buffer Parameters:
Distance:&nbsp;<input type="number" id="distance" size="5" value=1 />
<select id="unit" style="width:100px;">
  <option value="esriMiles">Miles</option>
  <option value="esriFeet">Feet</option>
  <option value="esriKilometers">Kilometres</option>
  <option value="esriMeters">Metres</option>
</select>

However, when testing this, it doesn't seem to be picking up the radius values, and won't draw the buffer.

Any ideas?


texts: Why using buffers in land use / land cover classification task?

texts:
I am following the tutorial How To: Land-Use-Land-Cover Prediction for Slovenia and I would like to understand why they used a positive buffer (expanding boundary limits) in the first part of the tutorial; and why they used a negative buffer (when they eroded the parcels: ground truth data) in the second part.
It is not the first time I observe buffers being used like this.


texts: Image classification history in ArcGIS Pro

texts:
I would like to examine the details for a supervised classification project performed previously in ArcGIS Pro. Specifically, I am interested in the segmentation characteristics used (spectral and spatial detail, etc.).
Is there a way to view this in history?
I can't find a record of it in the project raster function history pane.


texts: hasLayer and removeLayer methods not working on leaflet map as expected

texts:
I am trying check if Layer already exist on map...

IF(no)  THEN Addlayer on Map.
IF(yes) THEN removeLayer and AddLayer on map "because I dont want same Layer to be added again and again on map".

I think, I am close but I dont know why its not working... Help Please...

            // load GeoJSON from an external file
            $.getJSON(url,function(data){
                // add GeoJSON layer to the map once the file is loaded
                DistrictLayer = L.geoJson(data, {
                    style: {color: "#000", weight: 2, fillColor: "#03a6fa", fillOpacity: .6},
                    onEachFeature: function (feature, layer) {
                        layer.bindPopup("<strong>" + feature.properties['DISTRICT'])
                    }

                })
                if(map.hasLayer(DistrictLayer)){
                    map.removeLayer(DistrictLayer);
                    map.addLayer(DistrictLayer);
                    map.fitBounds(DistrictLayer.getBounds());
                }
                else{
                    map.addLayer(DistrictLayer);
                    map.fitBounds(DistrictLayer.getBounds());
                }

            });


texts: Removing T and Z in WFS timestamp o/p field by changing configurations in GeoServer

texts:
Geoserver is appending T and Z in the values of timestamp columns(postgres DB) in WFS o/p.
Current Output:2020-12-11T00:00:00Z
Expected Output:DEC-11-2020 00:00:00


texts: What does "Input geometry has unknown (-1) SRID" mean?

texts:
I am trying to run the following query to determine the area and number of points of a polygon:

Select osm_id, st_npoints(the_geom) as numPoints, st_area(transform(the_geom, 3395))/1000000 as areakm2 from water_table where st_isvalid(the_geom);

I receive the following error:

Input geometry has unknown (-1) SRID

Anyone know what the problem is?


texts: Right click on layer in TOCControl get a context menu

texts:
Is there a way to right click on a layer in the AxTOCControl and get a context menu?

I would like the options to be something like: {annotate layer, change annotation, move layer up or down the view stack}

Using ArcEngine 10 C# VS2010


texts: Loading .qml file for raster using Python

texts:
How can I modify this script in order to load .qml for all the raster (not vector in my case) files in the canvas?

layer = iface.activeLayer()

if layer.geometryType() == QGis.Point: # I need to change from vector to raster

    layer.loadNamedStyle('c:\\myQml.qml')
    layer.triggerRepaint()


texts: Showing geology data from Edina Digimap in QGIS

texts: I have downloaded Geology data from Edina Digimap and want to view this on QGIS. How do I do this once I have downloaded the geology data from Digimap?


texts: Translate sld:name and ogc:property name into GeoServer CSS

texts:
Given is following SLD rule:

<sld:Rule>
  <sld:Name>mytrack#</sld:Name>
  <ogc:Filter>
    <ogc:PropertyIsEqualTo>
      <ogc:PropertyName>jogging</ogc:PropertyName>
      <ogc:Literal>always#</ogc:Literal>
    </ogc:PropertyIsEqualTo>
  </ogc:Filter>
  <sld:MinScaleDenominator>100</sld:MinScaleDenominator>
  <sld:MaxScaleDenominator>1000</sld:MaxScaleDenominator>
  <sld:LineSymbolizer>
    <sld:Geometry>
      <ogc:PropertyName>alltracks</ogc:PropertyName>
    </sld:Geometry>
    <sld:Stroke>
      <sld:CssParameter name="stroke">#ff0000</sld:CssParameter>
      <sld:CssParameter name="stroke-width">1</sld:CssParameter>
    </sld:Stroke>
  </sld:LineSymbolizer>
</sld:Rule>

which I would write with css like:

[jogging='always#'][@scale>100][@scale<1000] {
  stroke: #FF0000;
  stroke-width: 1;
}

and what is translated by geoserver into following sld:

    <?xml version="1.0" encoding="UTF-8"?><sld:StyledLayerDescriptor xmlns="http://www.opengis.net/sld" xmlns:sld="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc" xmlns:gml="http://www.opengis.net/gml" version="1.0.0">
  <sld:NamedLayer>
    <sld:Name>Default Styler</sld:Name>
    <sld:UserStyle>
      <sld:Name>Default Styler</sld:Name>
      <sld:FeatureTypeStyle>
        <sld:Rule>
          <ogc:Filter>
            <ogc:PropertyIsEqualTo>
              <ogc:PropertyName>jogging</ogc:PropertyName>
              <ogc:Literal>always#</ogc:Literal>
            </ogc:PropertyIsEqualTo>
          </ogc:Filter>
          <sld:MinScaleDenominator>100.0</sld:MinScaleDenominator>
          <sld:MaxScaleDenominator>1000.0</sld:MaxScaleDenominator>
          <sld:LineSymbolizer>
            <sld:Stroke>
              <sld:CssParameter name="stroke">#0029FF</sld:CssParameter>
              <sld:CssParameter name="stroke-width">1.5</sld:CssParameter>
            </sld:Stroke>
          </sld:LineSymbolizer>
        </sld:Rule>
        <sld:VendorOption name="ruleEvaluation">first</sld:VendorOption>
      </sld:FeatureTypeStyle>
    </sld:UserStyle>
  </sld:NamedLayer>
</sld:StyledLayerDescriptor>

How can rename

<sld:Name>Default Styler</sld:Name>

to

<sld:Name>mytrack#</sld:Name>


texts: Reading multiple grid data with Python and ModelBuilder?

texts:
How can I read grid data (3 or 5 sources of grids) like slope, soil, cover, etc?

My goal is to build a simulator based on multiple grid data for harvesting decisions.


texts: Where to find the "contours from points" model in QGIS 2.6?

texts:
I can't find the example model "contours from points" in QGIS 2.6.
Is there another tool for creating contours from points, besides the QGIS Contour plugin?


texts: How to get a list of attributes in FME when there are no records in table?

texts:
How to get a list of attributes in FME when there are no records in the table?
I try use a PythonCaller but when there are no records, FME passes over this table.
In PythonCaller I write:

import fmeobjects
def FeatureProcessor(feature):
    print feature.getAllAttributeNames()

Is it possible to get attributes when there is no record?


texts: Snapping geometry/polygon to grid with minimum distortion

texts:
I am working on a project to transform polygons given in a pre-defined precision to a lesser precision, which I think is equivalent to snap them to a fixed spacing grid.
For example, in the example below, I have a polygon in yellow, I need to snap its vertices to the anchors of a grid (in green). I have tried to load my data to QGIS and use the snap to grid feature, which generates transformed polygon (in cyan). But you can see the transformed polygon has some large distortion, especially in the lower right corner. For this example, I would think the polygon (in red) in the second image is probably a more similar transformed version to the original polygon.
Is there any algorithm that works on this kind of problem, to snap a geometry to grid anchors with minimum distortion? Vertices of the transformed polygon do not need to be the closest to their correspondence on the original polygon. Adding or removing vertices are also allowed, but the resulting shape should be as similar as possible to the original polygon, and the resulting polygon vertices must be on the grid anchors.
I have not found a good metric to measure the similarity, but I am thinking the angle at vertices with two long sides shall be with minimum change. It sounds even like an optimization problem to me.


texts:
Having two points (longitude, latitude) on earth, if I am sure their distance is under 100 km; how much error will I have using tunnel distance approximation? (Defined at here)

Note: Currently I am using haversine formula, but since my points are normally in a 2 km range (at 99% of time and for other 1% it goes from 2 to 100 km) I think I do not need that complex calculation.


texts: How can one easily differentiate KML placemarks of polygon and icon type?

texts:
I'm using fastkml in python to parse a large kml file where Google Earth has made a new style for every placemark.

I'm trying to parse subfolders which contain a point placemark and/or a polygon placemark, whereby each placemark is assigned one of 5 icon styles or one of 5 polygon styles (respectively), then apply that style accordingly.

Eg.

root/fldr1/num1/pm1
if pm1 is Point:
    pm1.styleUrl = "#point_class2"
elif pm1 is Polygon:
    pm1.styleUrl = "#poly_class2"
else:
    raise Exception ("Nope")

Should/(can) I merge the styleUrls?
Is there a better way I should go about tidying 50000 lines of superfluous Google Earth-created styles?

NB. I asked a similar more specific question here


texts: Why is the 'straight line' path across continent so curved?

texts:
This is the result of mapping the straight line path from a point in US to Poland using Distance Measurement Tool.

Also, planes from Asia to US would travel almost over North Pole.

Why is the path so curved? I agree that this is a flat representation of a sphere, so I do expect some arc, but I don't think earth has this much curvature.

What am I missing here?


texts: Get MapBox Marker Coordinates Dynamically

texts:
I am new to MapBox and this may be a very simple solution, what I am trying to Achieve is to get the Coordinates of the marker which pops up when I search for Any Address/Location in the Search box. The Marker is a draggable, and it should fetch the updated coordinates when moved.
Please have a look at the below code:

<meta charset="utf-8" />
<title>Locate the user</title>
<meta name="viewport" content="initial-scale=1,maximum-scale=1,user-scalable=no" />

<link href="https://api.mapbox.com/mapbox-gl-js/v1.11.1/mapbox-gl.css" rel="stylesheet" />

    <link
    rel="stylesheet"
    href="https://api.mapbox.com/mapbox-gl-js/plugins/mapbox-gl-geocoder/v4.5.1/mapbox-gl-geocoder.css"
    type="text/css"
    />

Please add you Mapboxgl AccessToken.


texts: Error in projection or just a graphical error?

texts:
I've got a worldwide SRTM grid in WGS84 and want to reproject it to UTM coordinates.

If I do so in parts of the world are strange tiles of the map (especially in Eurasia), which are not part of the UTM stripe. You see a example in the map below for the UTM 5N in a Mercator projected map.
What are reasons for that?
As resamping technique I used bilinear.


texts: How to merge attributes fields and their values into a single table using a common key with FME

texts:
I am having difficulties with a relatively simple task. I have multiple transformers  which append an attribute to an "ElementID" coming from different transformers.
I would like to merge all outputs into a single table based on the ElementID field thus having a schema like "ElementID","Attribute01, Attribute02, Attribute03  instead of having 3 tables with ElementID, Attribute01 in table 1, and ElementID, Attribute02 in table 2 etc.

Not all ElementIDs have all three attributes..ie there will be missing attributes.

How do I achieve this with FME 2014.


texts: Problems with case statements in Leaflet's pointToLayer

texts:
I am trying to create a layer that has different types of business (separated by code properties in the geojson file) but when it executes the pointToLayer below I get an Uncaught TypeError: Cannot set property 'feature' of undefined. What am i doing wrong?

  var allbusinesses = L.geoJson(data, {
    style: allBusinessesStyle,
    onEachFeature: onEachFeature,
    pointToLayer: function(feature, latlng) {
      switch (feature.properties.code) {
        case "1":
          return L.circleMarker(latlng, {
            color: "#000"
          });
        ...
      }
    }
  });


texts: How to find X and Y of clustered features?

texts:
How can I find longitude latitude values of clustered features?

Please look the image. I hope you understand my question. I need to find out the X and Y of clustered features on every zoom levels.

Sample Code:

var defaultStyle = new OpenLayers.Style({
        strokeColor     : "#FFFFFF",
        strokeOpacity   : 1,
        strokeWidth     : 2,
        fillColor       : settings.color,
        fillOpacity     : 0.9,
        cursor          : "pointer",
        pointRadius     : "${point_radius}",
        pointerEvents   : "visiblePainted",
        label           : "${label}",
        fontColor       : "#000000",
        fontSize        : "12px",
        fontFamily      : "verdana",
        fontWeight      : "bold",
        labelAlign      : "center"
    }, {
        context: {
            label: function(feature) {
                var sumValue = 0;
                if(feature.attributes.count>1) {
                    for (var i=0; i<feature.cluster.length; i++) {
                        sumValue++;
                    }
                    return sumValue;
                }else{
                    return 1;
                }
            },
            point_radius: function(feature) {
                var pix = 10;
                if(feature.cluster) {
                    pix =Math.min(feature.attributes.count, 30) * 2.5 + 8;
                }
                return pix;
            }
        }
    });

    var onSelectStyle = new OpenLayers.Style({
        strokeColor     : "#000000",
        fontColor       : "#FFFFFF",
        pointerEvents   : "visiblePainted",
        fontSize        : "15px"
    });

    var strategy = new OpenLayers.Strategy.Cluster({distance: 30, threshold: 2});

    var ly=new OpenLayers.Layer.Vector(GSMapLayerName,{
            strategies: [strategy],
            styleMap:new OpenLayers.StyleMap({
                        "default":defaultStyle,
                        "select": onSelectStyle
                    }),
            displayInLayerSwitcher:false
        });
    map.addLayer(ly);

Note: I am using openlayers 2.11

Thanks for any reply


texts: ArcGIS Javascript API - Label renderers from Map Service for a FeatureLayer

texts:
I wonder if it is possible to use the labels renderers defined in a .mxd published as a Map Service (Feature Access enabled) for a FeatureLayer with the javascript API v4.3 just like it works with the geometry renderer?

I define a featureLayer like this:

var organismesLayer = new FeatureLayer({
  url: "https://infogeo.education.gouv.qc.ca/arcgis/rest/services/SandBox_Maxime/GDUNO/MapServer/1",
  labelsVisible: true
});

In the map, I see the layer style is the same than in the .mxd, but the labels are not visible. Do I have to define the LabelClass manually in the javascript code?

Is there a way to get the label renderers from the Map Service?


texts: GeoServer WMS layers not working on Azure app service

texts:
Tl;dr I'm hosting GeoServer (2.15.0) in an Azure app service Java container with Tomcat. I can access the admin panel, but none of the layer previews work, nor can I fetch them over WMS service. They just timeout after a while.

In July a couple of colleagues used GeoServer and leaflet.js to show some map data in a web app. Fast forward to today, they've both quit the company and I took over the project. Nothing was changed on the app service or the GeoServer since then (we have 2 separate instances of GeoServer) but they've both stopped working.
I've tried reinstalling GeoServer on one of the instances, changed different settings in Azure portal (Java & Tomcat versions), created a new app service and installed a fresh copy of GeoServer to no avail.
I've contacted my former colleagues and they say the last time it was a simple thing of unpacking the zip and it all worked (another colleague confirms it).

Before I try other options (using a virtual machine), has anyone encountered similar problems? Is there some trick to hosting GeoServer in an Azure app service?


texts: Automatically export numbered features to txt

texts: How to automatically count the numbers of the buildings in a layer (CM_Buldings layer in the pic) and output the number as a txt file in QGIS?


texts: Which topics are fundamental to ArcGIS for Desktop?

texts:
Similar to Which topics are essential for a crash course in QGIS?, which topics are fundamental to ArcGIS for Desktop?

Basic through advanced topics are welcome. Should be freely accessible to all.


texts: Creating donut holes without losing existing polygons using ArcGIS for Desktop?

texts:
Using ArcGIS 10.1 for Desktop is there a way to do a similar operation to creating a donut hole or island in a polygon without losing the existing features?

I have mapped small lakes as several different polygons based on vegetation type.

However, I want to map the vegetation surrounding the lakes, too - but if I create a circular polygon over the top of the lake and the surrounding area, it includes the area I've already mapped and given attributes to.

So I want to 'cut out' a hole in the middle of the new polygon without losing the polygons that are now covered by it.