texts: Why does Uncaught TypeError: Cannot read property 'trim' of undefined occur in esriRequest

texts:
I am trying to use esri/request to pull in some data from an API. I am have two weird problems. First, when I try to run this, it doesn't complete and throws the following error in the browser console: "Uncaught TypeError: Cannot read property 'trim' of undefined". I am not using trim() anywhere and am not sure why this error is being thrown (specifically on the line where I make my esriRequest). Second, if I try to test this in the browser by pasting in the function to the console, it says that esriRequest is not defined. I have checked the order of my required list and functions several times and they look like they are in the proper order. Am I missing something here?

require([
    "esri/map",
    "esri/basemaps",
    "esri/dijit/BasemapToggle",
    "esri/layers/ArcGISDynamicMapServiceLayer",
    "esri/tasks/ImageServiceIdentifyTask",
    "esri/tasks/ImageServiceIdentifyParameters",
    "esri/layers/WebTiledLayer",
    "esri/dijit/Search",
    "esri/config",
    "esri/request",
    "dojo/json",
    "dojo/promise/all",
    "dojo/domReady!"], function(
        Map,
        esriBasemaps,
        BasemapToggle,
        ArcGISDynamicMapServiceLayer,
        ImageServiceIdentifyTask,
        ImageServiceIdentifyParameters,
        WebTiledLayer,
        Search,
        esriConfig,
        esriRequest,
        JSON,
        all) {

        .....
       var url = "https://earthquake.usgs.gov/ws/designmaps/asce7-16.json?latitude=18&longitude=-66&riskCategory=I&siteClass=D&title=Default";
      esriRequest(url, {
        responseType: "json"
      }).then(function(response){
        // The requested data
        var seis_data = response.data;
      });


texts: How to plot distances distances in meters in QGIS?

texts:
My goal is to plot circles on QGIS canvas by taking distances in kilometers from MySQL base. This is code I use:

from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.core import *
from qgis.gui import *
from math import *

def funkcijaCrtanje(iface):

import MySQLdb
import sys

connection = MySQLdb.connect(host = "127.0.0.1", user = "root", passwd = "********", db = "mydb")
cursor = connection.cursor()
cursor.execute("select column from table")
data = cursor.fetchall()

edge = data[0]

def drawCircle(rb, center, edgePoint, N):
    r = sqrt(center.sqrDist(edgePoint))
    for itheta in range(N+1):
        theta = itheta*(2.0 * (pi)/N)
        # addPoint finkcija dodaje tocke QgsGeometry objektu i povezuje ih linijom
        rb.addPoint(QgsPoint(center.x()+r*(cos(theta)),center.y()+r*(sin(theta))))

class Kruznice(QgsMapTool):
    def __init__(self, canvas):
        QgsMapTool.__init__(self, canvas)
        self.canvas = canvas

        self.segment = 100

    def canvasPressEvent(self, e):
        point = self.toMapCoordinates(e.pos())
        n = 10
        color = [Qt.red, Qt.blue, Qt.green, Qt.yellow]
        for i, edge in enumerate(data):
            rb = QgsRubberBand(self.canvas)
            rb.setColor(color[i])
            rb.setWidth(2)
            ed = QgsPoint(point.x()+float(edge[0])*n*100, point.y()+1000)
            drawCircle(rb, point, ed, self.segment)

krug = Kruznice(iface.mapCanvas())
iface.mapCanvas().setMapTool(krug)

I'm not sure that I plot real distances on my map by using this. How can I solve this?


texts:
I'm attempting to use the split tool to divide a vector shapefile into two parts. I used a process that I've been successful with before:

Go to Edit > Split > Select my feature to split > Draw a line that goes all the way through the feature (past two edges) > then either click "finish" or double click to close my line

For some reason when I went through this process just now it split the feature along my line but then immediately erased the smaller polygon. Has anyone else had this problem and/or know how to fix it?


texts: QGIS doesn't read raster layer CRS properly - ESRI:54004 instead of EPSG:3857

texts: I've switched versions from 3.4 to 3.12 and now QGIS reads rasters created in older QGIS wrong. When I'm clicking "Set CRS - Set Project CRS from Layer" to check it out it shows ESRI:54004 (don't know what this is). When I set raster CRS manually to EPSG:3857 it's all fine. My question is why QGIS dosen't read EPSG:3857 from raster? Is it my fault, or is it bug? Should I re-save all my rasters?


texts: Creating a distance overlay that is computed by the distance of points/lines and recalculated on zooming

texts:
I created a fictional map and visualized it with d3. Now I would like to show a scale/rule so that the user can see which distance (kilometers,meters) he is currently using.

I was thinking about creating lines in QGIS (with which I created the map) and than setting the distance in km/m as attributes. After that using these information for the scale on the website.

Is that a good way or is there any easier and better way to do this.

The scale should also change the values when zooming for sure. So a picture is no solution.

Found this here: which is generating distances. But I'm not sure that a can adapt this for me: http://bl.ocks.org/mpmckenna8/def3735b8e1eef5128a7

This one is nearly what I'm looking for: http://www.qgistutorials.com/en/docs/making_a_map.html But I want to create that scalebar on my site and not just as an image. So how to transfer that information?


texts: Python error on previously working model

texts: I have a model built that's adding columns to a vector layer and creating a fixed width buffer layer.  Everything was working fine, but now when I go to edit the model it won't open and I'm receiving the following python error...


texts: Can declutter work with MultiPoint features?

texts:
I have a vector layer with Point data. There could easily be between 10K and 30K features to render, each feature with their respective label, so decluttering is needed.

My goal is to have the feature/"circle" of each point rendered (to visually see the data distribution), with the label/text being decluttered when needed. Like in this little screenshot where we see all features but only some of the labels:

A key point here is that normally to render Point vector data, one would make use of the ol.style.Circle class;

image = new ol.style.Circle({
radius: xxx,
fill: new ol.style.Fill({
    ...
}),
stroke: new ol.style.Stroke({
    ...
}),

And for labels we have to use the "text" element of the Style:

text: new ol.style.Text({
...
})

But the declutter mechanism is resolving conflicts between "image" and "text" elements. As from the API:

  Declutter images and text. Decluttering is applied to all image and text styles, and the priority is defined by the z-index of the style. Lower z-index means higher priority.

So this gives us a challenge rendering points using both the image and the text elements in the applied styles.

For geometry type 'Point' I have managed to find a way around this by using the "geometry" element of the Style, and I construct a new ol.geom.Circle allowing me to render points not defining the "image" element, but instead applying standard "fill" and "stroke" like we do for Polygons or Lines.

geometry: new ol.geom.Circle(feature.getGeometry().getCoordinates(), 2*resolution)

So the Style then becomes:

var featureStyle = new ol.style.Style({
geometry: new ol.geom.Circle(feature.getGeometry().getCoordinates(), 2*resolution),
fill: new ol.style.Fill({
  ...
}),
stroke: new ol.style.Stroke({
  ...
})

});

However, for 'MultiPoint' features, I see no way of working around the use of the "image" element.

Hence, the core question here is:
Is it not possible to render a MultiPoint with labels as long as the layer is in declutter mode? Have I missed some part of documentation stating this is not possible? Or do I still not understand all details of the decluttering mechanism?

If I have a MultiPoint feature using a Style with both image and text, and turn on declutter, it is not rendered at all on the map. See this Codepen as example: https://codepen.io/eivindronnevik/pen/GRKObvE

Currently the layer has declutter:true, set declutter to false and the MultiPoint feature will appear.

If I remove the text-element of the Style, then the features render. If I remove the image-element, the labels render. I would have understood if there were individual Points of the MultiPoint disappearing when decluttering comes into play, e.g. the "image" element of these two Points were too close - so only one of them are rendered. But in theory we can create a MultiPoint with 2 Points, one in Russia and one in the US, and still if declutter is on nothing will render.

Using Points, an image and text element of the same Style object does not "collide", declutter will either render both or none (in relation to other Points around). But it seems like for MultiPoint this is different, it is like the image and text element of each individual Point collide, hence nothing is rendered.

Any input or possible corrections of my understanding of decluttering is sought.


texts: What can I do in Mapbox Studio with the osm_id number for #buildings from Mapbox Streets?

texts:
My goal is to use OpenStreetMap buildings data to create map tiles that show buildings styled with information such as name, floors, and other information from the tags associated with the way/feature in the OSM database.

Is it possible to access the name of a #building through the osm_id in Mapbox Streets or by some other method?

What is the potential use of the osm_id in Mapbox Studio?

The specific use case here is a college campus map.


texts:
Following example through ESRI is to access the features of file geodatabase (GDB), I would like to know that is there a similar way to do this with personal geodatabase (MDB) using Java.

I am quite new to ArcObjects so if anybody could point me to right direction.


texts: Extracting lat/long from screen when using British National Grid projection in ArcGIS Desktop?

texts: Using ArcGIS Desktop 10.0, is it possible to click anywhere on the map (projected to British National Grid), but return the coordinates in WGS84 lat long?


texts: Selection feature acting weird on queried layer in QGIS

texts:
My select feature is acting weird. I am making a plugin that downloads an WFS layer and then filter out objects based on the features. But when i try to select the points by an area, the points that are selected seems very random. It appears that no matter where i make the selection, the the same points are selected. If I select two points, it's always the same two, if i select ten points, it's always the same ten.

I get the layer by this method:

        QgsMapLayerRegistry.instance().addMapLayer(WFSlayer)
        self.hideLayer(WFSlayer)
        self.iface.legendInterface().setLayerVisible(WFSlayer, False)
        layer_name = self.WFSlayer.name()
        query = "SELECT * FROM " + layer_name + " " + where #where statement created elswhere
        vLayer = QgsVectorLayer("?query=%s" % (query), "newLayer", "virtual" )
        QgsMapLayerRegistry.instance().addMapLayer(vLayer)
        self.canvas.setExtent(vLayer)

I have another method of adding layer in memory insetad, and this layer does not have the same problem.

I have no idea why this particular layer acts like this, while all others acts normally.

Area that is being selected:

The obects that gets selected


texts: Python Add-In Dialog with Textboxes in ArcPy?

texts:
With Python Add-In Wizard I created a ButtonTool with which I:

click in the map
get the coordinates of the mouseclick-location in a dialog with textboxes

I want  to receive the coordinates in a textbox to have the possibility to copy and paste the coordinate values (e.g. to use them  in a document).
I get a crash between TKinter and arcpy like I read in some posts before. Nevertheless I want to ask, if anybody has experience with this and has a suggestion.

Here is my code:

Import Tkinter
class ToolCoord(object):
    def __init__(self):
        self.enabled = True
        self.shape = "NONE"
        self.cursor = 3
    def onMouseDownMap(self, x, y, button, shift):
           pt = arcpy.PointGeometry(arcpy.Point(x,y))
           main = Tkinter.Tk()
           lb = Tkinter.Label(main, text = "UTM-Koordinaten:")
           lb.pack()
           t1 = Tkinter.Text(main, width=20, height=0)
           t1.insert("end", str(x))
           t1.pack()
           t2 = Tkinter.Text(main, width=20, height=0)
           t2.insert("end", str(y))
           t2.pack()
           main.mainloop()


texts: Find lines whose both start and end point ends in a polygon

texts:
I have a river dataset, which is dissolved, so between two lakes there will always be one single river line.

Now I want to find all rivers that starts in a lake, and ends in a lake. So if my rivers are stored as line features, and the lakes as polygon features, I  need to find all lines whose both start point and end point touches a polygon in the lake data set.

How do I do that in either ArcGIS or QGIS (preferrably without any programming).


texts: Adding field to table without path to table (python)

texts:
I want to add a field to my table in the FieldCalculator.

import arcpy
arcpy.AddField_management("veg8", "test", "TEXT")

An error message occured doing this and I figured out that I need to set a path to my table
Like this for instance

C:\\user\\...

But, since I work in the Modelbuilder and the Model will be running on other computers, the path to the table will be a different one and I dont want to update all my paths everytime I run the model on a different computer. Is there a way to tell python to use the table to which the FieldCalculator is "connected"?

Thanks
Chris


texts: Spatialite: 'no such function: GeomFromText' error

texts:
I've this PHP code on my Ubuntu 15.10 virtual machine

<?php
  $lon = 4.221542;
  $lat = 32.21115;

  $db_data_sessions = new SQLite3('./DataBase/MyDB.sqlite');

  # Loading SpatiaLite as an extension ...
  $db_data_sessions->loadExtension('mod_spatialite.so');

  $q="SELECT comune FROM comuni32632 WHERE
      WITHIN(Transform(GeomFromText('POINT(".$lon."
      ".$lat.")',4326),32632),comuni32632.Geometry)";

      echo "Query = " .$q;
      echo "\n\n";

      try {
        $stmt = $db_data_sessions->prepare($q);
        $results = $stmt->execute();

        if ($results->fetchColumn() > 0)
         {
           echo "OK";
           echo "\n\n";

           $theName = '';
           while($row=$results->fetchArray(SQLITE3_ASSOC))
             {
              $theName = $row['comune'];
             }
           //return $theName;
           echo "Result = " .$theName;
           echo "\n\n";
         }
         else {
            //return "NO DATA";
            echo "Result = NO DATA";
            echo "\n\n";
         }
     }
     catch(PDOException $e) {
        print "Something went wrong or Connection to database failed! ".$e->getMessage();
        //return "ERROR";
        echo "Result = ERROR";
        echo "\n\n";
   }
   $db_data_sessions = null;
?>

and when I try to execute I can see only the first echo output, so my current query text and not more. In the browser console I can see "Error 500"

When I try to execute from command line

php test.php

the result is

PHP Warning:  SQLite3::loadExtension(): SQLite Extension are disabled in /var/www/html/test.php on line 10
Query = SELECT comune FROM comuni32632 WHERE WITHIN(Transform(GeomFromText('POINT(4.221542 32.21115)',4326),32632),comuni32632.Geometry)

PHP Warning:  SQLite3::prepare(): Unable to prepare statement: 1, no such function: GeomFromText in /var/www/html/test.php on line 19
PHP Fatal error:  Call to a member function execute() on boolean in /var/www/html/test.php on line 20

Note that if I try to execute the query in Sqlite3 loading the Spatialite extension all works fine

Note that if I try to execute the PHP test code available here http://www.gaia-gis.it/gaia-sins/spatialite-cookbook/html/php.html all works fine so I think that PHP and Spatialite are working fine together

What is going wrong in my code?


texts: Flood Plain Mapping using HEC-GeoRAS and ArcMap?

texts:
How can I get complete bounding polygon of XS CutLines for flood plain mapping (Flood Inundation Mapping) using HEC-GeoRAS and ArcMap 10.3?

When I tried to perform this task, the bounding polygon formed was in an unexpected manner.

I want to get the bounding polygon without the gap that is visible around the river junction.


texts: Calculate if day/night time for point dataset

texts:
I have a point dataset (crime data for Stockholm) with timestamps and I want to calculate for each point if the certain crime happened during day or night time. So connecting the timestamp to sunrise and sunset. I know that there is a python tool called pyephem, but I don't know how I can use this tool to automatically calculate day/night time for my entire dataset.

So I would need a python code, which imports the timestamp column for each feature. In the next step it has to be calculated if it is day/night for every point and finally a new column should be added in shapefile with the calculated result (e.g. 'day' / 'night').


texts: Resizing selection of polygons to percentage of their initial size using ArcGIS Desktop?

texts:
How do I resize a selection of polygons to a percentage of their initial size (3000 polygons)?

I don't think the Arc "Scaling" tool will be useful as it scales the polygons around a point.  I would like to decrease their size all around their individual centroids.


texts: Why is my stacked area graph in ggplot2 returned as stacked lines?

texts:
I'm having a bit of a problem creating my stacked area chart in ggplot2 in R, for some reason it is returning as a stacked line;

e.g. data:

library(ggplot2)

Sector <- rep(c("S01","S02","S03","S04","S05","S06","S07"),times=7)
Year <- rep(c("1950","1960","1970","1980","1990","2000","2010"),each=7)
Value <- runif(49, 10, 100)
df <- data.frame(Sector,Year,Value)

and then using ggplot2 to create a traditional stacked area chart;

ggplot(df,aes(Year, Value,fill=Sector))+geom_area(aes(colour=Sector),position="stack")

for me, that returns a stacked line as per below. Is it my version of R? or my code? the code for ggplot2 looks like every other eg of a stacked area I can find online but I get odd results and I can't find any info.


texts: OpenLayers JavaScript bug when dragging outside the map element

texts:
I'm trying to learn OpenLayers. But the project I try to create has a bug when moving the mouse outside of the map HTML element when dragging the map. When you try to navigate by keeping the left mouse button down and releasing it outside of the map HTML element.

Here is a basic project from the OpenLayers website that has this exact problem: The Openlayers project

Here you can view the code of the project:

import "ol/ol.css";
import Map from "ol/Map";
import View from "ol/View";
import GeoJSON from "ol/format/GeoJSON";
import { Tile as TileLayer, Vector as VectorLayer } from "ol/layer";
import { bbox as bboxStrategy } from "ol/loadingstrategy";
import OSM from "ol/source/OSM";
import VectorSource from "ol/source/Vector";
import { Stroke, Style } from "ol/style";

var vectorSource = new VectorSource({
  format: new GeoJSON(),
  url: function(extent) {
    return (
      "https://ahocevar.com/geoserver/wfs?service=WFS&" +
      "version=1.1.0&request=GetFeature&typename=osm:water_areas&" +
      "outputFormat=application/json&srsname=EPSG:3857&" +
      "bbox=" +
      extent.join(",") +
      ",EPSG:3857"
    );
  },
  strategy: bboxStrategy
});

var vector = new VectorLayer({
  source: vectorSource,
  style: new Style({
    stroke: new Stroke({
      color: "rgba(0, 0, 255, 1.0)",
      width: 2
    })
  })
});

var raster = new TileLayer({
  source: new OSM()
});

var map = new Map({
  layers: [raster, vector],
  target: document.getElementById("map"),
  view: new View({
    center: [-8908887.277395891, 5381918.072437216],
    maxZoom: 19,
    zoom: 12
  })
});


texts: Will I lose project functionality on upgrading QGIS?

texts: This might be a stupid question but if I update to QGIS 3.2.1 will I have problems working on projects created in my current 2.18.20?


texts: Computing NDVI - Repeated calculation error message (Google Earth Engine)

texts:
I am doing  simple calculation for NDVI using Landsat 7 TOA Reflectance (Orthorectified). Landsat seven has the following bands: B1-Blue, B2-Green, B3-red and B4-nir. I used the following code:

//Compute NDVI

var red = filtered.select('B3');

var nir = filtered.select('B4');

var ndvi = nir.subtract(red).divide(nir.add(red));

Map.addLayer(ndvi, {min: 0, max: 1}, 'NDVI')

*"filtered" being my image that I am trying to calculate NDVI for.

When I run the code I get an error message in the Console that says "nir.subtract is not a function"....but it is. Has anyone experienced something similar, or is there a workaround for this problem?


texts: Encoding Greek alphabet with ArcGIS for Desktop?

texts:
I have a big type of Greek shapefiles, but a lot of time when I show the attribute table I see ?.????????? .

How can I encode them?

Is there any way to make ArcGIS 10.1 for Desktop encode them automatically?


texts: feature in the formOpen method of a custom form display two feature object

texts:
Qgis version : 2.18.1
I was having trouble when launching a custom edit form because of my use of feature in the formOpen method. I want to get an attribute from the current feature as this:

feature.attributes()[1]

I thought it was intended to reference the current feature that was edited and it kinda does it but when I print feature and it's id() it displays every time two feature objects one with id=0 and the second with the correct id. it's not really a big issue because it still works but I get a message each time I launch the form that the list index is out of range.
Does anyone know how to fix or at least circumvent this ?

Edit 1: I should precise that i have set a relation widget on this table whith referencing to other features from the same table.

Edit 2:
Alright so it seems the problem was deeper and my issue could possibly help others. I wanted to post another problem that arose from this custom form whith python logic but it seems that the two are linked.
I've set a widget to plot a graph when a button was clicked, it worked, but stopped after i had set up a relation manager as a widget too, with a relation set on the same layer. That's when the ghost non valid feature problem happened. The thing is that the widget used as the Qgis Relation Manager displays the same dialog so when i asked inside the formOpen(dialog, layer,feature) method to print the current feature id it printed the non valid feature from the Relation Manager (when there wasn't relations) or the children features and then the current feature of the original edit form.
The best solution seems to get rid of the Relation Manager when your relations are within the same layer.
I think i will replace the Relation Manager by customizing a QListWidget and a QLineEdit except if someone comes with a better idea.


texts: ImageCollection (Error) User memory limit exceeded. Mean annual land surface temperature landsat data

texts:
I am trying to calculate mean annual and monthly land surface temperature for three image collection of landsat8, Landsat5, and landsat7(image collection of landsat7 and landsat 5 are merged).
The problem is that there is an error:ImageCollection (Error)
User memory limit exceeded, when I am trying to join two image collection with different bands.

var filter = ee.Filter.equals({
  leftField: 'system:index',
  rightField: 'system:index'
});

// Create the join.
var simpleJoin = ee.Join.simple();

// Applt join
var thermalband8 = ee.ImageCollection(simpleJoin.apply(landsat8thermal, landsat8EM, filter));
var thermalband82 = ee.ImageCollection(simpleJoin.apply(landsat8EM, landsat8thermal, filter));

print('Joined', thermalband8, thermalband82);

var final_col8 = thermalband8.map(function(img){

  // Create a collection with 1 image
  var temp = ee.ImageCollection(ee.List([img]));

  // Apply join to collection 2
  // Resulting collection will have 1 image with exact same date as img
  var join = simpleJoin.apply(thermalband82, temp, filter);

  // Get resulting image
  var i2 = ee.Image(join.first())

  return img.addBands(i2)
});

print('band Names:',final_col8);

var LSt8function = function(image){
  var LST8EQ=image.expression(
    '(Tb/(1 + (0.001145* (Tb / 1.438))*log(Ep)))-273.15', {
      'Tb': image.select('Thermal8'),
      'Ep': image.select('EM')}).rename('LST8');
 return  image.addBands(LST8EQ)};

var LST8= final_col8.map(LSt8function);
print('band Names:',LST8);

Map.addLayer(LST8.select('LST8'),{min: -30, max: 32, palette: ['white','blue','green','yellow' ,'red']},'LST8');

I tried to work on it and the problem is not solved, I need to calculate the LST of each image to use it for estimating the average land surface temperature(annually, and seasonally)

code is here


texts: Add layer of a shapefile into ArcGIS map

texts:
I need to develop an Android application that shows a given shapefile within a map. The problem is that shapefile is only visible if I dont use any mapoption at all (in xml file I delete the mapoption field). Shapefile consists some points in the map and are not visible if mapoption is set to any option streets, hybrid.  https://developers.arcgis.com/android/10-2/api-reference/reference/com/esri/android/map/MapOptions.MapType.html

 try {
       shapefileFeatureTable = new ShapefileFeatureTable(Environment.getExternalStorageDirectory().getPath()+"/test.shp");

    } catch (FileNotFoundException ex) {
        Toast.makeText(this, "File: not found!", Toast.LENGTH_LONG).show();
        ex.printStackTrace();
        return;
    }
    FeatureLayer featureLayer = new FeatureLayer(shapefileFeatureTable);
    featureLayer.setRenderer(new SimpleRenderer(new SimpleFillSymbol(
            getResources().getColor(android.R.color.holo_green_dark),
            SimpleFillSymbol.STYLE.SOLID)));
    if (featureLayer.getSpatialReference()==mMapView.getSpatialReference())
        test.setText("Correct Spatial Reference");
    else
        test.setText("Incorrect Spatial Reference");
    mMapView.addLayer(featureLayer);


texts: Draggable CircleMarker error using GeoJSON

texts:
Draggable circlemarker works if the marker is defined. If I try to read in some GeoJSON data from a file I get the error:"Uncaught TypeError: mymarker.setLatLng is not a function" even if my GeoJSON is one point.

var mymarker = L.circleMarker([42.5, -73.5], {
    color: '#f03',
    fillColor: '#f03',
    fillOpacity: 0.80
}).addTo(map);

mymarker.on({
      mousedown: function () {
      map.dragging.disable();
        map.on('mousemove', function (e) {
          mymarker.setLatLng(e.latlng);
        });
      }
   });
   map.on('mouseup',function(e){
    map.dragging.enable();
    map.removeEventListener('mousemove');
   })

This code below fails: ( I pass it GeoJSON it draws but I can't move anything.)

var mymarker =L.geoJSON(data, {

onEachFeature: onEachFeature,

pointToLayer: function (feature, latlng) {
    return L.circleMarker(latlng, {
        radius: 8,
        fillColor: "#ff7800",
        color: "#000",
        weight: 1,
        opacity: 1,
        fillOpacity: 0.8
    });
}

}).addTo(map);

So as a defined point it works but as a point from GeoJSON it fails. My goal is to be able to move points to their proper location and save their location by updating their lat/Long. I can do this with draggable markers but they want circleMarkers. Any ideas?

https://jsfiddle.net/bchappell61/41nsc3xr/    defined inline (Works)
https://jsfiddle.net/bchappell61/ec34pj6m/  defined from GeoJSON(Fail).


texts: Is there any way to show progress bar while waiting that all layers are rendered on Leaflet map? In this specific case, while GeoJSON layers are rendered.


texts: How to check if vertex exists (on polyline) based on location of point layer using ArcPy?

texts:
Using ArcGIS for Desktop Advanced 10.1 (with all extensions)

I am building a centerline network that will be used by a MapObjects program to do routing based on "turns at any vertex" (as opposed to turns only at endpoints...)  What I am trying to do is identify overpasses/underpasses that have a shared vertex with the road that crosses it (but is not considered an intersection)  If there is a shared vertex, then the mapping program will (incorrectly) make a turn where there is not one physically possible.  What I've done so far is build a topology using the "Must not intersect (Line)" rule.  I export the errors to a feature class, and now I have a bunch of points at intersections that I need to check.  What I'm wondering is if there is a way to cycle through each point and check to see if a vertex exists where the centerline intersects the point.  Here's what I've come up with so far...

import arcpy

fcPoints = "C:/MyPointLayer..."
fcCenterlines = "C:/MyCenterlineLayer..."

rows = arcpy.UpdateCursor(fcPoints)
for row in rows:
    # For each point, select the centerlines that touch (intersect) that point
    arcpy.SelectLayerByAttribute_management(fcPoints, "NEW_SELECTION", "\"OBJECTID\" = " + str(row.getValue("OBJECTID")))
    arcpy.SelectLayerByLocation_management(fcCenterlines, "INTERSECT", fcPoints)
    # Need code here to cycle through each selected centerline and determine if a vertex exists at the point's SHAPE@XY location?
    #PSEUDO CODE:
    for each line in selection:
        vertexExists = False
        for each vertex in line:
            if vertex.xy == point.xy:
                vertexExists = True

        # If there is a matching vertice, then update a field in the point fc to "FLAG" it...
        if vertexExists == True:
            row.setValue("FLAG", "CHECK THIS INTERSECTION")
            cursor.updateRow(row)

Or.... is there an easier way...?  I would be looking at about 2000 points to check...


texts: Why ArcMap 10.2 crashes when running model?

texts: I`m running a simple model in ArcMap 10.2, which makes my ArcMap crashing (closing and then an error message pops in saying that ArcMap has encountered a serious error, and there is an option to report it to ESRI). Basically I want to rasterise some polygons, so I built a model: Iterate Feature Class> Project> Polygon To Raster. Does anyone have an idea why my ArcMap crashes when I run this particular model? (I do run a lot of different model, but have never encountered this issue before...)


texts: What do you think are the top current trends to watch in GIS Development?  Which capabilities/tools do you think will be most in demand to both developers and end users in the future?


texts: Does Shapely support serialization of GeoJSON with empty and non-empty coordinate sequences?

texts:
The shapely.geometry.shape method will fail with an IndexError if a Multi[Point|LineString|Polygon] has an empty coordinate sequence, and I am curious if this is the intended functionality, note the Shapely version is 1.7.1.  An MRE:
from shapely.geometry import multilinestring, shape

multipoint = {
    "type": "MultiPoint",
    "coordinates": [
        [0, 0],
        [],
    ],
}
shape(multipoint) # fails

multilinestring = {
    "type": "MultiLineString",
    "coordinates": [
        [[0, 0], [1, 0], [1, 1], [0, 1]],
        [],
    ],
}
shape(multilinestring) # fails

multipolygon = {
    "type": "MultiPolygon",
    "coordinates": [
        [[[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]]],
        [],
    ],
}
shape(multipolygon) # fails

This seems to only occur if some are empty, if all are empty, the method works fine:
multipoint = {
    "type": "MultiPoint",
    "coordinates": [
        [],
        [],
    ],
}
shape(multipoint) # succeeds

Apologies in advance if this is the wrong place to ask this question, but there seems to be no Shapely users group and the official documentation directs here.


texts: Loading Esri ArcGrid export file (e00) in R?

texts:
I am trying to load an ESRI ArcGrid export file into R. The file is located at- ftp://ftp.epa.gov/castnet/tdep/grids/n_tw/n_tw-2013.zip. The documentation for this file states, "Gridded data of the above variables are available in compressed ESRI ArcGRID export files" and not much else. Link to documentation  here. My major problem is that this is an e00 file, an old output format of Arc.

The most straightforward way to load a raster normally (if this wasn't an e00 file) would be:

require(raster)
require(rgdal)
test <- raster('/path/to/n_tw-2013.e00')

However, you will receive the error: n_tw-2013.e00 not recognised as a supported file format..

The RArcInfo package claims to be able to convert .e00 files into more useful Arc/Info binary coverage. To do this:

require(RArcInfo)
#first argument is the path to the e00 file, and the second argument is the new directory to create
e00toavc('/path/to/n_tw-2013.e00','/path/to/test1')

When I run this command it prints NULL, and creates two directories, test1, and info, however, both are empty. I'm not sure whats going wrong here.

How do I get this loaded into R as a raster so I can extract data to my specific lat/long?

Output of sessionInfo pasted here:

> sessionInfo()
R version 3.2.2 (2015-08-14)
Platform: x86_64-redhat-linux-gnu (64-bit)
Running under: CentOS release 6.7 (Final)

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
 [9] LC_ADDRESS=C               LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods
[7] base

other attached packages:
[1] rgdal_1.2-5  raster_2.5-8 sp_1.2-4

loaded via a namespace (and not attached):
[1] tools_3.2.2     Rcpp_0.12.3     grid_3.2.2      lattice_0.20-33


texts: Finding local maxima with a variable window size in Python

texts:
I am trying to find local maxima in a DEM (canopy height model) that represents tree tops. I can do this with a fixed window size. However I would rather have a variable window size based on the height of each cell with a function such as windowSize=x*0.1+1. I found this approach in a in R package ForestTools, based got the method from a paper by Popescu and Wynne from 2004.

At the moment I am trying to do the same in Python using the function peak_local_max from the Scikit package but this does not seem to have any options for variable windows. I have looked for other python packages/functions but so far without success. Ideally I would like to only use Python. The Forest tools has some issues with large files, I could also move to R if it has a better alternative.

Code so far:

chm_file = "my_chm.tif"
chm_dataset = gdal.Open(my_chm)
chm_raster = chm_dataset.GetRasterBand(1)
cols_chm = chm_dataset.RasterXSize
rows_chm = chm_dataset.RasterYSize
chm_array = chm_raster.ReadAsArray(0,0,cols_chm,rows_chm).astype(np.float)
local_maxi = peak_local_max(chm_array,indices=False, footprint=np.ones((10, 10)))

My code reads a raster as an array to find the local maxima. The parameter footprint gives the search window in this case. I have tried fitting a simple function there such as:

def varvalue(x):
y = x*0.1+1
return y

However I am not sure how to pass the cell value to this function and I wonder if this is even possible. I have looked at using other functions/packages but so far without any result. Most Python related answers I find are based on one dimensional data sets and GIS answers use QGIS or Arcmap.

Does anyone have an idea or hint how to approach this problem?


texts: I have just started developing with ArcGIS for Javascript and as I'm learning Dojo I'm seeing a reason to modularize development for a  more plug-n-play concept. I can't seem to find any good articles on Best Practices for Web Map development and am wondering if there are anything resembling industry standards for setting up new web-map applications in this environment.


texts:
Each time I run the a code, see below, the ArcGIS window updates the results and the process is extremely slow. Is there a way to let python know that it doesnot have to update the results each time it has calculated something?

Code:

linklayer = "NewUSConnections"
nodelayer = "NCNodes"
AllLinkIDs = []
for row in arcpy.SearchCursor(linklayer):
    AllLinkIDs.append(row.ID)
for LinkIDs in AllLinkIDs:
    ExtractedNodeIDs = []
    query = """ "ID" = %d"""%LinkIDs
    arcpy.SelectLayerByAttribute_management(linklayer, "NEW_SELECTION", query)
    arcpy.SelectLayerByLocation_management( nodelayer, "BOUNDARY_TOUCHES", linklayer , "", "NEW_SELECTION", "NOT_INVERT")

    for row in arcpy.SearchCursor(nodelayer):
        ExtractedNodeIDs.append(row.ID)
    try:
        with arcpy.da.UpdateCursor(linklayer, ["A_NODE","B_NODE"] ) as cursor:
            for row in cursor:
            row[0] = ExtractedNodeIDs[0]
            row[1] = ExtractedNodeIDs[1]
            cursor.updateRow(row)
    except:
        pass

PS: I can see the steps followed by ArcGIS.

Updating Results...
Adding Output to map...
Drawing Background...
Labeling...


texts: Labeling using 'wordwrap' and 'wrap on character' with mutiple fields in QGIS?

texts:
I'm having a problem with multi-field labeling both using  the 'wordwrap' function and the 'wrap on character' mode. In fact I would need to use the multiline mode only on the second of the two fields (setting the wrap_lenght to optimize the number of lines) but I have not succeeded in any way.

in the first image, on the left the correct formatting, on the right the formatting that I would like to correct but avoiding too many lines as on the right in the second image. The "number" is the first field while the "name" the second field.

I tried to solve the problem using the expression:

"FIELD_1" || '' || wordwrap ("FIELD_2", 12)

without changes in formatting,
and also using a rule based labeling with:

"FIELD_1" || ''
'' || wordwrap ("FIELD_2", 12)

but completely losing the formatting.

Is it possible to solve this problem?


texts:
I finally broke down and bought Neteler and Mitasova's book, Open Source GIS: A GRASS GIS Approach, because I need to learn command-line GRASS. On trying to work through the book's examples, I'm finding that there may be a problem with my GRASS installation, or more likely with the installation of the tk framework. Here's the code and the error message that I'm getting:

  gregory@Gregorys-MacBook> gis.m
  ~ gregory@Gregorys-MacBook> /usr/local/bin/wish: line 2:
  /usr/local/bin/../../../Library/Frameworks/Tk.framework/Versions/8.5/Resources/Wish.app/Contents/MacOS/Wish:
  No such file or directory

I installed GRASS and its dependencies using the binary installers from KyngChaos, so I'm pretty confident that the TK framework got installed correctly, but the error message seems to indicate otherwise. Also, I don't know how to interpret the section of the path in the error message that reads "/../../../".

Can anyone suggest a fix for this?


texts: GlobCover 2009 dataset - importing legend and colormap values into R

texts:
Does anyone know how to import a legend for a raster layer (loaded as a .tif) into R such that it can be attached to the raster file itself and used in plotting?

In this case, I am using the GLC2009 V2.3 data, available from: http://due.esrin.esa.int/page_globcover.php

The legend is stored separately as a .xls file, as well as .lyr and .dsr files. I believe the latter 2 extensions are for ArcINFO files but regardless I have never used either in R.

I can import the .xls legend as follows easily:

    GLC_Legend <- read.xlsx("./Globcover2009_V2.3/Globcover2009_Legend.xls",
            sheetIndex=1, header=TRUE)
    head(GLC_Legend)

    #   Value                                                     Label Red Green Blue
    # 1    11         Post-flooding or irrigated croplands (or aquatic) 170   240  240
    # 2    14                                         Rainfed croplands 255   255  100

But I am lost as to how I could utilize this data frame to re-create the color map, as well as to how I could use the .lyr or .dsr files in its stead. Any know-how on importing colormaps from files for mapping purposes in R would be wonderful!


texts: Downloading data from Nikon NPL 322 Total Station

texts:
What software do I need to be able to download the data that I have collected using the TS?
What are the requirements and what is the procedure of downloading data from the TS on to a new PC Windows 10 without any prior Survey software?


texts: spTransform doesn't seem to actually reproject my data

texts:
I'm trying to create spatial point data and lay it over a shapefile. I'm starting out with simple LatLons, converting to sp using Coordinates(), changing the CRS to the same as the shapefile, and then reprojecting using spTransform into the same projection as the shapefile. However, when I check the summary of my point data afterward, it still appears to be in LatLon form. E.g.:

# Before spTransform

summary(new)

Object of class SpatialPointsDataFrame
Coordinates:
                min       max
coords.x1  40.68967  40.73516
coords.x2 -73.96243 -73.90697
Is projected: TRUE
proj4string :
[+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333
+lat_0=40.16666666666666 +lon_0=-74 +x_0=300000.0000000001 +y_0=0
+datum=NAD83 +units=us-ft +no_defs +ellps=GRS80 +towgs84=0,0,0]
Number of points: 25

new<-spTransform(new,CRS(proj4string(bb)))

> summary(new)

Object of class SpatialPointsDataFrame
Coordinates:
                min       max
coords.x1  40.68967  40.73516
coords.x2 -73.96243 -73.90697
Is projected: TRUE
proj4string :
[+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333
+lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +datum=NAD83
+units=us-ft +no_defs +ellps=GRS80 +towgs84=0,0,0]
Number of points: 25

The shapefile I'm taking the projection from looks like this:

summary(bb)

Object of class SpatialPolygonsDataFrame
Coordinates:
       min       max
x 972619.7 1030453.5
y 146772.0  208564.1
Is projected: TRUE
proj4string :
[+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333
+lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +datum=NAD83
+units=us-ft +no_defs +ellps=GRS80 +towgs84=0,0,0]

I know that I've done this before, which is why I'm scratching my head at why spTransform leaves me with point data that is still in LatLon form?


texts: Modifying WebAppBuilder "LayerList' widget to use radio buttons instead of checkboxes?

texts:
I've downloaded the developer version of the Web app builder and installed in on my machine.
I can open up the widget '\LayerList\LayerListView.js', and I have been exploring the source code.

I would like to be able to modify the code in order to have a combination of both radio buttons, and check boxes.  I am fairly new to javascript programming, so my knowledge doesn't go much beyond being able to copy / paste and hack basic code.

Functionality needed:

By clicking a check box, it will enable or disable a map service
By clicking a radio button, it will enable a map service, while turning off other map services that are also exposed as radio buttons.


texts: Vehicle Routing Problem not utilizing all available routes

texts:
I have a Vehicle Routing Problem with:

319 orders, each with specific delivery quantities
8 depots
105 routes, all with the same capacity
840 route renewals, one for each route at each depot

Upon solving, only 273 orders are serviced by 47 routes utilizing 157 route renewals.

Can anyone explain why the remaining 58 routes were not used to serve the remaining 46 orders?


texts: bringToFront() method equivalent in OpenLayers

texts:
Does OpenLayers have a method similar to bringToFront in Leaflet?
In Leaflet this works as follows. If I have a layer called geometry this method will bring that layer to the top of the other layers so it won't be covered.
An exampe in Leaflet is:
var geometry=L.geoJSON(geojson).addTo(map);
geometry.bringToFront();

I would like to do something like this in OpenLayers.


texts: Collector for GIS will not sync offline maps

texts:
We are testing collector for GIS and its proving promising however there is one major flaw i cannot figure out.

When i take a map offline and add data to it and then bring it back to the office to sync it up it runs through the process and seems to complete however when i go and check the map on arcgis online nothing shows up. Also this is bi directional, when i add data to the arcgis online map it will not sync up with the offline map. So for what ever reason the offline map will not sync up with arcgis online even though it says it does.

Any one have any thoughts? or even better a way to get the data off the ipad with out sync'ing it.

Feature service is hosted on premise servers not on arcgis online.

ArcGIS 10.3.1, Server 10.3.1, Collector 10.3.3,  Apple Air iOS 8.4

I was able to save the data using a script found in this thhread (https://geonet.esri.com/thread/92080) but i am still really curios as to what was happening, ie what offline maps wouldn't not sync.


texts: How to specify which component of a file is imported when using a PluginDataSource?

texts:
In my application a plug-in datasource has been implemented to read specifically formatted files containing a mix of data types, including geometry in the form of polygons, lines, and points. Using the "Add Data" button in ArcScene/ArcMap, the file in question can be selected to reveal the different datasets (polygons or lines, usually) contained within the file, after which they can be added individually into the scene or map.

A functionality I'm looking to add involves automation of the add data process, so I need a way to specify which (one, some or all) dataset(s) within the file get added to the scene/map. However, what I have currently always adds the polygons (which are the default) and not line or point data.

Currently I'm creating feature layers as follows:

        // open workspace
        IWorkspace workspace = workspaceFactory.OpenFromFile(System.IO.Path.GetDirectoryName(filename), 0);

        // cast as feature workspace
        IFeatureWorkspace featureWorkspace = workspace as IFeatureWorkspace;

        // get a featureclass from the workspace
        IFeatureClass featureClass = featureWorkspace.OpenFeatureClass(System.IO.Path.GetFileName(filename));

        // create a new feature layer
        IFeatureLayer featureLayer = new FeatureLayerClass
        {
            Name = featureClass.AliasName,
            FeatureClass = featureClass
        };

The issue arises when getting a featureclass from the workspace, which is done by a string containing the name. It only works properly if I include the file extension in the argument to OpenFeatureClass(string), for example "3d_data.ext" instead of simply "3d_data". I find this confusing-- shouldn't the extension only be required when opening the workspace from the file?

When browsing the file via "Add data", the polygon data is called "3d_data", and lines is "3d_data-lines". Creation of the feature class fails if "3d_data-lines" is provided to OpenFeatureClass, just as it does if just "3d_data" is used. Do I have to Open a FeatureDataset from the file first, and then create individual featureclasses from the FeatureDataset?


texts:
I inserted a pdf (image) in my work.

But the pdf is rotated 180 degrees...
My objective is to georeferencing it.

How can I rotate it prior to georeferencing it?


texts: QGIS Python - 'NoneType' object has no attribute 'addFeatures'

texts:
Trying to convert a python console script to a standalone version of it. By following the general guidelines from QGIS-documentary https://docs.qgis.org/2.14/en/docs/pyqgis_developer_cookbook/intro.html

But unfortunately, I run in a "simple" problem, but unable to solve it...
This is my code:

import sys
sys.path.append(r"/Applications/QGIS.app/Contents/Resources/python")
sys.path.append(r"/Applications/QGIS.app/Contents/Resources/python/plugins")

from qgis.core import *

import json, requests
import processing

# supply path to qgis install location
QgsApplication.setPrefixPath(r"/Applications/QGIS.app/Contents/MacOS", False)
# create a reference to the QgsApplication, setting the # second argument to False disables the GUI
qgs = QgsApplication([], False)
# load providers
qgs.initQgis()

geometry_decode = [(0.0,0.0),(1.2,1.3)]

geometry = [QgsPoint(point[1],point[0]) for point in geometry_decode]

line = QgsGeometry.fromPolyline(geometry)
route = QgsVectorLayer("LineString?crs=epsg:4326&index=yes","route","memory")
pr = route.dataProvider()
print(route)
route.startEditing()

feature = QgsFeature()
feature.setGeometry(line)

pr.addFeatures([feature])

# Write your code here to load some layers, use processing algorithms, etc.
# When your script is complete, call exitQgis() to remove the provider and # layer registries from memory
qgs.exitQgis()

This is the error-message:

Error in sys.excepthook:
Traceback (most recent call last):
  File "/Applications/QGIS.app/Contents/Resources/python/qgis/utils.py", line 196, in qgis_excepthook
    showException(type, value, tb, None, messagebar=True)
  File "/Applications/QGIS.app/Contents/Resources/python/qgis/utils.py", line 107, in showException
    open_stack_dialog(type, value, tb, msg)
  File "/Applications/QGIS.app/Contents/Resources/python/qgis/utils.py", line 142, in open_stack_dialog
    iface.messageBar().popWidget()
AttributeError: 'NoneType' object has no attribute 'messageBar'

Original exception was:
Traceback (most recent call last):
  File "test_standalone_copy.py", line 33, in <module>
    pr.addFeatures([feature])
AttributeError: 'NoneType' object has no attribute 'addFeatures'


texts: Standalone Application QgsCoordinateReferenceSystem fails without OSGeo environment

texts:
I am writing a simple standalone QGis application. My code is like this:

class MainForm(QtGui.QMainWindow, Ui_MainWindow):
    def __init__(self):
        super(MainForm, self).__init__()
        self.setupUi(self)

        self.qgisMapCanvas.mapRenderer().setProjectionsEnabled(True)
        self.canvas_crs=core.QgsCoordinateReferenceSystem(3857)
        assert self.canvas_crs.isValid()
        self.qgisMapCanvas.mapRenderer().setDestinationCrs(self.canvas_crs)
...

The Python script works pretty well. Then I built .EXE package using pyinstaller, and the .exe package works well too within the OSGeo4W64 environment. However, once I move the package into another computer without OSGeo4W64 environment, it stops working. The Cmd error message shows that something seems to be wrong with QgsCoordinateReferenceSystem:

Qt: Untested Windows version 10.0 detected!
Running In An Application Bundle
Error in sys.excepthook:
Traceback (most recent call last):
  File "qgis\utils.py", line 196, in qgis_excepthook
  File "qgis\utils.py", line 107, in showException
  File "qgis\utils.py", line 142, in open_stack_dialog
AttributeError: 'NoneType' object has no attribute 'messageBar'

Original exception was:
Traceback (most recent call last):
  File "CallMainWin.py", line 324, in <module>
    win = MainForm()
  File "CallMainWin.py", line 49, in __init__
    assert self.canvas_crs.isValid()
AssertionError
[64268] Failed to execute script CallMainWin

If I remove the statement

assert self.canvas_crs.isValid()

from program, it will be working, the map will be displayed in the qgsMapCanvas, except that the CRS of the map is apparently wrong.

Could anyone help me out with this please?


texts:
At my job, most people have ArcReader and the 2 GIS staff will create PMF for people so they can use them for projects and reporting. Lately we've been having an issue where the ESRI basemaps will have chunks missing and a JPG logo (bottom middle) wont show up. These all happen when our colleagues are trying to print to PDF from ArcReader.

Does not happen in ArcMap, just Reader.

Any ideas as to what could be causing this?

This is an 11x17 map layout...logo issue doesn't happen in 8.5x10 layout.


texts: Calculating Specific Catchment Area in QGIS?

texts:
I have tried to calculate the Specific Catchment Area in QGis from Specific Catchment Area and Flow Width Tool available through Saga. However, I got just SCA product with valour 0. I got a reasonable product with Flow width.

Can anyone help me?


texts: Pass data through hyperlink to a separate sub form in Survey123

texts:
I have a parent form in survey123 that generates a unique ID when teh form is opened with the calculation: uuid() .
This ID is in a field called "inspectionid"
I have subsequent hyperlinks to "subforms" within this parent form. THey are distinctly separate forms that are more detailed inspections of parts of the parent form.  Basically the parent form say "is "X" present? Y/N"  If yes,here's a link to a sebform to describe "X". I've since been asked to link the subform inspection events to the parent form so I'd like to pass that auto generated "inspectionid" value to the subform into a field called "parentinspectionid"
The current hyperlinks are formatted like this and they launch fine (we may want to modify them later to open in the Survey123 app rather than the web app but for now, this is what I'm comparing against.
 Related Subfform
I've tried this:
 Subform Name
and other variations that were just totally misguided.
But the Survey123Connect preview does not show this as a clickable link but rather as static text so I can't tell if I'm even on the right track.
Can I format the hyperlink to pass the "inspectionid" from the main form to to the "parentinspectionid" of the sub forms?


texts: Creating button to clear all layers from data frame using ArcObjects and VB.NET?

texts:
I'm trying to create a button using VB.NET and ArcGIS 10 to clear the layers of my data frame.

I know how to add layers, but how do I delete/clear them?


texts:
I´m using entropy and GLCM function of Google Earth Engine (https://developers.google.com/earth-engine/guides/image_texture). I want to use NDVI band, obtained from landsat 8 image band algebra:
 //Agregar bandas
 var addVI = function(image) {
 var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI');
   return image.addBands(ndvi);
 };
 var withVI = filteredCollection.map(addVI);
 print ("conBandas", withVI);

//COMPOSICIONES
//Composicion de imagenes 2014
var cverano2014 = withVI.filterDate('2013-12-01', '2014-03-31');
print ("Coleccion2014", cverano2014);
//Buscar Imagen Limpia 2014
var listImages = cverano2014.toList(cverano2014.size());
var il2014 = ee.Image(listImages.get(0)).clipToCollection(table);
Map.addLayer(il2014, trueColour, "ImagenLimpia2014");

// Get the NIR band.
var iv = il2014.select('NDVI');
// Define a neighborhood with a kernel.
var square = ee.Kernel.square({radius: 3});
// Compute entropy and display.
var entropy = iv.entropy(square);
print('entropy', entropy);
Map.addLayer(entropy,
             {min: 1, max: 5, palette: ['0000CC', 'CC0000']},
             'entropy');

// Compute the gray-level co-occurrence matrix (GLCM), get contrast.
var glcm = iv.glcmTexture({size: 3});
print(glcm)
var contrast = glcm.select('NDVI_contrast');
print('constrast', contrast);
Map.addLayer(contrast,
             {min: 0, max: 1500, palette: ['0000CC', 'CC0000']},
             'contrast');

Well, when I run the code, console show me:
Image (Error)
Image.glcmTexture: Only 32-bit or smaller integer types are currently supported.
How do I reduce bit of my bands, or transform band to type integer?


texts: Table joined to a layer doesn't appear when uploaded to web map

texts:
I have a layer in my map to which I’ve joined an Excel spreadsheet that shows additional attributes.

Once the join was successful, I exported the joined layer to my map’s file geodatabase and added it back to my map.

From there I wanted to upload this joined layer to a web map (“share as a map service”) which should automatically update the web map application, but the extra fields from the join are not visible.

Is there a reason why this is occurring? Is there a suitable workaround?


texts: Save geodataframes from list to separate shapefiles

texts:
I have a list of geodataframes which I imported from a folder full of shapefiles:
import glob
import geopandas as gpd
shapefiles = glob.iglob('E:/folder/shapefiles/*.shp')

gdfs = [gpd.read_file(file) for file in shapefiles]

I then clipped them all:
clipped = [gpd.clip(shapefiles, boundary) for shapefiles in gdfs]

But now I need to save each geodataframe (list item) from the list to a shapefile of its correct name in an output folder. I used this code to create another list called filenames of the shapefiles names with the .shp at the end:
filenames = []
for names in items:
    if names.endswith(".shp"):
        filenames.append(names)
print(filenames)

But I dont know how to save each item in the clipped list of geodataframes using its corresponding name from the filenames list
Basically I want to perform this task:
clipped[0].to_file("E:/test_output/"+filenames[0])
clipped[1].to_file("E:/test_output/"+filenames[1])
clipped[2].to_file("E:/test_output/"+filenames[2])

etc etc
But to perform it iteratively and automatically. I need something like this:
for things in clipped:
    things.to_file("E:/test_output"+(for stuff in filenames return filenames[]))

But with correct syntax and that actually works, Thanks!


texts: Selecting features within Polygon from another layer using QGIS?

texts:
I'm using QGIS 1.8 with a PostGIS database and I have three PostGIS-layers in QGIS a lines layer, a points layer and a polygon layer.

What I'm trying to do is to select a polygon in the polygon layer and mark all points within this polygon.

For me the only way to do this is to use the button "Select by Polygon" but in this way, I cannot use an existing polygon for selection.


texts: Running Python outside of QGIS - (ImportError: No module named qgis.core)

texts:
I'm trying to run simple scripts on QGIS from outside the QGIS environment - as here

How to run a simple python script for QGIS from outside (e.g. Sublime Text)?

My code:

import sys
sys.path.append("C:\\Program Files (x86)\\QGISBrighton\\apps\\Python27")
import qgis.core

My output:

ImportError: No module named qgis.core

My qgis_core.dll is at

  C:\Program Files (x86)\QGISBrighton\apps\qgis\bin

QGIS install is here:

  C:\Program Files (x86)\QGISBrighton

In Control Panel > System > Advanced

PYTHONPATH set to

  C:\Program Files (x86)\QGISBrighton\apps\Python27

Path set to

  C:\Python27; C:\Program Files (x86)\QGISBrighton\apps\Python27; C:\Program Files (x86)\QGISBrighton\apps\qgis\bin

I've tried multiple permutations of setting my Path, PYTHONPATH and appending to sys.path on the fly (as detailed in many answers here).

I'm hoping by providing the above information someone can give me explicit instructions on how to get past

  ImportError: No module names qgis.core


texts: Calculating number of plants on orthophoto

texts:
I have orthophoto map of agricultural field and I need to count plants there. I have also multispectral map (NDVI) but I think it's not very usable here.

I know the best would be to have point layer on every single flower and then count these points. Please, is there any way how to create this point layer based on classification or something else, in QGIS or ArcMap? I know there are tools which are good for tree counting, but they are based on height canopy model and it's not very usable in my case as the height of plants is small. Below is the example of my orthophoto:


texts: Duplicate entries in table list when adding PostGIS layer in QGIS

texts:
When adding PostGIS layer in QGIS (version does not matter) i observe duplicate entries in the list as can be seen below (e.g. ax_gebaeude, ax_flurstueck, ax_gehoelz...), one with a valid geometry, the other(s) with a warning and a 'select geometry type' hint.

Adding the ones with valid geometry to the project causes no issue, but i am not sure what this means.

Can anyone clarify this behaviour?

Added/Edit:

geometry_columns view:

and everything's fine in browser panel:

and another behaviour in the DB Manager plugin. For geometry_types 'GEOMETRY' are marked with a '?' and it says 'no entry in geometry_columns' (wich is definitively not the case) but preview is no problem:


texts: Different extents of SOC and BDOD VRT files

texts:
I downloaded the SoilGrids layers from the WebDAV and reprojected them to EPSG:4326 using the following command:
gdalwarp -t_srs EPSG:4326 -of VRT ./sand_60-100cm_Q0.5.vrt ./sand_60-100cm_Q0.5_4326.vrt

Subsequently, I used GDAL to translate the VRT files to GeoTiffs:
gdal_translate -of GTIFF -a_srs '+proj=longlat +datum=WGS84 +no_defs' ./sand_60-100cm_Q0.5_4326.vrt ./sand_60-100cm_Q0.5_4326.tif

This worked out all fine, but I realised that SOC and BDOD are shifted by almost a whole pixel. Because I applied the mentioned GDAL commands to all files equally, I am sure that it must be an issue with the initial files. I went back and checked the VRTs and realised that the longitudinal extents of SOC and BDOD are three pixels smaller than all the other VRTs. This seems to have produced my issue.
Is this something you want to correct? Or can you suggest another workaround, so that all layers are projected equally?


texts: Get geographic coordinates based on pixel coordinates?

texts:
I'm working with Mask R-CNN to locate fires.

I have the following image. I have info about the geographic center coordinates of the image, and the orientation of the image in degrees let's suppose this image is oriented to north. I have also the center coordinates of the masks. How can I convert the center coordinates of the mask in pixels to geographic coordinates.


texts: Getting features in source for GeoJSON vector layer, and then altering them?

texts:
I have managed to do everything I need for getting and reading features using the technique from https://github.com/openlayers/ol3/blob/master/doc/faq.md#why-arent-there-any-features-in-my-source (Note that the link at the top of the FAQ page is not quite correct). (I am using OpenLayers 3.15)

The trouble comes when I try to actually alter a feature (eg by changing a property) then the .on('change') seems to fire again.  This means that if I have 6 features they are processed repeatedly, thus:  0, 01, 012, 0123, 01234, 012345. This can confuse things, repeat indefintely, or just take a long time.

I have tried several ways round this, starting with using .once rather than .on, but either the routines fire repeatedly (sometimes) or they do not fire at all (sometimes). I have also tried using ol.Observable.unByKey(key). You will realise that I am rather struggling in the dark.

You can see some of my attempts in my rough little demo at http://mapping4ops.org/ShowMapsDev/foreachfeature.html. What I have tried to do in this is to hide alternate features (by setting a 'hidden' property) in a vector layer containg 20 features, using 4 different techniques. The essence of the code is below, but the full code has console.log statements, loop limits etc. (Operation of this is best observed by seeing the console.log messages, using eg Chrome Inspector.)

Hide1) uses simple getSource and getFeatures:
[While this works well in some circumstances most times it just returns no features]

var source = vectorLayer.getSource();
if (source.getState() === 'ready') {
    var features = source.getFeatures();
    for(var i=0; i< features.length; i++) {DO THE HIDING}}

Hide2) uses simple getSource and forEachFeature
[Similar to Hide1, most times it just returns no features]

var source = vectorLayer.getSource();
if (source.getState() === 'ready') {
    source.forEachFeature(function(feature) {DO THE HIDING})}

Hide3) uses .getSource().once('change') and forEachFeature
[I thought using .once would get over the issue of repeated firing, but it does not seem to]

vectorLayer.getSource().once('change', function(evt){
    var source = evt.target;
    if (source.getState() === 'ready') {
        source.forEachFeature(function(feature) {DO THE HIDING})}})

Hide4) uses .getSource().on('change') and forEachFeature
[This is where I started from as this works fine when the fetures are not changed]

vectorLayer.getSource().on('change', function(evt){
    var source = evt.target;
    if (source.getState() === 'ready') {
        source.forEachFeature(function(feature) {DO THE HIDING})}})

If I try hiding (and unhiding) the features, what happens depends on the order I do things (reloading the page between each attempt). For example:

Case A) Using just Hide1 or Hide2, the hiding and unhiding happens OK - BUT I know I cannot use this in other circumstances because getsource is asynchronous.

Case B) 3 then 1: Hide3 first (nothing happens), and then Hide1 (hiding happens OK, but the Hide3 fires now and then the Hide1)

Case C) 1 then 3: Hide1 (hiding happens OK), Unhide1 (OK), and then Hide3 (nothing happens)

Case D) 4 then 2: Hide4 first (nothing is changed), and then Hide2 (the Hide4 fires now repeatedly, thus:  0, 01, 012, 0123, 01234 ..., and I abort it when it reaches a limit)

I did think I found the answer at https://freeopengis.blogspot.co.uk/2015/09/openlayers-3-load-start-load-end-events.html, but I do not understand it, and the fiddle does not work.

What I want is a way to access the features in the layer and manipulate them, that will work every time and not process the features more than once.


texts: QgsLayoutItemScaleBar.setUnits does not work

texts:
I'm writing a plugin that allows to automatically export a map with a predefined layout, based on the content of the canvas.
I need to have a scale bar in nautical miles, so in my code I add a QgsLayoutItemScaleBar (code below).
scaleBar = QgsLayoutItemScaleBar(c)
scaleBar.setLinkedMap(global_map)
scaleBar.setStyle('Line Ticks Up')
scaleBar.setUnits(QgsUnitTypes.DistanceNauticalMiles) # this line to set the unit of scale bar in Nautical Miles
scaleBar.setSegmentSizeMode(1)
scaleBar.setNumberOfSegmentsLeft(0)
scaleBar.setNumberOfSegments(2)
scaleBar.setMinimumBarWidth(15)
scaleBar.setMaximumBarWidth(40)
sbf = scaleBar.font()
sbf.setPointSize(9)
scaleBar.setUnitLabel("MN")
scaleBar.setHeight(2)
scaleBar.setPos(x + 2 * (w // 3) + 25, y + h - 20)
c.addItem(scaleBar)
QgsMessageLog.logMessage('SCALE BAR UNIT: {0}'.format(scaleBar.units())) # result is Nautical Miles

It seems that this line scaleBar.setUnits(QgsUnitTypes.DistanceNauticalMiles) doesn't work: when exporting my map, the scale bar is in map units. I tried to move this line of code, to use DistanceKilometers but still no change. When I ask scaleBar.units(), it's in Nautical miles so the information is understand but not apply.
Has anyone experienced this problem before? Could it be a bug?


texts: Calculate stats with GDAL Python or rasterio for each bands

texts:
I have MODIS EVI file stacking with 10bands.
I use Python 3 in Spyder.

I'm trying to rastersize shapefile and do mask to my MODIS EVI.

This is my script :

import os
import glob
from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt
import scipy.ndimage
import pandas as pd
import datetime as dt
import ogr
import rasterio as rio
import fiona
import rasterio
import rasterio.mask

inDir = 'D:/projetp/Doy_2000/'
driver = ogr.GetDriverByName('ESRI Shapefile')
os.chdir(inDir)

with fiona.open("C:/Doy2000/1996n.shp", "r") as shapefile:
    shapes = [feature["geometry"] for feature in shapefile]

with rasterio.open("D:/projetp/Doy_2000/MODIS_EVI_2000.tif") as src:
    out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)
    out_meta = src.meta

out_meta.update({"driver": "GTiff",
                 "height": out_image.shape[1],
                 "width": out_image.shape[2],
                 "transform": out_transform})

with rasterio.open("MODIS_MASK_2000.tif", "w", **out_meta) as dest:
    dest.write(out_image)

when i whrite out_image in spyder i hade this resultas

print(out_image)
[[[   0    0    0 ...    0    0    0]
  [   0    0 1682 ...    0    0    0]
  [   0 1841 1888 ...    0    0    0]
  ...
  [   0    0    0 ... 2180    0    0]
  [   0    0    0 ... 4033    0    0]
  [   0    0    0 ...    0    0    0]]

 [[   0    0    0 ...    0    0    0]
  [   0    0 1659 ...    0    0    0]
  [   0 1829 1830 ...    0    0    0]
  ...
  [   0    0    0 ... 3953    0    0]
  [   0    0    0 ... 4627    0    0]
  [   0    0    0 ...    0    0    0]]

 [[   0    0    0 ...    0    0    0]
  [   0    0 1706 ...    0    0    0]
  [   0 1547 1740 ...    0    0    0]
  ...
  [   0    0    0 ... 4419    0    0]
  [   0    0    0 ... 4564    0    0]
  [   0    0    0 ...    0    0    0]]

 ...

 [[   0    0    0 ...    0    0    0]
  [   0    0 1998 ...    0    0    0]
  [   0 2255 1841 ...    0    0    0]
  ...
  [   0    0    0 ... 3339    0    0]
  [   0    0    0 ... 3501    0    0]
  [   0    0    0 ...    0    0    0]]

 [[   0    0    0 ...    0    0    0]
  [   0    0 1992 ...    0    0    0]
  [   0 3337 3337 ...    0    0    0]
  ...
  [   0    0    0 ... 2058    0    0]
  [   0    0    0 ... 2080    0    0]
  [   0    0    0 ...    0    0    0]]

 [[   0    0    0 ...    0    0    0]
  [   0    0 7913 ...    0    0    0]
  [   0 7659 7294 ...    0    0    0]
  ...
  [   0    0    0 ... 5014    0    0]
  [   0    0    0 ... 4626    0    0]
  [   0    0    0 ...    0    0    0]]]

And when I visualize my out_image in QGIS I have this result

Do I need to delete the value 0 or not? And how can I show my out_images in Python?

The values are so big.  How can I apply a scale factor *1/10000 to have values (-1,1)?
- How can I calculate statistics ( mean of pixel of my out_image just inside the area not in black area  )! Because I need just mean of pixel for all bands to export after in dataframe!


texts: Splitting up polygon with grid using QGIS

texts:
I am using QGIS.

I have a county map as a single polygon, it has one FEATURE in which one FIELD is Area.
I have a grid which also a polygon that divides the county into squares and each grid square is a FEATURE which all have a FIELD called Area.

I want a new vector layer with each grid-square as a FEATURE in which the area of the country falling within that grid square is listed under the Area FIELD.
I keep trying to Clip, Extract and Union and I keep getting an "invalid geometry" error.  I know this is fairly straight-forward but I keep stumbling over this error time and time again.


texts:
I had a layer that I split by fields. The resulting files were given a name. I updated the names of the resulting layers in QGIS but have realized that the shapefiles have not been updated accordingly.

Is there a way to batch synchronize shapefile names with the layer names?

I know I can save as, but that is cumbersome with 75+ layers.


texts: Merging shapefiles from subfolders to corresponding new folder using ArcPy?

texts:
The dataset have shapefiles saved in different folders, as shown below.

The numbers of shapefiles in each folder are not equal. The content of one folder is shown below.

I need to merge all shapefiles in each folder and save the merge results in their corrsponding folder based on the FacIDxx, as shown below.

The following code works fine for the requirements mentioned above, but it takes lots of time when the dataset folder has many shapefiles to be merged.

import arcpy

#resetting the workspace
def find_all_fcs(workspace):
    prev_workspace = arcpy.env.workspace

    arcpy.env.workspace = workspace
    fcs = arcpy.ListFeatureClasses()

    arcpy.env.workspace = prev_workspace
    return fcs

filein0=r"C:\Users\TL\Desktop\812test\SplitAllRoutes"
arcpy.CreateFolder_management("C:/Users/TL/Desktop/812test", "Merge")
fileinIntersect = r"C:\Users\TL\Desktop\812test\Intersect"
fcs = find_all_fcs(filein0)
for fc in fcs:
    tempname_Intersect = fileinIntersect+'\\'+fc[:-4]
    fcs_Intersect = find_all_fcs(tempname_Intersect)
    arcpy.CreateFolder_management(r"C:\Users\TL\Desktop\812test\Merge", fc[:-4])
    input_Merge = []
    for fc_Intersect in fcs_Intersect:
        tempname_toMerge = tempname_Intersect+'\\'+fc_Intersect
        input_Merge.append(tempname_toMerge)
    output_Merge = r"C:\Users\TL\Desktop\812test\Merge"+'\\'+fc[:-4]+'\\'+fc[:-4]+'.shp'
    arcpy.Merge_management(input_Merge, output_Merge)

I wonder if someone can help me improve the efficiency of the code.


texts: How to make some labels BOLD using Expression Builder QGIS

texts:
I am using QGIS 2.14 and I am trying to make some of my labels bold to stand out. I'm not sure how to do this using the expression builder. I know this is wrong and I was wondering if someone could let me know how to correct it:

CASE WHEN
"Style" like '%2016%' then  Bold
ELSE Normal
END

"Style" is my layer of points, I need the 2016 labels in bold.

My only solution I have at the moment is making a 'Bold' column in my attribute table and using this field to select features to be bold. I think this clutters up my attribute table and different maps require different labels bold.


texts: Offset after clipping an image using rasterio.mask (illustrated with figures)

texts:
I want to clip an image named img_A of dimensions 14000x10000 using a polygon shape in order to get the same dimensions than another image named img_B of dimensions 1925x1372 (see figure 1 below). The polygon shape has the same extent than img_B. Pixel size is 0.1 m for img_A and img_B.

One important thing is that the pixels of img_A are not aligned with the pixels of img_B (figure 2). I mean there is a small offset between the pixels of both images (3 cm).

So I used rasterio.mask to clip img_A and obtain img_B extent using that code (figure 3):

import fiona
from rasterio.mask import mask

polygon = polygon.shp
raster = img_A.tif
dst = result.tif

with fiona.open(polygon, "r") as shapefile:
    geometry = [feature["geometry"] for feature in shapefile]

with rasterio.open(raster) as src:
    out_image, out_transform = mask(src, geometry, crop=True)
    out_meta = src.meta.copy()
    out_meta.update({"driver": "GTiff",
                     "height": out_image.shape[1],
                     "width": out_image.shape[2],
                     "transform": out_transform})

with rasterio.open(dst, "w", **out_meta) as dest:
    dest.write(out_image)

However, the final result is an image of dimensions 1926x1373. Moreover, the new image is slightly offset to the South-West, corresponding to the offset cited previously (3 cm), which means the image is outside of the polygon shape to the West and to the South (figure 4).

How can I obtain a well placed image with dimensions 1925x1372?

Figure 1 - Input dataset

Figure 2 - Pixel offset

Figure 3 - Overview of the result after clipping

Figure 4 - Problem


texts: Cannot open .zip file from browser in QGIS 3.6.3 in Windows 7?

texts: In QGIS 3.6 on Windows 7, I can not open the .zip file from Browser. But it can open  through Layer--Add layer--Add vector layer.


texts: Line features with equal parameters different visualisation

texts:
I have three different line features on a map all three with equal visualisation parameters (line width, line style, etc) but with one of them the dots/small lines that form the feature appear to be longer, so it seems to be a solid line when it is supposed to be a dotted line.

Some info that might be useful to solve this issue: every time I've had this problem with a file it's been created with the GPX segment importer tool and after I've saved them as a shapefile.


texts: Specifying min and max values for image while exporting from Google Earth Engine

texts:
I have this small code for getting nightlight image for a particular city using VIIRS Nighttime dataset.

var district = ee.FeatureCollection('ft:1PA2zwArj8EsplrX9eMxJ2H_TICyyx855KPnbJhC1','geometry')
.filter(ee.Filter.eq('name','Gurgaon'));
var nighttime = ee.Image(ee.ImageCollection('NOAA/VIIRS/DNB/MONTHLY_V1/VCMCFG')
              .select('avg_rad')
              .filterBounds(district)
              .filterDate('2013-05-01', '2017-05-31')
              .median())
var nighttimeVis = {min: 0.0, max: 60.0};
Map.addLayer(nighttime.clip(district), nighttimeVis, 'Nighttime');

Now I wish to export this image to my Google Drive. Usually I use the following script for exporting:

Export.image.toDrive({
  image: nighttime.clip(district),
  description: 'Nighttime_Image',
  maxPixels: 499295920080,
  scale: 30,
  region: district
});

But for this particular image, I need to specify the min and max values like I have done while using Map.addlayer (for viewing the image) by using nighttimeVis vector. How can I specify the same max and min before or in my Export command for this image?


texts: Extract LINESTRING ZM geometry (XYZM) in QGIS

texts:
I have a linear layer of roads (every segment divided into many parts) with linestring geometry. I would like to generate a WKT linestring geometry in attribute table - coordinate X Y together (!) with Z and M values (XYZM for each nodes)

Example of a correct result: LINESTRING ZM ((540762.573877557 519541.938412964 0 0, 540715.013082288 519508.218571241 0 58))

The main problem is how to generate a X and Y and Z and M values for each nodes with a single query (!) in field calculator in QGIS or at export to csv.
I need this as one field in the attribute table.


texts: Calculating surface using Lambert correction?

texts:
I'm using Lambert conformal conic, so we have distortions along any given parallel, I would like to get the area used to solve this distortion.

I have calculated the area using Shoelace formula, surprisingly it gave me the same area exactly as the area calculated by Arcgis, that means that Arcgis doesn't calculate Lambert Correction, this correction is necessary to correct distortions.

Please correct me, if this is wrong.

I would like to calculate the surface using Lambert correction, if you can help me to find the formula of Lambert Correction?

I'm using Arcgis 10.1, I would like to solve the problem using ArcPy.


texts: How can I fill the category field in a raster using GRASS

texts:
I have a tiff file with land use data. The land use class in in the value field, but the label field is empty.
Now I want to find out which is the area covered by each land use class, but GRASS does not find any category data.

So I wanted to copy the data in the value field to the category field, but I have not yet found any way of doing this.


texts: How do I update common fields in multiple feature classes with ArcPy?

texts:
I have a piece of code that I can use to update several fields in a specific feature class.  I would like to use this code to update all of the feature classes in the data set.  The data set would would be WATER and the feature classes would be WATERLINES, HYDRANTS, VALVES, AIR_RELEASE, ETC)  Each of the feature classes have the fields I would like to update:  LASER_LINK, ASSETID.

Here is the code i use to update a single feature class called WATERLINES:

# Python Script: updateAttributes.py
# This script autopopulates the LASER_LINK and ASSETID fields.

print "Start script: "

# Import the ArcPy package
import arcpy
import os

# feature class to be updated
fc = 'C:\Users\cbgibson\AppData\Roaming\ESRI\Desktop10.2\ArcCatalog\TEST_WSdata.sde\waterlines'

# set the current workspace
workspace = os.path.dirname(fc)

# start edit session
edit = arcpy.da.Editor(workspace) # Start an edit session. Must provide the worksapce.
edit.startEditing(False, True) # Edit session is started without an undo/redo stack for versioned data(for second argument, use False for unversioned data)
edit.startOperation()# Start an edit operation

updateFieldsList =  ["Laser_Link", "Project_nu", "IDNAME", "IDNUMBER", "IDNUMBER2", "ASSETID"]

with arcpy.da.UpdateCursor(fc, updateFieldsList) as cursor:
    for row in cursor:
         row[0] = "http://lf/weblink8WaterandSewer/search.aspx?dbid=0&searchcommand={[]:[WS Project Number]=" + '"' + str(row[1]) + '"}'
         row[5] = row[2] + str(row[4])
         cursor.updateRow(row)

# exit edit session
edit.stopOperation() # Stop the edit operation.
edit.stopEditing(True)# Stop the edit session and save the changes

print "Finished attribute update: "

I did figure out how to return all of the feature classes from the WATER data set by using this:

import arcpy
from arcpy import env

env.workspace = "C:\Users\cbgibson\AppData\Roaming\ESRI\Desktop10.2\ArcCatalog\TEST_WSdata.sde"
fc = arcpy.ListFeatureClasses("", "", "Water")

print fc

Im just not sure how to return all of the feature classes and run the update against every feature class.  Im sort of a newb at this and have been reading through quite a bit of material.  Enough to get me this far anyway.  Im just stuck!!!


texts: Shapefiles in geodatabase not visible in ArcGIS Pro but can open in QGIS

texts: I can open a shapefile from a geodatabase in QGIS but it isn't visible in ArcGIS Pro 2.5.0. Has anyone else had this problem.


texts:
I've done a layer stack of some Sentinel 2 data in QGIS. Is there a way of renaming the bands in the style menu (in order to reflect what the bands are showing, Red, Green, NIR etc)?

So far they're just called "Band 01" through to "Band 12".


texts: Assigning raster values to spatial point using R?

texts:
I have a SpatialPointsDataFrame and create a raster object. I assign each raster cell a value and like to extract that value for each spatial point. I.e. if a point lies within a raster, I'd like to end up with that raster cell number in a new column in my SpatialPointsDataFrame.

This is how I proceed:

library(raster)
library(rgeos)

I define a crs and this might be the problem (?).

crs <- "+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"

I've uploaded a subset of my SpatialPointsDataFrame here. To illustrate my problem I need a shapefile of South Africa:

za.spdf <- readRDS("pointsdf.Rds")

## Shapefile of south africa
za <- getData('GADM', country='ZAF', level=0)
za <- spTransform(za, crs)

In order to inspect, see this plot, everything seems fine:

plot(za)
plot(za.spdf, add=T, col="red", pch=20, cex=2)

I then create my raster object and assign values to it:

## Raster over south africa (100x100km)
r <- raster(extent(za))
res(r) <- c(100000, 100000)
values(r) <- 1:ncell(r)
proj4string(r) <- crs

I thought I could just transform the raster object into a shapefile and use over. However, this returns only NA values...

grid100 <- rasterToPolygons(r, n = 4, na.rm = TRUE)
spatial100 <- sp::over(za.spdf, grid100)
spatial100
          layer
 14777    NA
 14785    NA
 14793    NA
 14796    NA
 14800    NA

When plotting, the rastering and the map overlay nicely. However, where are my points?

plot(grid100)
plot(za, add=T)
plot(za.spdf, add=T, col="red", pch=20, cex=2)

I couldn't figure out what the problem is. I highly suspect that my SpatialPointsDataFrame does not correspond with the map and therefore with the raster.


texts: How to create a Digital Surface Model using ASTER level 1A/B data?

texts:
I am working on estimating Canopy height of mangrove. I want to employ ASTER level 1A or 1B data to generate Digital Surface Model (DSM). But after I searched it in Google, I cannot find an existed method to generate a DSM. So, how do I generate a DSM using ASTER level 1A/B data? Or, can the ASTER data (level 3 DEM) set be considered as a digital surface model (DSM) which provides heights of the surface objects?

Discussion:
After thinking for a while, I understand ASTER has theer subsystems: VNIR (visible and Inferred band), SWIR and TIR. The VNIR has the 3N and 3B band which can be applied to generate DEM. Moreover, DSM is a subset of DEM. Since the channels cannot penetrate canopy to reflect ground, the level 1 data of ASTER indicates the surface height information instead ground information. Thus, the raw data can generate digital surface model.


texts: Applying specific Rendering to Raster and export to JPEG?

texts:
We have a specific process which generate an interpolated single band 32-bit Float Raster in Erdas Imagine (.IMG) format. This process generates an image every hour.

We wish to apply a specific stretched renderer, which stretches from Green to Red, to the single band raster and then save it as a JPEG file.

I have found that it can be done through the UI in ArcMap by applying the .lyr file, and then exporting to JPEG 2000. The Jpeg 2000 is converted to the regular JPEG image by python.

We now want to automate this process and run it using some open source or free tool.

Is it possible to do something like this in GDAL or some other command line tool?


texts: Display only slopes facing west using ArcGIS for Desktop?

texts:
In ArcGIS for Desktop I have to Derive and symbolise areas considered susceptible to avalanche risk based on the following simple rules:

areas with a slope between 20 and 45 degrees,
facing towards the prevailing wind (west in this case)
not in forested areas.

By far I manage to display the Slopes and Aspects using tools from Spatial Analyst.

But how I can display the slopes in 3 different categories where they all facing west?

Photo attached shows what I manage to come up with


texts: Result of intersect between two shapefile is empty

texts:
I'm trying to intersect one shapefile with an other one in QGIS. The problem is, whatever I do, the result is always empty, (the files I'm working on are not empty).

I tried to change the scr, delete the files, change the order, but nothing, the result is always empty. I know it happened to some people, but none of the solution worked for me.

Do you have any idea of what might cause this issue?

I can give you my files, but they are just two basics shapefiles so I don't know if it's very useful.


texts: Attribute field name is not case-sensitive anymore after creating new layer

texts:
I joined a sqlite layer with a csv table (each attribute field name starts with a capital letter). When I create a new layer out of them ("save as", e.g. as a new sqlite layer with all data combined), the attribute field names in the new layer are changed to lower case. The entries in the table are not affected, though - only the attribute field names (= column title, to say it differently).

Does anyone have an idea how I can save a join of a spatialite layer and csv table and at the same time keep the original attribute field name that is case-sensitive? I also tried to use a different codification, but this did not help.


texts: Leaflet geojson styling leaves gaps between polygon

texts:
I'm styling a geojson polygon composed of multiple closely placed vertical segements which doesnot have a gap in between (i've checked the topology). I want to display just the fill and without stroke. But when doing so, using following code:

function getColor(d) {
    if (d<=40){return 'green';}
    if (d<=60 && d>40){return 'purple';}
    if (d>60){return 'orange';}
    }

    function style(feature) {
        return {
            weight: 0,
            opacity: 0,
            //color: 'white',
            //dashArray: '3',
            fillOpacity: 1,
            fillColor: getColor(feature.properties.id)
        };
    }

The result is as follows:

You can see the white lines in between polygons.


texts: Multiple rule based expressions with 2 layer styles in Atlas

texts:
I have 18 transects over a geographical extent and I am applying atlas to display a map for each transect.

All is going well and I'm venturing into expressions and rule based rendering $ID = @atlas_featureid . 'Transect' is my coverage layer displaying individual transects (successful so far) without other transects bleeding into each others layouts.

However the problem arises when I try to include another layer... a point shapefile with start and end positions for each transect line.

I've tried to apply rule based rendering on the point shapefile so that the atlas displays transect with corresponding start and end positions but other start and end positions are visible in transect features (see image below as an example).
I would also like to change the symbology of start and end points and label them too but I'm confusing myself with even more rules.

Any pointers to how I can navigate my current issues and achieve handling multiple layers with multiple rule based styles and labels while using atlas?


texts: Solving ‘Null geometry…’ error during distance matrix operation?

texts:
I have a tab delimited file with data about earthquakes epicenters and their intensities (signif.txt). I also have a shape file with point data representing populated places of the world (ne_10m_populated_places_simple.shp). I am trying to perform a ‘Distance matrix’ operation to find the nearest populated place for each point in the earthquake layer. However, performing ‘Vector> Analysis Tools> Distance Matrix…’ results following error message:

Traceback (most recent call last):
File "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\algs\qgis\PointDistance.py", line 141, in processAlgorithm
matType, nPoints, feedback)
File "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\algs\qgis\PointDistance.py", line 197, in linearMatrix
featList = index.nearestNeighbor(inGeom.asPoint(), nPoints)
ValueError: Null geometry cannot be converted to a point.

Execution failed after 17.02 seconds

I am following the steps mentioned here and then I have even tried to convert the delimited text layer ‘signif’ to shape file point data and then performing the same operation. However, the result is still the same.

Why is this happening and how do I remove this error?

I am using QGIS 3.6.0 Noosa 64bit (OSGeo4W Network Installer), on Windows 10 machine. Snapshots of the errors are also attached.


texts:
I'm trying to create a map with a time slider, basically duplicating this example: http://bl.ocks.org/javisantana/f9b073869fdd3b696f93. I must not be entering in my info correctly because when I copy the code and change it to my map, the time slider does not appear. In the script section of the code, I've changed the user from 'documentation' to 'phrusa' (my username), changed the column 'built' to 'date_stamp' (my column that specifies date), 'us_po_offices' to my table name, and I've put the link for my map's api, but I can't figure out what else I need to change to fit my table/data.
Do I need to change "Date" or "FullYear" or "Time" to have the slider interact with my data?

I've pasted the code I'm trying to use below:

<script>

  function changeLegend(start, end) {
      $('#legend').html(new Date(start).getFullYear() + " - " + new Date(end).getFullYear())
  }
  function addTimeSlider(sublayer) {
     var sql = cartodb.SQL({ user: 'phrusa' })
     // fetch time range
     sql.execute('select max(date_stamp), min(date_stamp) from tablename', function(data) {
       var range = data.rows[0];
       var max = new Date(range.max).getTime()
       var min = new Date(range.min).getTime()
       // update slider with range
       $("#slider").slider({
          range: true,
          min: min,
          max: max,
          values: [ min , (min + max)/2 ],
          change: function(event, ui) {
            // give feedback to the user on slide change
            changeLegend(ui.values[0], ui.values[1]);
          },
          stop: function( event, ui ) {
            // when user selects the dates, update the layer with the range
            var start = new Date(ui.values[0]).toISOString()
            var end = new Date(ui.values[1]).toISOString();

            // build sql
            sublayer.setSQL("select * from tablename where date_stamp >= '" + start + "' and date_stamp <= '" + end + "'");
          }
        });
        changeLegend(min, (min + max)/ 2);
     });
  }

  function main() {

    var map = new L.Map('map', {
      zoomControl: false,
      center: [35, 40],
      zoom: 3
    });

    L.tileLayer('http://{s}.api.cartocdn.com/base-light/{z}/{x}/{y}.png', {
      attribution: 'CartoDB · OSM data'
    }).addTo(map);

    cartodb.createLayer(map, 'http://phrusa.cartodb.com/api/v2/viz/...../viz.json', { legends: false })
     .addTo(map)
     .on('done', function(layer) {
       // add time slider on change
       addTimeSlider(layer.getSubLayer(0));
    }).on('error', function() {
      cartodb.log.log("some error occurred");
    });
  }

  // you could use $(window).load(main);
  window.onload = main;
</script>


texts: Exporting directions from Google to KML?

texts:
I used Google to draw directions from one location to another location, but can't find a way to export it as KML so I can import it into my Android smartphone.

Is there a way, or must directions be drawn in Google (Maps | MapsEngine)?


texts: How to Add image and hyperlink along with GetFeatureInfo Popup?

texts:
I want to cutomize the Geoserver GetFeatureInfo Popup. I want to add image against each parcel  I click. For this purpose I've read to change the content.ftl file of geoserver. But I found that file in wms.jar file. So, I extract it and after changing the content.ftl file, when i agin convert into .jar file, Geoserver doesn't work.
Is there a solution to it?
Any other suggestion to add the image along Featureinfo.


texts: How to find Entropy in land use (Voronoi polygons)?

texts: I have a shapefile of land use in Montreal. It is composed of polygons. There are seven land uses: Residential, Commercial, Industrial, Waterbody, Parks, Open Area and Government. I need to find Entropy from this layer. Can anyone suggest how to do that? I am using Arcmap 10.1


texts: Operations using JTS with CRS as WGS84 or similar geographic coordinate reference systems

texts:
I need to carry out the operations of the JTS, that is a API for processing geometry on the  2D dimensional plan. Operations as Intersection, Union, Buffer, Boundary... And the parameters of this operations would be in the same CRS, but this CRS could be WGS 84 or other similar geographic coordinate reference systems.
I could make this operations using JTS without transformations or I need use other tool to transform first to the Cartesian plane, and transform the result to the CRS of the source?


texts:
I am trying to make election map so I can draw chart and label with same chart color. So I have six election zone and running candidates may be up-to six or less. So if some election zone might have only two candidates so it should label only 2 not  the empty value. Labeling should be stacked and should match with the color code of the charts. Here is my code which worked if all the column has value. On candidate 4 which has only 4 won,t label it. All the document is attached.

def FindLabel ( [Sheet1$.Candiate1], [Sheet1$.Candiate2], [Sheet1$.Candiate3], [Sheet1$.Candiate4], [Sheet1$.Candiate5], [Sheet1$.Candiate6], [Sheet1$.VoteCount1],[Sheet1$.VoteCount2],[Sheet1$.VoteCount3],[Sheet1$.VoteCount4],[Sheet1$.VoteCount5],[Sheet1$.VoteCount6]   ):

if [Sheet1$.VoteCount1]> 0:
    [Sheet1$.Candiate1]= [Sheet1$.Candiate1]+" Vote Count= "+Sheet1$.VoteCount1]

if [Sheet1$.VoteCount2] > 0:
    [Sheet1$.Candiate2]= [Sheet1$.Candiate2]+" Vote Count= "+sheet1$.VoteCount2]

return "" +[Sheet1$.Candiate1]+"" +'\n'+""+[Sheet1$.Candiate2]+ ""


texts: Visualizing GPS points on interactive map

texts:
I am using version 3.0 of python, I want to visualize GPS points in my dataset on interactive map. Is there any python library that can help me to visualize my data on maps.

The input dataset in given in the link
https://drive.google.com/file/d/1JjvS7igTmrtLA4E5Rs5D6tsdAXqzpYqX/view

Three lines of the dataset:

lat,long,tep,ht,altitude,date,time,userid,date_time,Transportation_Mode,time_delta,trip_id,segmentid,Vincenty_distance,velocity,acceleration,jerk
39.97954667,116.3068133,0,0,40784.42587,8/29/2011,10:13:15,141,8/29/2011 10:13,walk,1,1,1,1.551676553,1.551676553,0.550163852,-1.017629555
39.97955833,116.3068233,0,0,40784.42588,8/29/2011,10:13:16,141,8/29/2011 10:13,walk,1,1,1,2.101840405,2.101840405,-0.467465704,0.732307066


texts:
I installed the new version of QGIS 2.16 and I want to import an image (SVG file) as symbol for my point, but the setting style supports only xml format in import. Before, with an older version of QGIS, importing SVG files was possible.

How can I do this in the new version?

if i check svg symbol, I don't have the windows svg groups and images, my screen shot it's different


texts: Select observations based on Quality Bands , Google Earth Engine

texts:
From the following NDVI collection, I want to keep only the observations with Good Quality (equal to zero).
 I try to do this with the SummaryQA but I get "Pattern 'SummaryQA' did not match any bands".
 My attempt:

var geometry = /* color: #0b4a8b */ee.Geometry.Polygon(
  [[[110.49279650854476, 39.617434795765476],
    [111.02014025854476, 36.57484614339866],
    [117.48010119604476, 36.92695433193529],
    [117.43615588354476, 40.35816406661026]]]),modis =ee.ImageCollection("MODIS/006/MOD13Q1")
    .filterBounds(geometry)
    .map(function(image){return image.clip(geometry)})
    .filter(ee.Filter.calendarRange(4,7,'month'));
    var modisNDVI = modis.select("NDVI");

    var maskQA = function(image) {
      return image.updateMask(image.select("SummaryQA").eq(0));
    };
    var best = modisNDVI.map(maskQA)
    print(best)


texts: Load a layer from SpatiaLite database, using Python

texts:
I'm new to SpatiaLite and PyQGIS.
I would like to upload a layer from SpatiaLite database. To do this, I tried with this code:

database = r"C:\sqlite\db.sqlite"
uri = QgsDataSourceUri()
uri.setDatabase(database)
schema = ''
table = 'table_name'
geom_column = 'point'
uri.setDataSource(schema, table, geom_column)
display_name = 'table_name'
vlayer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')
QgsProject.instance().addMapLayer(vlayer)

But this doesn't work, the program says "Unavailable layer! Layer datasource could not be found"
I think the problem may be the value of "geom_column".
The table has a column called "the_geo", where the row values are POINT(x_value, y_value).

I can't understand where I'm wrong, is there anyone who can kindly help me?

Thanks in advance.