texts:
I am trying to select a feature from the layer using my plugin. When working with my plugin in main Python function using QGIS, I receive an AttributeError when I call a class method from within the class itself.
I'm using QGIS 2.8.

Here is my code below:

def handle_mouse_down(point):
    """"Manage selection when clicking on the canvas."""

    layer = iface.activeLayer()
    canvas = iface.mapCanvas()
    layers = canvas.layers()
    clickTool = QgsMapToolEmitPoint(canvas)

    vl_to_query = [l for l in layers if l.name() != u'my_layer_name'][0]
    w = canvas.mapUnitsPerPixel() * 3
    rect = QgsRectangle(
        point.x() - w,
        point.y() - w,
        point.x() + w,
        point.y() + w
    )

    l_rect = canvas.mapSettings().mapToLayerCoordinates(vl_to_query, rect)
    vl_to_query.select(l_rect, False)
    feats = vl_to_query.selectedFeatures()
    vl_to_query.removeSelection()
    clickTool.canvasClicked.disconnect(self.handle_mouse_down)

When I run this function am getting an error like

  File "C:/Users/.qgis2/python/plugins\Myclass\Myclass.py", line 220, in
  handle_mouse_down point.x() - w,
  AttributeError: Myclass instance has no attribute 'x'


texts: Algorithm for smoothing DEM without losing max and min elevations

texts:
I am looking for a method/algorithm for smoothing a DEM raster similiar to something like lowpassing, however I want to keep the highest and lowest points, resp. peaks. Latter often are only represented by few cells and are therefore endangered to 'lose height', getting replaced by the mean of the neighbouring cells. Therefore, I think I'd need an algorithm that recognizes isolated top points and lowest points.

Any solution using Python, QGIS & corresponding software or even ArcGIS Desktop would be suitable. Open source is generally preferred..


texts: Computing MNDWI for Sentinel 2 Images in Google Earth Engine?

texts:
MNDWI is the normalized difference of Green Band and SWIR Band. In LANDSAT both are of a similar resolution i.e. 30m.

In the case of LANDSAT

var landsat8 = ee.Image("LANDSAT/LC08/C01/T1_RT_TOA/LC08_148044_20130425")
var ndwi = image.normalizedDifference(['B3', 'B6']).rename('MNDWI')

gives us MNDWI Band

whereas in the case of Sentinel 2 the resolutions of Green Band (B3 is 10m) and SWIR Band (B11 is 20m) are different.

How to change SWIR Band to more finer resolution so that I can apply normalized difference?


texts: static external HTML table renders smaller than expected in QGIS

