texts: spatial index intersection with qgsrectangle only computes polygones

texts:
I have the problem, that i only get polygon features inside my defined
rectangle. If i switch off my polygon layer(s) i get my lines.

TOC: Line-Layer on, Polygon-Layer on: i only get Polygons which are inside
TOC: Line-Layer on, Polygon-Layer off: i get my lines which are inside.

In the tutorials, they say, that i will get all features which intersects my
defined rectangle.
If I use the index.nearestNeighbor function I will get poylgones and lines together. What did I wrong?

from qgis.gui import QgsMapTool
from qgis.core import QgsMapLayer,QgsPoint,QgsMapToPixel,QgsRectangle,QgsFeature, QgsFeatureRequest, QgsGeometry, QgsSpatialIndex
from PyQt4.QtGui import QCursor, QPixmap
from PyQt4.QtCore import Qt, QVariant
from qgis.utils import iface
import itertools

class NearestFeatureMapTool(QgsMapTool):

def __init__(self, canvas):
    super(QgsMapTool, self).__init__(canvas)
    self.canvas = canvas
    self.cursor = QCursor(Qt.CrossCursor)

def activate(self):
    self.canvas.setCursor(self.cursor)

def canvasReleaseEvent(self, mouseEvent):
    feat=QgsFeature()
    for layer in self.canvas.layers():
        if layer.type() != QgsMapLayer.VectorLayer:
            # Ignore this layer as it's not a vector
            continue
        if layer.featureCount() == 0:
            # There are no features - skip
            continue

        layer.removeSelection()

        # Determine the location of the click in real-world coords
        layerPoint = self.toLayerCoordinates( layer, mouseEvent.pos() )

        feature_dict = {feature.id(): feature for (feature) in layer.getFeatures()}

index = QgsSpatialIndex()
map(index.insertFeature,feature_dict.values())

#define rectangles
displayrectangle = self.canvas.extent()

#print displayrectangle.toString()
ownrectangle = QgsRectangle(layerPoint.x(),layerPoint.y(),layerPoint.x()+200,layerPoint.y()+200)

#print ownrectangle.toString()
if displayrectangle.contains(ownrectangle):
    userectangle = ownrectangle
else:
    userectangle = displayrectangle

# Auswahl der Features innerhalb des Rechteckes
for i in index.intersects(userectangle):
    print feature_dict[i].attributes()

Update: Unfortunately the code structure made  problems and the PyQgis didn't care about this.

# -*- coding: utf-8 -*-

from qgis.gui import QgsMapTool
from qgis.core import QgsMapLayer,QgsPoint,QgsMapToPixel,QgsRectangle,QgsFeature, QgsFeatureRequest, QgsGeometry, QgsSpatialIndex
from PyQt4.QtGui import QCursor, QPixmap
from PyQt4.QtCore import Qt, QVariant
from qgis.utils import iface
import itertools

class NearestFeatureMapTool(QgsMapTool):
def __init__(self, canvas):

    super(QgsMapTool, self).__init__(canvas)
    self.canvas = canvas
    self.cursor = QCursor(Qt.CrossCursor)

def activate(self):
    self.canvas.setCursor(self.cursor)

def canvasReleaseEvent(self, mouseEvent):
     layerData = []
     feat=QgsFeature()
     index = QgsSpatialIndex()
     for layer in self.canvas.layers():
         if layer.type() != QgsMapLayer.VectorLayer:
            # Ignore this layer as it's not a vector
             continue

         if layer.featureCount() == 0:
    # There are no features - skip
             continue

         layer.removeSelection()

        # Determine the location of the click in real-world coords
         layerPoint = self.toLayerCoordinates( layer, mouseEvent.pos() )

         feature_dict = {feature.id(): feature for (feature) in layer.getFeatures()}

         map(index.insertFeature,feature_dict.values())

    #define rectangles
         displayrectangle = self.canvas.extent()

     #print displayrectangle.toString()
         ownrectangle = QgsRectangle(layerPoint.x(),layerPoint.y(),layerPoint.x()+200,layerPoint.y()+200)

     #print ownrectangle.toString()
         if displayrectangle.contains(ownrectangle):
             userectangle = ownrectangle
         else:
             userectangle = displayrectangle

     # Print features inside my rectangle
         for i in index.intersects(userectangle):
             print feature_dict[i].attributes()


texts: how to index measurements with spatial coordinates

texts:
I am beginner with the RTree I use libspatialindex, I will want to know how I insert a polygon / point with a measurement (statistic or measurement of weather)
do I have to change the structure ?

from rtree import index

idx = index.Index()
idx.insert(0, (0, 0, 1,1 ),my mesure)

yet I see in the class Node the attribute byte** m_pData; (it is he who contains the information of each leaf)


texts: Leaflet: popup showing information on the polygon

texts:
I have a Leaflet map. To the map I load two JavaScript files (A.js & B.js) which both consist of only polygons. Now my aim is, that when I click inside Leaflet on one of the polygons, that I get a popup showing the area of that polygon. The area is saved in the JavaScript file under Shape_Area.

"type": "Feature",
"properties": { "OBJECTID": 1, "Id": 1, "gridcode": 1, "Shape_Leng": 1805.0000000099999, "Shape_Area": 8781.24999991, "LandUse": "Urban" }

My code until now is:

function onEachFeature(feature, layer) {
                if (feature.properties && feature.properties.NAME) {
                        layer.bindPopup(feature.properties.NAME);
                    }
            }
var LayerA = L.geoJSON(UF,{style: style},{onEachFeature: function(feature, layer) {layer.bindPopup(feature.properties.type);}}).addTo(mymap);;
LayerA.eachLayer(function (layer) {layer.bindPopup(UF.Shape_Area);});

I get already my popup, but I do not know how to exactly get the Shape_Area of my JavaScript-file.

Anybody has a hint for that please?


texts: Using Trimble Geoexplorer with EM31 as external sensor?

texts:
I am using a Trimble Geoexplorer 6000 in conjunction with a Geonics EM-31. The EM-31 communicates conductivity data to the Trimble as an external sensor. It provides conductivity measurements as the Trimble collects location information.

I think that I've successfully connected the Trimble to the EM-31 via bluetooth. When I open the "Sensor" menu under the "Status" menu, it shows that data is being collected. However, when I download the data through pathfinder, I only receive the positions. There's no conductivity data associated with the points or line features that were collected.

Is there a specific property that needs to be altered to receive the external sensor data?

I found a similar question without an answer: how-do-i-export-sensor-records-from-a-trimble-geoexplorer-to-a-mapping-program


texts:
I'm looking for ways to customize built in widgets, but couldn't find anything in version 4.x.

In version 3.x you can do it in following way as mentioned in doc https://developers.arcgis.com/javascript/3/jshelp/localization.html:

var map;
require(["esri/map", "dojo/i18n!esri/nls/jsapi", "dojo/domReady!"],
  function(Map, bundle) {
    bundle.widgets.homeButton.home.title = 'Эзэмшлийн газар'
  }
);

But in version 4.x, there's only way to customize DOJO text, not arcgis:
https://developers.arcgis.com/javascript/latest/guide/localization/index.html


texts: How to display the map.InfoWindow outside of the map container?

texts:
I need to show an infoWindow outside the map (i.e. have it owned by the body of the HTML document, not the map div).
If i subclass InfoWindowBase I can do this myself but I really don't want to do that. Is there any way of doing this?

I have tried the following without success:-

     var infoWindow = new esri.dijit.InfoWindow({
      domNode: dojo.create("div", null, dojo.body())
    });
map.InfoWindow = infoWindow;

and
    map.infoWindow.place("", dojo.body());

thanks,

Phil


texts: Can we compare a Service in Arc GIS server with normal webservice

texts: In Arc GIS server there are many types of services are available (like feature service, Map Service, Geo-processing service etc.,). My question is in general web service whether it is ASP.NET web service or WCF service are like same or different with respect to map services in Arc GIS server.


texts: Getting values from adjacent polygons or nearest point

texts:
I have joined temperature values from the point layer to intersecting polygons but I have few polygons (yellow in the image) that don't contain any point data.
Can I somehow get average temperature values from adjacent polygons or get temperature value from the nearest data point? Now the temperature values in these polygons have NULL values because of lack of data points.


texts: AcrMap- How to use the dissolve tool with large amounts of polygons?

texts:
I am a college student working on an independent research project and I am trying to figure out how many points fall within a polygon.

I am examining the impacts of zoning in the Great Barrier Reef and I am currently focusing on coral bleaching occurrences.
I have already created a map in ArcMap with a layer for recorded coral bleaching (symbolized by levels of severity) and a layer file of the different zones (general use, conservation park, etc).

We already created a count field and a spatial join between the point shapefile (coral bleaching) and the polygon shapefile (zones). The problem we are having is that when we try to group or summarize the points to the zones the program is recognizing each area (over 100) as its own zone instead of grouping the zones together. For example if there are 45 specific areas that are classified within the general use zone the computer counts all 45 areas instead of lumping them together under general use. Right now it is stating that each zone only has 1 value within it.

Question- What tools in ArcMap would be able to group the polygons together based on a common field (zone_) even with 1000 different values.

The scenario is: Given a layer bleaching occurs and a layer of zone boundaries by zone-use type, calculate/summarize the number of bleaching occurrences in within each zone.


texts: pgrouting shortest path result in wrong order

texts:
I am using the Dijkstra shortest path algorithm to plot a route between two points using the latest version of pgrouting in postgresql 9.1 with postGIS 2.0.1.

I am able to get a correct result however the edges are in the wrong order so when I plot it in my app it does not draw correctly.

I have read several other posts with people having the same problem but none of the solutions provided have worked for me.

Here was the original query I was using:

SELET rt.gid, ST_AsGeoJson(rt.the_geom) AS coordinates, hh_2po_4pgr.gid
FROM hh_2po_4pgr, (SELECT gid, the_geom
               FROM dijkstra_sp_delta('hh_2po_4pgr', 336, 557, 0.1)) as rt
WHERE hh_2po_4pgr.gid=rt.gid;

I then read a post which suggested doing it this way to order by the id:

