texts: QGIS "Export/Add Geometry Columns..." not responding to redefining CRS

texts:
I created a vector layer using the EPSG:3857 / Pseudo-Mercator CRS. After randomly sampling points within it, I am now trying to export their lat/long values in a different projection (e.g., a UTM based projection). It seems that redefining the layer CRS via either the:

Rick-click on layer > Set CRS...

or

Vector > Data Management Tools > Define Current Projection...

does not induce changes in the coordinates that are output by the:

Vector > Geometry Tools > Export/Add Geometry Columns...

function.

It appears that only saving the vector layer as a new layer with a new CRS will induce changes in the coordinates that are output.

Is this the appropriate behavior? Are the Set CRS... and Define Current Projection... functions purely for visualization? Can others confirm that this is the behavior that occurs for them when they follow the same workflow?

Using QGIS 2.14.8-Essen.


texts: Convert from EPSG:27700 to latitude and longitude

texts:
I'm trying to transform EPSG:27700 coordinates to latitude and longitude. From my research I understand that WGS 84 (EPSG:4326) is the SRS that resembles latitude and longitude.

I'm using Proj4j to try and convert a test coordinate (those familiar with my previous question will recognise Kirkby In Ashfield station).

My Java code snippet:

String projection27700 = "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 "
    + "+x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs";
CRSFactory crsFactory = new CRSFactory();
CoordinateReferenceSystem sourceSRS = crsFactory.createFromParameters(null,
    projection27700);

String projection4326 = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs";
CoordinateReferenceSystem targetSRS = crsFactory.createFromParameters(null,
    projection4326);

CoordinateTransform transformation = new BasicCoordinateTransform(sourceSRS,
    targetSRS);

ProjCoordinate result = new ProjCoordinate();
ProjCoordinate input = new ProjCoordinate(450077, 356089);
transformation.transform(input, result);

System.out.println(result.x + ", " + result.y);

The result printed is:

-155.20355661333892, -35.91545136124713

This doesn't seem right - it neither makes sense as degrees or metres (to me at least). Have I missed out a conversion somewhere, either on my input data or in this output data I've produced?

I'm fairly happy my input data is right, based on this documentation, which confirms that +units=m will allow my input data to be in metres. It also suggests the output data should be in degrees, which is what I want.


texts: Improving Global Scene transition when changing location using bookmarks in ArcGIS Pro?

texts:
Using ArcGIS Pro 1.0.0 I have:

Created a new project using Global Scene.aptx
Switched the Basemap to National Geographic
Navigated to and created Bookmarks for some countries of different sizes: Australia and Singapore
Used the project's Options to set the Transition Time to 2.0 seconds on the Navigation tab
Used the Bookmarks to navigate alternately between Australia and Singapore

I expected that doing this would result in a display cache being created that would mean no re-rendering and a very smooth display during transition.

However, each time I use a bookmark to go to the other country there are a few seconds of blurry re-rendering that occur.

Are there settings available that should assist me to get the result I describe as expected above?


texts: I'm working with QGIS and I'm trying to geo-reference a tobin lease map.  The map datum: NAD27, projection: UTM Zone 14N.  The process to georeference is fairly straight forward but the results I'm getting is off.  The program is basically mirroring the image. I have tried it three times with the same results.  I'm not sure what I am doing wrong.


texts: Adding layer using QgsMapLayerRegistry within thread PyQGIS Plugin

texts:
I am developing Python QGIS plugin, where I am opening two threads in run function with one thread performing calculation while second thread is printing value looping and when first thread is completed i.e. when self.letsGo variable value changes to True, second thread should add the layer to the QGIS. Everything else is working perfectly fine i.e. it is printing the values in the log panel of QGIS, but when self.letsGo is changed to True, second thread detects it but does not add the layer to map. Is there any dependency of QgsMapLayerRegistry from threading and how can I fix it?

def addLayertoMap(self):
    #add an empty layer to map
    vectorLayer = QgsVectorLayer('Polygon?crs=epsg:3857', 'testLayer' , 'memory')

    while (True):
        QgsMessageLog.logMessage("Looping")
        if self.letsGo==True:
            QgsMessageLog.logMessage("Detected")

            # QgsMapLayerRegistry.instance().addMapLayer(vectorLayer)

            QgsMapLayerRegistry.instance().addMapLayer(vectorLayer)

            QgsMessageLog.logMessage(str(lyr2.name()))
            for lyr in self.iface.mapCanvas().layers():
                QgsMessageLog.logMessage(lyr.name())
            QgsMessageLog.logMessage("Removed")

            break;
        time.sleep(0.5);

def ouput(self):
    time.sleep(2)
    #connect to database and fetch results
    self.letsGo=True;

def run(self):
    t = threading.Thread(target=self.output)
    t2 = threading.Thread(target=self.addLayerToMap)
    t.start()
    t2.start()

UPDATE

I just figured out that QgsMapLayerRegistry  is adding layer to Project properties Identify Layer but not adding it to Layers Panel.


texts: ArcMap topology toolbar validate topology button is disabled?

texts:
I find a strange thing when I try to use validate topology button. Most of the buttons on topology toolbar are disabled.

Last month when I was using a trial license of ArcGIS, I can use these functions.

From the instruction I find these functions only not available for the basic license. But I am using a standard license. Cannot figure out the reason.

The following pictures are the license information and the topology toolbar during editing.


texts:
I have account in qgis cloud. When i login in qgis cloud from qgis, i received an error which is shown below

Can some give the solution for my problem ?


texts: How to create a map of Denmark with ggmap

texts:
I have installed the following two packages:

library(ggmap)
library(maps)

I have found code online to create a world map (see below). How do I zoom the plot on a country level? For example, Denmark.

ggplot()+
  borders("world", colour="gray50", fill="gray50")


texts: Packaging failed in ArcGIS Server 10.2.1

texts:
I use ArcGIS Server 10.2.1. on Windows 8 64 bit.

I'm getting the simple error: "Packaging failed" for each new service I try to publish. Any thoughts what might be causing this?

Thanks


texts: Seeking documentation on DEMProcessing from GDAL?

texts:
I am using python 3.X and I am transforming a DEM (from a raster GeoTIFF file) into a slop tif file.
I am using GDAL and rasterio to do that.

Specifically, I am using the method:

  gdal.DEMProcessing(output, input, 'slope')

My problem is I do not know what is the mathematics and what exact method GDAL uses to calculate the slope. I can't find the answer in the GDAL documentation.

Can anyone can point me to that?


texts: Changing ledger line color and size using QGIS

texts:
I am using QGIS 3.4.

I want to change the color and size of the ledger line when drawing an object in a line layer to make it more visible if the tracing-function works correct.

It's currently very thin and red. The dotted part is from the last set point to the cursor.


texts: File Geodatabase XY Resolution and XY Tolerance

texts:
Here - http://gisintelligence.wordpress.com/2011/04/19/limiting-xyz-geometry-drift/ I found the description of XY Resolution and XY Tolerance in "Arcgis World".