texts:
A colleague gave me a QGIS layout including a static html frame. The HTML is exported from MS Word (don't ask). It contains of a table (with shading), and some lines of text outside of the table.
It should have 17 cm width, i.e. both text and table should have 17 cm. Inserting the frame sets it to 18 cm, but the table width is downscaled to 16 cm. This also affects text size, which should be 10 pt.
[
I tried simplifying the Word-derived HTML file, to no effect: tables are always not using the full width of the HTML frame.
I have also looked to find some option in QGIS which could be forcing the mismatch/scaling, but wasn't able to find one. I think I might just be using the wrong search terms to find a solution on the web, including SE.
What should I try or search for?
Just FRT, the weird blank space in the text below the table are intended.


texts: Many to one relationship or join in QGIS

texts:
I have a road network layer that has a unique id "ID_2" for each road segment. I also have a table that contains the route for each bus network. Each row in this table contains the name of the bus and the road that the bus travels in the form of "ID_2".

I want to join these two in such a way that my road network layer must contain all the bus route numbers for each road segment.

Is there a way to do this in QGIS?

Edit: I have tried the following procedure already

Using relations: With parent as my road network layer and child as my bus route table and ID_2 as the relation. But I don't see it when I select the road network layer information.

Using join: This gives me only one route number for every road segment while actually each road segment could cater to multiple routes.


texts: How to programatically generate my own wms file?

texts:
I'm new to OpenLayers and as far as I understood layers are retrieved from an wms file. I know this is not precisely a development forum but I was just wondering if there is a way to fool this framework and make it read a dynamically generated wms with a programming language. I mean, instead of giving it the URL of the wms, make it load a PHP with certain parameters and make the map display the data I want. Would it be the same? Is there a way to build a wms file? Where can I find a good tutorial to do this?

My project involves a map like this one: http://maps.nrel.gov/biomass.

Look how different colors are painted over the map depending on the product I choose on the left. This is exactly what I'm trying to achieve. I have the data in a relational database, not in a wms that's why I'm asking if there's a way I can "make" my own wms or make OpenLayers framework read my data from it to display it on the map. Thanks


texts: javascript .match to check for right coordinate format

texts:
I am using javascript .match to check and correct right coordinate format after user writes it in textbox.

Working example:
input: lat=57.000000, long=24.000000
code: match(/(\d{2}).(\d{4,9})/gi)
returns: 57.000000,24.000000

How to write .match algoritme for MGRS coordinate:
input: 35VLD6787561669, mgrs=35VLD6787561669, etc.
code: ?
returns: 35VLD6787561669 (sequence: 2 numbers, 3 letters, 10 numbers)


texts: Set Metadata domain with GDAL (python bindings)

texts:
I'm trying to set geolocation metadata in a gdal vrt with python bindings.

ds = gdal.Open('foo')
ds.SetMetadata ({'X_BAND' : '1'})
gdal.BuildVRT('bar', ds)

Results in a vrt with the following metadata

<Metadata>
<MDI key="X_BAND">1</MDI>
</Metadata>

However, I need the metadata to read:

<Metadata domain="GEOLOCATION">
<MDI key="X_BAND">1</MDI>
</Metadata>

I'm able to accomplish this using ds.SetMetaDataItem("X_BAND", "1", "geolocation") but would like to know how to do it with gdal.SetMetadata

I know the problem is with my failure to understand how swig is translating the python to the C++ so hopefully someone can quickly fix my syntax.
From http://gdal.org/python/

SetMetadata(self, *args)
SetMetadata(self, char papszMetadata, char pszDomain = "") -> CPLErr
SetMetadata(self, char pszMetadataString, char pszDomain = "") -> CPLErr


texts: Why is there a black hole in the Heatmap?

texts:
After running the Heatmap tool in QGIS I like to color the geoTiff according to my wishes. Using Psuedo color works fine, but it leaves me with no options, as type of color scheme or transparancy on certain values, however chosing color-map instead of psuedo creates a black hole  hole right in the middle, where the intensity (heat) should be strongest.

Doesn't matter how I classify, 5, 10 or 20 classes...


texts: Fixing address on TeleAtlas, how long before Google Maps has this information?

texts:
I recently submitted information to TeleAtlas and they told me it was accepted and has been added.

When I search on Google Maps the address still does not appear.

Does anyone know how long Google Maps take to update with new information or it just a matter of TeleAtlas making the data available for Google?


texts:
I wanted to use the PostGIS functionality raster2pgsql, but I had to discover that the Python function was not installed with my PostGIS installation.

Do you know if raster2pgsql is generally part of the PostGIS installation or do I have to install a special package?

I installed PostGreSQL and PostGIS on my mac via homebrew.
> brew upgrade postgis tells my that the latest version is installed, which is obviously not the case (1.5.3 installed, but 1.5.5 available).

The complete installation is

select PostGIS_full_version();
postgis_full_version
---------------------------------------
POSTGIS="1.5.3" GEOS="3.3.2-CAPI-1.7.2" PROJ="Rel. 4.7.1, 23 September 2009" LIBXML="2.7.3" USE_STATS


texts:
How can I get two different raster layers in ArcMap, that have exact the same edge?

I tried clipping or extract by mask but the output looks exactly like the input.

I only want to have these values that have a cell in the other raster.

In other words: I have a red and a black raster. In the first image I want to get rid of the black cells of the underlying raster. In the second image I want to get rid of the red cells of the underlying raster.


texts: Calculating Angle between two GPS Coordinates in QGIS

texts:
i'm looking to create an angle attribute inside a layer of points that has a starting coordinate and a finishing coordinate. The angle needs to be 0-360.

Is there a field calculator command that I could use to calculate this?


texts: How is the data access cursor performance so enhanced compared to previous versions?

texts:
The data access module was introduced with ArcGIS version 10.1.  ESRI describes the data access module as follows (source):

  The data access module, arcpy.da, is a Python module for working with
  data. It allows control of the edit session, edit operation, improved
  cursor support (including faster performance), functions for
  converting tables and feature classes to and from NumPy arrays, and
  support for versioning, replicas, domains, and subtypes workflows.

However, there is very little information regarding why cursor performance is so improved over the previous generation of cursors.

The attached figure shows the results of a benchmark test on the new da method UpdateCursor versus the old UpdateCursor method.  Essentially, the script performs the following workflow:

Create random points (10, 100, 1000, 10000, 100000)
Randomly sample from a normal distribution and add value to a new
column in the random points attribute table with a cursor
Run 5 iterations of each random point scenario for both the new and old UpdateCursor methods and write the mean value to lists
Plot the results

What is going on behind the scenes with the da update cursor to improve the cursor performance to the degree shown in the figure?

import arcpy, os, numpy, time
arcpy.env.overwriteOutput = True

outws = r'C:\temp'
fc = os.path.join(outws, 'randomPoints.shp')

iterations = [10, 100, 1000, 10000, 100000]
old = []
new = []

meanOld = []
meanNew = []

for x in iterations:
    arcpy.CreateRandomPoints_management(outws, 'randomPoints', '', '', x)
    arcpy.AddField_management(fc, 'randFloat', 'FLOAT')

    for y in range(5):

        # Old method ArcGIS 10.0 and earlier
        start = time.clock()

        rows = arcpy.UpdateCursor(fc)

        for row in rows:
            # generate random float from normal distribution
            s = float(numpy.random.normal(100, 10, 1))
            row.randFloat = s
            rows.updateRow(row)

        del row, rows

        end = time.clock()
        total = end - start
        old.append(total)

        del start, end, total

        # New method 10.1 and later
        start = time.clock()

        with arcpy.da.UpdateCursor(fc, ['randFloat']) as cursor:
            for row in cursor:
                # generate random float from normal distribution
                s = float(numpy.random.normal(100, 10, 1))
                row[0] = s
                cursor.updateRow(row)

        end = time.clock()
        total = end - start
        new.append(total)
        del start, end, total
    meanOld.append(round(numpy.mean(old),4))
    meanNew.append(round(numpy.mean(new),4))

#######################
# plot the results

import matplotlib.pyplot as plt
plt.plot(iterations, meanNew, label = 'New (da)')
plt.plot(iterations, meanOld, label = 'Old')
plt.title('arcpy.da.UpdateCursor -vs- arcpy.UpdateCursor')
plt.xlabel('Random Points')
plt.ylabel('Time (minutes)')
plt.legend(loc = 2)
plt.show()


texts: Leaflet : Auto-creations of markers on Polyline

texts:
I am new to to Leaflet. I am drawing polylines on map but markers are pinned automatically (unwanted). Please suggest me that how can avoid auto-creation of marker ()

I was using chrome-Javascript debugger to find root-cause & I observed that when .addTo(map) is getting triggered, leaflet.js is placing markers by calling _initIcon()

My function

which is drawing Polyline on map

function Leaflet_Overpass(map,osmDataAsXml,filterpane)
    {
            var resultAsGeojson = osmtogeojson(osmDataAsXml);
            var resultLayer = L.geoJson(resultAsGeojson,
                    {
                        pane: filterpane,
                        style: function (feature)
                        {
                            switch(filterpane)
                            {
                            case 'entity1': return {color: "#0000ff", "weight": 10, "opacity": 0.8};break;
                            case 'entity2': return {color: "#ff0000", "weight": 7, "opacity": 0.7};break;
                            }

                        },
                        filter: function (feature, layer)
                        {
                            var isPolygon = (feature.geometry) && (feature.geometry.type !== undefined) && (feature.geometry.type === "Polygon");
                            if (isPolygon)
                            {
                                feature.geometry.type = "Point";
                                var polygonCenter = L.latLngBounds(feature.geometry.coordinates[0]).getCenter();
                                feature.geometry.coordinates = [ polygonCenter.lat, polygonCenter.lng ];
                            }
                            return true;
                        },
                        onEachFeature: function (feature, layer)
                        {
                            var popupContent = "
"+ feature.properties.tags.name + "
"+feature.properties.type + " : "+feature.properties.id+"";
                            layer.bindPopup(popupContent);
                        }
                    }).addTo(map);
            return;
    }

Leflet.js Code

which is placing marker

_initIcon: function() {
            var t = this.options
              , i = "leaflet-zoom-" + (this._zoomAnimated ? "animated" : "hide")
              , e = t.icon.createIcon(this._icon)
              , n = !1;
            e !== this._icon && (this._icon && this._removeIcon(),
            n = !0,
            t.title && (e.title = t.title),
            "IMG" === e.tagName && (e.alt = t.alt || "")),
            pt(e, i),
            t.keyboard && (e.tabIndex = "0"),
            this._icon = e,
            t.riseOnHover && this.on({
                mouseover: this._bringToFront,
                mouseout: this._resetZIndex
            });
            var o = t.icon.createShadow(this._shadow)
              , s = !1;
            o !== this._shadow && (this._removeShadow(),
            s = !0),
            o && (pt(o, i),
            o.alt = ""),
            this._shadow = o,
            t.opacity < 1 && this._updateOpacity(),
            n && this.getPane().appendChild(this._icon),
            this._initInteraction(),
            o && s && this.getPane("shadowPane").appendChild(this._shadow)
        },


texts: How to get schema locks in oracle spatial direct connect mode?

texts:
I have a geodatabase having relation classes, domains ,subtypes etc. When I am trying to open that featureclass in a multithreaded environment it often fails with error like
"Attempting to read proptected memory" or AccessViolationException

I wonder if this is because of geodatabase, may be it is getting hold of all the schema locks, I making a direct connect to oracle spatial 10g and no ArcSDE is installed.

Any Suggestions?


texts: Understanding GTFS static data relationships?

texts:
I'm trying to build a data store for the GFTS static data, but can't get my head around the relationships between Route, Stop, Stop_Times and Trips.

e.g
Route 2E runs between airport and city centre.
It has a fixed number of stops on this route.
An instance of the 2E is a trip which leaves the airport at a certain time (9:30) to travel to the city.  Its Stop times are the times that this trip specifically will be at designated stops.

However, this implies a link between Route and Stops.  Route 2E has 12 stops between Airport and city.  This never changes regardless of whether you catch the 9:30 trip or 11:30 trip.

So Stop Times are linked to Trip.  But stops should linked to Routes.  But every data model I have seen does not show the route stop relationship.

What am I missing ?


texts: What tools in Python are available for doing great circle distance + line creation?

texts:
I need to use Python to create a great circle distance -- both a number, and preferably some kind of 'curve' that I can use to draw in a client-side map. I don't care about the format of the curve -- be it WKT, or a set of pairs of coordinates -- but just want to get the data out.

What tools are there out there? What should I use?


texts: How do I set attributes when creating vectors using OpenLayers and WFS-T?

texts:
I have published a WFS-T layer (using geoserver) and am able to add, edit and delete features using OpenLayers and it's Vector layer.

But, I don't know how can I populate additional attributes using OpenLayers and WFS-T. Any pointers?

Let's say I have database structure:

ARMY (
ID NUMBER PRIMARY KEY,
NAME VARCHAR2(200));

INVADED_AREA (
ID NUMBER PRIMARY KEY,
INVADED_BY_ARMY_ID NUMBER NOT NULL,
AREA_GEOMETRY SDO_GEOMETRY,
CONSTRAINT ia_fk FOREIGN KEY(INVADED_BY_ARMY_ID) REFERENCES ARMY(ID));

and I have a web application where you first select an army (with e.g. ID 42) and then start drawing geometries of invaded areas using OpenLayers. How do I get OpenLayers to put value 42 into INVADED_AREA.INVADED_BY_ARMY_ID column?

EDIT:
I have found http://dev4.mapgears.com/bdga/bdgaWFS-T.html# which seems to do what I am looking for. Will update when I have investigated it more.

Also found bunch of related gis.stackexchange questions. I will look through them more carefully, at first glance none had an easy answer:

I want to edit corresponding postgresql data of a vector from openlayers- how to do it?
How to edit Attributes on WFS Layer through Popup?
How to edit feature attributes with openlayers?
Modify WFS with OpenLayers
adding non spatial attributes using geoserver
How to insert non-spatial data through WFS-T insert?

And from openlayers-dev: http://lists.osgeo.org/pipermail/openlayers-dev/2007-April/000520.html


texts: Layer Symbology Based on both Category and Quantity?

texts:
I am running ArcGIS 10.2.2 Version 10.2.2.3552, using ArcMap.

Here's my situation/what I would like to do:
I have a layer that has sewer lines. Some are main lines, some are laterals. I would like to display the pipe using a color ramp based on a Diameter column, but ONLY if it is a sewer main (there is a column that says if it is a lateral or main)

i.e.
If object is Sewer Main -> Use color Ramp depending on Diameter
If object is Sewer Lateral -> Use grayed out / Other symbology.

Is there any way to do this that avoids duplicating the objects of the layer that I want to display?


texts: Imitating symbology for DEM / Hillshading in ArcMap?

texts:
I am using ArcMap 10.2 and I am looking to copy the basemap styling from this map but I'm not sure how to go about doing that.

I have a DEM of the area and the hillshading layer used but I'm not sure how to get the DEM to be shaded from green to tan based on the elevation. It's not an option in my symbology for raster dater.


texts:
I've got an issue I am using GeoServer of boundlessSuite and when I tried to load a layer in OpenLayers 3 it doesn't show on the map. I found the problem but can't fix it.

On my GeoServer when I use the preview to see my layer this is the URL:

  http://localhost:8080/geoserver/test1/wms?service=WMS&version=1.1.0&request=GetMap&layers=test1:BoiteOptiqueTest&styles=&bbox=514092.0239487532,1977542.5898263564,522107.6474402136,1981434.1000406644&width=768&height=372&srs=EPSG:32620&format=text%2Fhtml%3B%20subtype%3Dopenlayers

but when I inspect element on Chrome this is what it tries to load:

  http://localhost:8080/geoserver/test1/wms?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&FORMAT=image%2Fpng&TRANSPARENT=true&LAYERS=test1%3ABoiteOptiqueTest&WIDTH=256&HEIGHT=256&CRS=EPSG%3A3857&STYLES=&BBOX=-10018754.171394622%2C2504688.542848654%2C-7514065.628545966%2C5009377.08569731

How can I fix that?

Here is my code:

  var layers = [
       new ol.layer.Tile({
         source: new ol.source.OSM()
       }),
       new ol.layer.Tile({
         extent: [-13884991, 2870341, -7455066, 6338219],
         source: new ol.source.TileWMS({
           url: 'http://localhost:8080/geoserver/test1/wms?',
           params: {'LAYERS': 'test1:BoiteOptiqueTest'},
           serverType: 'geoserver',
           isBaseLayer: false,
           projection: 'EPSG:4326'

         })
       })
    ];
  var map = new ol.Map({
    layers: layers,
    target: 'map',
    view: new ol.View({
      projection: 'EPSG:4326',
      center: [0, 0],
      zoom: 2
    })
  });

I tried using different projections like EPSG::32620, EPSG::900913 but nothing works.


texts: Error Running ArcPy Script- "Could not open the address table. Failed to execute (GeocodeAddresses). "

texts:
I have a script I am attempting to debug so I can automate this task I have to do a weekly basis every monday morning which is tedious. I am very new to python and coding in general.

This is a very simple script, I am given a excel sheet with data (addresses) and I have to geocode them. I then must clip out the points to only a specific area (polygon) and then I select a specific 'series' which is a field type for its attribute and have that become its own layer so I can manipulate its symbology uniquely. Below is my script.

All the paths are correct however the "Y" drive is a network drive shared between the entire organization. When I first came across this error my table was located there so I moved it locally on my desktop and the error persists.

The table is not open or being used by any other program by myself or anyone else so there should not be a file lock on it.

Anyone have any suggestions?

My code is below.

# Import arcpy module
import arcpy
import os

#Change user "jhead" to your local group username (amiller), this creates a temp folder on your desktop
#because the path for geocoding will not save outputbecause of spaces in path
newpath = r'C:\Users\jhead\Desktop\Temp'
if not os.path.exists(newpath):
    os.makedirs(newpath)

# Local variables:
Sheet1_ = "Sheet1$"
LakeWorth_Addressing_Parcels = "Y:\\Department for Community Sustainability\\GIS Division\\Internal Division Files\\GIS Addressing Tools\\LakeWorth_Addressing_Parcels"
table1_xlsx = "C:\Users\jhead\Desktop\table1.xlsx"
Geocode1_shp = "Y:\\Department for Community Sustainability\\GIS Division\\Internal Division Files\\Division Work Folders\\ABT\\New_project\\TempData\\Geocode1.shp"
LakeWorth_DBO_Municipal_Boundary_Poly_2016 = "LakeWorth.DBO.Municipal_Boundary_Poly_2016"
clip1_shp = "Y:\\Department for Community Sustainability\\GIS Division\\Internal Division Files\\Division Work Folders\\ABT\\New_project\\TempData\\clip1.shp"
liquor_select_shp = "Y:\\Department for Community Sustainability\\GIS Division\\Internal Division Files\\Division Work Folders\\ABT\\New_project\\TempData\\liquor_select.shp"

# Process: Geocode Addresses
arcpy.GeocodeAddresses_geocoding(Sheet1_, LakeWorth_Addressing_Parcels, "Street ADDRESS_LINE1 VISIBLE NONE;City City VISIBLE NONE;ZIP ZIP VISIBLE NONE", Geocode1_shp, "STATIC")

# Process: Clip
arcpy.Clip_analysis(Geocode1_shp, LakeWorth_DBO_Municipal_Boundary_Poly_2016, clip1_shp, "")

# Process: Select
arcpy.Select_analysis(clip1_shp, liquor_select_shp, "\"SERIES\" = '3PS'")

The error I am recieving is:
**

Runtime error  Traceback (most recent call last):   File "", line 29, in    File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\geocoding.py", line 212, in GeocodeAddresses     raise e ExecuteError: ERROR 000006: Could not open the address table. Failed to execute (GeocodeAddresses).

**
I would like to add that when I perform this tool manually (toolbox and select the tool and add the inputs manually) everything works perfectly. it grabs the tables and geocodes all the addresses no issues. So not sure why this is happening.

I think maybe its due to this geocode tool having multiple inputs, sheet, the table and the addressing locator, but in the script its missing the table inut and the model shows a dotted line connected rather solid line.


texts: show all points of carto-point-layer at low zoom-level

texts:
I have a Leaflet-map with a carto-layer(points), using the createLayer-method.  It seems that Carto is doing something internally, that not all points are shown at low zoom-levels. Each time zooming in, more points are shown.
How do i get control of this function, or, respectively how do i ask carto to show all points at low zoom-levels, regardless of collisions and overlapping?


texts: WGS84 to EPSG:20007, Pulkovo 1995/ Gauss Kruger zone 7

texts:
I have a shapefile drawn in Selected CRS EPSG 84, when I try to load it in my project in Pulkovo 1995/ Gauss Kruger zone 7, its missplaced, and when i change shapefiles crs in properties, still no change.

im using qgis 2.18.3


texts: Drawing border using R in shapefile based on dummy variable?

texts:
For an analysis in R I want to create a "new" border from a shapefile. The Shapefile comes with a dummy variable (zone_dummy) indicating the regions laying east and west from the border (see example below). Now I wanted to create a geomatric path along the border of the green and red municipalities:

Here is my shapefile and below an example code:

rm(list = ls(all = TRUE))

library(maptools)
library(RColorBrewer)

area <- readShapePoly("VZ250_GEM.shp")

kreise.data<-data.frame(area@data)

zone<-data.frame(c("Weiden i.d. OPf.", "Bad Kissingen", "Schweinfurt", "    Mellrichstadt", "Brückenau", "Rhön-Grabfeld",
               "Hofheim", "Ebern",  "Haßberge", "Coburg", "Hof", "Selb", "Kulmbach",
               "Marktredwitz", "Bayreuth", "Bamberg","Staffelstein", "Lichtenfels", "Kronach",
               "Stadtsteinach", "Naila", "Münchberg", "Rehau", "Wunsiedel", "Weiden",
               "Tirschenreuth", "Kemnath", "Neustadt a.d. Waldnaab", "Bad Neustadt a.d. Saale", "Neustadt b. Coburg", "Vohenstrauß", "Nabburg",
               "Overviechtach", "Waldmünchen", "Cham", "Roding",
               "Deggendorf", "Passau", "Kötzling" ,"Viechtach", "Regen", "Bogen", "Grafenau",
               "Wolfstein", "Wegscheid", "Wunsiedel i. Fichtelgebirge", "Schwandorf", "Freyung-Grafenau"
               ))
#create dummy to mark municipalities in zone
colnames(zone)<-c("GEN_K")
zone$zone_dummy<-1
#merge with original data
kreise.data.zone<-merge(kreise.data, zone, by="GEN_K", all=T)
## Replace NA's with 0
kreise.data.zone[is.na(kreise.data.zone)] <- 0
##Merge assembled Shape-file Data Frame to Original Shapefile
area.zone<-merge(area, kreise.data.zone)
#we only care about bavaria:
bayern<-subset(area.zone, area.zone@data$NUTS1_CODE=="DE2")
#plot a map of bavaria with zone:
my.palette <- brewer.pal(n = 3, name = "Set1")
spplot(bayern, "zone_dummy", cuts = 1 , col.regions = my.palette, col="white", colorkey=FALSE)

My question is now:
How can I create a border along the zone (=sharp border marking the path along the red and green municipalities)?

The idea is to (1) create this border in the shapefile (2) estimate the distance of each municipality in Bavaria to this border.


texts: In what units do area calculation in geographic coordinate system return?

texts:
Is it correct to call it degrees square?

Will it be the same if calculate the area using Plate Carrée projection?


texts: Dissolve (merge) polygons QGIS based on attribute values only when are adjacent

texts:
I am a fairly new user of QGIS having a question regarding dissolve functionality.

I am working with the NUTS-2 shapefiles provided by EU, which I have imported in QGIS. I have created a new layer with only Denmark and Netherlands as you can see in the picture below.

The corresponding attribute table is the following.

The fields describe the ID and name for each NUTS node, as well as the country code. I have introduced a new attribute field, called Merge_attr.

I have managed to dissolve the polygons with the same attribute value of Merge_attr. However, I would like to dissolve polygons based on the Merge_attr value only when they belong to the same country and potentially only when they are adjacent.

Can we constrain dissolve functionality to only polygons that are adjacent?

This means that the following should not take place (yellow polygons got merged even though not adjacent):


texts: Only select non null values using st_value on Postgis tiled raster

texts:
I have a table tmax_2016 that contains tiled raster data.

I believe the following query retrieves the raster value at a point on the given rast_date:

SELECT ST_Value(rast, ST_SetSRID(ST_MakePoint(-83.54005,31.59668),4269)) AS tmax FROM tmax_2016 WHERE rast_date='2016-09-28';

However since it's tiled, I always get 1 valid result and a bunch of nulls for the tiles that don't contain the asked for lat/long. How do I modify the query to not return the nulls?

The following pukes I believe because  ST_Value is an agregate function?

SELECT ST_Value(rast, ST_SetSRID(ST_MakePoint(-83.54005,31.59668),4269)) AS tmax FROM tmax_2016 WHERE rast_date='2016-09-28' WHERE tmax IS NOT NULL;

I could get rid of the nulls using a nested select but is there a better way?


texts: Batch processing 300 FBX files with textures using FME?

texts:
My goal is to provide a contractor with FBX files with textures with a file size that is roughly 10MB.

I have esri ASC files for my Area Of Interest(AOI) as well as a recent aerial capture. My AOI fits within 300 500m square tiles.

I had thought that the best way of doing this would be to use the RasterTiler and RasterPropertyExtractor transformers to give me files named with the bottom left of the tile so that the file name could be used to merge the imagery with the surface. Here is my file name expression e@Value(_lower_left_x)n@Value(_lower_left_y)

So I now have 300 files of each ASC, JPG and WLD files for the JPGs.

I can create the desired output using a TINGenerator and AppearanceSetter followed by a FBX writer.

I had to manually enter a matching JPG file to the JPG reader and fme_basename value seems to get lost passing through the TINGenerator. So it writes files with blank names ".fbx" and ".fbm" dir for texture.

I feel I should be able to use the FilenamePartExtractor transformer to pass _rootname to the second reader, but I haven't been able to make this happen. It might be easier to point the reader at the ECW mosaic of the entire AOI and pass it the geom of the footprint of the ASC tile to be clipped

I have found this article on batch processing which I intend to follow once I can get the above process working.


texts:
I have a point data as observer point with field attributes for Spot, offsetA, offsetB, Radius1, Radius2, Vert1, Vert2, Azimuth1 and Azimuth 2. For the raster file- a dsm with elevation data. They have the same projection and coordinate system

But when I run the viewshed analysis I get this error:

ERROR 010067: Error in executing grid expression. No valid Observation points present.
Failed to execute (Viewshed).

Any suggestions how to solve this?


texts: Simple javascript viewer for ArcGIS Server data

texts:
I would like a simple javascript viewer which will allow end users to do the following in their browser:

Add layers from arcgis server by clicking an "add data" button
See attributes with an info tool or table veiw
Measure distances and areas
Draw points, lines, polygons
Save the state of their map for viewing later

Is there a javascript template for ESRI javascript API that can do this?
Basically I want a local copy of the ArcGIS online viewer for users to access our internal ArcGIS Server data.
Everything I have found references ArcGIS Online, and not a local Arcgis Server.

Thanks!


texts:
I'm very interested by understanding how Contraction Hierarchies algorithm works. I have found this page : https://www.mjt.me.uk/posts/contraction-hierarchies/ and read a lot of stuff. I understood how the algorithm goes and work except for one part. I didn't grasp how the contraction strategy really work. (Important: I'm not fluent in Maths).

In the first example on the page given above, the contraction order is 6→8→3→0→5→7→4→1→2→9 and I can't understand why. The explanations are not clear to me and the original article contains too much maths in it.

Can someone, please, explain the strategy used to define the contraction order?


texts: How to change field defined as REQUIRED to NON_REQUIRED in ArcGIS Geodatabase

texts:
I have created some feature classes in a geodatabase and also some fields using this sample code:

arcpy.AddField_management("schools", "ref_ID", "LONG", 9, "", "", "refcode", "NULLABLE", "REQUIRED")

Since the fields where created with field_is_required's parameter as "REQUIRED", how to change it to "NON_REQUIRED" without re-create it?


texts: Calculating only values that are not null in field using ArcGIS Pro

texts:
I'm using ArcGIS Pro, the field calculator. I want to populate a hydrant layer with hydrant inspection data.
I have two fields from a joined table: field 1 and field 2.  I want to populate field 1 with the values from field 2, but only if the values ARE NOT Null. Basically I want to preserve the values in field 1 if the corresponding value in field 2 is Null.  I've been trying slight variations of the python code below.  But it just populates field 1 with all of field 2 whether there are nulls or not.
I'm a beginner when it comes to Python.
Codeblock:
def reclass(f1,f2):
        if f2 is not None:
            return f2
        else:
            pass

Expression:
reclass(!field1!,!field2!)


texts:
I have to tab some panels:

from PyQt4.QtGui import *
from PyQt4.QtCore import *

dws =  iface.mainWindow().findChildren(QDockWidget)
parent = iface.mainWindow().findChild(QDockWidget, 'Layers')
dockstate = iface.mainWindow().dockWidgetArea(parent)
for d in dws:
    if d is not parent:
        if iface.mainWindow().dockWidgetArea(d) == dockstate and d.isHidden()==False:
            print d.objectName()
            iface.mainWindow().tabifyDockWidget(parent,d)

How can I activate the parent?


texts: Remote sensing accuracy assessment - McNemar`s Test

texts:
I would like to know how you determine the f12, f21 parameters for use in the McNemar`s test when for example I have two confusion matrices such as below:

(where f12 = number of cells misclassified by A but not by B
And f21 = number of cells misclassified by B but not by A)

In the example with the two confusion matrices from above, I have the totals but am having difficulty in correctly determining the matrix elements f11, f12, f21, f22 below.


texts: I've imported a number of files that used to be in KML format into geoserver after putting everything into a PostGIS table. Now I would like to use CSS style information that was associated with the original KML files, but  Geoserver says the layers cannot be rendered because they are missing the attributes being called for. Is there any way I can add these attributes to the imported layers?


texts: Cannot bind Legend layer to Polygon layer in Leaflet/R

texts:
I try to plot a map of concentrations of chemicals. These have very different ranges, so a shared legend is not feasible. The ultimate goal is to have the different layers for the different chemicals selectable and show the respective legend.

From what I read in the documentation of addLegend(), exactly this should be handled by the group parameter that was set to the same name in the addPolygons() for that layer before. However, it just does not work. Does anyone have an idea for the reason? Here's what I tried:

library(spData)
library(leaflet)
library(sf)
library(tidyverse)

set.seed(1)

vietnam <- world %>%
  filter(name_long == "Vietnam")

grids <-
  st_make_grid(vietnam, square = FALSE)

grids_sf <- st_sf(
  id = 1:length(grids),
  geometry = grids,
  category = sample(letters[1:3], size = length(grids), replace = TRUE),
  value = rnorm(length(grids))
)

addPolygonsLayers <- function(map, data, layer_by, fill_by) {

  data_groups <- split(data, data[[layer_by]])
  data_groups_names <- names(data_groups)

  for (i in seq_along(data_groups)) {
    group_data <- data_groups[[i]]
    group_name <- data_groups_names[i]
    group_values <- pull(group_data, fill_by)

    pal <- colorNumeric(palette = "viridis", domain = group_values)

    map <-
      map %>%
      addPolygons(
        data = group_data,
        fillColor = ~ pal(group_values),
        fillOpacity = 1,
        group = group_name,
        popup = str_c(group_name, " </br> Value: ", group_values)
      ) %>%
      addLegend(
        group = group_name,
        opacity = 1,
        pal = pal,
        values = group_values,
        position = "bottomright"
      )
  }
  map <-
    map %>%
    addLayersControl(
      baseGroups = data_groups_names,
      options = layersControlOptions(collapsed = FALSE),
      position = "topright"
    )
  return(map)
}

leaflet() %>%
  addTiles() %>%
  addPolygonsLayers(
    data = grids_sf,
    layer_by = "category",
    fill_by = "value"
  )


texts: Clipping rasters to mask layer QGIS 2.18, 3.4, 3.6 not working

texts:
Having a problem clipping a singleband GEOTIFF of Canadian land uses to a vector Quebec provincial boundary.  When I do Raster -> Extraction -> Clip by mask layer, the program processes endlessly.  I never get an error message; it just doesn't stop. Out of curiosity, I let the process run for more than 24 hours in QGIS 3.4 just to see what would happen.  Never got a result.  I've also tried in QGIS 2.18 and 3.6. No luck. I assembled the raster coverage from five separate files that were each in their own UTM zone. I merged and reprojected to CRS 4269 (NAD 83). The resulting raster is quite large (~15 gigs), but I'm not sure that's the problem.  I made sure the vector is in the same CRS, so I don't believe that's the problem either.  I thought at first it might be some problem with my files, but I brought them to school and was able to clip the raster in ArcMAP.  It took a while to process -- maybe 20 minutes or so -- but that worked.  I've found a few mentions of similar-but not exactly the same-problems in other forums online.  Is this a known issue with QGIS or user error?  Is there a good work-around besides switching to ARC?

New:


texts: How to have map project East Asia connecting to Atlanta?

texts:
Is there a projected coordinate system that enables the eastern part of Asia to connect with Alaska?
my current setup is this:I wish for the entirety of east Asia and part of Alaska to be on the same side. So the end of the left-hand side to connect to the right-hand side in one map.


texts: Are these glitches in my map rendering caused by not projecting WGS84 coordinates to Cartesian?

texts:
I'm attempting to render continental geometric data stored within an R-Tree in a SQlite database. The data is store as long/lat WGS84 and I'm trying to render it to 256px PNGs for displaying on a slippy map. Keep in mind I have no training in GIS mapping and rendering and everything I know is self taught.

When I render it using OpenGL, I'm simply remapping the long/lat coordinates of each of the polygons returned from the SQLite query into a 0px to 256px x/y coordinate space for rendering on a slippy map tile using the following code:

X = (long - extent.X1) / (extent.X2 - extent.X1) * ImageSize;
Y = (((lat - extent.Y1) / (extent.Y2 - extent.Y1) * ImageSize) *-1) + ImageSize;

For the Y coordinate, I'm switching the axis so that 0:0 is the top left of the image.

The extent information in the above code is built by taking X/Y and Zoom data from a slippy map (leaflet.js in this case) and converted into long/lat bounds using the code provided by OSM here

Thats all I do before I send the X & Y data over to OpenGL.

At the higher zoom levels (5+). The polygons render okay, heres Australia at zoom level 5 (there are 1 or 2 exceptions around the equator where there are some glitches at zoom level 5, but they seem to completely disappear at zoom levels 6+):

However, at zoom levels below 5, things start to go hay wire, here's Australia and South Asia at zoom level 3:

There are significant glitches in the rendering, (repetition and misalignment with the rendered polygons on the tiles and stretching of the continents) And it only gets worse as I zoom out more. Here's part of the map at zoom level 2:

It seemed weird to me that everything would render well at the higher zoom levels but fail so badly at the lower ones. But then it occurred to me that the closer to earth I go, the flatter (and more 2D like) it becomes, so I'm now wondering if the reason for the major issues I'm having is because I'm not projecting to a Cartesian coordinate system first.

Can anyone confirm whether this is the case and if so, what should I be using to convert the coordinates from WGS84 to Cartesian X/Y?


texts: Why does geoserver layer disappear when baselayer is false?

texts:
I am trying to overlay a layer from Geoserver on an OpenStreetMap base layer.  If I load the Geoserver layer as a base layer, I can see it - but if I set baselayer to false it disappears.

My code is:

<!DOCTYPE HTML>

    <title>OpenLayers Demo</title>

        Shows the basic use of openlayers using a WMS layer

I've used firebug to investigate the calls being made to geoserver.  In the case where I add the WMS as a base layer I see the following:

whereas when baselayer is false it's

The first call returns a PNG showing my data, the second one is empty.

A final point is that my data in Geoserver has native SRS EPSG:4258 and declared SRS EPSG:4326 - I have set the "Reproject native to declared" option.


texts:
I have two tables of polygons: Table_1and Table_2. All the polygons of T1
are inside the polygons of T2.

I want to SELECT the id of T1 and the name of the polygon from T2 where the polygon is located.

I am using a spatial join :

SELECT T1.id, T2.name
FROM Table_1 AS T1
JOIN Table_2 AS T2
ON ST_Contains(T1.geom, T2.geom)

The trouble I have is that the query does not return any value. I don't understand why. I have already used this method for other queries and it worked. Having already done this doesn't help for this case.

UPDATE 1:

I tried with ST_Intersects and nothing was returned neither.

UPDATE 2:

I checked the SRID for the Table_2, none was defined. I corrected it but the results are still the same.

UPDATE 3:
Here is the geometry of the first row for Table_1

1; 01060000206A080000010000000103000000010000001F000000768C2BD8D7392441EBB09D7A7D1C5B410A497D15F03A24410BBC0F716C1C5B41366DF983803B24413246D86B5D1C5B4129CFBCC2DB3B2441D5616275491C5B416B1A9FE6863C244171AEA1B83D1C5B4157485EC7673D2441AF20131C3F1C5B41B76046F83A3E2441B02ED7F3391C5B413C64BBBEE93D244106E8ABEC1D1C5B414CF08D90B33D2441ECAE432D0A1C5B41C01F2375013E24414B069F9AF41B5B41B1E53CDA463D2441E057FA4EDF1B5B41987E1A55C83C244135ADDCA1DC1B5B417D513EA5E43A24410EF6856DBC1B5B41E6FEDBFB583A24414C69E48CB71B5B419B448F6F523A24417D139ABEB61B5B41D614CBE8DE382441DE529054A31B5B41FC1348DEF3362441A8D13B6B891B5B41CCC3BB8973352441194EE3937C1B5B412BBE50BBC9342441F232372E9C1B5B419782A0B82E3424419E9C1FE2D21B5B41F81E12BDAD33244166FE8D8F041C5B41FE3278D9AF3324412FDBB1FB231C5B4175D7441B2D332441878F829C391C5B4114C0D422DA3324415FA552FC491C5B412D03A336AE342441BEB77102521C5B41406F4356EC35244127E1BC3B5D1C5B410ACB77C47D362441BDC9A9145D1C5B4123FC9AA4A13724419EDFE6046A1C5B4134D3FF7E76382441B9F5A0997D1C5B41CA43370F16392441FCB951AD8A1C5B41768C2BD8D7392441EBB09D7A7D1C5B41

The geometries of my Table_2 have to many characters to fit in the update. The SRID is 2154.


texts: Implement Rtree index to custom objects list

texts:
I'm trying to improve a spatial calculation, and I think the most efficient way is to add a spatial index, but I'm not sure how to do it.

In one hand I have a list of custom objects, with a geom attribute load in memory from a PostGIS table.
In the other hand, a shapely.LineString.

What I would need to do is get the segments of the line that are inside each polygon.
What I'm doing now is iterate through the Polygons list, and check if the linestring intersects with each polygon, if does, add it to a list of segments with some polygon
for buff in polygons:
    intersection = multi_line.intersection(to_shape(buff.geom))
    if intersection:
        dict_segment = dict(segment_length = intersection.length, pol_name=buff.name)
        segments.append(dict_segment)

The main problem that I've seen is that I"m iterating through all polygons list instead of iterate only in a portion to save time.
For this I'd need a spatial index, but as far as I know, I only can set a spatial index to a list of shapes (polygons in this case) but not a list of custom objects, although it has the geom property.
Does anybody know how I could set the index to this object list?


texts: Increasing click selection tolerance in ArcGIS Online

texts:
Will the 2016 solution listed at Increase FeatureLayer click tolerance work in Web App Builder (ArcGIS Online) or perhaps a newer solution exists?
We have valid feedback:

I am having a little trouble finding the right place to click to get
the info box to pop up on some projects.  Can you add a selection
feature to make that easier?

Please see  I am seeing when attempting to use the suggested code in AGOL's (within WAB) Configure Popups (does a JS to Arcade converter exist?) and the code copy / paste below.
In 3.x, you can do this by overriding the "_calculateClickTolerance" function of the map's popupManager property. For example:
map.popupManager._calculateClickTolerance = function(graphicsLayers) {
    return 10;};

That is a very simplistic example, but you could implement whatever logic you want inside that function as long as you return an integer (preferably non-negative).
In ArcGIS Desktop
under Selection -> Selection Options -> Selection tolerance
exists and wondering where the similar option is in AGOL (WAB) for configure popups.
I've read through Your Arcade Questions Answered but have not found a solution yet.


texts: NDVI from ENVI to Arcmap stretched value

texts:
I calculated NDVI in ENVI by band math. When I drag and drop the image in ArcMap it becomes a black image with stretched value =0 for all pixels. However, by identity clicking each pixel has a value equal to NDVI.
Is there any solution to solve the problem?


texts: Stereographic projection of WGS84 ellipsoid

texts:
According to this source, an optimal projection for small, roughly square areas is the stereographic projection. The wikipedia page of the stereographic projection shows how to convert from spherical coordinates to polar or xy coordinates. However, I would like to convert WGS84 data (latitude and longitude) to xy coordinates and WGS84 uses an ellipsoidal model of Earth instead of a sphere.

My questions is:

What are the formulae for a stereographic projection on an ellipsoid centred at lat0,lon0 (point tangent to the plane) projected from the antipode of lat0,lon0?

Bonus questions:

1) What is the inverse transformation?

2) Is this projection still conformal?

I have read section 3.3.1.1 of the Geomatics Guidance Note 7, part 2 "Coordinate Conversions & Transformations including Formulas" (March 2019) and while I believe it may contain the answers to questions (1) and (3), I'm not sure what does it mean by "For the ellipsoid the  parameters defining the conformal sphere at the tangent point as origin are first derived". Does it first compute the parameters for a sphere tangent to lat0,lon0? Then uses the antipode of lat0,lon0 in that sphere instead of the ellipsoid? If so, why does it do it that way? To preserve conformality? If this is true then would there be any advantage in projecting stereographically from the antipode on the ellipsoid instead of projecting from the antipode in the conformal sphere?


texts: Identify zones intersected when a polygon is selected and generate the areas of the polygon parts that falls in different zones

texts: I am consuming two (2) datasets in a JavaScript application (field and zones). I want to be able to identify which zone/zones a field falls into. In many cases, some fields falls into 2 or more zones and I want to generate the areas of the field part that falls in different zones. I've configured a pop-up showing unique ID and Total Area of each field when selected, but this additional information is required for fields on zone boundaries and I'm thinking about what I can achieve using the ArcGIS geometryEngine and geometryService. How can I achieve this using JavaScript?


texts: How to union polygons and add attribute values of combined output features

texts:
I have a layer containing 10 polygons all overlapping at various locations.

Each polygon has a value of '1' in a column called "value"

How do I Union the polygons (rows) together such that the values in the "value" column are added for each output polygon?

What I'm after is the ability to add the values from the columns of the input features such that a combined overlay is produced. What I mean is, where 8 polygons overlap within the feature, the output union geometry has a value of "8".

And yes, this can be done using raster, but is there a way to do this using vector and SQL Spatial?

I would like to to this using SQL Spatial, ie. using PostGIS or SQLite.

There are discussions on this site on how to Union polygons together, but how to combine values from the inputs is something I cannot find.


texts: Problems Renaming Fields with Make Feature Layer

texts:
I am using Make Feature Layer as a tool in a model. You can see I have two instances of the tool. I'm having issues with the top one, Make Feature Layer (2).

I am using the tool's field renaming feature in order to avoid the complicated Add Field, Calculate Field, Delete Field process that is the other option for renaming. However, the tool is not seeing the correct field map from the input layer. It is unaware of the two required fields Shape_Length and Shape_Area. Therefore, when these two fields unexpectedly show up, the renaming is all offset by two fields. For instance, here is a table showing an example of my settings and desired result compared to my actual result.

Now to form this into some kind of question. Has anyone else experienced this? How can I get the tool to see the correct field map? Or is there another way around this?


texts: Removing duplicate features with same geometry in ArcMap?

texts:
Is it possible to uniquely identify each geometry in a feature layer?

Essentially I have a polygon feature layer with many rows of data represented by the same geometry. I would like 1 row for each unique geometry. So I'm looking for some out-of-box methodology to combine attributes based on the shape of the geometry.

I come from a Postgres world and can solve this problem there, but I don't know how to do it (or if it's possible) to do it non-programmatically with ESRI tools.