select rank, gid, ST_AsGeoJson(the_geom) FROM (
             select edge_id, ROW_NUMBER() OVER (PARTITION BY 1) as rank
             from shortest_path(
                 'SELECT gid as id, source::integer, target::integer, length::double precision as cost
                  FROM hh_2po_4pgr', 557, 350, false, false)
             )
as route join hh_2po_4pgr h on h.gid = edge_id order by rank;

This worked for one of my route finding queries but I have tried it with others and it doesn't work.

I read this: ftp://ftp.remotesensing.org/pgrouting/forum/pgrouting.postlbs.org/ticket/156.html

which says the wrappers return an id which can be used to preserve order so I then created this query:

SELECT id, gid, ST_AsGeoJson(the_geom) FROM dijkstra_sp('hh_2po_4pgr', 557, 350)

but the ordering is still not correct.

Is there any way to fix this?

Thanks


texts: How can I use measure tool with Openlayers 3?

texts:
I didn't see any mention for measure tool in OL 3 API.
Is there an alternative / a way to use a measure tool (line & area) with OL 3 ?


texts: Finding geometry anomalies using ArcGIS Pro

texts: Is there a simple way (could be hard as well)  to find geometry (polygons) anomalies. For example polygon should be nice rectangle, but doesn't always look so (like in example). I am talking about more than 100 000k objects and somehow I need to find all anomalies.


texts: Showing more than one style in GeoServer layer?

texts:
I am able to only display point type on the map when it is opened a layer in GeoServer layer preview. I want to display other shapes such as LineString.

How can I add this to Style Editor to show different shapes on one layer?

Below is my style script which I have copied from Style Editor in GeoServer.

Is it possible to show more than one shapes style in a layer?

If I ask in a different way, is  Symbolizer able to added in the script of XML to display when a layer is opened by OpenLayers in GeoServer?

<?xml version="1.0" encoding="UTF-8"?>
    <StyledLayerDescriptor version="1.0.0" xmlns="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc"
                           xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
                           xsi:schemaLocation="http://www.opengis.net/sld http://schemas.opengis.net/sld/1.0.0/StyledLayerDescriptor.xsd">
        <NamedLayer>
            <Name>M_DIZAYN_VANA</Name>
            <UserStyle>
                <Name>SERVIS_NOKTASI</Name>
                <Title>SERVIS_NOKTASI</Title>
                <Abstract>SERVIS_NOKTASI</Abstract>
                <FeatureTypeStyle>
                    <Rule>

                        <PointSymbolizer>
                            <Graphic>
                                <Mark>
                                    <WellKnownName>square</WellKnownName>
                                    <Fill>
                                        <CssParameter name="fill">#a43eb4</CssParameter>
                                        <CssParameter name="stroke">#000000</CssParameter>
                                        <CssParameter name="stroke-width">1</CssParameter>
                                    </Fill>
                                </Mark>
                                <Size>7</Size>
                            </Graphic>
                        </PointSymbolizer>
                    </Rule>
                </FeatureTypeStyle>
            </UserStyle>
        </NamedLayer>
    </StyledLayerDescriptor>


texts: Long Field Name in ArcGIS Online Hosted Table Strips Data

texts:
I've encountered an issue when uploading a CSV to ArcGIS Online as a hosted table where if the field name exceeds X number of characters (I have not discovered this number yet), all data in that field is stripped out of the table.

I have managed to work around this by shortening field names before the upload process, but I would like to understand this issue better. Is this a known limitation of ArcGIS Online, or is this a new bug?


texts:
I have got follow code:

for row in arcpy.da.SearchCursor(db_from, fields, where_clause="objectid < 11111"):
    value_from = list(row)
    for rowto in arcpy.da.SearchCursor(db_to, fields, where_clause="uid = '" + str(row[0]) + "'"):
        value_to = list(rowto)
        if(set(value_from) == set(value_to)): # data is same
            same +=1
            print "data is same"
        else:
            with arcpy.da.UpdateCursor(db_to, fields, where_clause="uid = '" + str(row[0]) + "'") as cursor:
                for x in cursor:
                    cursor.updateRow(value_from)
                    updated +=1
                    print "data update"
    else:
        to_cursor.insertRow(value_from)
        print "Inserted"

print "Same Data: ", str(same)
print "Updated Data: ", str(updated)

I am getting next error:

data update
data update
data update
data update
data update
Inserted
Traceback (most recent call last):
  File "C:\code\Python\work\dbsync.py", line 69, in <module>
    for x in cursor:
RuntimeError: workspace already in transaction mode

Where I should to close transaction? Am I right understand that issue is rise when I do insert when transaction (cursor?) of update is still active?

I am using ArсPy


texts:
I have a problem with projection.

When I try to georeference an image in QGIS, I set up CRS WGS 84 Pseudo Mercator (EPSG: 3857).

That info appears in georeference settings and it is also my default projection system.

When I start georeferencing, the final raster comes with a USER 100004 CRS.

I tried to change it but when I open tiff file again in properties is still USER CRS.


texts: wcs Time Dimension Always Returning Last Available Date

texts:
I currently have a geoserver instance with a working wms time dimension layer. I am able to ask for the raster geotiff for a specific day (2016-1-14) with the following request and geoserver responds correctly.

http://myserver/geoserver/workspace/wms?service=WMS&request=GetMap&layers=layername&time=2016-1-14&bbox=-125.020833333333,24.0625,-66.479166666662,49.937500000002&width=1600&height=800&srs=EPSG:4269&format=image/geotiff

Now I want to get the same functionality from wcs but geoserver always returns the raster geotiff for the last available date which is currently the 24th instead of the date that I pass in via the time param which is the 14th. My wcs request looks like this (again I want the 14th, but geoserver responds with the 24th):

http://myserver/geoserver/gdd/wcs?service=WCS&version=2.0.1&request=GetCoverage&CoverageId=mycoverage&time=2016-1-14&format=image/geotiff

If I have the coverage described via

http://myserver/geoserver/ows?service=WCS&version=2.0.1&request=DescribeCoverage&coverageId=mycoverage

I do see the different days listed:

<gmlcov:metadata>
<gmlcov:Extension>
<wcsgs:TimeDomain default="2016-01-24T00:00:00.000Z">
<gml:TimeInstant gml:id="gdd__agdd_anomaly_td_0">
<gml:timePosition>2015-01-01T00:00:00.000Z</gml:timePosition>
</gml:TimeInstant>
<gml:TimeInstant gml:id="gdd__agdd_anomaly_td_1">
<gml:timePosition>2015-01-02T00:00:00.000Z</gml:timePosition>
</gml:TimeInstant>
<gml:TimeInstant gml:id="gdd__agdd_anomaly_td_2">
<gml:timePosition>2015-01-03T00:00:00.000Z</gml:timePosition>
</gml:TimeInstant>
<gml:TimeInstant gml:id="gdd__agdd_anomaly_td_3">
<gml:timePosition>2015-01-04T00:00:00.000Z</gml:timePosition>
</gml:TimeInstant>
<gml:TimeInstant gml:id="gdd__agdd_anomaly_td_4">
<gml:timePosition>2015-01-05T00:00:00.000Z</gml:timePosition>
</gml:TimeInstant>
<gml:TimeInstant gml:id="gdd__agdd_anomaly_td_5">
<gml:timePosition>2015-01-06T00:00:00.000Z</gml:timePosition>
</gml:TimeInstant>
<gml:TimeInstant gml:id="gdd__agdd_anomaly_td_6">
<gml:timePosition>2015-01-07T00:00:00.000Z</gml:timePosition>
</gml:TimeInstant> ... up to 2016-01-24T00:00:00.000Z

I've also tried setting &TIME=2016-01-14T00:00:00.000Z  but still no dice.


texts: Intersecting multiple classified rasters

texts:
I have 3 raster files:

A transform aspect (TRASP) layer (classified into 2 strata-(0-0.5, 0.5-1))
A canopy cover layer (classified into quantiles-5 breaks)
A canopy height layer(classified into quantiles-5 breaks)

The purpose is to take these three layers and do 2x5x5 to get 50 strata,
How can I do this in the raster calculator or any other tool in ArcGIS Desktop?
I tried using the simple "TRASP"*"canopy_cover"*"canopy_height" in the raster calculator but I got a big maximum value, so I don't know if it's correct or not. I also tried using Con(((IsNull(TRASP))|(IsNull(canopy_cover))|(IsNull(canopy_height))),0,1) suggested in another question, but that returned an error.


texts: How to get an ESRI spatial reference in XML for Pathfinder Export?

texts:
I am trying to export from Trimble Pathfinder into an ESRI Geodatabase, in UTM10N (EPSG:26910), however I do not have a spatial reference XML setup in the export preferences.  I looked into how to get this and found instructions that said I could use ArcCatalog 10.0 on a standard license to do it.  I have access to a standard license but it is 10.3.  I tried using the spatial reference exported from that but Pathfinder crashes on export when using it.

I have looked online and tried some XML variants posted at epsg.io and other websites, but nothing is working, I suspect because the ESRI XML variant to be used in a file geodatabase is somewhat customized (?).

Does anyone have leads on how to get this XML snippet so I can do this export to file geodatabase?


texts: Clearing data frame projection of current mxd using ArcPy?

texts:
I am required to open a bunch (50 per week) of mxd's and 'Clear' the Coordinate system from the data-frame to make sure my students have all of their data projected in the same (correct) Coordinate System, to see if their data 'lines up'.

I do this by hand, over and over and over, by using the steps below:

-Right Clicking the Data Frame

-Go to the Coordinate System Tab of 'Data Frame Properties'

-Click on the Globe dropdown arrow

-Scroll to 'Clear'

-Click 'OK' or 'Apply' as I watch to see if any of the data disappears.

There has GOT to be a way I can copy/paste a short line of Python into the Python Window to automate this. I know I could write a complete script that would check the coordinate system of each layer in the map without actually opening the map, but I'm just looking for a quick fix that I can do while I also evaluate their cartography.

I am using ArcGIS 10.2.1 for Desktop with an Advanced level license.


texts: SQL Query for percent of buffer filled with objects

texts:
I am trying to write a query to find the percentage of a buffer polygon (Called 250_Buffer) that is filled with or contains other polygons in another layer (Called Development) any help that you guys could give me would be awesome.

I am hoping for an answer output as a percentage of 100 but I could also get something in square meters/acres/miles and calculate the percentage myself.


texts: How to add a table from ArcGIS Online to ArcMap

texts:
When you're signed into ArcGIS Online within an ArcMap session, there is an option Add Data From ArcGIS Online:

I have a service in ArcGIS Online which consists of 2 feature layers and 1 table:

When I add this to ArcMap using the above option, the table isn't included. Is this a bug, or is there a workaround?

Thanks


texts: Improving search results in qgis2web delimited text layer

texts:
I'm using QGIS 3.10.11 to generate a webmap using qgis2web. In the 'Export to webmap' (in export as OpenLayers) settings window, I am adding 'layer search' and selecting a csv text field for the search. The export runs OK.
When I open the map in a browser (MS Edge) and type a search term, the list of results returned includes many unwanted hits. This is because the search is finding any record which contains the letters in the search term. so a search for 'box' will return 'Box farm' but also returns 'Roxby farm Cottage'. I have tried entering the term in single and double quotes, but this returns  no hits. Clearly, the more text that is entered the better the results. However, some users will only enter 4 or 5 letters in the search box which causes the problem. I have about 2000 rows (records) in the csv file.
How can I force only a specific string of text to be returned in search results?


texts: Error message: Invalid geometry Refactor Fields QGIS

texts:
I'm trying to do a simple modification of changing field type to vector table by using the tool Refactor Fields. However, the output log returns the following error message, every time:

Traceback (most recent call last):
File "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\algs\qgis\FieldsMapper.py", line 132, in processAlgorithm
return super().processAlgorithm(parameters, context, feeback)
_core.QgsProcessingException: Feature (18) has invalid geometry. Please fix the geometry or change the Processing setting to the "Ignore invalid input features" option.

Execution failed after 0.17 seconds

It would be OK for me to ignore invalid input features, but there is no such option in the processing settings of this tool (see image). I'm using QGIS 3.6.1. Any suggestions?


texts: QGIS: Add attribute from polyline layer to point layer

texts:
I want to add an attribute from a vector line layer to a vector point layer.

The vectorline is railwaynetwork. It has been divided into segments automatically. To process the data further I had to convert the lines into points. I wanted defined distances in between so I chose to use the plugin QChaninage. Now I got points in my defined distance. But in the attribute table there aren't any of the original attributes that where in the vectorlines table. How can I connect those attributes with my points?

I already tried to use tools from the toolbox (e.g Join attributes by location tool) but until now I could not solve it.


texts: Displaying feature points from an imported csv file (GEE)

texts:
I have imported the following CSV file into GEE (https://docs.google.com/spreadsheets/d/14XYsKQYylJ77uhNVhiLqUYwyGgBhE0N1h4F5QTn2Cvk/edit#gid=832877374).
However, even though I can print the data, I wasn't able to display it using the usual procedure. (No error Messages)

Asset: https://code.earthengine.google.com/?asset=users/farhatlokmen/ecoModeling/Falco

Code: (https://code.earthengine.google.com/?scriptPath=users%2Ffarhatlokmen%2Ftrials%3Atest)

var falcoCherrug = ee.FeatureCollection("users/farhatlokmen/ecoModeling/Falco");

print(falcoCherrug);
Map.addLayer(falcoCherrug.draw({color: 'ff0000', pointRadius: 5}), {}, 'Falco Cherrug');


texts: Rasterize a shapefile generate a black TIFF

texts:
I run the following codes in jupyter notebook to rasterize the shapefile: building.shp, but it generates a black TIFF, which in QGIS shows only 0 and nan values.

Here are the source codes

Could anyone help me to find what the problem is?

import gdal
from osgeo import osr
from osgeo import ogr

raster_path = 'C:\\exercise_data\\training_manual_data\\epsg4326\\buildings.tif'
shapefile = 'C:\\exercise_data\\training_manual_data\\epsg4326\\buildings.shp'

# 1) opening the shapefile
source_ds = ogr.Open(shapefile)
source_layer = source_ds.GetLayer()

# 2) Creating the destination raster data source

pixelWidth = pixelHeight = 0.001 # depending how fine you want your raster
x_min, x_max, y_min, y_max = source_layer.GetExtent()
cols = int((x_max - x_min) / pixelHeight)
rows = int((y_max - y_min) / pixelWidth)

target_ds = gdal.GetDriverByName('GTiff').Create(raster_path, cols, rows, 1, gdal.GDT_Byte)
target_ds.SetGeoTransform((x_min, pixelWidth, 0, y_min, 0, pixelHeight))
band = target_ds.GetRasterBand(1)
band.FlushCache()

# 4) Instead of setting a general burn_value, use optionsand set it to the attribute that contains the relevant unique value ["ATTRIBUTE=ID"]
gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[1])

# 5) Adding a spatial reference
target_dsSRS = osr.SpatialReference()
target_dsSRS.ImportFromEPSG(4326)
target_ds.SetProjection(target_dsSRS.ExportToWkt())


texts: Creating drive time polygons using Open Source tools?

texts:
Is anyone using open source tools and/or data to generate drive time polygons that radiate out from a given location?  For example, I have used ESRI products like Network Analyst and Business Analyst to generate polygon boundaries that represent the maximum distance one can travel along a road network from a given point within a user defined period of time (say 5, 10, and/or 15 minutes).  Are there any similar open source or web-based alternatives?

I saw this post, but it really doesn't give alternative options:
Can you make a travel region polygon with with Google Maps API?

I was particularly curious to know if anyone is (or can) use OpenStreetMap to do drive time analysis?

Also, is there any implementation of the pgRouting utilities that could be used to generate polygons?
http://pgrouting.postlbs.org/


texts: Calculating length (in meters) of GPS track in ArcGIS Desktop?

texts:
I'm trying to calculate the length of a GPS track in ArcGIS 10.3.1. What I've done so far is:

Converted track using the "Points to Line" tool
Opened attribute table
Made new "Length" column
Opened "Field Calculator"
Entered in Python mode, "!Shape!.getLength("GEODESIC","METERS")"

Will this get me the length of the line in meters? It seems like it does but I wanted to double check.


texts: Wrong boundary when plotting Canada's shapefile in Python

texts:
I want to plot the Canadian shapefile with province boundaries highlighted and hollow filling.

I have downloaded the shapefile from here.

This is the script I use to visualize the shapefile:

import shapefile
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
#   -- input --
sf = shapefile.Reader("myfolder\gpr_000a11a_e.shp")
recs    = sf.records()
shapes  = sf.shapes()
Nshp    = len(shapes)
cns     = []
for nshp in xrange(Nshp):
    cns.append(recs[nshp][1])
cns = array(cns)
cm    = get_cmap('Dark2')
cccol = cm(1.*arange(Nshp)/Nshp)
#   -- plot --
fig     = plt.figure()
ax      = fig.add_subplot(111)
for nshp in xrange(Nshp):
    ptchs   = []
    pts     = array(shapes[nshp].points)
    prt     = shapes[nshp].parts
    par     = list(prt) + [pts.shape[0]]
    for pij in xrange(len(prt)):
     ptchs.append(Polygon(pts[par[pij]:par[pij+1]]))
    ax.add_collection(PatchCollection(ptchs,facecolor=cccol[nshp,:],edgecolor='k', linewidths=.1))
ax.set_xlim(-160,-40)
ax.set_ylim(40,90)

The output is confusing, in that the Northern border is completely wrong and it does not resemble the real one. I think there are unwanted lines that connect points far away from each other. As a result, the islands in Northern Nunavut have been merged in a single polygon. How to obtain a neat border?


texts: Modify polygons so that they don't overlap and area stays the same

texts:
I have a set of polygons and they can overlap with each other, like this:

I want to modify them in such a way that they don't overlap and the resulting surface area stays the same. Something like this:

It is okay if the shape or the position changes. The main thing is that they should not overlap with each other and the area should not change much (I know the area changed a little in the second image but I drew it manually thus let's just assume that the areas did not change).

I am trying to do it programmatically with the help of Python. Basically I stored polygons in a PostGIS database and with the help of a script I want to retrieve them and modify them.

I am very new to GIS and thus this seems like a difficult task.

What is the correct way of doing it? Is there an algorithm that solves this kind of problems?


texts:
I want to run this script, but it always gives me this error message:

Traceback (most recent call last):
File "C:\PROGRA~1\QGIS21~1.1\bin\pct2rgb.py", line 121, in
ct_size = ct.GetCount()
AttributeError: 'NoneType' object has no attribute 'GetCount'

I tried it with different images - all the same problem...

Can anybody explain what tcp means?


texts:
I am wondering if it is possible to set a custom transformation in Qfield. We are doing surveing in austria and unfortunately the most common Coordinate-System is not ETRS89 based.
In QGIS it is possible to copy a NTV2-Grid into the \proj folder and to select this as the standard transformation. Is this also possible in Qfield?
I have seen that with the latest version it is possible to select a grid for vertical shift. It would be great if this is also possible for horizontal coordinates.


texts: Autonumbering using ArcGIS Web AppBuilder?

texts:
I work for a water utility division and we are using ArcMap 10.1 with attribute assistant. I used the AutoNumber method to auto increment and populate a field with the next highest integer. Unfortunately this method does not seem to work in ESRI Web AppBuilder.
Anytime I edit my layer and create new fields using the Editor widget provided by Web AppBuilder, the value of the AutoNumber field is set as 0. However, it works perfectly through the Editor session in ArcMap.
Is there a way to make it work through the Web AppBuilder ? Or is there another option to auto populate a field when using Web AppBuilder and considering daily updates ?


texts: Libraries and/or ideas for navigating WMS time-parameter

texts:
I'm currently working on a project where we need to incorporate the TIME-dimension of WMS-layers (see http://docs.geoserver.org/latest/en/user/services/wms/time.html#wms-time) in our web map.

The technicalities of the time-dimension is rather straight-forward, but in terms of creating a good user-interface things are more difficult. What I'm wondering is if there exists some (open-source) javascript library for handling a "timeslider" of sorts, where the user can skip through time and update the map, preferrably also with a "play"-button to play back an animation.

My searches haven't turned up anything, so I don't expect such a library to exist. What I do suspect exists are several implementations of similar functionality. If you know of any such implementations (be it web, desktop or mobile) I'm interested in hearing about them, to get ideas to how this can be solved.

So, specifically, what I'm asking is:

Do you know of any open source openlayers-components/plugins/examples of advanced time-navigation?
Do you have any input on how the GUI for such a component could look like? Screenshots, links to programs/apps or websites.

(A bit more background: I'm working on a web map that uses OpenLayers 2.x, and we currently have 20-something map layers that supports the TIME-dimension. These ranges from 48 hours into the future to 100 years back in time (with different resolutions of course). The available timestamps for each layer can be parsed from the GetCapabilities-document. This means that I need to create (or find) a rather generic component with no hard-coded dates or layers.

I do realize that creating a timeslider-component that handles both the year-by-year and hour-by-hour aspect will be difficult. Some compromise could be done here.

Technically I'm thinking of basing a solution on OpenLayers, moment.js, Backbone.js and bootstrap, but I'm open to suggestions here as well)


texts: How big does a circular geofence need to be for the curvature of the earth to matter?

texts:
I'm working on adding geofences to a mobile application. I know the proper way of calculating distance between two points on a sphere is using the great-circle or Haversine formula. I don't know if using those formulas are needed for relatively small scale applications. This leads to my question.
At what point does the distance given by a great-circle calculation diverge significantly from a planar calculation?

For the purposes of this discussion assume greater than 5% error is significant.

The geofence is simply a center point and a radius. If the point is within the radius, by distance, it is considered within the boundaries of the fence.

I am not sure if it matters, but I'm just worried about circular geofences at this time with a user defined radius between 50 and 200 meters.


texts: How to load Landsat imagery over ROI defined as a polygon rather than a point? GEE

texts:
In google earth explorer you can load Landsat image around a roi. You can then conduct analysis on this image. I want to expand my region to cover a whole country or continent rather than just the small box defined by the roi. I've tried changing my roi to a polygon instead but this still only returns one image.

// filtering the image collection
var image = ee.Image(ee.ImageCollection('LANDSAT/LC8_L1T_TOA')
  .filterBounds(roi)
  .filterDate('2015-01-01', '2015-12-31')
  .sort('CLOUD_COVER')
  .first());
Map.addLayer(image, {bands: ['B4', 'B3', 'B2'], max: 0.3}, 'landsat');

You can see in the image I have my region of interest on the coastline of Bangladesh. I need to expand the Landsat image to cover a larger area.


texts: Can a custom javascript web application be uploaded and run on ArcGIS online?

texts:
I'm researching arcgis online and have a subscription account. From what I can make out from the documentation and some experimentation, I've found the following:

You can create a feature service or tiled map service and have ArcGIS online host it (but not a dynamic map service). You can reference from anywhere using the REST service.
You can create a "web map" that gives you a viewer and simple widgets, scale bar, zoom, etc. Using the ID of the web map, you can reference anywhere either arcgis online or your own site.
You can create a "web application" from one of a couple dozen template that gives you a few more options for customization such as popups and even editing widgets. You can run it on ArcGIS online. You can also download the template, customize if desired and run on your server.
You can upload a reference to a "web application", basically a bit of publicity and add it to your to your gallery. However, the web application only runs on your own webserver.

Question: Can you upload your custom web mapping application to ArcGIS online and have them host it? I'm understanding that this is not possible. Can someone confirm or deny? Thanks!


texts: QGIS 2.4.0-Chugiak: How to produce raster colour map

texts:
By carefully deriving the NDVI values from Landsat Band 3 and Band 4 after doing atmospheric correction by applying DOS1 method in semi-automatic classification Plugin QGIS, I was to calculate the Green Vegetation Fraction by using equation as I mentioned in my previous post post1 (also mentioned below, here), in which @WhiteboxDev. nicely explained the problem that might would occur. I had put the wrong NDVImin and NDVImax range and this was happened due to not seeing in the metadata rather by seeing visually in QGIS. In my second post post2, @SaultDon. pointed me the that mistake that actually, I was doing.

GVF = (NDVI - NDVImin) / (NDVImax - NDVImin)

After using the minimum and maximum values shown in the metadata, I recalculated the Green Vegetation Fraction by using the same equation that i have mentioned above, now the ranges starts form 0.12-0.90 in metadata, not visually.

My question is that what is the difference between the values ranges that have shown in visually and metadata.
And also I want to make the color map that shows the ranges between 0-0.20, 0.20-0.40, 0.40-0.60, 0.60-0.80, 0.80-1.0 in a index form.
I have done this in QGIS by clicking the layer properties in style menu and by selecting the singleband psedocolor in renter type and to generate new colour map I have selected continuous mode and 0 as minimum and 1.0 as maximum values. I have changed these values manually, did not take those values that already shown in it (the same values that visually displayed and these are not the actual values). Am I doing right or how can I do it otherwise to make the such figures that Hu and Jia (2010) made in figure 4 in attached paper, I have also attached band 3 and band 4 data of landsat paper, band 3 and band 4.
Briefly explain about how to make the colour maps in style menu of layer property dialogue box by using mode "continuous/equal interval" with minimum and maximum values. If One would manually change these values then what happened with the colour map that will display.


texts: QGIS calculate radius of an arc to attributes

texts:
I have QGIS 2.18.16 installed.
My problem is, how to calculate radius of an arc / curve feature (line)? The lines are originated from .dgn file, which I read with FME and wrote to PostGis DB. I have created primary keys etc, and the table is fully editable. I calculated lengths for the arcs successfully, but cant figure out how to calculate radius for these curves. This ( https://www.mathopenref.com/arcradius.html ) might give some idea on the math side, though I can't figure out how to use it in QGIS calculator. So I need to calculate the "R" for all the lines i have in the DB.

Bellow is an example. My data has Arcs / Curves on different table than "normal" polylines.


texts:
I am newbie in postgis.

I have a raster landsat images in table "ras" and a polygon shapefile image in table "shp".

Now I want to clip raster with polygon using postgis, but I have a problem

The raster and polygon have the same SRID(wgs84), but when using st_intersects(), the result table are all FALSE

Can anybody explain that?


texts: TMS or WMS source(s) for DEM data tiles?

texts:
I'm spinning my wheels trying to find an existing, free source for TMS or WMS elevation data that doesn't involve a long-winded clickabout or an emailed url. This question has a similar purpose, but the answers offered aren't compatible with what I want to accomplish.

My goal isn't to show tiles in a map, rather I want to automatically calculate some details about a user-selected drainage corridor (stretches people want to canoe/float, to be specific). For example, once the user sets an area (+/- 10 miles ^2) and choses their endpoints, I'd like to derive information like drop in elevation between endpoints, overall slope of the run, and graph of the longitudinal profile for the reach would be a nice start. ..my thinking is, maybe I can implement difficulty metrics.

Ideal data traits:

30-meter pixels,
png or tiff, tiles to mosaic or a custom BBox area is fine
source offers full coverage of the US (at least the 48 contiguous)
actual elevation values (i.e. I'm not looking for hillshades or textures.)

Example bounding area:

bbox=-93.424564,36.993675,-93.331009,37.056170

Alternatively, any approach conducive to auto-download of DEM tiles or clips based on a fixed grid or path/row combination would work, so long as there is either a preexisting grid to work with, or enough theory to pre-generate a grid. Of course, key to this is access to a repository of data that allows direct url or ftp access.

I'm very curious to see what people think.

http://geobrain.laits.gmu.edu/wcs4dem.htm It seems the makers of DEM Explorer (above) provide a WCS gateway that serves geotiff DEM clips for BBox regions. Thanks to @BradHards for finding this. Here's an example call to the WCS service against the example bounding region I noted above:

http://geobrain.laits.gmu.edu/cgi-bin/gbwcs-dem?service=wcs&version=1.0.0&request=getcoverage&coverage=SRTM_30m_USA&bbox=-93.424564,36.993675,-93.331009,37.056170&crs=epsg:4326&format=image/geotiff&store=true


texts: Openlayers 4 multiple geometries with different styles

texts:
Is it possible to set feature geometry to multiple polygons and style them separately?

Something like this ... imagine polygons is array with 4 polygons

var PolygonFeature = new ol.Feature({
    geometry: new ol.geom.Polygon(polygons)
});

or

new ol.geom.GeometryCollection([
    new ol.geom.Polygon(polygons[0]),
    new ol.geom.Polygon(polygons[1])
]);

and then

PolygonFeature.getGeometry()[1].setStyle(customStyle);


texts: network coverage for mobile phones - europe

texts:
I am looking for GIS data with the network coverage (mobile phones). Something like the layer shown in the cell mapper website where I can select different countries and all the providers to export a GIS layer with the towers heatmap would be great.

Although a nice raster with coverage given by all the providers would be perfect, I am guessing that maybe something like this is not public, I can also manage to work with presence/absence of network coverage by administrative regions in Europe.

Any help is appreciated, it is the first time I am trying to find & use this type of data, so I might have failed to find any obvious websites, or what I am looking for is not possible to obtain.


texts: Can mosaic dataset be used as input to ArcGIS Pro's Image Classification Wizard?

texts:
I am currently not licensed for the Image Analyst or Spatial Analyst extensions of ArcGIS Pro.  However, I am contemplating trialling/purchasing Spatial Analyst to perform image classification on a mosaic dataset.

I have been reviewing ArcGIS Pro Help pages like:

Overview of image classification;
The Image Classification Wizard; and
Mosaic datasets

but have not been able to find anything that says whether a mosaic dataset can be used as input to ArcGIS Pro's Image Classification Wizard.

If you know that it can, then can you point me at an ArcGIS Pro Help page that confirms it?


texts: Edit an old georeference raster projection in QGIS

texts:
I have created a map from a picture. And now I have some more detailed pictures and georeferenced them. Is there any way to quickly access them later on.

I mean everything is fine with georeferencing. The picture is projected right where it shall be at the end. But I would like to know if there is a way to add some reference points or even delete some.

I know I could delete the picture and load the raster again and start from a new georeferencing.

So is there a possibility to edit georeferenced pictures?


texts:
Kyngchaos installer is not finding python3. I installed python3 using HomeBrew and is not located at /usr/local/Cellar/python/3.6.4_4/Frameworks/Python.framework/Versions/3.6/bin/python3.6

How do I tell the QGIS installer to recognize that python3 is installed?

I had previously installed GDAL using HomeBrew, so tried reinstalling (or installing another version wherever it might end up) thinking it might help with the situation, but the QGIS3 installer still can't find python 3.

The installation error is a popup that comes up as soon as the installer is open and says "python3 QGIS requires Python 3.6."

For non Mac users, most seem to rely on Kyngchaos to provide an installer for QGIS and he has just released the installer for v3.

Update: answered my own questions below. Don't use HomeBrew unless you know enough to point the QGIS installer/QGIS to the right version (which I certainly don't), but use the python3 installer version linked to in the ReadMe.


texts: Difference between Landsat 7 and Landsat 8 TOA reflectance computation

texts:
In section 3.1 of this article it is explained how the Top of Atmosphere Reflectance for Landsat 7 and Landsat 8 are computed. The equation for Landsat 8 is not mentioned explicitly, but can be found here on page 61 (PDF page 69).

What is the actual difference? (Why) are the products still comparable?


texts: How to find out if a WFS instance supports filtering via CQL?

texts:
I am considering using Common Query Language (CQL) to add filter parameter to a HTTP curl query which runs against a WFS server. I found out that GeoServer supports CQL via WFS.

Can I tell from a GetCapabilities document if the WFS is served by GeoServer?
How can I test if CQL is supported by the WFS?


texts: Webservices for GPS (Upload on DB and display it)

texts:
i'm doing an application where my users would be able to upload their gps tracks. My idea is to put their tracks on mysql database and display it on appropriate map.
Question:How can i do that? My first idea is use garmin webservice to upload data, and then, show that on a map using openlayers lib. This is a good approach? What direction should i take?

thanks in advance


texts: Sequence of Iterate Field Values in ModelBuilder?

texts:
I'm pretty new to ModelBuilder and have been working on a tool that separates an exotic vegetation point feature class by species, then separates each species by treatment action, and then date.

I made a model using Iterate Field Values to create a new output feature class for each species in the input feature class.  This portion worked great, but now I want to be able to go through each output feature class and iterate for another field value: treatment action.  I have been researching nested models/integrating a model within a model but have not been successful.  When I try to insert a new iterator, all the options are grayed out.

Here is a screenshot of the model that has worked so far and produced new feature classes for each species (taxon):

Is there any way to make the output values act as an input for the next iteration?


texts: Overlay single, non georeferenced image in Leaflet?

texts:
I'm wondering whether it would be possible to somehow overlay a 3rd party, changing non-georeferenced image such as:

If the image was static, I'd use mapwarper.net. But it's not. So is there a way to add coordinates to the image's corners, then bend it on the fly, allowing it to be overlaid?

UPDATE

Confirming that ImageOverlay works:


texts: OpenLayers 5: this.source.loadFeatures is not a function with Cluster

texts:
I'm trying to change the example of the Clustered Feature using my points from geojson layer but there is something wrong.

var random_points = new ol.layer.Vector({
  title: 'Random Points',
  source: new ol.source.Vector({
    url: '{% url 'random_points_geojson' %}',
    format: new ol.format.GeoJSON(),
  }),
});

var clusterSource = new ol.source.Cluster({
  source: random_points,
});

var clusters = new ol.layer.Vector({
  source: clusterSource,
  style: function(feature) {
      style = new ol.style.Style({
        image: new ol.style.CircleStyle({
          radius: 10,
          stroke: new ol.style.Stroke({
            color: '#fff'
          }),
          fill: new ol.style.Fill({
            color: '#3399CC'
          })
        }),
        text: new ol.style.Text({
          text: size.toString(),
          fill: new ol.style.Fill({
            color: '#fff'
          })
        })
      });
    }
});

Withouth the cluster function I can see the random points vector on the map but if I use the code above I see in console this error:

Uncaught TypeError: this.source.loadFeatures is not a function
    at e.loadFeatures (VM215 ol.js:7)
    at e.prepareFrame (VM215 ol.js:7)
    at e.renderFrame (VM215 ol.js:7)
    at e.renderFrame_ (VM215 ol.js:7)
    at e.<anonymous> (VM215 ol.js:7)
Cluster.js:117


texts: Attribute information not showing in Silverlight viewer popup

texts:
I have feature service which I am accessing through a Silverlight viewer. The popups contain all the available fields (about 40), most of which have coded value domains. Looking at the popup, these fields appear blank. If I click on the Edit Values button, the fields' existing attribute information appears (which I can see in ArcMap).

This might be an issue with version 3.1 of the viewer; similar issues are mentioned on the ESRI forums. I have not had this problem before, and I could see the values in the viewer three days ago. The only thing that's different is I had to republish the service after making cosmetic changes to the map. The data, field aliases etc all stayed the same.


texts:
On calling API http://osrm.stackbox.xyz/table/v1/driving/78.4444076,11.6000205;78.5929274,11.5958152, I am getting a 2x2 matrix.

When I put the same coordinates with some more co-ordinates like this:
http://osrm.stackbox.xyz/table/v1/driving/78.4444076,11.5870854;78.5929274,11.5958941;78.5938641,11.6000205;78.5961901,11.5958152, it gives me expected 4x4 matrix.

The issue is it's giving different distance(durations) matrix for points 78.4444076,11.6000205;78.5929274,11.5958152


texts: What cartographic line symbol could make a coastline represent water like this old map?

texts:
I live in a coastal area with lots of islands. Sometimes with a black and white map it's hard to distinguish land from sea. So I'm trying to replicate this old map's style.

I'm trying to do this in QGIS but if it's easier in ArcGIS I'd be happy to hear how to do it in there as well. At first I tried to make multiple simple lines and offset the duplicates but it offsets them inside as well as outside.  So then I made an SVG with offset lines and I got closer to what I'm aiming for but since the SVG is a rectangle shape there are gaps and overlaps at every angle.

Does anyone have any tips?


texts: Multiple leader lines on one annotation (non feature-linked)

texts:
Is it possible to link multiple leaders lines to one piece of annotation?

I'm working in ArcGIS 10.2 Desktop Advanced with non-feature linked geodatabase annotation. I have hundreds of numbers labeling various polygons, but some polygons share a number. Instead of creating multiple numbers/annotations for each polygon I'd like to have one number with multiple lines pointing to each polygon.

So, in the image below, the number 5 corresponds to multiple polygons. Here I've moved 5 from it's ideal location to show that the second leader line is not connected to the number annotation. My workaround has been to create another piece of geodatabase annotation with a blank text string to serve as a placeholder for a leader line. The other option was to use map/graphic annotation to create the second leader line, but with the multiple disadvantages of not being in the same annotation feature class as the number annotation.

My goal would be for each piece of number annotation to be connected to multiple leader lines like appears below.


texts: English translations of QGIS commands, algorithms and menu items

texts:
For posting here at Gis.StackExchange, I often have to find the english equivalent of commands, functions, algorithms, menu items etc. of my german language version of QGIS. Sometimes, a word-by-word translation results in the exact equivalent in an english-language QGIS-installation, but sometimes they are named quite differently. So it's an effort to always find the correct description to point to the relevant menu, algorithm etc. Is there somewhere a kind of lexicon/dictionary with translations? Something like:
german: Eigenschaften, englisch: properties
german: Werkzeugkiste, englisch: processing toolbox
etc.
To have a translation of the official QGIS terminology would help a lot.


texts: Django ORM geo function ST_AsMVTGeom always returns type "bytea" not "geometry"

texts:
I want to use POSTGIS function ST_AsMVTGeom with django ORM. Therefore I create following geo function:
from django.contrib.gis.db.models.functions import Transform, GeoFunc

class AsMVTGeom(GeoFunc):
    function = "ST_AsMVTGeom"
    geom_param_pos = (0, )
    output_field = GeometryField(srid=4326)

and call it with following query (example bounding box):
mvt_geom = (
        municipality.objects
        .annotate(
            mvt_geom=AsMVTGeom(
                Transform("geom", 3857),
                Value("ST_MakeBox2D(ST_Point(0, 0), ST_Point(4096, 4096))"),
                4096,
                0,
                False,
            )
        )
    )
    print(str(mvt_geom.query))

This gives me something like this (for simple model containing a geom column):
SELECT "map_municipality"."id", "map_municipality"."ags", "map_municipality"."nuts", "map_municipality"."name", "map_municipality"."type", ST_AsMVTGeom(ST_Transform("map_municipality"."geom", 3857), ST_MakeBox2D(ST_Point(0, 0), ST_Point(4096, 4096)), 4096, 0, False)::bytea AS "mvt_geom" FROM "map_municipality"

But as I want to use outcome in further queries, I need type of mvt_geom to be geometry not bytea. If I remove casting to bytea, everything is fine.
I somehow cannot change it to type geometry, trying for hours...
Even explicit casting to GeometryField does not help.
This user seems to have a similar problem:
https://stackoverflow.com/q/56475858/5804947


texts:
I have several tools working with QGIS rubberband
before activating my new map tool i want to make sure there is nothing on rubberband
previously up to QGIS 3.2 I was able to clean the rubberband via the following code
rbs = [ i for i in canvas.scene().items()
                    if issubclass(type(i), qgis._gui.QgsRubberBand) ]
if len(rbs)>0:
                    for rb in rbs:
                        if rb in canvas.scene().items():
                            canvas.scene().removeItem(rb)

This code no longer works(as QGIS adds some default rubberband and when I try this QGIS crashes
now the question is how to have similar functionality (cleaning map canvas from all rubberbands before starting)


texts:
I am working through this tutorial that I think was intended for 2.8x -  https://www.qgistutorials.com/en/docs/making_a_map.html

Step 24 advises on how to add a grid on a map in Layout view. But there is no longer a CRS menu in Item Properties. In the tutorial, the grid appears when the coordinate reference system is set to WGS84 in the CRS Selector dialog box. Then the interval is set a 5 degrees. I have set the CRS to WGS84, but has not appeared on the map, and I don't see how to set the grid interval to 5 degrees.

How do I do this in version 3.0?

Here is the map with the version 3.0 Item Properties:


texts:
I have an ArcPy script that does a whole bunch of processes with ArcGIS 10.1 to write the output data to a SDE geodatabase. What I want to achieve is to enable an online user to call the script and then do the job on the backend automatically.

I figured out CGI and made a mapping to py on IIS.

However, everything I called the script from the web. It threw out the following error:

  HTTP Error 502.2 - Bad Gateway The specified CGI application
  misbehaved by not returning a complete set of HTTP headers. The
  headers it did return are "Creating catid files
  ESRI\Server10.1\Configuration\CATID\arcgis-response-streamer.catid.ecfg
  Create JARID File...
  outfile=ESRI\Server10.1\Configuration\JARID\arcgis-response-streamer.jarid.ecfg
  Creating catid files
  ESRI\Server10.1\Configuration\CATID\arcgis_indexer.catid.ecfg Create
  JARID File...
  outfile=ESRI\Server10.1\Configuration\JARID\arcgis_indexer.jarid.ecfg
  Content-Type: text/plain;charset=utf-8 Hello World! starts Hello
  World! ends ".

I am just wondering if ArcPy script can be exposed to CGI like a normal python script?


texts:
I ran into this problem in PostGIS:

raster2pgsql -c -I -C -s 4326  -N -9999  -M -Y  -F ./test2/*.tif  rdem2 > rdem2.sql

  Processing 84/84: ./test2/dem2_4326_06_14.tif

  ERROR: diff_rastinfo: Could not run raster alignment test

Does someone have any ideas about possible reasons for this problem ?


texts: How to set visibility scales for ArcGIS for Server image service?

texts: I have created a mosaic dataset in ArcGIS 10.3.1 and from it I created an image service which works okay. However, how would I set the visibility scale range on this service so that, for example, the service only displays between 1:10,000 and 1:50,000? I could of course load the layer into an mxd and set the visibility scale here and then publish as a map service, but then that wouldn´t be an image service any more.


texts: QGIS - Average a DEM into a single, flat plane

texts:
I have a XYZ point cloud that I am interpolating into a DEM. The data represents a surface that is relatively flat. I am hoping to be able to create a DEM of a single flat plane that best represents the input DEM.

Any pointers?


texts: Missing labels when using arcpy.mapping.ExportToPDF

texts:
I have a script that is run on a daily basis to export a folder of MXDs to PDFs and I have noticed some peculiar behavior. Many, but not all of the PDFs will be missing some labels or their label orientation gets changed. These errors do not occur when I manually export from arcmap. In an attempt to isolate the problem, I wrote a code snippet to simply export a PDF and it behaved normally, no errors.

Is there a reason that stepping through a folder of MXDs would cause this error?

I have included the code below.

    import arcpy
    import os
    arcpy.env.workspace = r"\\fs1\shared\GIS\Map_Documents\NightlyUpdates"
    ws = arcpy.env.workspace
    outPDF = r"\\fs1\shared\GIS\Map_Documents\NightlyUpdates\PDF"
    MXDlist = arcpy.ListFiles("*.mxd")
    for MXD in MXDlist:
        root = os.path.splitext(MXD)[0]
        title = '\\'+ root + ".pdf"
        tempMXD = arcpy.mapping.MapDocument(os.path.join(ws,MXD))
        arcpy.mapping.ExportToPDF(tempMXD, outPDF+title, image_quality="BETTER")
        del tempMXD


texts: ST_GeomFromGeoJSON strips coordinates, is it a bug?

texts:
I'm storing some geo paths in PostGIS using GeoJSON and the ST_GeomFromGeoJSON()-function but I just found an issue with it:

SELECT ST_NumPoints(ST_GeomFromGeoJSON('{"type":"LineString","coordinates":[[-122.023468,37.330232,0],[-122.02345,37.330233,0],[-122.02345,37.330233,0],[-122.023405,37.330236,0]]}')) = 4; -- FALSE

SELECT ST_NumPoints(ST_GeomFromText('LINESTRING (-122.023468 37.330232 0, -122.02345 37.330233 0, -122.02345 37.330233 0, -122.023405 37.330236 0)')) = 4; -- TRUE

It seems to strip the two 'middle' coordinates because they are the same. And I can buy that it would do that but it doesn't do it consistently, like when I use ST_GeomFromText() with the same list of coordinates. It also does not strip the coordinates when using "type": "MultiPoints" instead, but it would be the wrong semantics.

It's a problem for me because I need the points to match a sequence of timestamps as well or data will be lost. And the reason I'm using GeoJSON is that the app server is written in Javascript so it just seems natural.

I'm running Postgres 9.4.1 and SELECT PostGIS_Full_Version() gives me:

POSTGIS="2.1.5 r13152" GEOS="3.4.2-CAPI-1.8.2 r3921" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.11.1, released 2014/09/24" LIBXML="2.9.1" LIBJSON="UNKNOWN" RASTER


texts: Choosing the range of the correlogram with an irregularly shaped area

texts:
I am doing correlograms (Moran's I) between 125 plots. Those plots are located in large patches in my study area (see figure).

The entire study area is approximately 2000 x 2000 m. I know that the usual range for a correlogram is 1/2 the value of the length of the study area. But in my case, there is large section of the 2000 x 2000 area that are not sampled.

How to choose the range of my correlogram?


texts: Remove color table from TIFF in QGIS 3.8 because virtual raster does not work due to color table issue

texts:
I have a set of 14 rasters that I would like to create a virtual raster (vrt) from. When I create the vrt one of the rasters does not appear in the new file. The raster have identical formats except the missing one displays  in the layer styling box and the others show . It seems like this is due to the missing one not having a color table or palette where as the others do. I have gotten around this by creating a vrt of the tif without the color palette, opening the vrt with a text editor  and changing <ColorInterp>Gray</ColorInterp>  to <ColorInterp>Palette</ColorInterp> and then creating a vrt with the other 13 and the hacked vrt.

Is there a way I could instead convert the file to have a color table or better remove the color table from the others?
From my tests warp does not strip the color table.


texts:
According to the QgsDistanceArea docs, an azimuth for computeSpheroidProject should be

  in radians, clockwise from North

So, to build a new point with computeSpheroidProject i use math.radians(azimuth)

pt1 = QgsPointXY(30.30, 60.60)
distance = 100
azimuth = 45

da = QgsDistanceArea()
da.setEllipsoid('WGS84')
pt2 = da.computeSpheroidProject(pt1, distance, math.radians(azimuth))
print(pt1, pt2)
>>><QgsPointXY: POINT(30.30000000000000071 60.60000000000000142)>
<QgsPointXY: POINT(30.30129068361398836 60.6006346117680792)>

But when needed to find an azimuth between pt1 and pt2 it doesn't return the original value (45):

azimuth = pt1.azimuth(pt2)
print(azimuth)
>>>63.81727521235025

Looks like I haven't quite figured out azimuth calculations. Why results are different?


texts: KML raster not showing images when using another users machine

texts:
I have converted some geotiff's into raster KML file packages, (which consist of a KML file and two folders - images and files). I have created these KML's as superoverlay KML files.
I am able to successfully load the KML into Google Earth EC on my computer and see the output.
I have directed a different user on another PC to load the KML file - which they can, however it only takes them to the extent of the data and does not show the images associated with the KML. The tree structure on the right also does not show the associated files as per on my machine.
Additionally, the user in question that is having the problem has full read/write access to the folder in which the KML exists and also the images + files subfolders.
Why do the images not show for the other user, but they do for myself?
Sample KML contents:
<Document>
<name>e_100y_max_d</name>
<visibility>1</visibility>
<NetworkLink id="kml_1">
<name>e_100y_max_d1</name>
<Region>
<LatLonAltBox>
<north>-21.73210771550464</north>
<south>-21.96438026963675</south>
<east>116.50377196418752</east>
<west>116.25510594725066</west>
</LatLonAltBox>
<Lod>
<minLodPixels>256</minLodPixels>
<maxLodPixels>-1</maxLodPixels>
</Lod>
</Region>
<Link>
<href>files/e_100y_max_d1.kml</href>
<viewRefreshMode>onRegion</viewRefreshMode>
</Link>
</NetworkLink>
</Document></kml>

Addtionally, the KML files in the 'files' subfolder, contain references to the images as per here:
<href>../images/e_100y_max_d1_2_1_0.tif</href>


texts: Get position/rotation of free labels in QGIS

texts:
I want to export the position and rotation of labels with placement "Free" as I need them in Mapbox, which has inferior label placement compared to QGIS.

I've tried with various polygon-to-point operations (centroid, pole of inacessibility) and they don't give me the same result and in any case don't contain rotation.

Essentially, I want to copy this to Mapbox:

Similar to Turn line labels in QGIS to a new point shape layer - but I have a polygon layer and that method won't work.


texts: Export long lat from polyline as csv file in QGIS

texts:
I want to draw a line in QGIS and then, export the long-lat of each point into a Comma Separated Value (csv) file. I tried doing this by creating a new shapefile layer and then using the Add Feature tool to draw my line. I then tried to export my line via Save As on the layer and exporting as CSV.

Unfortunately, the only thing that gets saved to the CSV file is the id of the layer. I've been unable to figure out how to export the coordinates. Some help please.


texts: What are the differences between formats for storing styles?

texts: I'm looking for a format for store information about data styles for my web application. There are different formats for describing vector data styles. For example, SLD, Mapnik XML, and there may be others. What are the strengths and weaknesses of these formats? What are the other alternatives?


texts:
I am wanting to come up with a way to speed up my Imagery Index layer.
To give you an overview I have about 500 imagery files. They cover extents all over Australia.
I am wanting to create a shapefile with the outlines of the imagery extents and populate it with the date the imagery was taken and region its in. I am then wanting to turn this into a layer so I can use it to see what areas our imagery covers and it makes it easy to find which image I am after.
Currently I am just dragging in each image, drawing a box around it and then populating it with year and area.
I was told there is was way I can point my Python script to a folder and it will
draw the imagery extents as in that folder and add source location to the metadata.
Any ideas?
Below shows shapefile extent over imagery.


texts: How to let public query Enterprise Data in ESRI Webapp

texts:
I've got a webapp set up to display the information on backflow preventer devices in our city, but would like the public to be able to track their own information. Currently, they have to enter their address into the search bar to find their house, then click on the backflow point to find the information we have regarding when the last test was performed and the next test is due by. I would like to automate this somewhat so that all they have to do is enter their address and the information on test dates performed and due is displayed.
The attribute table has one field for address number, one for street name, and one each for date of last test and date the next test is due.


texts: Vegetation detection of anthropic crops with Google Earth Engine

texts: I am very new using Google Engine and I need to know how to detect only agricultural crops on slopes without them mixing with that natural vegetation on hills. I have tried different indexes (like NDVI in different time periods) and Qgis SCP plugin, but I still have a margin of error and would like to fix it with GEE if that is possible. Thanks a lot


texts: Adding a latitude/longitude graticule in CartoDB (e.g. WMS)

texts:
I'd like to add a basemap that has a graticule in latitude/ longitude to the Map View of a CartoDB project.  This could be in the WMS format for example, but I haven't been able to find one that works.

For my ArcMAP work, I first tried this NOAA WMS link, but got an "unspecified error" message. When I tried it in CartoBD, I get options for 1, 5 and 10 degree, but these end up with a blank background.

http://maps.ngdc.noaa.gov/arcgis/services/graticule/MapServer/WMSServer?request=GetCapabilities&service=WMS

I can email NOAA about the error, but ideally, I would prefer a satellite or terrain basemap that has lat/longs included in the map.  So, I thought I'd ask the forum first.

For example, I have an ArcGIS server link from NOAA that works great in ArcMAP:  GIS Servers\arcgis on maps.ngdc.noaa.gov\arctic with many layers of interest.


texts: Editor tracking on geometric network feature classes?

texts:
Editor tracking is not working on my geometric network feature classes.

Although I enabled the editor tracking - dates and times are not recorded in layers while editing.

Is it impossible? or a local problem with my geodatabase?


texts: ArcCatalog 10.1 tree doesn't refresh after adding new custom IGxObject at runtime

texts:
I am using ArcObjects 10.1 (.NET) to customize ArcCatalog: adding new nodes in catalog tree like in the Esri's example. But when I add new child node to already expanded parent, it doesn't show up in tree. Actually, catalog tree node never updates again after being expanded.

Any solutions?

Tried this solution, but it doesn't work.


texts: st_transform reprojection issue with R and sf

texts:
According to https://geocompr.robinlovelace.net/reproj-geo-data.html i tried to reproject world map to Lambert with st_transform from sf package

library(spData)
library(sf)
library(tmap)

# test
plot(world["continent"])

# reprojection Lambert centered on Atlantic ocean

world2 <- st_transform(world["continent"],  crs = "+proj=laea +x_0=0 +y_0=0 +lon_0=-30 +lat_0=10")

strange mapn(yellow)

# plot
 plot(world2)

 # plotting with tmap (what I want but with an issue)
 tm_shape(world2) +
 tm_fill()

 tmap_mode("plot")

 error : The shape world2 is invalid. See sf::st_is_valid


texts: Using ST_Difference to remove overlapping features?

texts:
I am trying to use ST_Difference to create a set of polygons (processing.trimmedparcelsnew) that do not contain any of the area covered by another set of polygons (test.single_geometry_1) using PostGis 2.1 (and Postgres SQL 9.3). Here is my query:

CREATE TABLE processing.trimmedparcelsnew AS
SELECT
    orig.id, ST_Difference(orig.geom, cont.geom) AS difference
FROM
    test.single_geometry_1 cont,
    test.multi_geometry_1 orig;

But the resulting polygons have not been trimmed, instead they seem to have been split where they intersect with the other layer. I have tried just running the select without putting the result in a table and everything else I can think of, but I can't seem to get this function to work.

I have attached a picture of the result

After comments, I have tried adding a WHERE clause. I want the parcels that have no intersections, and the intersecting areas of the other parcels removed (the layer test.single_geometry represents contamination I want removed from my parcels). I tried an intersect but of course I actually want the non intersections so I am now trying a disjoint. I have also tried adding the orig to my table but the documentation for ST_Difference (http://postgis.net/docs/ST_Difference.html) does say it returns the exact geometry I need (a geometry that represents that part of geometry A that does not intersect with geometry B), so I am confused as to why I would want the original polygon in my table instead. Anyway, here is my modified code:

CREATE TABLE processing.trimmedparcelsnew AS
SELECT
    orig.id, ST_Difference(orig.geom, cont.geom) AS difference, orig.geom AS geom
FROM
    test.single_geometry_1 cont,
    test.multi_geometry_1 orig
WHERE ST_Disjoint(orig.geom, cont.geom);

Following from dbaston's answer I have now tried:

CREATE TABLE processing.parcels_trimmed AS
SELECT id, COALESCE(ST_Difference(geom, (SELECT ST_Union(b.geom)
                                         FROM test.single_geometry_1 b
                                         WHERE ST_Intersects(a.geom, b.geom)
                                         AND a.id != b.id)), a.geom)
FROM test.multi_geometry_1 a;

The result of this is just a copy of test.multi_geometry_1. Though now the splitting is no longer occurring.

I tried the earlier version, but again just get a copy of test.multi_geometry_1:

CREATE TABLE processing.parcels_trimmed_no_coalesce AS
SELECT id, COALESCE(ST_Difference(geom, (SELECT ST_Union(b.geom)
                                         FROM test.single_geometry_1 b
                                         WHERE ST_Intersects(a.geom, b.geom)
                                         AND a.id != b.id)), a.geom)
FROM test.multi_geometry_1 a;

I am starting to wonder if there is something else I am doing wrong? The proceeding statement is:

DROP TABLE IF EXISTS processing.parcels_trimmed_no_coalesce;

And I am running the queries from the PostgreSQL SQL query window and Openjump.

The statement I use to see the table is:

SELECT * FROM processing.parcels_trimmed_no_coalesce;

In the interest of simplification I have now reduced this query down to just:

SELECT id, COALESCE(ST_Difference(geom, (SELECT ST_Union(b.geom)
                                         FROM test.geometriestocutagainst b
                                         WHERE ST_Intersects(a.geom, b.geom)
                                         AND a.id != b.id)), a.geom)
FROM test.geometriestocut a;

This still results in just the original polygons (test.geometriestocut) when the desired result is the original trimmed against test.geometriestocutagainst.


texts: How to align Aerial images from MassGIS with OSM using python?

texts:
MassGIS Dataset: I'm using the MassGIS Data: USGS Color Ortho Imagery (2013/2014) aerial images of Boston. Basically, the dataset is divided into small patches of jpeg2000 images, where each of the images have a corresponding xml file. One important piece of data in the xml file is the bound of that image in latitude and longitude. They are using EPSG:26986 while reporting this. The images are ortho-rectified.

OSM Dataset: I donwloaded just the building polygons of Boston using the Overpass Turbo API. As far as I understand, they are denoted using WGS84 datum and the corresponding EPSG is 4326.

What I want to do: For some image from the MassGIS dataset, I want to draw the building polygons on top of it using python.

What I tried: For the given image, I retrieve the bounding latitudes and longitudes of the image from the corresponding xml file. EPSG:26986 and EPSG:4326 seems similar enough, so I don't convert one to another. Then for all of OSM building polygons that lie inside the MassGIS image bound, I find the corresponding pixel position of each node of the building polygons using a linear transformation-

x = (longitude - longitude_bound_min) / (longitude_bound_max - longitude_bound_min)
                                    * image_pixel_width
y = (latitude - latitude_bound_min) / (latitude_bound_min - latitude_bound_max)
                                    * image_pixel_height

Doing the above, I am getting a mis-alignment between MassGIS images and OSM building polygons, as you can find here. My understanding was that since the images are ortho-rectified, this transformation should be enough. I would be really glad if anyone gave me any hint about what I am misunderstanding here.

Edit: For one of the MassGIS aerial images, I get the following denoted as a bound for the image in the xml file-

<bounding>
  <westbc>-71.240667</westbc>
  <eastbc>-71.222039</eastbc>
  <northbc>+42.196636</northbc>
  <southbc>+42.182783</southbc>
</bounding>

There is also a corresponding .aux file for image, running `gdalinfo' on which gives me-

    Driver: HFA/Erdas Imagine Images (.img)
Files: 19TCG150725.aux
Size is 5000, 5000
Coordinate System is:
PROJCS["NAD_1983_UTM_Zone_19N",
    GEOGCS["GCS_North_American_1983",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS_1980",6378137,298.257222101]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-69],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
Origin = (315000.000000000000000,4674000.000000000000000)
Pixel Size = (0.300000000000000,-0.300000000000000)
Corner Coordinates:
Upper Left  (  315000.000, 4674000.000) ( 71d14'26.40"W, 42d11'46.62"N)
Lower Left  (  315000.000, 4672500.000) ( 71d14'24.68"W, 42d10'58.02"N)
Upper Right (  316500.000, 4674000.000) ( 71d13'21.04"W, 42d11'47.89"N)
Lower Right (  316500.000, 4672500.000) ( 71d13'19.34"W, 42d10'59.29"N)
Center      (  315750.000, 4673250.000) ( 71d13'52.87"W, 42d11'22.96"N)
Band 1 Block=64x64 Type=Byte, ColorInterp=Undefined
  Description = Layer_1
  Metadata:
    LAYER_TYPE=athematic

Band 2 Block=64x64 Type=Byte, ColorInterp=Undefined
  Description = Layer_2
  Metadata:
    LAYER_TYPE=athematic

Band 3 Block=64x64 Type=Byte, ColorInterp=Undefined
  Description = Layer_3
  Metadata:
    LAYER_TYPE=athematic

Band 4 Block=64x64 Type=Byte, ColorInterp=Undefined
  Description = Layer_4
  Metadata:
    LAYER_TYPE=athematic


texts: Field Type Date invalid when run in Python, but not in Model?

texts:
Using ArcGIS Desktop Advanced 10.1 and PyScripter

When running the arcpy.FeatureClassToFeatureClass_conversion, I get the error that the field type of Date is not valid.  However, when running this same procedure in the model where the python was exported from, it runs successfully.
I have verified that the field type is in fact Date for the field and that the field length is sufficient.
This seeems to be an issue after moving to 10.1

Error Message:

ExecuteError: ERROR 000278: Field Mapping error: Field type invalid for output field xxx and input field xxx

Any ideas?


texts:
We are looking for a way using SQL Server spatial queries (or spatial library if needed) to do the following:

We have our own set of spatial data for custom roads (single line) and many point along that road. Some are checkpoints and some are just pass throughs.

We need a way to query what the next checkpoint is heading to a destination based on a current location on that line.

Is this possible with SQL Server spatial?


texts: ArcGIS Pro SDK OpenRelationshipClass does not return results

texts:
I have been struggling with this issue for hours now and am wondering if anyone else has encountered this.

I have a feature class (ADDRESS_MANAGER.GIS_INF_ADDRESS_P) and a table (ADDRESS_MANAGER.GIS_TBL_ADDRESSCARD_B).

Using ArcCatalog, I have verified the IDs of the feature class is 987 and the address card table is 981.  I have verified the relationship exists from the feature class to the address card table (Relationship = ADDRESS_MANAGER.AddressCardHasAddress / Label = AddressCard / Origin=ADDRESS_MANAGER.GIS_INF_ADDRESS_P).

However, no relationship classes are returned when I run this code:

var sourceId = sourceTable.GetID().ToString(); // returns "987"
var destId = destTable.GetID().ToString(); // returns "981"
IReadOnlyList<RelationshipClass> relationshipClasses = geodb.OpenRelationshipClass(sourceId, destId);

I have verified the geodb object does indeed point to the geodatabase that holds both of these layers -- geodb is set directly from the feature class.  Both layer IDs are correct, and there is definitely a relationship class between them.  I am porting old ArcObject code over that did something similar, and I cannot figure out how to actually get this RelationshipClass anymore without hardcoding a reference to it.

This is ran in QueuedTask.Run per the documentation.  It doesn't seem to matter if I set up geodb by opening the appropriate SDE file, or if I set it from the feature class (the geodb is an versioned enterprise geodabatase -- Oracle).  Note that I can open the relationship class by using (RelationshipClass relationshipClass = geodb.OpenDataset<RelationshipClass>("ADDRESS_MANAGER.AddressCardHasAddress"))

This is in ArcGIS Pro 2.2 - using C# / .NET 4.6.2


texts:
I have a function who uses a layer already styled using another function and saves it as a gpkg and I want to save the style in that gpkg:
def saving_gpkg(styled_layer):
    options = QgsVectorFileWriter.SaveVectorOptions()
    options.driverName = "GPKG"
    QgsVectorFileWriter.writeAsVectorFormatV2(styled_layer, "nice_filename", QgsCoordinateTransformContext(), options)
    styled_layer.saveStyleToDatabase("polygon", "style automatically created", True, "")

The gpkg is created, but if I open the created file, it doesn't have the defined style, just a random singlesymbol. I'm missing something?


texts: Conversion from Raster to Vector in Google Earth Engine?

texts:
I have applied unsupervised classification on an image in Google Earth Engine using k means clustering algorithm. I have made 15 classes and I have visualized those 15 classes by adding a map layer. It is perfect. All 15 classes are represented by a separate color using 'randomVisualizer'.

How can I convert this raster image into a vector feature collection with 15 objects?

All the polygons of class 1 to become a single object. So my feature collection would have 15 objects and I can select the each class with its id.If possible please demonstrate this via dummy example using any image and unsupervised classification.


texts: How to add attributes from polygon layer to point layer when layers have different projections

texts:
I have two layers in different coordinate systems:

polygon layer with catchment areas and names of the catchments (EPSG 21781)
point layer with coordinates and flood flows (EPSG 2056)

They of course overlay nicely in QGIS because I have set their CRS correctly.

Goal:
Add the attribute "catchment name" to the point layer

What I have tried so far:

Join attributes by location is no option because it is a large dataset and crashed when I tried it.
I have tried to use SAGA Add attribute to point but it needs both layers in the same CRS and fails to add the attributes to the point layer.

Is there another tool to add attributes from polygon to point layer?
I am using QGIS 2.14.6 on Win7


texts: PostGIS nearest point with LATERAL JOIN in PostgreSQL 9.3+

texts:
I need to obtain on each element on one table the closest point of another table. The first table contains traffic signs and the second one the Entrance Halls of the town.
The thing is that I can't use ST_ClosestPoint function and I have to use ST_Distance function and get the min(ST_distance) record but I am quite stuck building the query.

CREATE TABLE traffic_signs
(
  id numeric(8,0) ),
  "GEOMETRY" geometry,
  CONSTRAINT traffic_signs_pkey PRIMARY KEY (id),
  CONSTRAINT traffic_signs_id_key UNIQUE (id)
)
WITH (
  OIDS=TRUE
);

CREATE TABLE entrance_halls
(
  id numeric(8,0) ),
  "GEOMETRY" geometry,
  CONSTRAINT entrance_halls_pkey PRIMARY KEY (id),
  CONSTRAINT entrance_halls_id_key UNIQUE (id)
)
WITH (
  OIDS=TRUE
);

I need to obtain the id of the closest entrance_hall of every traffic_sign.

I already know how to do it by joining the whole two tables. See question and answers here: PostGIS nearest point with ST_Distance

In this mentioned thread it is said that a new cool feature is out in PostgreSQL 9.3+. I would like to use this cool feature in PostgreSQL 9.3+ called LATERAL JOIN to achieve this same thing, but there is not much info about this feature and I don't figure out how it would be in this case to get the closest point of another table.

Hope someone can clear this out.


texts: Technical and usability differences between OpenLayers and Google Maps API V3

texts:
I know this is supposed to be a wiki question, but the feature has been taken of the site

What are the main differences between OpenLayers and Google Maps api?
What are the strengths and weaknesses of each system?
Where would you prefer one system over the other?
etc?

These are just a few questions, feel free to answer more.


texts:
I added a shapefile to the PostGIS using,

ogr2ogr -f PostgreSQL PG:"host=server_ip user=username dbname=dbname password=password" yourshapefile.shp;

Now I want to remove it. What is the command line code for that?

The table list is as follows \d,

List of relations
 Schema |         Name         |   Type   |  Owner
--------+----------------------+----------+---------
 public | pimapine             | table    | smaranh
 public | pimapine_ogc_fid_seq | sequence | smaranh
(2 rows)

And when I execute DROP TABLE public.pimapine there is no change to the output of \d.


texts: Join attributes to nearest without repeating feature

texts:
I have two layers, one of point and one of polygons, each with approximately 1000 features. I need to join the attributes of the points to the nearest polygons.

The problem is that point attributes cannot be repeated in more than one polygon, for example if point 1 has already been assigned to any polygon it cannot be assigned to another polygon.

Point -> Polygon

1     ->   A

2     ->   B

3     ->   C

4     ->   D

5     ->   E

I tried "nnjoin and join by nearest" but it doesn't allow iterating that way.
Does anyone know how it is possible to do this? I use QGIS.


texts: Anyone have any idea how to create (or where to find) the Waterman "Butterfly" projection in ArcMap?


texts: GeoPandas: how to merge together adjacent polygons?

texts:
I have a GeoPandas dataframe like the following one
f,ax=plt.subplots()
gdf.boundary.plot(ax=ax)

I would like to merge/dissolve the adjacent polygons.
gdf2 = gpd.geoseries.GeoSeries([geom for geom in gdf.unary_union.geoms])

but I get an error
GEOSException: IllegalArgumentException: Argument must be Polygonal or LinearRing

If I check the validity of the geometries I get:
for i in gdf.index:
    print(gdf['geometry'][i].exterior.type)

LinearRing
LinearRing
LinearRing
LinearRing
LinearRing

and:
all(gdf.geometry.is_valid)
True


texts:
I am trying to compute my own metric using LidR. I would like to compute the canopy cover metric as
(# of first returns > 2m) / (total number of first returns)

The package suggested this function from the Examples to compute this metric

myMetrics = function(z,rn){
  first  = rn == 1L
  zfirst = z[first]
  nfirst = length(zfirst)
  above2 = sum(z> 2)
  above2
  x =(above2/nfirst)*100
  x
  # User's metrics
  metrics = list(
    above2aboven1st = x, # Num of returns above 2 divided by num of 1st returns
    zsqmean = sqrt(mean(z^2))  # Quadratic mean of z
  )
  metrics
  # Combined with standard metrics
  return( c(stdmetrics_z(z),metrics))
}

metrics = grid_metrics(las, ~myMetrics(Z, rn=ReturnNumber))

However, the results of the metric using this function does not work correctly, since my values range between 0 and 300

it would be possible to calculate a density metric in height intervals as in FUSION or Lastools?

e.g number of all returns between 0.2 and 3 meters/number of total returns in 20 meters cell


texts:
I have a javascript code that worked well ,
but when convert to python it's doesn't work and raise invalid geojson error .
my code in java script is :
var point = ee.Geometry.Point([52.55845586988192,29.853906673142543])
var dataset = ee.Image("KNTU/LiDARLab/IranLandCover/V1")
var collection = dataset.sample({
 region:point,
 geometries: true,
});
print(dataset2.getInfo());

My code in python is :
point = ee.Geometry.Point([52.55845586988192,29.853906673142543])
dataset = ee.Image("KNTU/LiDARLab/IranLandCover/V1")

collection = dataset.sample({'region':point,'geometries': True})
print(dataset.getInfo())


texts: List of possible events from event handler

texts:
This may be an elementary question but I cannot figure out how to word it so that I can google it properly.

map.on('click', function(evt) {
    var coordinates = evt.coordinate;
    alert(coordinates);
});

What I'm trying to find is the list of possibilities inside of "evt". i.e. - evt.pixel, evt.coordinate, evt.dragging, etc. Is there a list somewhere?


texts: Is exporting to feature class the simplest way to make joins permanent?

texts:
I am a new user doing simple tasks with ArcMap 10.1

I have a shapefile and an excel file, and I want the join to be a permanent part of the shapefile. The way I've been going about this is to export the shapefile into a feature class within a database. Is this the most efficient way to go about saving the joined attributes?


texts: DLTM : Error:Only lon/lat coordinate systems are supported in geography

texts: How can I store and retrieve position in DLTM projection in PostGIS geography.


texts:
I need to display many images (icons) on an OpenLayers map. I'm using v4.6.5

Previously, I was simply displaying 1 icon, and I was using PHP on the backend to determine a proper scale factor. This is no longer possible because I'm displaying up to 250 images, and it takes ~16 seconds to have PHP read the height of each image to adjust the scale factor for all of them, which is far too long of a load time for my webpage.

Users are inputting the images, so it's not possible for me to set them all as a fixed size. Here's a screenshot example of what I'm dealing with currently:

Whereas what I need is for all of those flag images to be set at the same size. I have not been able to figure out through Google or the OL documentation how to set a fixed or absolute size for an image like this. Nor have I been able to find a workaround solution for this instance.

EDIT: Here is a sample excerpt from my code:

       new ol.layer.Vector({
            source: vectorSource = new ol.source.Vector({
                features: [
                    new ol.Feature({
                        geometry: new ol.geom.Point(ol.proj.fromLonLat([<?php echo $row['longitude']; ?>,<?php echo $row['latitude']; ?>])),
                        name: '<?php echo PWBASEURL . "nation/id=".$row['id']; ?>'
                    })
                ]
            }),
            style: new ol.style.Style({
                image: new ol.style.Icon(/** @type {olx.style.IconOptions} */ ({
                    anchor: [1, 1],
                    anchorXUnits: 'fraction',
                    anchorYUnits: 'fraction',
                    src: '<?php echo $row['flag']; ?>',
                    scale: 0.05
                }))
            })
        })