I have a question - my cadastral data in shapefile format have accuracy 1 cm. So, if I want to work with my data in Arcgis File Geodatabase I should use XY Tolerance 0,01 m and XY resolution 0,001 m ( http://goo.gl/uqDrG

is that ok? I mean - that values are proper for topology matters and so on?


texts: How to to fix writeRaster error ("trellis")?

texts:
How do I export a .tiff file for further editing in QGIS
I used the following code to produce the added map.
I get the added error when trying to export using writeRaster?
library(RColorBrewer)
my.palette <- brewer.pal(n = 7, name = "YlGn")
spplot(AGCmap1, "AGCmap_1", main = "Carbon density per ha", zlim=c(0,80),
       col.regions = my.palette, cuts = 6, col = "transparent")

writeRaster(AGCmap1, "AGCmap1.tif")
Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘writeRaster’ for signature ‘"trellis", "character"’


texts:
I'm trying to port 3D rendering into QGIS through VTK/pyvista renderer.
Everything works fine ... untill I quit QGIS. There vtk QGIS crash, but I can't figure out why ...
I've issued the problem at pyvista but still nothing...
How to reproduce:
how to reproduce (inside QGIS Python console):
import subprocess
from qgis.PyQt.QtCore import QStandardPaths

interpreter = QStandardPaths.findExecutable("python")
command = (
    f"{interpreter} -m pip install pyvista --user"  # --user so it works even without admin rights (default windows install)
)
# run qgis-python pip
subprocess.check_call(command, shell=True)

import pyvista as pv

then quit QGIS, it should crash ...


texts: How to remove spikes in polygons with ArcGIS 10/Python?

texts:
In the picture below you can see a polygon layer, representing steep slopes, extracted from a TIN. The polygons have a lot of spikes in them. I would like to remove these spikes while the rest of the borders stay on the exact same location. I would like to do this in a Python script, using arcpy or Python functions.

Further info: ArcGIS 10.0, Python 2.6.5, polygon layer is in a GDB


texts: Python script to implement a batch import of CAD files into Geodatabase?

texts:
Does anyone know where I can get a python script that can do a batch import of CAD files into a Geodatabase.

I have one I wrote here, but it is only able to import one CAD file at a time. It imports the CAD files and split it into feature classes and stores it in the geodatabase via the feature dataset.
It works perfectly well but I need to loop the script to search through a given folder for all the cardfiles in it and import them under different feature datasets and put all the datasets into one geodatabase.

# Name: CadtoGeodatabase.py
# Description: Create a feature dataset
# Author: Irene
# Import system modules
import arcpy
from arcpy import env

# Set workspace
env.workspace = "C:/data1"

# Set local variables
input_cad_dataset = C:\Users\iegbulefu\Documents\info\91036c01.dwg" output_gdb_path = "c:/data/cadfile.gdb"
output_dataset_name = "cadresults"
reference_scale = "1500"
spatial_reference = "Nad_1983_10TM"

# Create a FileGDB for the fds
arcpy.CreateFileGDB_management("C:/data1", "cadfile.gdb")

# Execute CreateFeaturedataset
arcpy.CreateFileGDB_management("C:/data", "cadfile.gdb")
arcpy.CADToGeodatabase_conversion(input_cad_dataset, output_gdb_path, output_dataset_name, reference_scale)


texts:
I am trying to acquire a method of detecting sub cm elevation without manually measuring from the field.  Can LIDAR data produce this for me?  Example:  A 1500 foot long pipe conveying water from one point to another.  I need to know the rise and falls of the terrain where the pipe will be lying.  Can LIDAR produce this for me?  I am needing coverage for agricultural land throughout the Delta.

Thanks.


texts: I have a point cloud and  need to process a very detailed boundary around the outer points. Is there a tool in QGIS? Convex hull is not detailed enough.


texts: Leaflet MultiLineString Popup and Tooltip Display Issue

texts:
I have a geoJSON layer of watercourses for which I want the multi line strings and also tooltips and popups to appear in a Leaflet map on instantiation.

The below code will not display popups or tooltips, but the watercourses do appear:

    function waterPops(feature, layer){
      layer.bindPopup(feature.properties.NAME)
      layer.bindTooltip(feature.properties.NAME)
    };

    var lineStyle = {
        "color": "#00b3ff",
        "weight": 1,
        "opacity": 1
    };

    var water = L.geoJSON(waterExtract, {
      style: lineStyle,
      onEachFeature: waterPops,
    }).addTo(mymap);

Something I tried with a strange result: If I load the above water variable into a layer group WITHOUT addTo(mymap) for var water, obviously the watercourses layer doesn't load on instantiation, but when I turn the layer on in the layer group the popups and tooltips will work, along with the watercourses. Similarly, if I load the water variable into a layer group WITH addTo(mymap) for var water, on instantiation the watercourses appear, but tooltips and popups do not work. If I turn the layer off then on again, the tooltips and popups will work and watercourses will appear.

Popups, tooltips, and watercourses will work for the below code if I turn the layer on in the layer control.

    function waterPops(feature, layer){
      layer.bindPopup(feature.properties.NAME)
      layer.bindTooltip(feature.properties.NAME)
    };

    var lineStyle = {
        "color": "#00b3ff",
        "weight": 1,
        "opacity": 1
    };

    var water = L.geoJSON(waterExtract, {
      style: lineStyle,
      onEachFeature: waterPops,
    });

    var baseLayers = {
      "Map": streets,
      "Satellite": googleSat
    };

    var overlayMaps = {
      "Thing": thing,
      "Watercourses": water,

    };

    L.control.layers(baseLayers, overlayMaps).addTo(mymap);

For the code below, watercourses appear, but popups and tooltips will not work on instantiation. If I turn the layer off then on in the layer control, then everything works:

    function waterPops(feature, layer){
      layer.bindPopup(feature.properties.NAME)
      layer.bindTooltip(feature.properties.NAME)
    };

    var lineStyle = {
        "color": "#00b3ff",
        "weight": 1,
        "opacity": 1
    };

    var water = L.geoJSON(waterExtract, {
      style: lineStyle,
      onEachFeature: waterPops,
    }).addTo(mymap);

    var baseLayers = {
      "Map": streets,
      "Satellite": googleSat
    };

    var overlayMaps = {
      "Thing": thing,
      "Watercourses": water,

    };

    L.control.layers(baseLayers, overlayMaps).addTo(mymap);

Ideally I would like to have the watercourses, popups, and tooltips working on instantiation, whether in a layer control group or alone.


texts: Python Console XML creation - Not releasing memory?

texts:
I have created a module to prepare an XML file. Since xml.etree does not output a pretty XML I use dom.

****Module:****
    from xml.etree import ElementTree
    from xml.dom import minidom
    import gc

 def prettify(elem):
        """Return a pretty-printed XML string for the Element.
        """
        rough_string = ElementTree.tostring(elem,"utf-8")
        reparsed = minidom.parseString(rough_string)
        tempString = str(reparsed.toprettyxml(indent="  "))
        reparsed = None
        return tempString

import xml.etree.ElementTree as ET

try:
    class source():
        def __init__(self):
            self.top=ET.Element("nrml",{"xmlns:gml":"http://www.opengis.net/gml","xmlns":"http://openquake.org/xmlns/nrml/0.4"})
            self.subtop=ET.SubElement(self.top,"sourceModel")
        def clean(self):
            self.top=None
            self.subtop=None

    class simpleFaultSource():
        sourceInstance = source()

        def __init__(self, **kwargs):
            for k, v in kwargs.items():
                setattr(self, k, v)
            self.simpleFaultSourceAttributes={}
            self.simpleFaultSourceAttributes["id"]= self.id
            self.simpleFaultSourceAttributes["name"]=self.name
            self.simpleFaultSourceAttributes["tectonicRegion"]="ActiveShallowCrust"
            projectName = self.filename
            simpleFaultSource.sourceInstance.subtop.set("name", projectName+":"+self.layername )
            self.createTags()
            self.addText(**kwargs)

        def createTags(self):
            self.simpleFaultSource= ET.SubElement(simpleFaultSource.sourceInstance.subtop, \
      "simpleFaultSource",self.simpleFaultSourceAttributes)

    except Exception as e:
        print("Printing at module!")
        print(e.message)
        print(e.__class__.__name__)
    finally:
        gc.collect()

Once I execute the main file once, the output file is correct.

****Main:****
from qgis.core import *
from PyQt4.QtCore import *
import processing
from qgis.utils import iface

import getpass
username = getpass.getuser()
import sys
exePath=r'C:\Users\10011\Desktop\Control\qgis\02_JG_Summit\02_scripts'
sys.path.append(exePath)
import oq_classes

try:

    layer = iface.activeLayer()
    #layer = processing.getObject(input_layer)
    featureList=[]
    prjPath=QgsProject.instance().fileName()
    splitPath=prjPath.split("/")
    filename=splitPath[-1]
    path="/".join(splitPath[:-2])
    savePath=path+"/05_output/"

    if layer.geometryType() == QGis.Line:
            layerType = 1
            LT = "Line"
        for f in layer.getFeatures():
            g = f.geometry()
            l = g.asPolyline()
            temp=""

            for e in l:
                temp = temp + "\n"+str(e[0])+" "+str(e[1])
            temp=temp+"\n"
            kwargs={"name": f['name'], "id" : f['id'], "pos" : temp,"layername":layer.name(),"filename":filename}
            tempInstance=oq_classes.simpleFaultSource(**kwargs)
            featureList.append(tempInstance)
        newPath=savePath+oq_classes.simpleFaultSource.__name__+".xml"
        string1=oq_classes.prettify(oq_classes.simpleFaultSource.sourceInstance.top)

        with open(newPath,"w") as f:
            f.write(string1)
except Exception as e :
    print(e.message)
    print(e.__class__.__name__)

finally:
    gc.collect()

Every time I run the script again, the children node (simpleFaultSource) gets duplicated under the parent node (sourceModel) unless I close QGis and open it again. Summarized output is shown below, everything below the red line is duplicated within the same parent node when executing the script again .

My understanding is that the import, classes and instances should not be in memory after the script has finished execution. Is this an issues with my implementation, Python or with QGis (version 2.14)?


texts: Finding point of intersection between: a line segment of two LatLngs, and a line with initial point and bearing

texts:
I have a line segment AB made up of two LatLngs, A and B.

I have another LatLng C and a bearing x degrees.

How can I calculate the point D where a line from C with bearing x degrees meets the line segment AB, if at all it meets?

All the lines and line segments are geodesic and are not straight.


texts: With just a web browser we want our field staff to collaboratively make complex maps with secured private data--looking for options?

texts:
Simple Question

With just a web browser we want our field staff to collaboratively make complex maps with secured private data.

What are our options?

Note: Please provide specific product examples to help me and others start the research process and ensure we dont overlook a great solution.

I will give a concrete example in response to the comments ( @brad and @iant ) that my earlier version of the question above is too general.

Concrete Real-World Scenario Example

For example, consider there are 100 local people spread around a province in Afghanistan and they only have a web browser. The objective is to have them collaboratively making a map of the safety and danger zones in their province.

Requirements

They do NOT require advanced spatial analysis such a routing etc.

Nevertheless, the final map-data product needs to be structured (geoJSON, shapefile, etc) so some one can later use the data for more complex spatial analysis.

The key functional requirement is that these people need to upload to the collaboratively edited base map the following new features with associated labels:

point locations of violent incidents with a labels
polygon zones of dangerous and safe zones with labels
line roads with labels

More elaborate question

What are the specific products I should look into where the data will be secure, away from the public eyes --- ideally something that connects with Google's Datastore, which is approved for SBU data (sensitive but unclassified data) --- most likely some product where spatial data is stored and served internally ?


texts: Field calculating a unique date and time

texts:
I want to field calculate a field to be have a have the date and time (down to the second) and two letters after (to be assigned by user).

For example, if the record was created on September 8th, 2016 at 7:15 am and 03 seconds it would look like this: 090816071503XX

I also need this to be unique though, in some cases this field calculate will be performed on over 5000 records. By doing a complete field calculate in real-time, this means all 5000 records would show the same time. So I need a code that will field calculate a "fake", but realistic time (i.e. can't have hour 25). I don't care that this time is incorrect, but I do want the date to be correct.

I imagine this is something that can be scripted in Python?

I have a VBA (.cal) script I can load into the field calculator which will calculate the date and time, but as I said, I need it to return unique times.

'enter your initials
INIT="XX"

'month string - add zero
If LEN(DatePart ("m", Now (  ) ))=1 then
    MO = "0"&DatePart ("m", Now (  ) )
Else
    MO = DatePart ("m", Now (  ) )
end if

'day string - add zero
If LEN(DatePart ("d", Now (  ) ))=1 then
    DA = "0"&DatePart ("d", Now (  ) )
Else
    DA = DatePart ("d", Now (  ) )
end if

'year string - two digits
YR = RIght(DatePart ("YYYY", Now (  ) ),2)

'HR string - two digits
If LEN (DatePart ("h", Now (  ) ))=1 then
    HR = "0"&DatePart ("h", Now (  ) )
Else
    HR = DatePart ("h", Now (  ) )
End if

'MN string - two digits
If LEN (DatePart ("n", Now (  ) ))=1 then
    MN = "0"&DatePart ("n", Now (  ) )
Else
    MN = DatePart ("n", Now (  ) )
End if

'SS string - two digits
If LEN (DatePart ("s", Now (  ) ))=1 then
    SS = "0"&DatePart ("s", Now (  ) )
Else
    SS = DatePart ("s", Now (  ) )
End if


texts: Again Band swap when reading created image with Gdal

texts:
As previously seen here I created a function to save a raster with gdalon python:

def array2raster2(fname,  rasterOrigin, pixelWidth, pixelHeight,array,  proj, shape):

    originX = rasterOrigin[0]
    originY = rasterOrigin[1]
    drv = gdal.GetDriverByName("GTiff")
    dst_ds = drv.Create(fname, shape[1], shape[0], 4)#, gdal.GDT_UInt16)
    dst_ds.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    dst_ds.SetProjection(proj)
    dst_ds.GetRasterBand(1).WriteArray(array[:, :, 0])  # write r-band to the    raster
    dst_ds.GetRasterBand(2).WriteArray(array[:, :, 1])  # write g-band to the raster
    dst_ds.GetRasterBand(3).WriteArray(array[:, :, 2])
    dst_ds.GetRasterBand(4).WriteArray(array[:, :, 3])

    dst_ds.GetRasterBand(1).SetRasterColorInterpretation(gdal.GCI_RedBand)
    dst_ds.GetRasterBand(2).SetRasterColorInterpretation(gdal.GCI_GreenBand)
    dst_ds.GetRasterBand(3).SetRasterColorInterpretation(gdal.GCI_BlueBand)

    dst_ds.FlushCache()
    dst_ds=None

And added the last 3 lines for SetRasterColorInterpretation()and worked great when using 3 channels.

But now I am using now 4 channels, 4Th channel for NIR band, as on SPOT satelite images or Pleiades, but when opening again on ENVI or QGIS the 4th band is shown as the first, what do I need to do to fix that? I checked on gdal.GCI_values and only shows for 3 bands, I understood that gdal.GCI_AlphaBandis for transparency band so I don't need that...yet.

And second issue, I have the borders of the image as 0, how can I set that value to Nodata so I don't have borders visualized on ENVI or any other GIS?


texts:
I am wondering if it is possible in ArcMap to isolate a particular subtype within a feature class to be visible in the map.  I have a feature class with several subtypes based on an attribute field.  I would like to make it so that I can turn these subtypes on or off with a checkbox in the TOC the way that you can with layers or groups.

This would essentially be the next level down in the tree:
Layers->Groups->Feature Classes->Subtypes


texts:
I recently completed a geostatistics course, and now trying to apply my knewly acquired knowledge to a real case.
I have a dataset with 2 122 groundwater level measurement points. I would like to extrapolate them into a raster map. The previous group that worked on this matter used some inverse distance weighting. But I thought kriging might make more sense. Which is why I decided to try it out and perform some LOOCV to compare methods.
I do not have any NA in the original dataset, but when using function krige.cv(), I get 85 NAs. This messes up with computing zscore statistics.
I would like to know why I have these NA values, and whether it is fine to simply remove them.
This is my call of the function:
GLGkrige1 <- krige.cv(GLG~1, ground_water_points, model=vgmGLG, nmax=100)

And this is the summary of the result:
Object of class SpatialPointsDataFrame
Coordinates:
                min       max
coords.x1  4.749619  5.017867
coords.x2 52.289492 52.425470
Is projected: NA
proj4string : [NA]
Number of points: 2122
Data attributes:
   var1.pred          var1.var         observed         residual            zscore               fold
 Min.   :-5.2735   Min.   :0.0000   Min.   :-5.630   Min.   :-1.77453   Min.   :-25624692   Min.   :   1.0
 1st Qu.:-1.9587   1st Qu.:0.2195   1st Qu.:-2.040   1st Qu.:-0.12146   1st Qu.:        0   1st Qu.: 531.2
 Median :-0.8083   Median :0.2301   Median :-0.810   Median : 0.00069   Median :        0   Median :1061.5
 Mean   :-1.3896   Mean   :0.2390   Mean   :-1.463   Mean   : 0.00135   Mean   :        0   Mean   :1061.5
 3rd Qu.:-0.4698   3rd Qu.:0.2467   3rd Qu.:-0.460   3rd Qu.: 0.12386   3rd Qu.:        0   3rd Qu.:1591.8
 Max.   : 0.7157   Max.   :1.7875   Max.   : 1.590   Max.   : 3.55142   Max.   : 25624692   Max.   :2122.0
 NA's   :85        NA's   :85                        NA's   :85         NA's   :85

I just noticed that my first run on the krige.cv() where I forgot to set a nmax has 1488 NAs, so I suspect this has to do with the nmax.
Could these be points that do not have enough neighbours?

The problem came from two points sharing the same coordinates (and different values). Now everything is running smoothly again.


texts: How can Flocking Algorithms be used in GIS?

texts:
The first flocking algorithm was written by Craig Reynolds in 1986. A recent browser implementation can be found here.

The notes refer to "The use of Flocks to drive a Geographic Analysis Machine" (1998) by James Macgill and Stan Openshaw which

  uses a flock model, with communication
  between boids, to better search for
  clusters in spacial datasets, just as
  a natural flock provides better
  foraging than individual birds could
  manage.

What other real world examples of this algorithm have been used in GIS, and do any GIS applications implement this functionality?


texts: How to count the number of entries a user has selected in a multivalue field of a tool dialogue using Python and ArcGIS 10.2.1

texts:
I am currently creating a script that will read data in one (to many) Excel spreadsheet(s) and would like to know how to get a count of the number of spreadsheets entered by the user. The script will only read the data in the GISData$ sheet of each spreadsheet file. So far I have the below code but need something to replace 'rowCount' as this doesn't seem to be working. I would like to just get a count of the number of spreadsheets selected.

spreadsheets = arcpy.GetParameter(1)

for idx in range(0, spreadsheets.rowCount):
    input_spreadsheet = str(spreadsheets.getValue(idx,0))
        input_spreadsheet += r"\GISData$

...do something...

Thanks in advance.


texts: ArcGIS REST API ElevationSync result is empty?

texts:
I would like to use the ArcGIS REST elevation service.

Providing a list of locations I would like to receive a list of elevations relative to these locations.

After some struggling with the syntax I was able to create a working (synchronous) request:

https://elevation.arcgis.com/arcgis/rest/services/Tools/ElevationSync/GPServer/Profile/execute?InputLineFeatures={
   "displayFieldName":"",
   "geometryType":"esriGeometryPolyline",
   "spatialReference":{
      "wkid":4326
   },
   "fields":[
      {
         "name":"OID",
         "type":"esriFieldTypeOID",
         "alias":"OID"
      },
      {
         "name":"Shape_Length",
         "type":"esriFieldTypeDouble",
         "alias":"Shape_Length"
      }
   ],
    "features": [
    {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "type": "MultiPoint",
        "coordinates": [
          [
            52.50,
            13.30
          ],
          [
            52.60,
            13.40
          ],
          [
            52.70,
            13.50
          ],
          [
            52.80,
            13.60
          ]
        ]
      }
    }
  ],
   "exceededTransferLimit":false
}&ProfileIDField=&DEMResolution=&MaximumSampleDistance=&MaximumSampleDistanceUnits=&env:outSR=&env:processSR=&returnZ=true&returnM=true&f=json