texts: Searching pre-1998 Indonesia Kabupaten (Level 2) shapefile

texts:
I am trying to find some "old" (pre-1998) Indonesia boundary shapefiles at the Kabupaten and province levels.
It is important that they are pre-1998, or at least dating before the peak of the secession of the regencies. Does anybody know where I can find them?


texts: Re-compiling Python Addin after editing script?

texts:
I have been developing a few python addins and I'm just wondering what is the best way to re-compile the addin anytime I edit the script.

Currently, if I make an edit to the script, it seems like I have to

Close ArcMap
Go to My addin folder and double click the 'makeaddin' python file
Double click the 'Python Addin File
Restart ArcMap for the changes to come into effect

If I don't do these steps, then the edits aren't recognized. This gets pretty annoying and time consuming when I want to place maybe a simple print line somewhere or I do a small edit just to test a result.

Is there a better way of doing this?


texts:
I'm trying to improve performance for the query below. No matter how I write the query (subquery in FROM clause, subquery in WHERE clause) postgres insists on running all ~570K rows through the expensive ST_DWITHIN function even though there are only 60 rows where county=24. How can I get postgres to filter on county=24 BEFORE running through the postgis func which seems to me would be much faster and far more efficient? 700ms isnt cause for too much concern but as this table grows to 10M+ I'm concerned about performance.