texts:
I am looking for a solution to a problem I have using ArcGIS Online.  I am currently creating 5 Tiled Map Packages which I tile on my ArcGIS for Desktop then Publish to my ArcGIS Online Account to save using any credits.  Once the Tiled Map Package is Published I can then email my colleague the URL (REST) which he then feeds into his ArcGIS Server.

He has asked me if I can group the Services together into one URL so that he can display the 5 Published Services as a Group instead of individual Layers.

Is this possible?


texts:
Once a new QGIS plugin is selected and loaded from the "manage and install plugins" toolbar selection, the plugin automatically appears somewhere on QGIS' toolbar panels.

Often, the new plugin icon is grouped with other icons and cannot be separated out from its assigned group, which may not be ideal.

If the default placement position is not the one desired, how is it moved to another location?


texts: Is city-supplied GIS data sufficient for transit operators or is own set of GIS data needed?

texts:
City GIS data is usually built in a format where streets/avenues/boulevards/etc are modeled as lines that are geographically located in the middle of the streets/avenues/boulevards/etc.  Is that sufficiently precise for big transit operators to manage their buses (real-time bus management and real-time information to users) or do they create their own set of GIS data that gives a more accurate model of bus routes?

If public (or city) GIS data is not sufficient, and given that more accuracy creates larger data sets, where do you draw the line, figuratively and literally?  Do you take the 'middle line' from the public GIS data and just apply an offset on each side to represent both directions of a bus line? And if so, is that offset variable between different segments (since an avenue is wider than a street)? Or do you create 'from scratch' a line that is as close to the actual bus route (and one line per direction)?