But unfortunately the response result is empty

{
   "results":[
      {
         "paramName":"OutputProfile",
         "dataType":"GPFeatureRecordSetLayer",
         "value":{

         }
      }
   ],
   "messages":[

   ]
}

What am I doing wrong?


texts: Get Result of ICalculator without actually populating a field (ArcObjects)

texts:
I have two feature classes: FeatureClassA and FeatureClassB.  I am running a cursor on FeatureClassA and I have some logic statements (spatial and attribute logic) that compares the current feature in FeatureClassA to FeatureClassB and returns one matching record in FeatureClassB.  I want to then do a field calculation on a field in FeatureClassA using the fields in the returned record from FeatureClassB.

To achieve this, I have the user build a field calculator expression using the ICalculatorUI2 interface.  This dialog is familiar to ArcGIS users.

private void btn_Click(object sender, EventArgs e)
{
    ICalculatorUI2 calcUI = new CalculatorUIClass();
    calcUI.Table = FeatureClassB as ITable;
    calcUI.DoModal(ArcMap.Application.hWnd);
    string expression = calcUI.Expression;
    calcUI = null;
}

If I wanted to calculate a field in FeatureClassB, I could just use ICalculator.

ICalculator calculator = new Calculator();
calculator.Expression = calc.Expression;
calculator.Field = "MyField";
calculator.Cursor = calcUI.Table.Search(null, false);
calculator.Calculate();

The problem is, I don't want to actually calculate a field in FeatureClassB.  I want the result of the field calculation to put into a field in FeatureClassA.  Is there anyway to get the result of the field calculation?


texts: Many of my Fusion Tables are missing off my Google Drive account, has anyone else experienced, and hopefully resolved, this issue?


texts: Highlighting feature by Longitude and Latitude in OpenLayers

texts:
I used Qgis2Web and I'm really a beginner. Is there a way can I highlight a feature using their coordinates?
I'm currently developing a website which suggests location based from pre-defined inputs. And the only thing i have is the lat and lon. I hope these helps :)


texts:
I am trying to evaluate an equation in the raster calculator and I need latitude values. I have the DEM for the region and am wondering if it is possible to create a latitude raster based on the DEM?

Thanks a lot!


texts: How to construct a WFS 2.0 GetCapabilities XML request to Geoserver?

texts:
I am trying to use WFS 2.0 specifically because it supports temporal filtering. My GetCapabilities.xml file is

<GetCapabilities
version="2.0"
service="WFS"
xmlns="http://www.opengis.net/wfs"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:schemaLocation="http://www.opengis.net/wfs
http://schemas.opengis.net/wfs/2.0/wfs.xsd"/>

I use curl to make the request i.e.

curl -X POST -H "Content-type: application/gml+xml;version=3.2" -d @GetCapabilities.xml http://ws-beta.int.unavco.org:9090/geoserver/wfs

But I receive the error:

<ows:ExceptionReport xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:ows="http://www.opengis.net/ows/1.1" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="2.0.0" xsi:schemaLocation="http://www.opengis.net/ows/1.1 http://ws-beta.int.unavco.org:9090/geoserver/schemas/ows/1.1.0/owsAll.xsd">
  <ows:Exception exceptionCode="NoApplicableCode">
    <ows:ExceptionText>java.lang.IllegalArgumentException: java.lang.ClassCastException@44895b09
java.lang.ClassCastException@44895b09</ows:ExceptionText>
  </ows:Exception>
</ows:ExceptionReport>

I know the problem is almost certainly in my construction of the XML file because it works over HTTP:

http://my.geoserver.url/geoserver/wfs?service=wfs&version=2.0.0&request=GetCapabilities

I seek a working version of GetCapabilities.xml for WFS 2.0.

I have also tried using GET rather than POST, i.e.

curl -X GET -H "Content-type: application/gml+xml;version=3.2" -d @GetCapabilities.xml http://ws-beta.int.unavco.org:9090/geoserver/wfs

I get a different error in this case:

<ows:ExceptionReport xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:ows="http://www.opengis.net/ows/1.1" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="2.0.0" xsi:schemaLocation="http://www.opengis.net/ows/1.1 http://ws-beta.int.unavco.org:9090/geoserver/schemas/ows/1.1.0/owsAll.xsd">
  <ows:Exception exceptionCode="MissingParameterValue" locator="request">
    <ows:ExceptionText>Could not determine geoserver request from http request org.geoserver.platform.AdvancedDispatchFilter$AdvancedDispatchHttpRequest@66acc6c5</ows:ExceptionText>
  </ows:Exception>
</ows:ExceptionReport>


texts: How can I export symbology and labels with QGIS Leaflet to webmap?

texts: If I make a map in QGIS export it with Leaflet the symbology and labels don't appear on the webmap. Is it possible to export them? If yes, how to export the all map?


texts: Is there a way to bring ArcGIS 10.3 tools/toolbars into ArcGIS Pro?

texts: I am currently working on a project that involves using Network Analyst, and I really dislike the way that it is set up in ArcGIS Pro vs. the toolbar from ArcGIS 10.3. I was wondering if there is any way that I can incorporate the ArcGIS 10.3 Network Analyst Toolbar, as well as any other tools/toolbars into ArcGIS Pro. If this requires using Python I am a little new to it, but exposed to say the least and willing to use it as a solution. Any ideas?


texts: Leaflet mouseover popup makes whole map move

texts:
I was following a couple of examples to get my leaflet map to show the popup for the marker when I hover over it.

This example to show popup on hover.

I also used example to do leaflet marker cluster groups very heavily.

For some reason, when I hover over the marker to show the popup, the popup box is sized to hold the popup data and that box size seems to make the whole map move, instead of positioning the popup info within the map box.  How do I keep the popup box from moving the whole map?

This is the code to show the marker popup info, which is very close to leaf-demo from the leaflet marker cluster demo mentioned above, and GitHub info here:

var markerClusters = L.markerClusterGroup();

for ( var i = 0; i < markers.length; ++i )
{
  var popup = markers[i].name +
              '
' + markers[i].city +
              '
IATA/FAA: ' + markers[i].iata_faa +
              '
ICAO: ' + markers[i].icao +
              '
Altitude: ' + Math.round( markers[i].alt * 0.3048 ) + ' m' +
              '
Timezone: ' + markers[i].tz;

  var m = L.marker( [markers[i].lat, markers[i].lng], {icon: myIcon} )
                  .bindPopup( popup );

  m.on('mouseover', function (e) {
     this.openPopup();
  });
  m.on('mouseout', function (e) {
     this.closePopup();
  });

  markerClusters.addLayer( m );
}