Also to note, p.id is a primary key, p.zipcode is an fk index, z.county is a fk index, and p.geom has a GiST index.

Query:

EXPLAIN ANALYZE
  SELECT count(p.id)
  FROM point AS p
  LEFT JOIN zipcode AS z
    ON p.zipcode = z.zipcode
  WHERE z.county = 24
    AND ST_DWithin(
      p.geom,
      ST_SetSRID(ST_Point(-121.479756008715,38.563236291512),4269),
      16090.0,
      false
    )

EXPLAIN ANALYZE:

Aggregate  (cost=250851.91..250851.92 rows=1 width=4) (actual time=724.007..724.007 rows=1 loops=1)
  ->  Hash Join  (cost=152.05..250851.34 rows=228 width=4) (actual time=0.359..723.996 rows=51 loops=1)
        Hash Cond: ((p.zipcode)::text = (z.zipcode)::text)
        ->  Seq Scan on point p  (cost=0.00..250669.12 rows=7437 width=10) (actual time=0.258..723.867 rows=63 loops=1)
              Filter: (((geom)::geography && '0101000020AD10000063DF8B52B45E5EC070FB752018484340'::geography) AND ('0101000020AD10000063DF8B52B45E5EC070FB752018484340'::geography && _st_expand((geom)::geography, 16090::double precision)) AND _st_dwithin((g (...)
              Rows Removed by Filter: 557731
        ->  Hash  (cost=151.38..151.38 rows=54 width=6) (actual time=0.095..0.095 rows=54 loops=1)
              Buckets: 1024  Batches: 1  Memory Usage: 3kB
              ->  Bitmap Heap Scan on zipcode z  (cost=4.70..151.38 rows=54 width=6) (actual time=0.023..0.079 rows=54 loops=1)
                    Recheck Cond: (county = 24)
                    Heap Blocks: exact=39
                    ->  Bitmap Index Scan on fki_zipcode_county_foreign_key  (cost=0.00..4.68 rows=54 width=0) (actual time=0.016..0.016 rows=54 loops=1)
                          Index Cond: (county = 24)
Planning time: 0.504 ms
Execution time: 724.064 ms


texts:
I am going to build my on On-The-Fly projection system, just like QGIS's. A re there any C++ Libraries that will allow this?

QGIS does not to resample the raster while on-the-fly reprojection, but to distort the image.

Didn't figure out how QGIS do this after read the source code.


texts: Clip and spatial join of hydrology polygon layer with district layer using ArcPy?

texts:
I have two shapefiles I am attempting to join spatially. One is a county polygon file and the other is a hydrology polygon layer spread throughout various districts within it.

I'd like to output a shapefile and attribute table that breaks the hydrology file by district so pieces of the hydrology file are separated and can be associated with respective borders that they fall within.

What are the steps to achieve this?

I am running ArcGIS desktop 10.4. Right now I'm trying to clip the hydrology layer from the district layer. Spatial join just doesn't break the hydrology file.it gives most of the hydrology parcels to one district and although I can see which districts it covers in a list the polygon itself is unchanged. So I couldn't show which part of the lake is contained within a particular district. I'm trying to select a district, clip the hydrology file based on that one district then output a shapefile. I'll do this for all. Trying to do it with python.


texts: Split a shapefile of lines into even spaced segments

texts:
I am using ArcGIS ArcMap 10.2.

I have a shapeless of polylines and want to split them each line into 100 m segments. I can see under the editor toolbar the split tool and that can divide each line into equal parts and accepts an input distance.

My issue is that my shapefile contains 1000's of polyline, (road map of a state). Is there a way I can do each line of my shapefile in one go? I can't merge the roads together as I still require the associated attributes.


texts: Wrong Scale after Reprojecting to a Projected Coordinate System

texts:
I have some data that was recorded by our team in the field in Trimble TerraFlex.  I've been able to export that data to Shapefiles and imported them in QGIS (2.18.0) in WGS 84 (4326) and are spatially correct when I pull in a Google Satellite Image through Openlayers WGS 84/Pseudo Mercator (3857).

What I'm trying to twofold:

A) Convert this data to a properly scaled DXF.

B) Get areas of the polygon layers and lengths of the line layers.

If I just export the data as is, it comes into AutoCAD at a really small scale (like any sort of area or length measurement starts with 0.000xxx).  My understanding after doing a fair amount of research is that this is due to the data being in a GCS rather than a PCS.  QGIS is trying to export using DMS instead of meters and is giving me really tiny numbers.

So, I need to reproject to a PCS, correct?  The area in question is just south of Branson, Missouri and thus is in UTM 15N.  If I change the layer CRS to this, the data moves on screen even if I have "On the Fly Transformation" selected and I can no longer use the Google Satellite imagery to spot check my locations.  This is okay.

Unfortunately, when I export these layers to DXF with this UTM 15N PCS, the scale is still really, really tiny.  Does anyone see any glaring issues with my process or know how to troubleshoot this?  It's frustrating because I've done a good deal of work trying to understand how projections work but I still can't solve what seems like a relatively simple problem.

I can upload my project to a host and let someone work with it if that would help.


texts: Details on Topcon theodolite data format (<3595959+3595959+****d094)

texts:
I am interfacing to a Topcon DT-102 (circa 1995) digital theodolite and wonder if anyone can help with the transmission format.  It is perhaps shared with Sokkia or Leica/Wild.

Background

The unit has only the ability to measure horizontal and vertical angles.
By trial and error I have found it is transmitting at 1200 baud 7 data even parity, and continuously sends messages twice a second.

When the LDC display shows 0° 53' 10" vertical and 7° 46' 40" horizontal measurements, the unit transmits:

<0005310+0074640+****d094
 |      |       |    |
 |      |       |    + checksum?
 |      |       + distance?
 |      + horizontal measurement sdddhhmm
 + vertical measurement dddhhmm

(each message is then followed by ETX, CR, LF, ie ASCII values 3, 13, 10).

Question

Does anybody have any details on this format?

Specifically

What is this format called?
Do you have a reference for it and a detailed specification?
What does +**** represent?  Perhaps a distance for a similar model with a rangefinder?
What is d094? A checksum?  Any details?


texts: Interpreting Cost Distance Tool (ArcMap) results?

texts:
I have a raster data where pixel value of each cell represents the potential minimum travel time across the cell in old world. I need to find the minimum travel time of some points from Addis Ababa (Ethiopia) in ArcMap.

So I use Cost Distance to generate the raster with the least accumulative travel time for each cell. However, the generated raster doesn't make sense. For example, in the original raster, the source cell (i.e. including Addis Ababa) has the value of 0.263932, and its neighbor cell has the value of 0.194596, but the pixel value of the neighbor cell in the generated raster is 0.003143 !!!

I am using ArcMap 10.6. The coordinate system of the raster data is WGS_1984.

What am I doing wrong? Here is a screen shot of the tool.


texts: How do you use the Raster Calculator in ArcGIS to combine rasters using the Con statement

texts:
So, I'm busy doing an assignment where we have to use a Weighted overlay to generate a raster of suitable places to plant crops. I've put all the factors in the weighted overlay and generated my output but now I need to consider the constraints. These are built up areas and conservation areas that need to be classified as 'not suitable'. The lecturer says that we need to use a CON statement in the raster calculator to get the '0' values from the constraint's raster put onto the weighted overlay raster.

But, I have no idea how to do that. All the resources I have found use the format CON (in raster = 0, x, y) but i don't know how to fill the Y value with an entire raster... if that makes sense?


texts:
I have to style some Rasters with SLD and get the legend from Geoserver to show it in my webapp.

By default, labels are placed right to the Color Ramp.

Is there any option, by using SLD or maybe any parameters I can pass to the request for the legend, so that I can place labels at some other position?

So far I tried with LEGEND_OPTIONS, but it seems there's no parameter that fits my needs.

.../geoserver/temper/wms?REQUEST=GetLegendGraphic&VERSION=1.0.0&FORMAT=image/png&LAYER=temper:DICIEMBRE2017&WIDTH=20&HEIGHT=5&transparent=true&LEGEND_OPTIONS=fontSize:12;fontStyle:bold;fontColor:0xFFFFFF;dpi:180;


texts:
I have gdal installed and (I though) working fine. I use a lot rasterio and sometimes gdal directly, but today I wanted to do a mosaic and I wanted to try gdal_merge.py, but when I try to launch the tool I get this error message:

Traceback (most recent call last):
  File "C:\Anaconda3\Scripts\gdal_merge.py", line 4, in <module>
    __import__('pkg_resources').run_script('GDAL==2.1.1', 'gdal_merge.py')
  File "C:\pyqgis\lib\site-packages\pkg_resources\__init__.py", line 2927, in <module>
    @_call_aside
  File "C:\pyqgis\lib\site-packages\pkg_resources\__init__.py", line 2913, in _call_aside
    f(*args, **kwargs)
  File "C:\pyqgis\lib\site-packages\pkg_resources\__init__.py", line 2940, in _initialize_master_working_set
    working_set = WorkingSet._build_master()
  File "C:\pyqgis\lib\site-packages\pkg_resources\__init__.py", line 635, in _build_master
    ws.require(__requires__)
  File "C:\pyqgis\lib\site-packages\pkg_resources\__init__.py", line 943, in require
    needed = self.resolve(parse_requirements(requirements))
  File "C:\pyqgis\lib\site-packages\pkg_resources\__init__.py", line 829, in resolve
    raise DistributionNotFound(req, requirers)
pkg_resources.DistributionNotFound: The 'GDAL==2.1.1' distribution was not found and is required by the application

I don't know what happens here. I try to run directly the gdal command line and I get one message saying I don't have osgeo installed (I actually have it). So I am pretty lost here. I have been reading a lot of post about gdal installation but I do the same I always see, and the proof is that I can work with gdal (from python or the cmd), but I don't realize where is my mistake.


texts: leaflet-kml Lat Lng changes and zoom don't work with KML

texts:
I am trying Leaflet KML plugin as below:
https://github.com/windycom/leaflet-kml
const map = new L.Map('map', { center: new L.LatLng(28.412, 77.317), zoom: 11 });

It works, loads the KML over OSM, running as localhost on server but if am changing lat longs or zoom level, it still shows the same location area always. It seems like these settings are not giving any effect. I want to point it to the specific lat long of the user who views it and with specific zoom level.
Here is full code:

<meta name="viewport" content="width=device-width, initial-scale=1.0,  maximum-scale=1.0, user-scalable=no" />
        <link rel="stylesheet" href="http://unpkg.com/leaflet@1.4.0/dist/leaflet.css" />

If I remove the KML file, then lat long pointing and zoom works fine.
It would be great if you clarify one more thing. I have KML file of my whole country that size around 15-20 MB. If I load entire KML then the map will download and load  full KML of 15-20 MB on user side or the KML loads in chunks as and when user drags to other locations of country.


texts:
In my data I have many occurrences of 2 lines that follow the same path, but in opposing directions. Applying LineMerger to this data results in the lines being joined together, but not merged into each other.

e.g. Given two lines with the same nodes

--->--->--->---

---<---<---<---

I want either:

--->--->--->---

or

---<---<---<---

But instead I get

--->--->--->--+
              |
---<---<---<--+

or

+-->--->--->---
|
+--<---<---<---

Is there a way to obtain my desired result in JTS?


texts: Join linestrings in shapefile of watercourses that have same NAME and PERENNIALITY

texts:
I have a shapefile (link below) of only watercourses comprised of many short linestrings.

OpenJUMP displays this about each segment

LINESTRING (
    153.05422 -26.44881350000003,
    153.05448 -26.448845000000006,
    153.05506349999996 -26.448870499999998,
    153.06856249999998 -26.448498500000028
)

I wish to join all the separate linestrings segments that have the same NAME and PERENNIALI (Perennial or Non Perenial)

For example: I wish to have all of the Noosa River that is perennial as one selectable piece and another one for Noosa River that is non perennial, and so on for all the watercourses (if there is a gap in the data then I would get a separate line for each piece)

I have used Dissolve from the processing toolbox and this works fine in QGIS but when I save and open another program like JOSM editor the segments are not joined and I do wish to use in JOSM.

I tried 'Vector Geometry>Multipart to singleparts' but cannot work out how to join on the NAME and PERENNIALI attributes.

I notice that the OGR Dissolve only allows one dissolve field option in the Processing Toolbox screen but at the bottom displays the GDAL/OGR console call as ...
ogr2ogr /var/folders/z5/.../OUTPUT.shp "/Volumes/.../Noosa_Hydrolines_named.shp" -dialect sqlite -sql "SELECT ST_Union(geometry) AS geometry, NAME FROM 'Noosa_Hydrolines_named' GROUP BY NAME" -f "ESRI Shapefile"
So, can I alter that to use both the NAME and PERENNIALI and run from the Terminal console?
If so, what changes do I make?

After saving as a shapefile I wish to open in JOSM and select the joined bits a one segment.

I am using MacOS Mojave 10.14.2, QGIS 3.4.1-Madeira (Kyngchaos)

I am a beginner with QGIS or in fact any GIS so simple instructions are preferred.

https://www.dropbox.com/s/y77zrwpsrh9qmmk/Noosa_Hydrolines_named.shp?dl=0


texts: Extracting all points within a certain distance and height difference from a Raster

texts:
I have a raster containing tens of millions of points all with X, Y and Z data defined. I need to figure out a simple way to extract and create a new raster or vector with ALL points at a 720 m distance from each other (X,Y) and at 120 m elevation difference (Z) from each other.

I have ZERO knowledge of SQL and Python. I have been trying to do this on VBA and came up with a couple of algorithms but the processing time is unreasonable and unrealistic. I am sure there must be a simple GIS approach to accomplish this but cannot seem to find it.

I am using ArcMap.


texts:
I'm adding layers to a GeoExt map panel map dynamically in response to a user's selection. I have a corresponding FeatureStore definition from which I create the layers. When a new layer is added, I'd like it to be the top most layer.

var vecLayer = new OpenLayers.Layer.Vector(natureOfAttack, {
            styleMap: new OpenLayers.StyleMap({
                'default': style
            }),
            protocol: new OpenLayers.Protocol.HTTP({
                url: "./server/maps/incidents.php",
                params: {natureOfAttack: natureOfAttack},
                format: new OpenLayers.Format.GeoJSON()
            }),
            strategies: [new OpenLayers.Strategy.Fixed()]
        });
        layers.push(vecLayer);

        // manually bind store to layer
        me.getIncidentMapStore().bind(vecLayer);

        mapPanel.map.addLayers(layers);

        // some more controls
        var selectCtrl = new OpenLayers.Control.SelectFeature(vecLayer);
        mapPanel.map.addControl(selectCtrl);

        vecLayer.events.on({
            'featureselected': function(e) {
                me.showIncidentDetailPopup(e.feature);
            }
        });

        selectCtrl.activate();

        // for dev purpose
        map = mapPanel.map;
        mapPanel = mapPanel;

I've then added:

var cLayer = map.getLayersByName(natureOfAttack)[0];
    map.setLayerIndex(cLayer, map.layers.length+1);

in an attempt to make the layer the top most layer, but this doesn't seem to work.

Having thought about it for a while, I came to the conclusion that I should probably have a listener in the store definition as shown below:

Ext.define('CMD.store.SampleStore', {
    extend: 'GeoExt.data.FeatureStore',
    model: 'CMD.model.Summit',
    autoLoad: false,

    listeners: {
        add: function(thisStore, records, index) {
            console.log('FeatureStore: A new layer added...');
        }
    }
});

but the event in the listener is not executed. Any idea on how I can achieve this will be appreciated.

UPDATE:

I've added:

vecLayer.events.on(
        'loadend': function(layer) {
            map.setLayerIndex(layer, map.getLayerIndex(layer)+5);
        }
    });