Not exactly an expert in this field, so I hope my explanation is clear!

@ Jason Baker and dassouki :

I agree that centerlines are sufficient when you are planning a bus route, my question is more from the operating side. More and more transit operators (and in all but the smallest towns, the operator is a corporation separate from the city itself, albeit with close ties to the city administration) are using real-time bus management systems.  Buses are equipped with systems to continuously broadcast their position to the control center which takes that information and projects it on a city map.  Controllers can then follow every bus in service and make adjustments in real-time as they are needed.  The management system can also send to users real-time information on ETAs at every bus stop in the city.

From my understanding of those types of bus management, they take the coordinates of the bus and they match it to the GIS data they have.  Obviously, buses don't follow centerlines.  On narrow streets, the delta between the centerline and the actual bus position is probably not significant, but most buses take wider boulevards or avenues, where the difference is significant, or at least I expect it is. I'm expecting that difference could create various accuracy issues. For example, when the management system tries to determine if the bus is on-time, late or early, or what information it should send the users about ETA at a bus stop.  I'm guessing that transit operators refine the city GIS data to create a new dataset that is more accurate and better answers their needs.  I see two ways of refining the data (see original post), I'm wondering which one is most widely used, or if there are other ways to address the accuracy issue of city GIS data.

Do not hesitate to point out something that isn't quite clear in my explanation, English is not my native tongue.  ;-)