map.addLayer( markerClusters );


texts:
I am trying to write the following SQL command using case clause.

with x as
(select the_geom as geom from table1
where st_dwithin(the_geom, st_setsrid(st_makepoint(28.959495,41.019913), 4326), 0.00012) limit 1)
select case count(geom)
when 0
    then st_setsrid(st_makepoint(28.959495,41.019913), 4326)
when 1
    then st_closestpoint(x.geom, st_setsrid(st_makepoint(28.959495,41.019913), 4326))
end as the_geom
from x group by geom

As can be seen under the with clause, there could be either one row or no row in the record x. Code is working fine where there is one row in x but when the x is null (means count(geom) = 0) then the when 0 operation is not getting executed. I am getting the following output.

Plus, after modifying my SQL code a bit like the following (only the when 2 is changed!)

with x as
(select the_geom as geom from table1
where st_dwithin(the_geom, st_setsrid(st_makepoint(28.959495,41.019913), 4326), 0.00012) limit 1)
select case count(geom)
when 0
then st_setsrid(st_makepoint(28.959495,41.019913), 4326)
when 2
then st_closestpoint(x.geom, st_setsrid(st_makepoint(28.959495,41.019913), 4326))
end as the_geom
from x group by geom

The output is like the following screenshot when the count(geom) = 1.

Queries:

How to run the when 0 condition when the record x is empty?
What is the difference between the two SQL outputs (as shown in the screenshots)?


texts: How get the attributes list from a feature in OpenLayers 3

texts:
I'm getting a vector layer from a WMS server, but I don't know which fields are included in it.  In previous version feature object had a attributes collection but I can't find it now.

Is there a way of searching for feature attributes?

Right now I can use Feature.get('Name') and can get the value, but only if I know 'Name' is an attribute. Now imagine I do find a [Name] field or something similar.

Is possible do a  getFeaturesByAttribute('Name', strName);

Also available in previous version. But can't find it in OpenLayers 3.


texts: Multiply whole column with the field calculator in QGIS

texts:
I have a layer with the field size. I did some mistake with the value and now I need the values multiplied by 16.

But instead of going through 500 lines and recalculating them I would like QGIS to do that.

So I tried to mark all features and then I hit the field calculator. Now I selected to take the values and chosen the size column. In the expression I had $Id * 16.

But the result was that I got:
0
16
36

instead of e.g.:
47
38
67

So how can I take a whole column with simple operations?


texts:
It seems like old versions of OpenLayers supported SLD format. There was a construct like:

var format = new OpenLayers.Format.SLD();

But I cannot find it in modern OpenLayers library. I checked ol.format and see that there is no ol.format.SLD().  Why is that? Was it replaced by another format (some which is better than SLD)?


texts: Polygonize.py outputs into a single file or simple merge?

texts:
Given an unique raster topographic file at start, I want to get vector polygons for elevations thresholds 1m, 200m, 500m, 1000m, 2000m into a single shp file. I first slice my raster into 5 (or more) rasters, then polygonize each, and then have to merge them.

How to get these outputs into a single shape file ? Either by sending them as layers directly into the same file, or by creating individual shp files to merge into one final shp.

Current workflow details: It seems I must do this because poligonize works on pixels with a same value (elevation) x.

1. Generate raster slices:
I used gdal_calc.py to generate kinds of vertical slices, with one .tif for each level :

gdal_calc.py -A crop.tif --outfile=level001.tif --calc="1*(A>0)"     --NoDataValue=0
gdal_calc.py -A crop.tif --outfile=level200.tif --calc="200*(A>200)" --NoDataValue=0
...
...
...
...

All conserved pixels in level001.tif have an elevation of 1.  All conserved pixels in level200.tif have an elevation of 200m. Opened together in QGis it give slices like this:

2. Polygonize my rasters slices (files):

I run polygonize.py separatedly on each of these several raster files to get my vector layers:

gdal_polygonize.py level001.tif -f "ESRI Shapefile" levels001.shp levels001 elev
gdal_polygonize.py level200.tif -f "ESRI Shapefile" levels200.shp levels200 elev
...
...
...
...

3. Merge into one shp file.

  How to merge these shp outputs ?

Alternatively, I'am also interested by smoother workflow (loop).


texts: QGIS 3.4.10 for Mac users...Digitzing polygons

texts:
I have just installed the LTR version of QGIS (3.4.10) and am having trouble telling QGIS when I have finished a polygon. In fact, I can't.

It used to be in earlier versions that to end a polygon you would type ctrl-click. This is the equivalent to right-click in PCs.

It doesn't work at all in QGIS3.4.10.

Has anybody else encountered the same problem? What is the solution?