to the vector layer. This makes the created vector layer the top most layer. However, the error TypeError: a.setZIndex is not a function is generated.

UPDATE:

My environment: GeoExt 2, ExtJS 4, OpenLayers 2.12.


texts:
I am working on a QGIS project that involves me using the Google Hybrid Map Plug-in. I am due to spend time travelling over the next few days and won't have an internet connection. I was hoping someone could let me know if it's possible to continue using the plug-in without internet? For example by saving the relatively small I need for offline use.

Thanks
Sarah


texts: How remove geometry vertexes on hover OpenLayers 3+

texts:
I am using modify interaction and would like to be able to delete vertexes by either clicking delete + mouse on the vertex or delete + hover vertex. delete + mouse works, delete + hover does not as this.dragSegments_ is empty when accessed in removeVertex_ therefore it results in point I am hovering over not being deleted.

Taken from Modify.js

removeVertex_() {
    this.gdFeaturesNotEdited = [];
    const dragSegments = this.dragSegments_;
    const segmentsByFeature = {};
    let deleted = false;
    ...

Only places I see this.dragSegments_ being populated are handleDownEvent and insertVertex_. From the naming of the variable and functions it seems to me that vertex deletion on hover is not something that was designed to happen in OpenLayers 3+.

var modifyInteraction = new ol.interaction.Modify({
    features: selectCtrl.getFeatures(),
    deleteCondition: function(event) {
        var deletePressed = keys['46']; // 46 = delete key
        var isDeleteAction = event.type == 'pointermove' || event.type == 'pointerdown';

        return deletePressed && isDeleteAction ;
    }
});

Code that needs to be modified in TomazicM's answer to work with linear ring start/end coordinate deletion:

...
if (dist <= 8) {
    selectedCoord = coord[k];
    coord.splice(k, 1);

    // if deleting a coordinate at index 0 then coordinate at index coord.length - 1 will also be deleted because this is a coordinate that "opens" and "closes" linear ring
    // need to handle this by creating a new "closing" coordinate
    if (k == 0) {
        coord.splice(coord.length - 1, 1);
        coord.push(coord[0]); // creating new "closing" coordinate by duplicating first coordinate in array
    }

    setTimeout(updateGeometry, 150, feature, coord);
}
...


texts: How to initialize a PostGIS database in PostgreSQL single-user mode?

texts:
Update: Reposted this question on StackExchange / Database Administrators: https://dba.stackexchange.com/q/152496/57043

I'm using PostgreSQL 9.2.15 with PostGIS 2.0.7 through the standard packages available on CentOS 7 with EPEL.

I have a server script to automatically setup a PostGIS database, by first executing initdb, then using postgres --single to create a new database with PostGIS extension. This works roughly like this:

initdb -D ...

postgres --single -D ... postgres <<'EOF'
    CREATE DATABASE mydb;
EOF

postgres --single -D ... mydb <<'EOF'
    CREATE EXTENSION postgis;
EOF

However, the last command produces the following strange error message:

ERROR:  type reclassarg[] does not exist
STATEMENT:  CREATE EXTENSION postgis;

If I instead run the database server and execute CREATE EXTENSION postgis via sql, everything works fine:

postgres -D ... &

# Wait until service it up and running

psql mydb <<'EOF'
    CREATE EXTENSION postgis;
EOF

But that is not suitable for a database initialization script.

So how can I initialize a PostGIS database in PostgreSQL single-user mode?

(Also, is this question appropriate here, or should it go to another StackExchange community?)


texts: Why do I get the error "OGR error setting feature:5" while changing the attribute of the layer?

texts:
Getting the error while changing the name of one attribute in the point layer.

Errors: ERROR: 1 attribute value change(s) not applied.
Provider errors: OGR error setting feature 5:


texts: Getting map document name/path from published map service using ArcPy/ArcObjects?

texts:
I have several map services published to ArcGIS Server 10.0 SP2 and I have layer metadata files exported to HTML that I want to link to the layers in the map service.

Currently the metadata HTML files are named based on the data source (eg 'Zip_Code_Poly.htm') but I need to correlate the layer name/ID from the map service to that file. I was thinking that the easiest way to do it would be to run an arcpy script to get all the layer data sources.

I've run into an issue though that the client (this is a web app) needs to look up these HTML files by service name and layer name/ID but without hard coding a mapping from service -> MXD I can't figure out what map document to use in arcpy.

Is there a way using arcpy or something else that I can look at a published map service and get the map document it was created from?


texts: How to split a single OSM layer of a Key (eg highway) containing multiple Values into separate layers of the values in QGIS?

texts:
I am using Quick OSM in QGIS 3.12.3-București. When I 'Run query' for a specific geocodeArea, eg London, and I specify the 'Key', eg Highway, but I leave the 'Value' blank; (img 1). The result is one layer in QGIS with all of the Values of highway on it; (img 2) (well its 3 layers - Nodes, lines and hatches. but for the sake of this question, let's just say it's one layer).
What I would like to achieve is to have all of the Values for the Key (eg highway) to be on individual layers in QGIS. So I can avoid having to 'run query' for each and every value of 'highway'. As this is a laborious process and there must be a faster way.
I have looked into the 'Overpass query' method, but iv'e only found it gives me the same result of everything in one layer.
Is there a piece of code I can add into the Overpass query to tell it to load the values on to separate layers? Or is there a way to split the layer after it has been loaded, into separate layers of the Values of the Key?


texts: How to store the Shapefile M attribute in Geopandas/Shapely?

texts:
Let's say I  have a geopandas.GeoDataFrame with a shapely.geometry.polygon.Polygon as geometry type.

The problem is that I have an attribute per vertex of the polygon, which could theoretically be stored in the shapefile's M slot. The only problem is that I could not find an option to include this M into the polygon. The shapely documentation says:

  The Polygon constructor takes two positional parameters. The first is an ordered sequence of (x, y[, z]) point tuples [...]

Further, when using the geopandas.GeoDataFrame.to_file method, I couldn't find an option to insert this kind of information.

Can anyone help?


texts:
I updated from Geoserver 2.1.3 to 2.8.0 and I copied my SLD styles.

I had all the images in a folder name "symbols" and my SLD was like

<PointSymbolizer>
  <Graphic>
    <ExternalGraphic>
      <OnlineResource          xlink:type="simple"xlink:href="symbols/aerodrome.p.16.png" />
       <!-- href has a "symbols/" prefix to access the folder -->
      <Format>image/png</Format>
    </ExternalGraphic>
    <Size>32</Size>
  </Graphic>
</PointSymbolizer>

I also copied the "sympbols" folder to the data_dir\styles folder.

Images are now in data_dir\styles\symbols

Used to work in 2.1.3, but not in 2.8.0, I cannot see the images in the preview.

What am I missing here?


texts:
I recently started using JavaScript API in attempt to construct a web app. After taking a recent ESRI online course, I thought I had some base knowledge on how to do this.

This geoprocessing service I'm using from our development server is supposed to take an input of a ZipCode and in turn display delivery routes associated with said ZipCode.

I've looked at some sample code on the ArcGIS Javascript API help site and cannot figure out how to execute this with its parameters.

If anyone has any insight or where I can look for documentation on how I can initiliaze this task and make it appear on my WebMap, I'd be forever indebted.

Thanks!

Below is the REST services page for the geoprocessor.

Execution Type: esriExecutionTypeSynchronous

Parameters:
Parameter: ZIP
Data Type: GPString
Display Name ZIP
Description: ZIP code or list of comma, space, or tab delimited 5 digit ZIP codes
Direction: esriGPParameterDirectionInput
Default Value: 38138
Parameter Type: esriGPParameterTypeRequired
Category:

Parameter: Rte_Box
Data Type: GPString
Display Name Rte_Box
Description: R = carrier route, B = PO Box
Direction: esriGPParameterDirectionInput
Default Value: R
Parameter Type: esriGPParameterTypeOptional
Category:

Parameter: Output_Feature_Class
Data Type: GPFeatureRecordSetLayer
Display Name Output_Feature_Class
Description: Geometry and delivery stats for routes that meet search criteria
Direction: esriGPParameterDirectionOutput
Default Value:
Geometry Type: esriGeometryPolyline
HasZ: false
HasM: false
Spatial Reference: 102100  (3857)

Fields:
OBJECTID ( type: esriFieldTypeOID , alias: OBJECTID )
ZIP_CRID ( type: esriFieldTypeString , alias: ZIP_CRID , length: 13 )
ZIP_CODE ( type: esriFieldTypeString , alias: ZIP_CODE , length: 5 )
CRID_ID ( type: esriFieldTypeString , alias: CRID_ID , length: 4 )
CRID_IND ( type: esriFieldTypeString , alias: CRID_IND , length: 1 )
CRID_NBR ( type: esriFieldTypeString , alias: CRID_NBR , length: 3 )
CITY_STATE ( type: esriFieldTypeString , alias: CITY_STATE , length: 32 )
BUS_CNT ( type: esriFieldTypeDouble , alias: BUS_CNT )
TOT_CNT ( type: esriFieldTypeDouble , alias: TOT_CNT )
DROP_CNT ( type: esriFieldTypeDouble , alias: DROP_CNT )
RES_CNT ( type: esriFieldTypeDouble , alias: RES_CNT )
DROPSHIP_KEY ( type: esriFieldTypeString , alias: DROPSHIP_KEY , length: 6 )
FACILITY_NAME ( type: esriFieldTypeString , alias: FACILITY_NAME , length: 35 )
LESS_THAN200_IND ( type: esriFieldTypeString , alias: LESS_THAN200_IND , length: 1 )
FAUX_IND ( type: esriFieldTypeString , alias: FAUX_IND , length: 1 )
LT_200_IND ( type: esriFieldTypeString , alias: LT_200_IND , length: 1 )
DS_KEY ( type: esriFieldTypeString , alias: DS_KEY , length: 7 )
FAC_NAME ( type: esriFieldTypeString , alias: FAC_NAME , length: 35 )
Shape_Length ( type: esriFieldTypeDouble , alias: Shape_Length )
Features: None.

Parameter Type: esriGPParameterTypeRequired
Category:


texts:
I am looking for alternative (more efficient) methods for clipping large feature class datasets. When using the 'Clip' tool in ArcGIS, often times the tool will crash after many hours of processing. Without getting into specifics regarding file size or computation capabilities, I was wondering if there were better methods for clipping large datasets. Looking around on the Web, I've seen GDAL Ogr2Ogr mentioned several times. I was hoping someone could explain what exactly ogr2ogr is and if it truly is a more efficient way to clip data vs. using tools in ArcGIS.

I'm pretty new to the Geospatial field so my technical jargon isn't very vast yet.

The .gdb in question is 4.8 GB with a feature count of 21,858,032 features. ArcMap 10.2 is being used. As for computer specs: Intel Xeon @ 3.50GHz, 32.0 GB of RAM and it's a 64 bit Operating System.


texts:
I want to build a algorithm (Matlab), which detects if a 3D flight path intersects with clouds (cloud=cube at first). Flight path waypoints are on a great circle and between waypoints I move on rhumb lines. I am struggling a bit with the use of the coordinate frames. Within the geodetic frame (Lat/Lon/Alt) I can define position well, but I can´t use vector calculation. With in the ECFS I can use vectors but I get issues with the
curvature of the earth on long distances. ENU does´t work because the clouds might be bigger than 100 Nm.
I would love to have a system where I can enter (Lat/Lon/Altitude) and can calculate with vectors. But how do I manage the difference between rhumb line and great circle.

Does anyone have an idea?


texts: Pgrouting functions and geoms type not found. Install failed?

texts:
I have installed a postgresql 9.1 and postgis 2.0 from source.

I couldn't launch this :

# Add pgRouting launchpad repository ("stable" or "unstable")
sudo add-apt-repository ppa:georepublic/pgrouting[-unstable]
sudo apt-get update

# Install pgRouting packages
sudo apt-get install postgresql-9.1-pgrouting

So I compiled and installed pgrouting 2.0 (after some hours searching for dependencies).
I created the extension on my database in postgresql.

I included the function from the sql files pgrouting.sql but the functions I need are in pgrouting_legacy.sql and pgrouting_dd_legacy.sql

When I try to load them, the error I get is : psql:/usr/share/postgresql/9.1/contrib/pgrouting-2.0/pgrouting_legacy.sql:299: ERROR: type "geoms" does not exist

Postgresql and Postgis are working fine...

What did I wrong ? Maybe I forgot something or the install failed ?

I followed this documentation : http://pgrouting.org/docs/1.x/install.html

http://www.bostongis.com/PrinterFriendly.aspx?content_name=pgrouting_osm2po_1


texts: Adding new polygon attribute column based on location withing another polygon in R

texts:
I have here an example with two layers - one containing two large polygons Pol_large, and second containing 4 small polygons Pol_small. For Pol_small
I want now to find the location of each of the small polygons within the large polygons creating a new attribute column location and transporting the attributes from the Pol_large$id column. It would look following:

Pol_small$location
[1] "A1" "A1" "A2" "A2"

Not sure which funciton in sf or sp would be best to do it.

library(sp)

create_poly <- function(coordin,id_number){
  coordinates = matrix(c(coordin[1], coordin[3],
                         coordin[1], coordin[4],
                         coordin[2], coordin[4],
                         coordin[2], coordin[3],
                         coordin[1], coordin[3]),
                       ncol = 2, byrow = TRUE)
  P_1 = Polygon(coordinates)
  Ps_1 = SpatialPolygons(list(Polygons(list(P_1), ID = id_number)), proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
  ID_1 <- sapply(slot(Ps_1, "polygons"), function(x) slot(x, "ID"))
  df <- data.frame(rep(0, length(ID_1)), row.names=ID_1)
  Ps1_df <- SpatialPolygonsDataFrame(Ps_1, df)
  Ps1_df
}

# Create shapefile with polygon boundaries of interest:
Pol_1 <- create_poly(c(50,100,10,50),1)
Pol_2 <- create_poly(c(110,170,0,50),1)
Pol_3 <- create_poly(c(55,65,40,45),1)
Pol_4 <- create_poly(c(70,90,30,50),1)
Pol_5 <- create_poly(c(145,160,20,40),1)
Pol_6 <- create_poly(c(120,135,10,30),1)

# Join the polygons into one layer and assign attributes:
Pol_large <- do.call(rbind, list(Pol_1, Pol_2))
Pol_large$id <- c("A1","A2")
Pol_small <- do.call(rbind, list(Pol_3, Pol_4, Pol_5, Pol_6))
Pol_small$id <- c("b1","b2","b3","b4")

# Plot
plot(Pol_large, axes = TRUE)
plot(Pol_small, add=T, col = "blue")

Is this solution optimal? For some reason, in my real large example, this does not seem to work exactly as it should - specifically some smaller polygons with boundaries touching the edges of the large polygons do not get any assignment, they return integer(0)

library(sf)
a <- unlist(st_within(st_as_sf(Pol_small),st_as_sf(Pol_large)))

for (i in 1:length(a)){
  g <- a[i]
  Pol_small$location[i] <- Pol_large$id[g]
}


texts: Fixing ArcGIS Pro crash with opening notebooks after installing Anaconda

texts:
I have had problems with opening notebooks in ArcGIS Pro 2.5 since installing 64-bit Anaconda 3 on my PC (Windows 10 OS). When I create a project in ArcGIS Pro and try to open a notebook within project (or open Python window or Analysis tab), software fatally crashes.
I uninstalled Anaconda from my PC and hoped it would help ArcGIS Pro begin behaving normally again, but I still have same error. Based on what I've read here, it seems that some people have had success with installing a 32-bit version of Anaconda that runs concurrently on same machine as ArcGIS Pro.
Does anyone have suggestions of how I might fix ArcGIS Pro issue with opening notebooks?


texts: How to add feature class to MXD with ArcPy (Python)?

texts:
I am using ArcGIS 10.1 and I am trying to add a feature class that I create to an MXD, so that I can join to it.

This all needs to be done using Python.

Is this possible?


texts:
Let us say, there are two tables, table1 and table2. Both have one column (geometry_1 and geometry_2 respectively) of text type.

    table1

    geometry_1
    text
1   LINESTRING(37.00 59.00, 37.05 59.32)
....

    table2

    geometry_2
    text
1   LINESTRING(37.89 59.55, 38.05 60.3)
....

Both tables have more rows, but for now I want to focus on the following question. Is there any way to evaluate if there are "similar" geometries between the tables. By "similar"? I mean that if the distance between the respective points of the lines is no higher than a certain value (eg. 100m), the two lines are considered identical.

I have tried this:

select geometry_1 from table1 t1
    where exists(
             select 1 from table2 t2  where
             (select ST_Equals(ST_GeomFromText(t1.geometry_1),ST_GeomFromText(t2.geometry_2)) = true));

However the query takes too long. I m not sure if what I wrote is correct.

There are only linestrings with two or more vertices
Linestings with reversed directions are considered identical
If lines have more vertices and vertex B is in the middle of the line so that it does not bend the line,  A-B-C, A-C, A-B, B-A and C-B are identical
SRID=4326

Update1:

create index geom_1_idx ON table1 USING GIST (geometry_1);
select table1_id, table2_id  from table1 t1, table2 t2
              where (select st_DWithin(t1.geometry_2 , t2.geometry_2, units_of_srid) = true)

Consider it as initial approach to the problem regardless the above conditions


texts:
I'm trying to see how good a custom topographic map could be using only free data

For the moment I'm just playing with the data so of course the colours are not pretty.
I did some classification and used some land cover data anyway I have a raster layer as background with the main land cover basically swamp forest and no vegetation.
Since the result of a vectorization was awful and hard to fix I'm going to stay with the raster as background, the question is:

How to avoid the squared shape on my map? I would like to smooth the area to have a prettier map. I tried the bilinear interpolation but it creates white boundaries around the area and with the Nearest neighbour it's not pretty.

I tried some generalization tools as well but they are way too strong!

So what would you do?

Here is what I have done so far:


texts:
I have MCD43A3 product and I am using that to get the albedo values(Black sky and White Sky). I have some in-situ snow depth data, corresponding to which I want to get these albedo values.

For most of the in-situ data I'm getting a value of 32767, which is a fill value for the product. For some data I'm getting, a value like 908,889,976 which is multiplied by scale factor of 0.0010 to get the albedo values like .908,.889,.976.

My question is what should I do about the 32767 values which I am getting for most part of my data? Does getting this type of value mean that the image has noise or error? I have little experience with MCD products, any kind of help would be appreciated.
Thanks,
GBh


texts:
I would like to use QField and do structural measurement. Tablet includes the compass as well as a clinometer to quickly measure lineation on the rock.
Can we add these field in QField?


texts: I am using ENVI 5.1 and when I use the Data Manager to load several images and if I want to load two adjacent images, why would one of them load and the other doesn't?. Am I missing something here ?


texts: Read Sentinel-2 zip-file into memory using gdal in python

texts:
I was wondering if there's a way to read a Sentinel-2 zipfile directly into memory through python without extracting the whole thing or selected contents.

I figured since GDAL does have a driver for the S2 format, it should be somehow possible to open the zipfile in python and then read the bands into arrays separately.

Another approach I was thinking about is the zipfile module in python. I came as far as accessing the single files in the archive, but I haven't managed yet to read them into memory or pass them over to GDAL.

Edit:

Sorry, I forgot to mention I'm generally looking for a platform-independent solution

Edit2:

Also, I intend to stay away from ESA's snappy

Solution:

Alright, after the helpful tips I got here I finally implemented the following solution. I ended up installing gdal 2.1.3 from conda-forge inside a new environment to bypass any dependency issues. There's still an issue with the UTF encoding using vsizip which can be fixed setting the corresponding environment variable as discribed in  this issue on github.

#!/usr/bin/env python

import gdal
from zipfile import ZipFile
import re
import os

## vsizip bugfix
os.environ['CPL_ZIP_ENCODING'] = 'UTF-8'

#choose tile
tile = 'T36PXV'

tileREXP = re.compile(r'.*(%s).*B(02|03|04|05|06|07|08|8A|11|12).jp2$' %tile)

# sort bands 2 - 12 with 8A after 8
bndIX = [1,2,3,4,5,6,7,9,10,8]

with open('S2_FILELIST.txt') as flist:
    for filename in flist:
        with ZipFile(filename,'r') as zfile:
            bands = [x for x in zfile.namelist() if re.match(tileREXP,x)]
            bands.sort()
            bands =  [x for (y,x) in zip(bndIX,bands)]
            for b in bands:
                ds = gdal.Open('/vsizip/%s/%s' %(filename,b))
                arr = ds.ReadAsArray()

                # DO STUFF

                ds = None


texts: Determining the Coordinate System of shapefiles

texts:
I have 8 shapefiles which have a common coordinate system but none of them have it defined i.e. "unknown" is written when I look in the properties dialog box. However, I can see the extent which I am sharing here:

Top: 42.349145
Bottom: 41.710962
Right: -82.430463
Left: -83.126193

I need to define the projection of these shapefiles so that I could do further analysis. I know that the features these shapefiles represent are located in Windsor, Ontario i.e.they should have a projection of WGS 84 UTM Zone 17N. Can I directly define this projection without knowing about the geographic coordinate system?


texts: Problem with Orfeo Toolbox integration in QGIS

texts:
I'm getting this error:

OTB execution console output
'"C:\OSGeo4W\bin\otbcli_MeanShiftSegmentation"' is not recognized as an internal or external command,
operable program or batch file.


texts:
I am simply trying to delete a TIF that I created with GDAL. I have tried, simply, to use shutil which errors with:

NotADirectoryError: [WinError 267] The directory name is invalid:
'T:/p003M01p0001/extent_files/t215_E.tif'

My first attempt was to use:
filename = out_loc + '/' + temp_tif
driver = ogr.GetDriverByName("GTiff")
if os.path.exists(filename):
    driver.DeleteDataSource(filename)

but this just resulted in a None driver for some reason (AttributeError: 'NoneType' object has no attribute 'DeleteDataSource'). This proves os recognizes the file exists (so why does shutil not work?). I have confirmed the existence of the GTiff driver using:
driver_list = []
for i in range(gdal.GetDriverCount()):
    driver = gdal.GetDriver(i)
    driver_list.append(driver.GetDescription())
    print(gdal.GetDriver(i).ShortName)

I am a bit perplexed that ogr couldn't delete it, but even more perplexed on why shutil doesn't recognize that it exists. There must be something I am doing wrong creating the tiff.
The overall purpose of the following code is to get the boundary of valid TIF values.
Current code:
from osgeo import gdal, osr
from osgeo.gdalconst import GA_Update
import ogr
from core.gdal_functions import GDALFunc

input_tif = 'T:/p003M01p0001/t215.tif'
out_loc = 'T:/p003M01p0001/extent_files'
copy_tif = out_loc + '/t215_E.tif'
temp_tif = out_loc + '/shaded.tif'
shp = out_loc + '/t215_Ex.shp'

# get copy of tif
src_ds = gdal.Open(input_tif, gdal.GA_ReadOnly)
srs = osr.SpatialReference()
srs.ImportFromWkt(src_ds.GetProjectionRef())
dst_ds = src_ds.GetDriver().CreateCopy(copy_tif, src_ds, 0)
dst_ds = None
src_ds = None

# Edit data to remove 0 to turn it into nodata
nodata = 0
ras = gdal.Open(copy_tif, GA_Update)
for i in range(1, ras.RasterCount + 1):
    ras.GetRasterBand(i).SetNoDataValue(nodata)
ras = None

# create temp tif of where actual values exist
drv = ogr.GetDriverByName('ESRI Shapefile')
outfile = drv.CreateDataSource(shp)

outlayer = outfile.CreateLayer('polygonized raster', srs=srs)
newField = ogr.FieldDefn('DN', ogr.OFTReal)
outlayer.CreateField(newField)
command = 'python C:/Users/david/AppData/Local/Programs/Python/Python38/Lib/site-packages/osgeo/scripts/gdal_calc.py ' \
          '-A {infile} --A_band=1 --outfile {outfile} --calc="(A>0)*1" --quiet'.format(infile=copy_tif, outfile=temp_tif)

process = subprocess.Popen(command, shell=True)
while process.poll() is None:
    time.sleep(0.5)

#create polygon extent
tif = gdal.Open(temp_tif)
band = tif.GetRasterBand(1)
gdal.Polygonize(band, None, outlayer, 0)
outfile = None
tif = None

shutil.rmtree(copy_tif)


texts: The underlying connection was closed: The connection was closed unexpectedly. AOInitiliaze

texts:
I create a WCF Service.

When I try to Initialize ArcGIS License with AOInitialize I get error:

  The underlying connection was closed: The connection was closed
  unexpectedly.

public void initializeArcGISLicenses()
    {

            RuntimeManager.Bind(ProductCode.EngineOrDesktop);
            AoInitialize ao = new AoInitialize();
            if (ao.IsProductCodeAvailable(esriLicenseProductCode.esriLicenseProductCodeAdvanced) == esriLicenseStatus.esriLicenseAvailable) //Error is here
                ao.Initialize(esriLicenseProductCode.esriLicenseProductCodeAdvanced);
           /* else if (ao.IsProductCodeAvailable(esriLicenseProductCode.esriLicenseProductCodeBasic) == esriLicenseStatus.esriLicenseAvailable)
                ao.Initialize(esriLicenseProductCode.esriLicenseProductCodeBasic);
            else if (ao.IsProductCodeAvailable(esriLicenseProductCode.esriLicenseProductCodeStandard) == esriLicenseStatus.esriLicenseAvailable)
                ao.Initialize(esriLicenseProductCode.esriLicenseProductCodeStandard);*/
    }

What's the problem?
Sorry for my English.


texts: Splitting RGB image into separate bands using Drone (DJI Phantom imagery)

texts:
Previously I have made marine habitat maps using multiband imagery and remote sensing techniques (e.g blue, green and red bands from satellite imagery).

I am now attempting the same process however using drone imagery (DJI Phantom 3 pro). Is it possible to obtain separate bands from the RGB images that are captured and if so how can I do this?


texts:
I have sfc LINESTRING and I would like to remove overlapping lines. Object data  was created from an algorithm that joins lines with minimum distance. Using unique will not remove the initial lines.
In the example you can see that line 2 and 3 + other lines were the initial lines, then they were used to produce line 1, other lines were removed but 2 and 3 remained. and length of lines proves that!. Extracting the line with the maximum length will not work in my case because in some cases I have a non-overlaping line which I would like to keep!
> class(data)
[1] "sfc_LINESTRING" "sfc"
> length(data)
[1] 3
> data
Geometry set for 3 features
geometry type:  LINESTRING
dimension:      XY
bbox:           xmin: 181.4494 ymin: 4.441853 xmax: 278.9707 ymax: 14.72429
CRS:            +proj=longlat +datum=WGS84 +no_defs
LINESTRING (181.4494 4.462799, 181.4722 4.45583...
LINESTRING (181.4494 4.462799, 181.4722 4.45583...
LINESTRING (181.4494 4.462799, 181.4722 4.45583...

> st_overlaps(data,data)
although coordinates are longitude/latitude, st_overlaps assumes that they are planar
Sparse geometry binary predicate list of length 3, where the predicate was `overlaps'
 1: 2, 3
 2: 1, 3
 3: 1, 2

> st_length(data)
Units: [m]
[1] 11628206  7143693  9899478

Criteria to remove:

The line should be overlapping with other lines
The line is not the longest line.


texts: Create Point Feature with XY and Generate Multiple Buffers (ArcPy Tool)

texts:
I am trying to set up a custom tool via a python script in ArcMap. I want the tool to be able to create a point feature with X/Y input and then pass that point into the buffers.

My thinking is that someone could use the tool to find what is inside the radius from their given point that they input themselves without having to create a feature class first (I hope that makes it clearer).

Am I doing this right by creating the point feature and using the GetParameterAstext string? What exactly would go in the multiple ring buffer at the bottom?

import arcpy
import os

# Create a Point object from coordinates then set up buffer rings
pnt.X = arcpy.GetParameterAsText(0)
pnt.Y = arcpy.GetParameterAsText(1)
arcpy.Point(pnt.X, pnt.Y) = arcpy.GetParameterAsText(2)

BuffRing1 = arcpy.GetParameterAsText(3)
BuffRing2 = arcpy.GetParameterAsText(4)

distances = [BuffRing1, BuffRing2]
measure = ['miles']

arcpy.MultipleRingBuffer_analysis()


texts: Point of intersection issues with PostgreSQL/PostGIS?

texts:
I have a table containing a set of points of intersections.

Most of the roads intersects whith just one road BUT i have some roads which intersects with two, therefore three different point of intersection has been created. BUT, i want only one point of intersection for all the three distinct roads.

Any suggestion how I can merge the 3 points to just one where there is such situation.

I am using postgresql with postgis

The following is the table where i have the points of intersection

The following is a graphical screenshot of the situation.