@ Brad Nesom :

You're correct as far as displaying the information.  I'm looking at it from the management system's algorithm level: any information the system calculates and/or provides logically stems from the delta between the actual bus position, given by an onboard GPS, and the expected position. Now, if the expected bus position (the planned bus route if you will) is created by using basic city GIS data (i.e. the planned bus route follows the centerlines of streets/boulevards/avenues/etc) then the delta with the actual bus location can be quite large, which would, I expect, create accuracy issues with the information the management system calculates and/or provides.  If I'm correct, how are those accuracy issues addressed?  Do transit operator/agencies refine the city GIS data to create a more precise dataset that then allows them to have planned bus routes that are closer to the actual path a real bus would follow in the real world?

The more I explain myself (and thanks to all for making me do so), the more I get the feeling that only a GIS expert within a transit agency could give me the answer I'm looking for.


texts: How do I use Django to execute PostGIS query to get polygons within four points?

texts:
Here's a sample PostGIS query I use to get geometries within four points:

SELECT *
FROM myTable
WHERE ST_MakeEnvelope(-97.82381347656252, 30.250444940663296, -97.65901855468752, 30.29595835209862, 4326) && ST_Transform(myTable.geom,4326);

With this query, I can get all rows within those four points. How do I execute this query or similar queries in Django or GeoDjango?