[I have tried ctrl-click in other apps such as word, and it functions as it should. I have also completely rebooted the computer but that doesn't help. I am running OS Mojave--10.14.6]


texts: Postgis st_intersects unexpected results

texts:
Trying to update one table of parcels with st_intersects from another.
Results always update all table records and don't discriminate using st_intersects.
Aim of query is to update road_class value with 1 for parcels closer than 100 meters to roads and with 2 for parcels between 100 and 800 meters of roads.
Distance selections like st_dwithin gives very slow queries.

Results of this queries are always to all parcels get same value and results doesn't discriminate for being overlayed or not.

What's wrong?

Both geometries uses same srid, 4326

And Postgresql versions are:
-PostgreSQL 10.0, compiled by Visual C++ build 1800, 64-bit
-POSTGIS="2.4.0 r15853" PGSQL="100" GEOS="3.6.2-CAPI-1.10.2 4d2925d" PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 2.2.2, released 2017/09/15" LIBXML="2.7.8" LIBJSON="0.12" LIBPROTOBUF="1.2.1" RASTER

UPDATE parcels
SET
road_class=2
FROM parcels p INNER JOIN road2 r ON
st_intersects(p.geom,r.geom);

UPDATE parcels
SET
road_class=1
FROM parcels p INNER JOIN mi.r1 r ON
st_intersects(p.geom,r.geom);

Same results with:
FROM parcels p,r1 r
WHERE
st_intersects(p.geom,r.geom);

Screenshot:
parcels in brown.

r1 is diagonal polygon -road buffer of 100 m.

r2 is disgonal cross pattern -road buffer since 800 meters-.

Selected parcel should get value 2, rest value 1


texts:
I just upgraded to QGIS 2.18 Las Palmas. I exported an SVG file from Print Composer, as usual, but the SVG file crashes Adobe Illustrator. In fact, it locks up my entire system.

The SVG file is only 11.5 MB. This has never happened to me before, and it has brought my professional workflow to a standstill.

UPDATE:

I tried it in the old 2.8 Wien version and got the same problem. It appears the offending layer was a shapefile for water polygons. Too many nodes I guess. The SVG file imports fine without that layer.

I even tried exporting the water polygon layer by itself, but Illustrator still could not open it. That's pretty sad. I'm just surprised that this is the first time it's happened to me.


texts:
I have a saved webmap containing layers pointing at an ArcGIS Server which is now switched off. I need to update webmap, changing the layer URLs to a new ArcGIS Server URL. Layers were added straight to the map using Add Layer From Web. I need the old webmap so I retain styling and popup configuration. My web map ID is also referenced in a JS app. Apparently I can edit the URLs for layers and apps, but not maps. Why? Anyone know a way to edit the layer URLs?

From esri help:
https://doc.arcgis.com/en/arcgis-online/share-maps/item-details.htm

"URL (Apps and ArcGIS Server layers only)—Changing the URL of an ArcGIS Server layer is only recommended when you move the service to a different server. If you've changed the service—for example, added or removed layers—you should add the service as a new item. Any map that uses the moved layer will appear broken because it references the old URL. To fix the map, delete the layer, add it again with the new URL reference, and save the map. Any layer properties you saved to the item are still available. If you edit the URL of a changed service, the map may not display correctly, even after you readd the layer. You cannot edit the URL of a hosted feature or tile layer.


texts: Mac OS X 10.6.8 QGIS not using user-installed Python 2.7

texts:
My QGIS does not recognize the recently installed python 2.7 on my mac osx 10.6.8. it still uses the python version 2.6 that is automatically installed by apple.

Opening qgis I always receive an error message telling me that the module importlib is not found and prevents the processing plugin from working.

How can I tell QGIS to use the installed python version 2.7 (In /Library/Python) instead of 2.6 (In /System/Library/Python)?

QGIS was installed from Dakota cartography nightly builds.

This is different from OS X Python vs QGIS Python as it does not ask for the reason for the problem but for a solution.


texts: Creating a raster mosaic from multiple un-referenced jpegs

texts: I have multiple historical aerial photographs that I would like to use to compare coastal change over time in relation to urban development. The aerials I have obtained come in tiles that are not georeferenced but they are also not big enough to get a good RMSE when georeferencing them individually. So to fix this I would like to turn each of the individual tiles into one large mosaic to use for georeferencing. However, when adding them into Arcmap they each just layer on top of each other which does not help me at all. Is there anyway that I can manually (or otherwise) arrange them into their correct order so that I can then combine them into a single mosaic and subsequently accurately georeference?


texts: Does any body knows if there are a ArcSDM version or port for QGIS?

texts:
ArcSDM is a series of python analysis tools for ArcGIS developed by Don Sawatzky, Gary Raines y Graeme Boham-Carter in 2010. At present is maitained by the University of Campinhas, Brazil and Geological Survey of Finland.

ArcSDM or SDM (Spatial Data Modeller) "is a collection of tools for adding categorical maps with interval, ordinal or ratio sacle maps to produce a predictive map of where something of intereset is likely to occur." It includes tools for performing spatial analysis using Weights of Evidence, WofE and Logistic Regression, etc.

The original code can by seen at:

http://www.ige.unicamp.br/sdm/

or

http://newprojects.gtk.fi/mpm/ArcSDM5/

These series of tools are used mainly for geosciences, mining exploration, natural hazards, etc.

Does anybody knows if there are a project or port to QGIS?


texts: Why can't a user with db_owner role version data?

texts:
I've made an enterprise geodatabase and added a user with db_owner permissions.

When I add a database connection in ArcCatalog and connect to this geodatabase using the user I created, I'm not able to version the data.

When I right click on a table or dataset in the enterprise geodatabase in ArcCatalog the menu is all grayed out except for the "Locks..." selection.

I've added the user to other enterprise geodatabases with the same permissions and everything worked fine. The user was able to version data.

What could be wrong here?

If I sign into the geodatabase that isn't working with the SDE login, I am able to version the data.


texts: Setting lower and upper limits to classes in kriging interpolation, ArcGIS

texts:
I am using the kriging interpolation tool. I included two dummy values within my dataset:- an upper limit and a lower limit.

I ran the tool in arc map 10.1 software. However, I want to change the upper and lower limits of the kriging interpolation classes to the limits I set in the dataset. I will then   manually classify the kriging interpolation into specific classes.

How do I set the upper and lower limits for the interpolated classes?


texts: Why does raster still have pixels around edge when clipped?

texts:
I am trying to clip a raster in ArcMap 10.6 (derived from 'make nedcdf raster') to the shape of the Ireland polygon and I keep getting back pixels around the edge:

I have tried 'Clip' in 'Data management' (I have ticked 'Use input features for clipping geomotry' AND 'Maintain clipping extent').
I have also tried 'extract by mask' in Spatial analyst. None of them cut the raster with the Ireland Shapefile. I created this netcdf in R.
Could that be an issue?
Also, I'm wondering if it's because the raster doesn't have coverage for the entire shape of Ireland?


texts: How to switch rule-based styles in loop?

texts:
I'm a working on a QGIS script that could activate/deactivate rule-based styles of the vector layer in a loop (one active layer per iteration). On the every iteration this scenario should:

Activate one of the prearranged style rules for a given vector layer
and deactivate all  other;
Save current view as JPEG screenshot;

In my particular case the vector layer is a PostGIS table that contains the contours derived from some raster layer. Rules are based on the distinct values of the isolines. I take screenshot of the map for the every activated rule and finally put them all into the animated gif.

Since list of the rule-based styles can be very long, I'm looking torward the PyQGIS code snippet that could go through the one in a loop.


texts: ArcGIS Javascript API - Dynamic label placement in feature layer

texts:
Is it possible to use dynamic label placement with a feature layer in the javascript API v3.21 ? For instance, if two point features are superposed, I would like to display their labels at different placement.

The LabelClass documentation does not specify this possibility. Is there any solution to achieve this?

Otherwise, is it possible to set a priority to features for labeling. In case two points features are superposed, is it possible to define which feature's label will be displayer on the map?


texts: Prevent user from sliding the map to the right or left

texts: Is there anyway to prevent user from sliding the map to the extreme right or left?  something like setting a maximum/minimum zoom level in arcgis maps.


texts: Odd, old coordinate formatting from 1997 GPS

texts:
A typical coordinate is given as 68 35 067 N 123 54 817 W.

I thought maybe some old timers could help me out. The author says in his report of the time that they were using a GPS, but I contacted him and he doesn’t remember any of the details. The points can be located on a map using the first two pairs of two digits, so I’m pretty sure it’s degrees minutes, but the last three digits don’t make sense. Since some are greater than 600 it’s not seconds with another significant figure, i.e, a missing decimal point.

Any ideas? I’m thinking there was some “convention” early in the personal GPS receiver days. The original information is scattered throughout the text at http://traditionalarctickayaks.com/horton_trip_journal.htm

I'd like to use his info but without knowing how to interpret the last three digits I'm left with a ±2km precision.


texts:
How to disable anti-alias in QGIS using PyQGIS? I need to create output rasters (tiff / png) which has not used antialised using QGIS Composer / Atlas. Do you have some PyQGIS commands for novice (not programmer) which command need to use?

Indeed unclicking "Enable antialising" will create raster without antialising except texts and symbols, which will appear still antialised. I need to create whole raster without antialising. Also texts and symbols need to force to use without antialising.


texts:
I am trying to extract values from my rasters (slope, aspect, landcover etc.) for my animal group locations. I used Spatial Analyst > Extract Multi-values to Points, but when I look in my attribute table after the process the results are the same for all the points and make no sense. For example 0 for a landcover, -3x10-38 for a slope, etc.

Do you have an idea what I am doing wrong?

I am working with ArcGIS 10.


texts: Has anyone experienced duplication issues when exporting large datasets to .CSV on ArcGIS Pro?

texts:
I am working with a shapefile in ArcGIS Pro that has 3.2 million observations and about 20 columns. The .shp file is 87MB and the .dbf is just over 2GB.  I am attempting to accomplish an export.

When I export the data via the standard ArcGIS Pro instructions (http://pro.arcgis.com/en/pro-app/help/data/tables/export-tables.htm), the table successfully exports as a .csv but at around row 1.6 million all the observations begin to repeat the same observation. No observations are dropped or added, and the original shapefile remains intact.

Has anyone else experienced this?

Edit: after replicating the issue in QGIS, it seems that it has trouble getting past the 1.6millionth observation, and then copies that same observation down. In QGIS, it simply halts at that 1.6millionth observation (i.e. feature count directly differs between shapefile and the resulting table).


texts: Local Cartesian (NED) to Geodetic (Lat,Lon,Alt) Conversion - Max Distance from Reference Point

texts:
Looking at this:

https://www.mathworks.com/help/map/ref/ned2geodetic.html:

[lat,lon,h] = ned2geodetic(xNorth,yEast,zDown,lat0,lon0,h0,spheroid))

The obvious question is by how much the current location on local Cartesian NED (xNorth,yEast,zDown) can be far away from the reference point (lat0,lon0,h0)?

For example if the reference point is in Australia and the current location is in the US, does those "mathematical tricks" still will be valid?

I assume that there is an intermediate step of-

Local Cartesian NED to geocentric ECEF: https://www.mathworks.com/help/map/ref/ned2ecef.html

Since they also using Geodetic reference point - probably first converting the Geodetic reference point to ECEF, and then there is the questioned "mathematical tricks" to "add" current local Cartesian NED to that reference point (already represented in ECEF) so the current local Cartesian NED will be in ECEF - Thus the "mathematical tricks" are some sort of Cartesian/Vectorial manipulations, I suspect how long they holding in the context of maps, geodesy etc?

Regarding converting Geodetic reference point into ECEF in the intermediate step mentioned above, and in the final step converting back current Cartesian point (already in ECEF) to Geodetic, they use:

Transform geocentric (ECEF) to geodetic coordinates: https://www.mathworks.com/help/map/ref/ecef2geodetic.html

Assuming this is pretty much straight forward formulas with known limitations (Characteristics of different methods to convert ECEF to LLA) So I don't see how they interfere here...


texts: How to remove circles from new polygons in QGIS?

texts: I've just started using QGIS (2.6.0-Brighton).  Every time I create a new polygon, a circle appears in the center (see below).  I thought this might just be a feature of QGIS editing that helps you to distinguish between new polygons, but the circles stay there even in the print composer and exported maps.  Any ideas on how to remove them?


texts: can i merge resolution with float single data type

texts:
Can I merge resolution with both image is float single data type?
I have done this, but the output image become white.
Do the data type must be in unsigned-8 bit for the merge to be success?


texts:
Where can I find a shapefile from a specific city?

I am looking from a shapefile from Valencia (Spain). I only find files from the entire country. I want a specific shapefile with some polygons where I can filter a given point from the city and discover if it is in that polygon.

Can I download it from Google Maps?


texts: Creating unique buffer for each set of points separately and in one procedure using QGIS?

texts:
How could I create a buffer for each set of points (colour group) separately and in one procedure in QGIS 3?

The image below shows a subset of points for buffering.


texts: Leaflet routing control change marker icon

texts:
Im looking for a way to show different markers for the two waypoint that are initialized by leaflet.routing.control.

This is what i tried:

var controlWalk = L.Routing.control({
        lineOptions: {
           styles: [
               {color: 'white', opacity: 0.9, weight: 9},
               {color: '#FC8428', opacity: 1, weight: 3}
           ]
        },
        waypoints: [
            start,
            goal
        ],
        createMarker: function (i, start, n){
            //for (i = 0; waypoint.length; i++){
            return L.marker (start.latLng, {
                    draggable: true,
                    bounceOnAdd: false,
                    bounceOnAddOptions: {
                        duration: 1000,
                        height: 800,
                        function(){
                            (bindPopup(myPopup).openOn(map))
                        }
                    },
                    icon: WalkMarker
            }).bindPopup(myPopup).openPopup();

The result of this are two markers with the same icon.
Any idead how to change the function to solve the problem?

Theres also another behavior that i dont understand: when the icon appear, the popups don't. Only by clicking the marker they appear.


texts: Set transparency on an extent indicator in ArcMap

texts:
When creating an extent indicator in ArcMap 10.3, is it possible to set the transparency of the background colour of the extent indicator frame? I want the extent rectangle to have a colour fill, but not a solid colour as below:

In the frame properties dialogue box, I only seem to have the option to select a solid colour:


texts: Does saving multiple MXDs negatively affect ArcMap performance?

texts:
I'm using ArcMap 10.6. I have a "Master" MXD that has about 30 layers. My project is to make ~30 exhibit maps at a county scale for a legal issue. After I create the ~exhibit maps, the law team will review and ask for revisions to each exhibit.
My workflow is to use the master to create "Exhibit 1", export a PDF and save the MXD with the same name. Rinse and repeat for the 30 maps. This will allow quick revisions. My problem is that I've noticed ArcMap has gotten substantially slower after creating about 10 of these MXDs.
Does saving multiple copies of MXDs negatively affect ArcMap performance?
Is there a solution to streamline the program?


texts: Export shapefile to filegdb with specific projection

texts:
I tried to export a simple shapefile with polygons via "save as" > Format "ESRI fileGDB" and choose "DHDN / Gauss-Kruger zone 3". So there is the error message :
"Error: Failed at creating table for \project (General function failure.)"

With choosing the projection "WGS84" it works proper, but the data are needed in "DHDN / Gauss-Kruger zone 3". Any suggestions about that? Or solutions?

Thanks!


texts: Grouping Points by Location using ArcGIS Network Analyst

texts:
What I am trying to achieve is pictured below, but programmatically. I have a table of 400 locations from which I eyed out an equal amount of points and divided them into 3 routes with Zip Code Tabulation Areas as my zones.

Are there any toolsets that group points by location, up until a certain number of records?

In my example I have a ceiling of 140 locations per route.

I am using ArcMap 10.2 Desktop Advanced


texts: ArcGIS tool which shows all cells a point flows into?

texts: Is there a tool that is effectively the inverse of flow accumulation of a pour point - as in, if I have a specific point that is the site of potential pollution, is there a tool which will show all cells that are downstream of it?


texts: Create a fishnet in a 'for' loop with Arcpy

texts:
When running the following script, I get the error saying that I need a value for “origin_coord” and for “y_axis_coord”.
So is it not enough to define a template dataset? But how can I fill the “origin_coord” and “y_axis_coord” parameters with Arcpy working with ArcMap 10.0? Because in the loop are a lot of shapefiles from different cities, so I need a loop.

import arcpy
import os
from arcpy import env

env.overwriteOutput = True
env.workspace = r"D:\Users\julia\out_09_09\urbanA"
env.qualifiedFieldNames = False
#make a list with input cities as shapefiles.

fcList = arcpy.ListFeatureClasses()

for shpFile in fcList:
    shpFileName= os.path.splitext (shpFile) [0]
    print shpFileName # works
    # Process: Make Feature Layer
    arcpy.MakeFeatureLayer_management(env.workspace + "\\" + shpFile, "shpFile_ly")
    # join the city with the reclasstable
    jointable = r"D:\Users\julia\out_09_09\reclass.dbf"
    arcpy.AddJoin_management("shpFile_ly", "CODE", jointable, "joinCODE", "KEEP_ALL")
    # Process: Copy Features
    outFeatureClass = shpFileName + "_join.shp"
    arcpy.CopyFeatures_management("shpFile_ly", outFeatureClass)
    # Process: Polygon to Raster
    outraster = shpFileName + "_join.img"
    cellSize = 200
    arcpy.PolygonToRaster_conversion(outFeatureClass, "CODE", outraster, "CELL_CENTER", "NONE", cellSize)
    # Process: Make Raster Layer
    arcpy.MakeRasterLayer_management(outraster, "outraster_ly")#works
    # Process: Add Join
    arcpy.AddJoin_management("outraster_ly", "Value", jointable, "code", "KEEP_ALL")
    # Process: Copy Raster
    outraster2 = shpFileName + "_join2.img"
    arcpy.CopyRaster_management("outraster_ly",outraster2) #works
    # Process: Raster to Polygon
    outPolygons = shpFileName + "reclass.shp"
    arcpy.RasterToPolygon_conversion(outraster2, outPolygons, "NO_SIMPLIFY", "Value")#works

    # Process: Create Fishnet
    cellSizeWidth = "200"
    cellSizeHeight = "200"
    outFishnet = shpFileName + "_net.shp"
    arcpy.CreateFishnet_management(outFishnet, "", "", cellSizeWidth, cellSizeHeight, "0", "0", "", "NO_LABELS", outPolygons, "POLYGON")

Is there someone who can help me? I’m trying it without the Iteration but it is not working even if I add a layer before? Why?

import arcpy
from arcpy import env
env.overwriteOutput = True
env.workspace = r"D:\Users\julia\erste_aufg"

# Process: Make Feature Layer
#arcpy.MakeFeatureLayer_management(r"D:\Users\julia\erste_aufg\de013l_hannover    \de013l_hannover\de013l_hannover.shp", "hannover_ly")

#Process: Create Fishnet
outFeatureClass = r"D:\Users\ju\ers\result\hannover_fisch.shp"
cellSizeWidth = '200'
cellSizeHeight = '200'
templateExtent = r"D:\Users\ju\ers\de013l_hannover\de013l_hannover\de013l_hannover.shp"

arcpy.CreateFishnet_management(outFeatureClass, "", "", cellSizeWidth, cellSizeHeight, '0', '0', "", "NO_LABELS", templateExtent, "POLYGON")


texts: Is it possible to update layer source in QGIS

texts:
I've loaded a few shapefiles as layers in QGIS, and I'd like to make a duplicate of the project and all associated shapefiles into a single directory. However, all of the layer sources are still pointing to the original shapefiles.

Is there an easy way to update the paths to the duplicated shapefiles?


texts: Specifying XY extent when rasterizing using ArcGIS Desktop?

texts:
I'm making some maps for landscape genetics in ArcGIS Desktop 10.3 - raster maps of various variables (habitat type, river network, road network etc). I'm having trouble rasterizing a polyline road network to the specific size I want. The vector rasterizes fine but due to the extent of the road network being smaller than other data types (habitat, river network etc) that I have the resulting raster file has different min-max XY coords. I need all of the maps to be exactly the same size (and clipping to the road network raster would mean loosing data in other layers). I'm currently using the polyline to raster tool in the arc toolbox>conversion tools>polyline to raster.

Does anyone know if there is another tool in ArcGIS Desktop that I could use for this?


texts: How tool to split polylines into two-part lines at vertices in ArcGIS for Desktop 10.1?

texts:
I'm using ArcGIS 10.1 SP1 at and we have a problem. There is a serious topology bug that I've mentioned in another question that is affecting our workflow. At present, Esri has stated the fix will be in 10.2, but my government is not able to upgrade to 10.2 so easily. Polylines flagged by topology as overlapping or not being covered by something are normally easily subtracted. With the bug, the subtraction affect the entire polyline, beyond the topologically incorrect area.

My workaround is to use the Split Into COGO Lines tool, which is on the COGO toolbar, to split the offending polyline at all nodes, including where the subtraction should end (checking for topology errors generates nodes). Then, when I do the subtract, it only wipes out the offending lines. However, it also adds COGO attributes to all the other segments along the former polyline. This attribute update behavior is undesirable.

Is there a simple tool that would split a line at all nodes or vertices indiscriminately--and either leave COGO attributes in place, or wipe them?


texts:
Are there any recent comparisons/studies of ArcGIS Server vs. open source GIS servers in common tasks like providing a WMS/WFS, filtering data/search performance, and so on?

We received a lot of critic about featuring this thesis on our blog as the thesis is outdated and a little bit picky on serving WFS only. Citation from reddit: "This is like comparing two cars based on cup holder design."

There is an older presentation given here but like written: It's from 2010


texts:
I am trying to get from a particular point more informations about the surroundings (which are also points). So I wrote this:

SELECT *
   FROM osm_point
   WHERE ST_DWithin(way,(ST_SetSRID(ST_MakePoint(12.2385578688105, 51.8435615931295), 4326, 1000);

But there is an error

  ERROR: FEHLER:  LINE 3: ...T_MakePoint(12.2385578688105,
  51.8435615931295), 4326, 1000) ^ SQL state: 42601 Character: 133

What is wrong?


texts:
I am trying to create with PROJ a custom CRS in QGIS 3.2.0 - but I fail.

I have a local, cartesian coordinate system and a transformation report (2D-Helmert) with the following parameters: translation of the x-axis, translation of the y-axis, scale factor in ppm and rotation angles. With this Helmert-transformation the local coordinates should be transformed into cartesian WGS84 coordinates.

So I am trying something like this:
proj=helmert
    x= 2956.5331  y= 5917.9548  z= -28.7120
    rx= -0° 00' 09.19941" ry= -0° 00' 02.81345" rz= -1° 39' 10.90398"  s= 1 + 0.1347 ppm

I get an error message - am I doing something completely wrong?

Glad about any answers...

Best
Susanne


texts: Calculating quartiles of raster data within polygons using v.rast.stats

texts:
I'm trying to calculate the first and third quartiles of raster data within multiple polygons. The data spans the US and Canada, and there are over 2000 polygons. I've been using QGIS 3.14.0, but I also have access to ArcGIS if necessary.
My initial go to is using the QGIS GRASS tool v.rast.stats since it has the quartile option (unlike Zonal Statistics), but the program stops responding after I press run. I don't know if my inputs are wrong or if the tool is broken. The version of GRASS I'm using is 7.8.3. I've done a test run with much smaller data, but the same thing happened.
If the tool is broken, then what other tools or methods could help with this calculation? I'm not that familiar with coding, but I'm willing to try use R or Python.
Thanks in advance.

Here are my inputs. The two statistics I used are First Quartile and Third Quartile.

This is the program log. I can't copy the text because it freezes when I click on it, but the red text says: "WARNING: Concurrent mapset locking is not supported on Windows"


texts: geoserver returns incorrect number of points when cql used

texts:
I'm using Geoserver 2.2.4 to support WMS queries. All testing is done via the "Layer Preview" page.

If I send a request with no CQL a map full of points is returned as expected.

If I add a CQL filter the map that is returned has the wrong number of points. That is, I get about 50 points on the map when the underlying SQL query (pulled from the log and executed separately) returns over 1200 points.

In both cases the same style is used. The only difference is the use of a CQL filter.

How can I get all 1200 points displayed on the map?


texts: Adding drag interaction breaks ability to move map location in OpenLayers 3?

texts:
I am working with a version of OpenLayers 3.  I am attempting to add the ability to drag a feature on the map.  The problem I am running into is when adding a drag interaction to the map, I break the base functionality of moving around the map location.

I created a jsfiddle to demonstrate.  When opening, try changing location of the map.

http://jsfiddle.net/mcroteau/w8bj1htb/

Javascript

var mapTile = new ol.layer.Tile({
    source: new ol.source.OSM()
});

var lineStyle = new ol.style.Style({
    stroke: new ol.style.Stroke({
        color: 'blue',
        width: 3
    }),
    fill: new ol.style.Fill({
        color: 'rgba(0, 0, 255, 0.1)'
    })
})

var verticeStyle = new ol.style.Style({
    image: new ol.style.Circle({
        radius: 6,
        stroke: new ol.style.Stroke({
            color: 'rgba(0, 0, 255, 0.8)',
            width: 2
        }),
        fill: new ol.style.Fill({
            color: 'rgba(0, 0, 255, 0.5)'
        })
    }),
    geometry: function(feature) {
        var coordinates = feature.getGeometry().getCoordinates()[0];
        return new ol.geom.MultiPoint(coordinates);
    }
})

var mapView =  new ol.View({
    center: [-11000000, 4600000],
    zoom: 3
});

var map = new ol.Map({
    layers: [mapTile],
    target: 'map',
    view: mapView
});

/**
* TODO : Issue 1 Drag Interaction
*
* adding drag interaction to map breaks ability to move map location
*
**/

var dragInteraction = new ol.interaction.Pointer({
    handleDownEvent : function(event){
        var feature = map.forEachFeatureAtPixel(event.pixel,
            function(feature, layer) {
              return feature;
            });

        if (feature) {
          this.coordinate_ = event.coordinate;
          this.feature_ = feature;
        }

        return !!feature;
    },
    handleDragEvent : function(event){

        var feature = map.forEachFeatureAtPixel(event.pixel,
            function(feature, layer) {
              return feature;
            });

        var deltaX = event.coordinate[0] - this.coordinate_[0];
        var deltaY = event.coordinate[1] - this.coordinate_[1];

        if(feature){
            var geometry = feature.getGeometry();
            geometry.translate(deltaX, deltaY);

            this.coordinate_[0] = event.coordinate[0];
            this.coordinate_[1] = event.coordinate[1];
        }

    },
    handleMoveEvent : function(){
        //console.log('handle move')
    },
    handleUpEvent : function(event){
        //console.log('handle up event')
    }
})
map.addInteraction(dragInteraction);

var features = new ol.Collection();
var vector_source = new ol.source.Vector({
    features: features,
});

var featureOverlay = new ol.layer.Vector({
    source: vector_source,
    style: [ lineStyle, verticeStyle ]
})
featureOverlay.setMap(map);

var geometryFunction = function(coordinates, geometry) {

    var start = coordinates[0];
    var end = coordinates[1];

    if (!geometry) {
        geometry = new ol.geom.Polygon(null);
    }

    geometry.setCoordinates([
        [start, [start[0], end[1]], end, [end[0], start[1]], start]
    ]);

    return geometry;
}

var draw = new ol.interaction.Draw({
    condition: function(event) {
        return ol.events.condition.shiftKeyOnly(event)
    },
    features: features,
    type: 'LineString',
    geometryFunction: geometryFunction,
    maxPoints: 2
});

map.addInteraction(draw);


texts: Loading Topology to Parcel Fabric gives Must Contain Features Error?

texts:
I've been working on migrating our parcels to the parcel fabric with hopes of migrating to the LGIM.  I imported several CAD files, cleaned a large test area, and validated the topology (0 errors).  Everything gets loaded into the staging geodatabase that is provided by ESRI, and the only two feature datasets I'm working in are Type6LotUnits and Type7TaxParcels.  Type6 and type7 data loads just fine within the staging area.  Type6 data will load into fabric with no issues and display my lots as it should, but whenever I try to load type7 (the actual parcels) I get an error saying that the source feature class must contain features.

\USER\xxxx\Final Parcel

Fabric\Staging_Geodatabase.gdb\Type7_TaxParcels\Type7Parcels" # # "0.1 Meters"

UNJOINED_GROUP DEGREES_MINUTES_SECONDS QUADRANT_BEARING COMPUTE_AREA ACRES "0.5 Meters" 3_1908_TO_1980

Start Time: Thu Jul 14 16:17:50 2016

Verifying input topology dataset...

Validating input dataset coordinate system...

Importing topology.

Loading topology lines...

The source feature class must contain features.
Failed to execute (LoadTopologyToParcelFabric).

Failed at Thu Jul 14 16:17:57 2016 (Elapsed Time: 6.85 seconds)

I'm literally stumped at this point.  I've tried to load the data on different machines, locally vs. over the network.  I've loaded and reloaded everything, including moving the data to new staging databases/fabrics.  Nothing I do seems to work.  Has anyone encountered this issue?  If so, how was it fixed?  As you can see, no error code is mentioned, and the internet seems to not yield anything for an answer.


texts: I want to make the cells which values "1" to null in QGIS with Grass or Saga tools (exactly Saga).


texts: DB Manager - Oracle DB - Importing vector layer fails prompting "ORA-01741: illegal zero-length identifier" error

texts:
When trying to import a vector layer, to an existing Oracle DB connection, using DB Manager core plugin, the operation fails. An error message prompts with "ORA-01741: illegal zero-length identifier":

This happens due to a missing geometry field name parameter in the SQL statement as can be seen in the two adjacent double quotes before the MDSYS part:
CREATE TABLE "GIS_WAT_RH"."test"("" MDSYS.SDO_GEOMETRY)

This is puzzling because in the Import vector layer wizard, the geometry column is populated with a geometry field name though grayed out:

Why does the geometry field isn't being populated in the DB Manager SQL statement?
How to Reproduce

Make an Oracle DB connection via Browser or DB manager
Open DB Manager.
Navigate to the desired schema inside the previously created Oracle
connection
Import Layer/File...
The Import vector layer wizard appears. Input data from a valid
datasource.
Press OK


texts:
I need a little help/ ideas on how to convert tab files.

The problem: I have many TAB files produced daily by 3 party app, I need to batch import the files into postgres.

Most files import fine, however some files has mixed content like linestrings and multilinestrings in one file.

I can't import that type of files.

I can open the files in qQGIS, and save as a shapefile and then import to the database, but it needs to be done in a batch job.

This is what I'm using and what is working for files that only has linestrings.

C:\Program Files (x86)\FWTools2.4.7\bin>ogr2ogr -skipfailures -overwrite -t_srs "EPSG:4326" -s_srs "EPSG:4326" -a_srs "EPSG:4326" -f "PostgreSQL" PG:"host=localhost user=zzz dbname=geo password=zzz" C:\gisdata\Cables_1.tab  -nln tmpimport.Cables_1

This is the error I get when importing a mixed content file

Warning 1: Geometry to be inserted is of type Multi Line String, whereas the layer geometry type is Line String.
Insertion is likely to fail
ERROR 1: INSERT command for new feature failed.
ERROR:  Geometry type (MultiLineString) does not match column type (LineString)

Anyone has any advice?


texts:
I have a road network of center lines of MultiLinestrings. They have nodes at every intersection. I try to understand how I can merge roads together that have the same category as  attribute (there are about 80 different categories). Roads with the same category are not always connected, so I guess I need to have the result back as a MultiLineString again but grouped on attribute value.

Is there a way to do this?

Edit:
I played around some more and almost got the result I wanted.
Simply using

SELECT  "category", ((ST_Collect(geom))) as geom FROM s.road GROUP BY "category"

But this seem to be returned as a Geometrycollection of Multilinestrings and cannot be loaded into QGIS.
Any advice?

Edit 28.07

Final solution:

SELECT  "category", ((ST_Union(geom))) as geom FROM s.road GROUP BY "category"

Using ST_Union instead of ST_Collect gave me the results I wanted. For beginners in SQL like myself can I explain why. ST_Union dissolve overlapping and overlapping regions whereas Collect does not. The consequences using ST_collect in this case with MultiLineString is that the result is a Geometry collection that is not a supported type for QGIS to load.


texts:
(This is my first attempt to use R, so I apologize for any stupidity. Also, please let me know if this is better posted in Stack Overflow.)

I'm trying to use the Microsoft building footprint dataset to make a map of buildings in my city (Cambridge, MA). I'm following this gist as a guide. These lines of code are giving me problems:

options(tigris_class = "sf", tigris_use_cache = TRUE)
zctas_national <- zctas(cb = TRUE)

selected_zip <- zctas_national %>%
filter(ZCTA5CE10 %in% c(02139))

When I try to print(selected_zip), I get nothing back. Or rather, I get a table that looks like this:

[1] ZCTA5CE10  AFFGEOID10 GEOID10    ALAND10    AWATER10   geometry
<0 rows> (or 0-length row.names)

If I replace 02139 with 85308 (the zip from the gist), everything is fine:

  ZCTA5CE10     AFFGEOID10 GEOID10  ALAND10 AWATER10
1     85308 8600000US85308   85308 44254346   872904

Anyone know what's going on?


texts: Seeking appropriate CRS for Philippines while learning QGIS?

texts:
What is the appropriate CRS for Philippines?

I'm trying to figure out what is the best CRS for Philippines while learning QGIS.


texts: Attaching ArcGIS Online attachments in bulk?

texts:
I have an ArcGIS for Desktop Advanced license.  I have an attachment enabled feature service in ArcGIS Online.  I have thousands of files I wish to attach to it.  By current process convention, these files are named according to 2 fields in the attribute table of the feature service (for current and future scripting needs).

What I CAN do:

I can attach documents one by one to a local copy of the feature service using attachment manager and successfully sync those attachments.
I can create a script that creates a match table for the AddAttachment tool with a field for the Object ID in the feature service and a field storing the file path.
I can use this match table to attach all of these files to the local copy of the feature service at once.

What I CANNOT do:

Successfully sync this local copy and all of its new attachments to ArcGIS Online.

The Add Attachment tool and its associated arcpy function require the file geodatabase to be outside of an edit session to run.  My first guess as to why, after running this tool successfully on a local copy of an AGOL layer, the results can't be sync'ed has to do with it being run outside of an edit session.  I'm guessing (but don't actually know) that Desktop -> AGOL sync process sends over the instruction set performed during edit sessions (the attachment manager process for adding attachments requires an edit session).  The actions performed outside of an edit session are simply not sent through the sync process, meaning any tools that can't be run in an edit session, like Add Attachments, can't be sync'ed to AGOL.

Is there a known way to attach to an ArcGIS Online feature service in bulk, whether through arcpy, REST service, or some other tool I don't know of?


texts: Programmatically extract edge points from many polygons?

texts:
Objective: Extract the four Latitude, Longitude points which correspond to the corner edges of each polygon in the Military Grid Reference System (MGRS) shapefile.

How can this be done in QGIS?

In this picture the red arrows point to the places I would like to know the Latitude, Longitude:

Here is what the MGRS shapefiles look like in QGIS 2.8.16, notice the corners look good:

When I used:
y_min(bounds($geometry))
x_min(bounds($geometry))
y_max(bounds($geometry))
x_max(bounds($geometry))

And then plotted the data in Tableau Desktop, the result is some of the polygons are overlapping slightly and the corners aren't matching up exactly.

Here are the points which make up 14SPC:

GZD,LATITUDE,LONGITUDE,MGRS,PointID
14S,32.519328,-97.93512,14SPB,1
14S,33.434729,-99,14SNC,1
14S,33.420739,-97.924236,14SPC,1
14S,33.403103,-96.848967,14SQC,1
14S,34.322005,-97.912856,14SPD,1
14S,33.434729,-97.93512,14SPB,2
14S,34.341303,-99,14SNC,2
14S,34.336476,-97.924236,14SPC,2
14S,34.322005,-96.848967,14SQC,2
14S,35.238088,-97.912856,14SPD,2
14S,33.434729,-96.848967,14SPB,3
14S,34.341303,-97.912856,14SNC,3
14S,34.336476,-96.826229,14SPC,3
14S,34.322005,-96,14SQC,3
14S,35.238088,-96.802458,14SPD,3
14S,32.519328,-96.848967,14SPB,4
14S,33.434729,-97.912856,14SNC,4
14S,33.420739,-96.826229,14SPC,4
14S,33.403103,-96,14SQC,4
14S,34.322005,-96.802458,14SPD,4


texts:
I have the following set of points (there are far more). Originally they were given to me without any decimals. I was able to correctly plot them into arcmap however they don't line up with my other shapes, they float off to no where. I've tried re projecting and define projection, but they still go no where. I can only conclude they are incorrect co-ordinates, as when I try to google maps a point (-33.726298, 127.575079) google maps can't find it.

Points in the area I am looking at should look like -37.837266, 145.040987. I cannot understand how I was able to create shapes from these points that look correct, but don't line up with my layers, other than the distance between each point is still correct. They are meant to be in GDA94, but I still can't locate them in the Melbourne area where they should be.

 x               y
3372629.81  1275750.79
3370659.42  1276207.50
3370415.87  1275014.01
3368757.31  1275354.90
3368493.70  1274241.98


texts: Scaling surface reflectance product in ENVI?

texts:
I am using ENVI 5.4, specifically 'Data Stretch' tool in a bid to scale down my Landsat 8 OLI surface reflectance product values. I have populated the dialogue box (as below) and compared the histogram of the resultant raster with the original one.

However, there is anomaly noticed in its values (see histogram). I anticipate same histogram stretch however this is not the case.

Why there is such discrepancy in the results?


texts: Given a line/shape how to calculate the area that resemble how the US Census Geography's Inland Water is calculated

texts:
I am not a US resident so not quite comfortable asking the US Census directly about this, but I'm just curious how the US Census's Inland water geography is derived, i.e. the methodology of what geoprocessing operation are used with what parameter.
I want to kind of apply the same definition of Inland Water used by the US Census (see GARM Ch15) onto other line/shape I have, but I can't figure out how to do it; particularly it's how to derive that (or those?) 1 nautical mile line between two points on a given continuous line (I guess basically have the program figure out where should be an inlet and ensure that inlet has the maximum area possible, considering for the entire length of the inlet there could be multiple location where it's narrow to 1 nml and you'd want the one furthest out).
either ArcMap or QGIS references are fine, though QGIS is somewhat more preferred as that's free-to-use; also tbh I'm not that skilled with both and my GIS skill are not that great nowadays--been away from this field for while, so this is just something I wondered for some time.
New Edit: after some thought the best I can come up so far is to guess that for starter the coastline should probably be a polygon object (fake the other line segments?) since there might be a problem with line that whatever solution might end up working the wrong side of the line so maybe a polygon would preempt that...


texts: Open layers supercluster integration - is it possible to fetch underlying features for each cluster?

texts:
I am using the open layers 3 with supercluster at the moment, extending the example provided here : https://jsfiddle.net/d5gwvz2y/2/

My question is: How does one fetch all the underlying features behind a cluster if clicked? In open layers cluster implementation you could use cluster forEachFeature, or just getFeatures.

Here is a sample of the clusters working using default cluster implementation http://www.wwarn.org/surveyor/NMFI/ when a cluster is clicked, an array of features are returned, which can then be iterated through.


texts: Mouse click event used gis application with much data

texts:
I will create an application that using my water pipes, electric cables, gas pipes service in a city. All of streets and roads including water pipes. My service is publishing pipes. My service like following.

Pipe Service

domain.com/arcgis/rest/services/Water/MapServer
domain.com/arcgis/rest/services/Electric/MapServer
domain.com/arcgis/rest/services/Gas/MapServer

Feature Layers

Pipe Lines (Polyline): domain.com/arcgis/rest/services/Water/MapServer/0
Water Valves (Point): domain.com/arcgis/rest/services/Water/MapServer/1

I wanna show pipes on map and users can click on pipes and get features about pipe. Which way can be make sense for this solution?
I think use FeatureLayers on map and on click method. How does it affect the performance. If I use ArcgisDynamicMapServiceLayer, I can not click pipes.


texts: Calculating distance from points to different class in raster using ArcGIS Spatial Analyst?

texts:
I am looking for some advice on how to calculate the distance from a point to a different class of a raster.

Specifically, I have observations (points) of animal use overlaid on land cover data.

For each point I want to calculate the distance to the nearest cell of a different cover type.

For instance, if a point is located in a grassland patch, what is the distance to the nearest cell that is not grass.

I am using ArcGIS 10 and have a Spatial Analyst license.

I think raster calculator might be the best approach for this but could anyone provide an equation that would accomplish this?


texts: Comparing raster layers and finding the combined minimum

texts:
I have 5 Raster files each depicting a soil property, all of the same area. I want to find the location with the "worst" combined soil properties. It's also important to note that the values of the layers don't have the same values. I tried r.series but I'm not sure if it is the right way to do it. From my understanding r.series compares all rasters at once and picks the min (if selected) for each location?

I also heard that I can use the raster calculator to get what I want, but I have no idea how to form an equation that does what I want.

Another thing I thought about was to convert the layers to vector data and try working with them, but since the layers are really large I don't want to wait 5 hours only to find it wouldn't work.

I know this may not be a straight forward question, but I've tried a lot of things and nothing works and at this point I'm kind of lost.


texts:
I am testing the trial version of FME 2014 in order to make sure I can build some workflows on it. If its possible then we plan to proceed with buying a license.

Right now I have build a workflow which partially works. I have some questions on how to accomplish something and I was wondering if you can help me with.

Right now my workflow is separated into two workspaces. The first one takes an Excel (applies some attribute filter) and imports everything in a database in postgresql (with sql executor). This workflow doesn't include any spatial data.

The second one takes a dwg file, applies some transformations and it outputs a shapefile with point geometry plus some other attributes (id etc.).

Now, the table from the first procedure and the shapefile from the second procedure have in common this id.

What I want to do is to complete this workflow by merging them into a big table in PostgreSQL (with PostGIS extension). The crucial part here is that I want this table to get updated later with more values (inserting more values).

What is the procedure I have to use?

Right now I have tried to use SQL Executor in order to pass the data from the ShapeFile into the same database but it doesn't work cause there is geometry user attribute available to choose. Why does this happen?


texts:
Every year I spent lot of time manually adjusting annotation labels for our 2018 stations (250+ points & labels in a A4 page). To save effort, I want to reuse most of these labels for 2019.

Can I join 2019 station point features and copy/save the 2018 annotation labels as 2019 labels and add/drop labels; few stations get removed or few new ones added?


texts: .sjoin Geopandas returns too many results

texts:
I have two dataframes, both containing geometry columns. The first dataframe contains MULTIPOLYGONs while the second one contains POINTs. My aim is to join the dataframes so each POINT will be assigned to corresponding POLYGON. DF_polygons has 2480 rows and DF_points has aroud 150.000 rows (10 tractors with GPS driving over fields, data coming in every few seconds).

I used this code to join dataframes (also tried within and contains):

points_in_polygons = gpd.sjoin(DF_points, DF_polygons, op='intersects', how='inner)

The problem is, it returns a result with approximately 30.000.000 rows in all 3 cases. The desirable result would have 100.000 rows - so each unique point would be assigned to one of 2480 multipolygons.

Is it possible that the problem is in overlapping polygons? How can I deal with it? Which packages should I use? I am completely new to GeoPandas.


texts:
I am trying to follow a simple tutorial on Dotspatial web page "Sample code that demonstrates how to create a Multi Polygon Feature Set" I add the DotSpatial.Data.dll to the project references but I can not find the "DotSpatial.Geometries.dll"!

Can you please let me know where can I get it?


texts: What is Google Map’s knowledge of historic street names?

texts:
I am geocoding historic addresses in Ireland using Google Maps (free access) but in some cases Google Maps returns an address with a different street name. For example, when I enter the address of the Forthriver Primary School that was recorded on Ballygomartin Road, Northern Belfast in 1990 it is geocoded to a school ('point_of_interest', 'school') on Cairnmartin Road.

Is Google Maps utilizing knowledge of historic street names but providing the updated street name? How do I know it is locating the correct address, especially in cases where there is no point of interest to identify (just an old street name)?


texts: ListFeatureClasses isn't returning anything in GDB

texts:
I am writing a script that performs a series of operations on line, point, and polygon files in my file geodatabase matched through a unique number in their filenames. My geodatabase looks like this:

sd_junc_1.shp (point file)
sd_junc_2.shp (point file)
sd_test_1.shp (polygon file)
sd_test_2.shp (polygon file)
sd_test_streets_1.shp (line file)
sd_test_streets_2.shp (line file)

I'd like to be able to find each file that ends with the number 1, perform the operations, then do the same operations on files that end with 2, and so on. Here's the code I have so far:

#import os
import os

# Local environment
arcpy.env.workspace = r"E:\LA 221\Subdivisions\Wake County\wakecounty.gdb"

# Initialize counter
counter = 1

# Loop

while counter < 6:
    counter_string = str(counter)
    # Create variables from feature classes in the GDB
    Clipped_streets = [fc for fc in arcpy.ListFeatureClasses(feature_type = 'Line') if os.path.splitext(fc)[0].endswith(counter_string)]
    Clipped_streets_Junc = [fc for fc in arcpy.ListFeatureClasses(feature_type = 'Point') if os.path.splitext(fc)[0].endswith(counter_string)]
    Subdivision_boundary = [fc for fc in arcpy.ListFeatureClasses(feature_type = 'Polygon') if os.path.splitext(fc)[0].endswith(counter_string)]

    # Process: Select access points
    selected_access = arcpy.SelectLayerByLocation_management(Clipped_streets_Junc, "BOUNDARY_TOUCHES", Subdivision_boundary, "", "NEW_SELECTION", "NOT_INVERT")

This code fails at the last line (selected_access = ...) with the error "RuntimeError: Object: Error in executing tool" and I'm guessing because there is nothing in the "Clipped_streets_Junc" and "Subdivision_boundary" variables.


texts:
I am working on trying to create a tool which would allow me to type in any city name (in this cause European cities), and zoom into that city. My script will run but nothing is selected. I don't know  what I'm missing.

import arcpy
from arcpy import mapping
#set workspace
mxd = arcpy.mapping.MapDocument("CURRENT")
#set the dataframe
fc = ("N:/data2016/Lab13/cities.shp")
cities = arcpy.mapping.ListLayers(mxd, "Cities")[0]
#Set local variables
cityName = arcpy.GetParameterAsText(0)

cities_layer = arcpy.MakeFeatureLayer_management(fc, "cities_lyr")
arcpy.SelectLayerByAttribute_management("cities_layer", "NEW_SELECTION", CITY_NAME = cityName)

arcpy.mapping.ListDataFrames(mxd)[0].zoomToSelectedFeatures()
arcpy.RefreshActiveView()


texts: How can I join two data bases keeping all records for both DB?

texts:
I'd like to know if there's one way to join two data bases but keeping all the records for both DB? I have a point layer of 300000 locations and I want to assign them the attributes from a .csv table that contains info about people with health problems.

The .cvs file contains several individuals for the same location and once I've made the join with the point layer, only the first record is joined.


texts: Python conditional statement for range of numbers isn't working

texts:
Using ArcGIS desktop 10.7. The field Shape_area_Hect is float and the field Core_patch is string. Error for following syntax:
def Reclass(Shape_area_Hect):
  if (Shape_area_Hect >= 12150):
    return "Core"
  elif (Shape_area_Hect >= 486 and Shape_area_Hect < 12150):
    return "Patch"
  elif (Shape_area_Hect < 486):
    return “Less than patch”

Core_patch = Reclass(!Shape_area_Hect!)


texts:
I'm very new to GIS and have recently installed QGIS on my mac desktop computer, running Yosemite 10.10.

I'm working with DEM data for stream network extraction, but the majority of grass raster commands are not functional. For example, when attempting a stream network extraction, I get this output:

It seems to me that the first major problem is that QGIS cannot find the command r.stream.extract. However, it CAN find some other (simpler) commands, such as r.aspect. I am completely unfamiliar with digging around in the guts of these types of problems and would really like to know how to remedy this problem.

Is it as simple as somehow completing my collection of r.* commands, and if so, how?

Or is there something more subtle going on here?


texts:
I'm using pyrtree to create an index for a point in polygon operation. I am currently trying to insert 1.8 million rows into the index and it takes more than an hour to do so. I am unable to use rtree in pyspark and hence chose to use pyrtree. I read about bulk inserting into rtree but I did not find anything similar in pyrtree. Below is the script I am using for the insert (single insert)

from shapely.geometry import Point, shape
from pyproj import Proj, transform
from pyrtree import RTree,Rect
places_idx_new=RTree()
for index,row in places_data.iterrows():
    pt=Point(transform(Proj(init='EPSG:4326'),Proj(init='EPSG:3857'),row['longitude'],row['latitude']))
    pt_buff=pt.buffer(places_data['category_buffer'][index])
    places_idx_new.insert(index,Rect(*pt_buff.bounds))

Please let me know if there is a way to bulk load into pyrtree or if there is a faster way to do this and If I can provide any further information on this


texts: Empty error when running Table to Excel tool

texts:
The Table to Excel tool in the Conversion toolbox does not work for some reason. When I open it, it shows a red cross in the top left corner. When I try to run it an error pops up that says nothing. I've tried different input tables, different output locations and tried to do it it in python, which gives the "serious application error" and quits ArcMap. The result does not even show up in the geoprocessing window.

I have ArcMap 10.5.1.


texts: Updating URL in GeoJSONLayer using ArcGIS JS API

texts:
I am working on a simple node/react application that allows for the the request of data using the URL parameter in GeoJSONLayer. The way that I set up the API I need to be able to change the url to get a new set of data. I have set up a useEffect() to load the map that has all of the map elements and have told it to listen for changes on the treeURL.
  const [treeURL, setTreeURL]  = TreeURL();
  useEffect(async()=>{
    loadMap(treeURL)

  },[treeURL]);

The TreeURL function is just a simple useState() to default to getAll which is the default way of getting all of the data back. This part is working fine the map renders with the data that is fetched.
const TreeURL = () => {
    const getAll = 'getAll'
    const [treeURL, setTreeUrl] = useState(getAll)
    return [treeURL, setTreeUrl]
}

When I go to update the data using the setTreeURL() the function runs. treeURL prints out into the console, but useEffect() is not called and the map does not update.
props.list.forEach((item) => {
        elements.push(
            <li class="mr-3">
                <button
                    key={item}
                    className="inline-block border border-blue-500 rounded py-1 px-3 bg-blue-500 text-white"
                    onClick={()=>{
                        setTreeURL(`getByParams?CONDITION=${item}`)
                        console.log(treeURL)
                    }
                    }>{item}
                </button>
            </li>)
    })