EDIT: I must also note that my geom column is a MultiPolygon type.


texts: What is the OSGeo4W-equivalent for Linux?

texts:
OSGeo4W is a windows based stack of GeoFOSS which is based upon cygwin.

I know there is a package listing describing what is inside the OSGeo4W binary package and I could fetch and install every package on my own.

But it would be nice to install all packages at once: Is there anything similar for Linux available?

By the way, I'm using Sabayon/Gentoo Linux based on the portage package manager which is not compatible with ELGIS.

Update concerning the reasons:

Do you really need to install all the packages at once? Well: yes. I keep de- and re-installing packages because of certain options, dependencies, etc. I'm currently trying to get file.gdb [this and that] to work with qgis. But I need to de-install GDAL from package manager (which will break certainly a lot of dependencies) and re-compile it with fgdb support. I was looking for some way to install all-at-once. As stated in this question, it took me like one week to get QGIS running with full GRASS- and plugin-support.
Well we have to use OSGeo4W because Windows, unlike every other OS, doesn't have its package manager. Our lives would be a hell of lot easier if it did. Understood. But as I'm not using Ubuntu/Debian I only have some of the important packages in my official repositories, which is quite annoying. That's why I was looking for some meta-package like OSGeo4W.

The OSGeo-Live-distribution looks pretty good though.


texts:
I created a mosaic from multiple overlapping GeoTif files in ArcMap 10.2. Than I created seamlines between the images to define which part of the scene is taken from which original image. Now, I want to export them as one georeferenced tif such that the seamlines are regarded.
Does anyone know how?