texts: Joining by 'almost' identical attributes using QGIS?

texts:
So I have two point shapefiles - their attribute table consist of names, lat, lon, and others.

I would like to join the attributes of shapefile 1 to shapefile 2. A normal join between the names would to the job.

The problem is that the two datasets have slight differences in the 'name' field as well as in the lat/lon. For example, if I have one point in The Hague, Shapefile 1 will calle it 'Den Hague' whereas Shapefile 2 would calle it 'T. Hague'. Also spatially they are not identical (they do not overlap) so the location between both of them would differ by few kms in spite of being the same point. I have over 7000 points which make it tough for manual corrections. When some points are indeed correct.

Is there a way to do a matching which takes into account ~80% of similarity in the names?

I am using QGIS.


texts:
I'm trying to follow this tutorial so I have a method for creating contours in QGIS because I find the contours generated by the contour plugin a bit strange (very sharp edges and whatnot).

I can't get past the first stage of creating the raster using IDW interpolation.

I get the following error when I "Run as batch process":

  NameError: name 'QgsDoubleSpinBox' is not defined

I am using QGIS 3.6 Noosa

I'm wondering if anyone can successfully create the raster and how to do it? Can anyone successfully do the whole tutorial in 3.6 or whatever version they're using?

There is downloadable input files with the tutorial and so I downloaded them and I add the layer as a vector layer. I have tried using the .dbf, .shx and .shp but they all produced the same error for me.

I'm new to QGIS so I might be doing something wrong.

Points to Raster to Contour Lines Tutorial


texts: Using 'sp' method to subset SpatialPointsDataFrame with SpatialPolygon returns error in parLapply()

texts:
I have numerous very large data sets of ship tracking data (Automated Information System; AIS). These are originally in a file geodatabase and can be downloaded here (this is a smaller file for testing purposes).  After reading in the file using readOGR, I subset the data by attributes and then into latitudinal bands like this (warning - file is large and takes a while to download and read in. You will also need ReadFGDB driver support in your version of rgdal):

  # Libraries
  library(rgdal)
  library(sp)
  library(parallel)

  download.file("https://coast.noaa.gov/htdata/CMSP/AISDataHandler/2014/01/Zone9_2014_01.zip", "temp_file.zip")
  unzip("temp_file.zip", exdir="temp")
  broadcast <- readOGR("temp/Zone9_2014_01.gdb", layer = "Zone9_2014_01_Broadcast")

  # Remove unneeded fields
  broadcast@data <- broadcast@data[,c("SOG", "BaseDateTime", "Status", "VoyageID", "MMSI")]

  # Subset the data
  broadcast <- broadcast[broadcast@data$SOG>0.5,] # remove all points w speed <0.5
  broadcast <- broadcast[broadcast@data$Status==0 | broadcast@data$Status==7 | broadcast@data$Status==8,] # remove all points at anchor

  # Split data set into 4 regions
  coords <- coordinates(broadcast)[,2]

  chnks <- list()
  chnks[[1]] <- broadcast[which(coords<=37),]
  chnks[[2]] <- broadcast[which(coords>37 & coords<=40),]
  chnks[[3]] <- broadcast[which(coords>40 & coords<=44),]
  chnks[[4]] <- broadcast[which(coords>44),]

Now chnks is a list with four SPtsDFs which I want to clip with a polygon, EEZ.  (EEZ is the US West Coast Exclusive Economic zone).  I'll substitute a simple polygon for EEZ here:

EEZ <- SpatialPolygons(list(Polygons(list(Polygon(cbind(c(-117, -122, -129, -129, -124, -124, -117),c(32, 30, 39, 47, 49, 39, 32)), hole=F)), "EEZ")))

If I use lapply to do spatial subset on chnks, it works fine:

chunked <- lapply(chnks, FUN = function(x){
# Clip the AIS points to EEZ
return(x[EEZ, ])
})

But if I try to use parLapply, it fails with an error:

  # create cluster
  cluster <- makeCluster(4, outfile = "")
  clusterExport(cl=cluster, varlist = "EEZ")

  chunked <- parLapply(cl = cluster, chnks, fun = function(x){
    # Clip the AIS points to EEZ
    return(x[EEZ, ])
  })

  4 nodes produced errors; first error: object of type 'S4' is not
  subsettable

It seems like the sp method for [ (which references points[!is.na(over(points, geometry(polygons))),] according to this post) is not being loaded to the workers.

Does anyone have an idea how to avoid this error or explicitly export the sp method to the parallel clusters?

Addition

I also tried using plyr and llply() with .parallel = T and got the same failure message.  I would also note that the sp package seems to load fine on all of the workers, which I would think should allow the [ subsetting method for spatial objects.


texts: Add an individual layer from ArcGISDynamicMapServiceLayer to the map

texts:
I have a ArcGISDynamicMapServiceLayer containing many layers. I'd need to add only individual layers to the map, like layers 2,5 and 10.

The following code would add all the layers to the map:

url = "http://...MapServer";
mapServiceLayer = new ArcGISDynamicMapServiceLayer(url, {
     className: layer,
     opacity: opacity
});
map.addLayer(mapServiceLayer);

How can I add only i.e. layer 5 to the map? (.../MapServer/5)

There is a similar question here but answers seem not to be helpful for the front end solution.


texts: Split string of coordinates and put it in array


texts: Computing intersection of distance circles around fixed points

texts:
(to simplify calculations, I assume a spherical earth. All distances are assumed measured "as the crow flies")

Suppose we have the lon/lat coordinates for two points (p_1,p_2) less than Nkm apart.

Then, around point p_1 I draw a  circular polygon c_1 of the set of all points at distance exactly Nkm from p_1 and likewise for p_2 and c_2.

Now c_1 and c_2 will cross at two points. Is there a formula to find the  coordinates of these two crossings? If such a thing exists, a pointer to an R/python code doing this would be just as good as a formula.


texts:
I'm looking for a tool that would be able to automatically snap features closer than a given distance. Is there such a tool/plugin for QGIS/GRASS/OpenJump? Would it be possible to achieve that with PostGIS?

I have a river network where some of tributaries don't touch the river they flow in and I need to get them snapped.


texts: How to use Landsat surface reflectance pixel quality band?

texts:
I'm using TIMESAT for time-series analysis and Landsat 7-8 Collection 1 images. Within this process I would like to use Landsat pixel quality band for assigning reliability degrees to the pixels using the pixel QA band.

When using CFMask Band the range is easily defined (0 clear; 1 water;2 cloud shadow; 3 snow; 4 cloud). But when I try to use the Pixel QA band the pixel values were different than the values for interpretation given at https://landsat.usgs.gov/landsat-surface-reflectance-quality-assessment.

I have the opportunity to use QGIS, ArcGIS and ENVI for image processing.

How I can use or process the pixel QA band?


texts: Openlayers vector layer markers disappear near dateline

texts:
See demonstration of the bug here: http://jsfiddle.net/3ECSb/4/

I have seen and solved similar problems before with OpenLayers but the usual "three world" concept isn't fixing it this time so I'm putting this out there in case the slight differences can help someone identify the problem.

I have a bunch of points I am adding to my vector layer.

var point = new OpenLayers.Geometry.Point(result['longitude'], result['latitude']);
point.transform(map.displayProjection, map.getProjectionObject());

var feature = new OpenLayers.Feature.Vector(
        point,
        {lakeData:result}
    );
vectorLayer.addFeatures(feature);

The map is centered on New Zealand and the points appear. I zoom in one step and they are gone. Then I pan the map so NZ is almost off the map on the right (the important part is the antimeridian is gone) and the points reappear. If I click and drag the map around the points stay on the screen. However if I release the mouse when the antimerdian is still to the right of center the points disappear. I can then center the map on the Chatham Islands (-176 Longitude) and the points reappear and will stay on screen even if I go back to the previous position when they were not visible.

So to clarify. If I am centered at 179 latitude the points are only visible if I was previously centered at a negative latitude since the last time the antimeridian was brought on to screen from the right.

I have sphericalMercator and wrapDateLine set on every layer, I'm using Google layers as my basemaps, I have tried adding copies of the points at +/- 360 from their original places. But nothing has worked.


texts: Positioning labels in center of clipped polygons?

texts: I am using ArcMap 10.1  and would like to label a polygon layer that i have clipped.  It appears that the label is set for the geometric center of the polygon that falls outside of the clipping.


texts: Is it possible to use NamedStyle inside of UserStyle with MAPSERVER SLD?

texts:
I want to use a NamedStyle inside a UserStyle for the ElseFilter of my SLD-Definition:

<NamedLayer>
<Name>MyLayer</Name>
<UserStyle>
<FeatureTypeStyle>

<Rule>
<Name>Selected</Name>
<ogc:Filter>
<ogc:PropertyIsEqualTo><ogc:PropertyName>my_id</ogc:PropertyName><ogc:Literal>31117</ogc:Literal></ogc:PropertyIsEqualTo>
</ogc:Filter>
<PolygonSymbolizer>
<Fill>
<CssParameter name="fill">#008000</CssParameter>
</Fill>
<Stroke>
<CssParameter name="stroke">#BD0026</CssParameter>
<CssParameter name="width">30</CssParameter>
</Stroke>
</PolygonSymbolizer>
</Rule>

<Rule>
<Name>Others</Name>
<ElseFilter/>
<PolygonSymbolizer>
<Stroke>
<CssParameter name="stroke">#000000</CssParameter>
</Stroke>
</PolygonSymbolizer>
</Rule>

</FeatureTypeStyle>
</UserStyle>

</NamedLayer>

Something like this:

<Rule>
<Name>Others</Name>
<ElseFilter/>
<NamedStyle>
    <Name>the_style</Name>
</NamedStyle>
</Rule>


texts:
I have some data on frequency of different road segments, and need to visualise it as a heatmap. Is there some utility/plugin for Mapbox or Leaflet which does this kind of visualisation. I am only able to find Leaflet plugins for visualising points as heatmap.

My current approach is to convert the frequency into a colour range, and then put it in GeoJSON. But it suffers from the problem where I can't assign proper weights to the polyline. I need the line to exactly overlap the road (proper width and following the road curve), so that it seems like highlighting the road.

Another approach could be breaking the line into multiple points, encoding frequency in them and making a heatmap out of that. But that would be too artificial, and I think that there should be a better elegant way out there.


texts: QuickOSM using wild cards on lines (rivers)

texts:
I'm trying to find all rivers with a name containing specific word (or part of it). F. ex. all 'Danube-' rivers, or 'Danu-', 'Dan-' and so. In another topic I've found useful query, which after example modification looks like this:
<osm-script output="xml" timeout="250">
<union>
    <query type="nwr">
        <has-kv k="river"/>
        <has-kv k="name" regv="Danube(.*)"/>
        <bbox-query {{bbox}}/>
    </query>
</union>
<print mode="body"/>

and it finds nothing.
I tried also this:
<osm-script output="xml" timeout="250">
<union>
    <query type="nwr">

        <has-kv k="name" regv="Danube(.*)"/>
        <bbox-query {{bbox}}/>
    </query>
</union>
<print mode="body"/>

and it will find all 'Danube' restaurants and hotels, but not the river itself.
Can someone explain to me how to find also rivers (lines)? I need it for my historical research (toponymy). I might add I'm not great at coding.
Edit: Danube was just an example. I don't need Danube itself, but other rivers and wild card is a must.
I'll give you better example - I'm searching for the 'Stob-' rivers and I need query which will find rivers like 'Stob-nica', 'Stob-ianka' and 'Stob-niczanka' and others 'Stob-' rivers whose suffixes I don't know. Can't do that without wild card.


texts:
I am creating a python toolbox within ArcGIS 10.3 and I am having issues opening it up in ArcMap. I have been editing it, saving it, and running it within ArcMap this whole time to test it little by little, and now when I double click on the tool within the toolbox it won't respond. I checked to make sure that I do not have multiple scripts running, or that I am not hiding the pop up window for the tool behind other things open on my screen. This leads me to believe that there is something wrong with my python script (below).

import arcpy

class Toolbox(object):
    def __init__(self):
        """Define the toolbox (the name of the toolbox is the name of the
        .pyt file)."""
        self.label = "Toolbox Label Property"
        self.alias = "Toolbox Alias Property"

        # List of tool classes associated with this toolbox
        self.tools = [WiFi_Buffer]

class WiFi_Buffer(object):
    def __init__(self):
        """Define the tool (tool name is the name of the class)."""
        self.label = "Buffer WAP"
        self.description = "This python script tool will buffer WAPs"
        self.canRunInBackground = False

    def getParameterInfo(self):
        """Define parameter definitions"""
        param0 = arcpy.Parameter(
        displayName="Wireless Access Points",
        name="wap_pts",
        datatype="GPFeatureLayer",
        parameterType="Required",
        direction="Input")

        param1 = arcpy.Parameter(
        displayName="Obstruction Layer",
        name="obstructions",
        datatype="GPFeatureLayer",
        parameterType="Required",
        direction="Input")

        params = [param0, param1]
        return [params]

##    def isLicensed(self):
##        """Set whether tool is licensed to execute."""
##        return True
##
##    def updateParameters(self, parameters):
##        """Modify the values and properties of parameters before internal
##        validation is performed.  This method is called whenever a parameter
##        has been changed."""
##        return
##
##    def updateMessages(self, parameters):
##        """Modify the messages created by internal validation for each tool
##        parameter.  This method is called after internal validation."""
##        return
##
##    def execute(self, parameters, messages):
##        """The source code of the tool."""
##        return

I have left the parts that I have not reached yet commented out, and it has been running the first parameter just fine, it wasn't until trying to create the second parameter (param1) that the tool wouldn't run. I don't see anything wrong with the code, but I am new to building python tools so there could be something obviously wrong with it. Any ideas on how I can get the tool to run again?


texts: Loop through directories and mosaic all raster TIFFs in each dir

texts:
I have 200+ directories, each directory contains x amount of TIFFs which I want to mosaic.
I have a snip of code that will list the TIFF contents of all dirs under a parent directory, now I just need to parse the list (per dir), pass it to the mosaic program and run all TIFFs in the first dir, then loop to the next dir in the list and repeat the mosaic operation.
import arcpy
import os

#Set the workspace enviorment setting
arcpy.env.workspace = r"C:/Work/SpaceNet/Test2"

walk = arcpy.da.Walk(datatype="RasterDataset")

for dir_path, dir_names, file_names in walk:
    for filename in file_names:
        print(os.path.join(dir_path, filename))

Returns (from 2 dirs in this example, 00114 and 00628):
C:\Work\SpaceNet\Test2\00114 SN6_Train_AOI_11_Rotterdam_SAR-Intensity_20190822100114_20190822100413_tile_678.tif
C:\Work\SpaceNet\Test2\00114 SN6_Train_AOI_11_Rotterdam_SAR-Intensity_20190822100114_20190822100413_tile_679.tif
C:\Work\SpaceNet\Test2\00114 SN6_Train_AOI_11_Rotterdam_SAR-Intensity_20190822100114_20190822100413_tile_680.tif
C:\Work\SpaceNet\Test2\00628 SN6_Train_AOI_11_Rotterdam_SAR-Intensity_20190823100628_20190823100943_tile_5430.tif
C:\Work\SpaceNet\Test2\00628 SN6_Train_AOI_11_Rotterdam_SAR-Intensity_20190823100628_20190823100943_tile_5431.tif
C:\Work\SpaceNet\Test2\00628 SN6_Train_AOI_11_Rotterdam_SAR-Intensity_20190823100628_20190823100943_tile_5432.tif


texts:
I am designing a water distribution system for a rural village in Guatemala.  I have created contours for the area.  I need to lay out the pipe routes and need to know how to get the elevation of the lines to determine pressure drops in the system.  I assume I create a TIN surface and project my lines onto the surface.

How do I do that?

Go from a contour shapefile to a polygon surface and then draw the lines and project them onto that surface.


texts:
I’m trying to create a torque viz on the CartoDB dashboard following this StackExchange post: Using different colors on Torque CartoDB

I followed that answer, but my marker colors are still all the same.

SQL and CartoCSS are below, and here is a link to the viz: https://pgrasmussen.cartodb.com/viz/c3b82e78-8f4f-11e3-8156-0e625a1c94a6/map

Any ideas?

SQL

SELECT *,
CASE
WHEN type='Theft' THEN 1
WHEN type='Accident' THEN 3
WHEN type='Vandalism' THEN 7
END as type_number
FROM damage_nigeria

CartoCSS

/** torque visualization */
Map {
-torque-frame-count:512;
-torque-animation-duration:30;
-torque-time-attribute:"cartodb_id";
-torque-aggregation-function:"count(cartodb_id)";
-torque-resolution:2;
-torque-data-aggregation:linear;
}

#damage_nigeria_viz{
  comp-op: lighter;
  marker-opacity: 0.9;
  marker-line-color: #FFF;
  marker-line-width: 1.5;
  marker-line-opacity: 1;
  marker-type: ellipse;
  marker-width: 6;
  marker-fill: green;

}
#damage_nigeria_viz[frame-offset=1] {
  marker-width:8;
  marker-opacity:0.45;
}
#damage_nigeria_viz[frame-offset=2] {
  marker-width:10;
  marker-opacity:0.225;
}
#damage_nigeria_viz[value=3] {
  marker-fill: red;
}
#damage_nigeria_viz[value=7] {
  marker-fill: blue;
}


texts: Long description layer in legend of ArcGIS for Desktop?

texts:
I work with arcview 10.3 and i have layer with long description (green),hence the legend spread throw out the map:

I have to use a singel symbol for the layer

and not use the unique value as describe in
How to allow for multiple lines in legend?

I wonder if there a way to shorten the description line by sliding the line into another line under ?


texts: What PostGIS [topology] provided tools may be used for data quality assurance?

texts:
I'm learning about techniques to ensure GIS data quality. Googling
around, I've found a few papers about theory and best practices (here and here), and a couple of short presentations about the topic using ESRI
tools (here and here)

For example, the concepts of "tolerance" and "precision" are globally
defined in ArcGIS (I didn't know it), and I think Oracle globally
manages the concept of "tolerance" too. On the other hand, as far as I
know, PostGIS doesn't expose the user-defined precision model
supported by GEOS, and simply uses the full precision of floating
point numbers.

So, my question is: what PostGIS [topology] provided tools may be used
for data quality assurance? Has anyone experience with this topic? Or
maybe I'm just selecting the wrong tool for this.


texts:
I am working in ArcMap 10.

I am trying to divide a polygon with another polygon. I don´t want to clip any part away, but only to divide the polygon?

The polygons are stored in two different layers.

I would like to do it using the ModelBuilder tool.

I have tried the split function but it is not working.
In this tool I have to choose an input and a field that the division is based on, but none of the fields are working.


texts: Zip Code Boundary Maps for Commercial Use?

texts:
I need images of cities with features, streets, zip code names and zip code boundaries.  For example, take this screen shot:

This is very close to what I need, except not every zip code name is listed in each boundary.  If I could get city images like this, with complete zip code names and that allows commercial use, that would be perfect.

I don't need zooming or tiling.  Static images like the above are perfectly fine.  I don't even mind capturing the images myself if there's a service that will let me get this data, provided I can use the images commercially.

Do you know of any sources? I googled, but what I found was either lacking some features (as in the above) or did not allow commercial use.


texts: How to hide map to print only the legend on a PageLayout?

texts: If I'm printing a Pagelayout, is there any API for hiding the MapFrame or stop it from getting draw on the printout? I would like to only print the Legend on the paper.


texts: Unrestricted codec for creating raster files with interior pyramids files

texts:
I'm looking for a tool that will allow me to mosaic a collection of images into a file format with performance similar to MrSID or ECW, but without the licensing restrictions that those formats come with. (I believe the .sid limit is 50 mb with the ECW limit is 500 mb).  Unfortunately, the tool that I am working with doesn't recognize external pyramid files (.rrd) or .tif with internal pyramids.

Does that leave me any good options?


texts: GDAL installation on MAC OSX via conda fails, linked due to Libpoppler and font configuration

texts:
When installing GDAL via Conda, install fails, on OSX Mojave and Python 3.7:

dlopen(/Users/tomkom/anaconda3/lib/python3.7/site-packages/osgeo/_gdal.cpython-37m-darwin.so, 2): Library not loaded: @rpath/libfontconfig.1.dylib
  Referenced from: /Users/tomkom/anaconda3/lib/libpoppler.78.dylib
  Reason: Incompatible library version: libpoppler.78.dylib requires version 14.0.0 or later, but libfontconfig.1.dylib provides version 13.0.0

I am unable to locate the cause. Installed using:

 conda install -c conda-forge gdal

As per standard.

The error seems to indicate issues with libpoppler and a font library (libfontconfig). I am unsure about the role of these libraries and how to update them.


texts: In developing a PyQGIS plugin, I would like to be able to print out variable values to ensure that my program has the desired functionality. I am aware that I can print messages with something like QgsMessageLog.logMessage("hello"), but this does not work for variable values.


texts: GeoServer SLD rounding values after division

texts:
I am doing  normalisation of two fields on GeoServer.

     <ogc:Div>
       <ogc:PropertyName>PERSONS</ogc:PropertyName>
       <ogc:PropertyName>LAND_KM</ogc:PropertyName>
     </ogc:Div>

Now I want to round the result of Division to n decimal places.


texts: How to render points from CartoDB with Leaflet and Leaflet Vector Layers?

texts:
I have a table in CartoDB containing some points. I want to render these points dynamically with Leaflet. There is cartodb.js for such tasks, but it does not support full customization of Tooltips, as I have seen. Instead I found a Leaflet-Plugin for that.

In my HTML I have the following JS-Code, but it does not work except showing the Base Map. What am I doing wrong? Are there other solutions for Leaflet?

var map = new L.Map('map').setView([51.0, 10.0], 6);
var defaultLayer = L.TileLayer.provider('MapQuestOpen.OSM');
map.addLayer(defaultLayer);

var baseLayers = {
"Mapquest OSM": defaultLayer,
"OpenStreetMap German Style": L.TileLayer.provider('OpenStreetMap.DE'),
    };

// Use Leaflet Vector Layers plugin to show data from CartoDB

var cartodblayer = new lvector.CartoDB({
    user: "user",
    table: "example_data",
    scaleRange: [6, 20],
    symbology: {
        vectorOptions: {
            icon: new customIcon({
                iconUrl: "/icon/icon.png"
            })
        }
    },
    popupTemplate: '<p>{homepage}</p>',
    singlePopup: false
});

cartodblayer.setMap(map);

map.addControl(new L.Control.Layers(baseLayers, '', {collapsed: true}));


texts:
I have a customized tool developed and running on ArcGIS server 10.0. using Web ADF
The tool is used to draw a rectangle on map and mark selection to the selected area in a Graphics layer.
This tool works fine in IE (obviously for IE11, I have to set the compatibility mode on and set to the browse page).
But when I try the same tool in Chrome, the rectangle drawn enlarges to some extent on the map and selects unwanted area too.

Is there any map rendering option available in ArcObjects 10.0 to resolve the issue.


texts: Generate polygon separating date and values using qualitymosaic in GEE

texts:
I am adapting a series of scripts that I found on the internet, in particular, some of Rodrigo E. Principe. One of them links the pixels with the dates, which in theory would allow a mapping using these two parameters.
I am trying to generate polygons from burnt areas, however, those that are spatially contiguous end up becoming one polygon, even if some of them have different dates.
I managed to generate a mapping with the polygons, however, I need them to be separated coherently with the dates (similar to a featurecollection).
If anyone can help me, I will be grateful.
That's how far I got. Here is the code link:https://code.earthengine.google.com/?scriptPath=users%2Fmaximoicmbio%2Fdefault%3ADate_QualityMosaic


texts: How to run Processing algorithms in QGIS 3 Plugins?

texts:
In QGIS 3 I am trying to run Processing algorithms in a worker thread in a QGIS plugin (not a Processing plugin) using processing.run, and experience crashes.

I guess I need to provide a QgsProcessingContext when running algorithms, but have not found a way to get a proper one.

mycontext = QgsProcessingContext()
inpbuff = processing.run("native:buffer", {'INPUT': self.inpvl, 'DISTANCE': radius, 'SEGMENTS': 10, 'END_CAP_STYLE': 0, 'JOIN_STYLE': 0, 'MITER_LIMIT': 1, 'DISSOLVE': True, 'OUTPUT': 'memory:'}, feedback=None, context=mycontext)

Both with and without the context parameter, I get many identical warnings before QGIS crashes:

Warning: QObject::setParent: Cannot set parent, new parent is in a different thread

What is the recommended way of running Processing algorithms in a worker thread in a Plugin in QGIS 3?


texts: Apply cloud mask to sentinel 2 image in SNAP

texts:
I want to apply one of the cloud masks in the sentinel 2 file to the sentinel image.

The sentinel 2 image I have downloaded features files concerning cirrus and opaque clouds at 10m, 20m and 60m.

how do I apply one of these masks to the RGB image and cancel out the clouds for future processing?

Is there a specific tool? or equation?

I am using SNAP software.


texts: Windows 7 BigSQL packaged PostGIS ST_Transformation queries not working

texts:
we've recently migrated from the EnterpriseDB packaged PostgreSQL 9.5/PostGIS 2.3 to the BigSQL verision (9.6/2.3), PostGIS was installed via PGC.  Platform is Windows 7 64bit and data migration was via a database dump.

Unfortunately this has resulted in the PostGIS ST_Transform queries no longer functioning. Previously this would work -

"UPDATE secure.lp_long_mx SET nztm_geom = ST_Transform(nzmg_geom, 2193) WHERE lp_long_mx.cc_x >=2000000 AND lp_long_mx.nztm_geom IS NULL;"

Now it throws the error of -

  NOTICE:  PostGIS was unable to transform the point because either no
  grid shift files were found, or the point does not lie within the
  range for which the grid shift is defined. Refer to the ST_Transform()
  section of the PostGIS manual for details on how to config

  ERROR:  transform: couldn't project point (2.67867e+006 6.45903e+006
  0): failed to load datum shift file (-38)

If I run the same query on a Linux box where I dumped the same database it runs no problem.

"SELECT PostGIS_Full_Version();" results are these

Windows 7 64bit -

"POSTGIS="2.3.2 r15302" GEOS="3.5.0-CAPI-1.9.0 r4084" PROJ="Rel.
4.9.2, 08 September 2015" GDAL="GDAL 2.0.3, released 2016/07/01" LIBXML="2.9.4" RASTER"

Ubuntu 16.04 64bit -

"POSTGIS="2.3.2 r15302" GEOS="3.5.1-CAPI-1.9.1 r4246" PROJ="Rel. 4.9.2, 08 September 2015" GDAL="GDAL 1.11.3, released 2015/09/16" LIBXML="2.9.3" LIBJSON="0.11.99" RASTER"

Windows spatial_ref_sys on database also seems fine when I run

"select * from spatial_ref_sys where srid in (2193,27200);"

I had a hunt around for other similar problems and saw some references to file locations and checked PostGIS has correctly set the Windows GDAL_DATA Environment Variable (C:\PostgreSQL\gdal-data).

Anyone have any ideas on how to resolve or better diagnose - I'm more a GIS than database persons so am clueless on the debugging side. I'm also stuck on Windows 7 so the Linux box isn't a solution. Seems (naively) to me an issue with the BigSQL PostGIS installation.


texts:
Just a conceptual question... when writing a plugin using Plugin Builder 1.8.3, where does the script get connected to the application instance objects?

As an example (taken from eventlayerplugin.py):

class EventLayerPlugin:
    def __init__(self,  iface):
    self.mIface = iface
    ....

def run(self):
    dialog = EventLayerDialog( self.mIface )
    if dialog.exec_() == QDialog.Accepted:
        lineLayer = QgsMapLayerRegistry.instance().mapLayers()[dialog.lineLayer()]
        eventLayer = QgsMapLayerRegistry.instance().mapLayers()[dialog.eventLayer()]
        ....

where does iface and QgsMapLayerRegistry get assigned to the actual objects in the running instance of QGIS? From my beginner's perspective they seem to 'magically' refer to those objects already existing in the running instance, yet nowhere do I see them assigned to represent the instance objects. I.e. where in the plugin script are the QGIS objects made known to the script?

Clear as mud? If so, let me know and I'll try to figure out a better way to phrase this.


texts: What sub-meter geolocation solutions exist and how can I integrate them into a web and mobile mapping workflow?

texts:
I am the manager of seven water systems. What I am trying to do is map out every meter and valve on to a map and then my workers can have the maps on their iPhones so when they have to go fix a water leak or something they can just look on their map and know where the meter and the valve is.

I've plotted a few meters on my phone and then got the longitude/latitude and transferred that into Google Earth, but there are a few problems with that:

My phone is usually 10 or so yards off on my "current location". I need to be within a few feet to plot water meters.
It takes forever to take the points that I have plotted on my phone and key them into Google Earth.
Once they are on Google Earth my workers still do not have them on their phones so it does them no good while they are on the job.

So does anyone know of a much easier way to do this? I just need to plot the points very accurately and then be able to access the points from different iPhones and computers.


texts: How to select random objects in eCognition?

texts: I classified a Landsat image in eCognition and now I want to select random objects so that I can compare their class to ground-truth data and perform the accuracy assessment. How can I do this in eCognition?


texts: Not plotting NA-Values in raster with python

texts:
I have a raster that I cropped beforehand based on a shapefile with

gdalwarp -srcnodata '-1000 -255' -cutline {shape} {vrt_path} {tiff_path}

The thing is, that there are many NAs that get simply plotted as black pixels. The shapefile is ways smaller than the raster and when I plot it with rasterio, I barely can see anything of the raster. How could I mask the NA-Values out, so that they aren't part anymore of the raster?


texts: Creating a Goolge Fusion Table with GeoJSON data

texts:
I'm trying to import a csv file of data into Google Fusion Tables. Included in the data is column of data that contains a polygon for each row which is defined as GeoJSON. All the data imports without issue and I can copy the GeoJSON into geojsonlint.com to confirm it's valid.

However, if I select an individual row and click the "geocode..." link under the geometry field it always reports "Geocoding failed".

Is it actually possible to import GeoJSON polygons into Fusion Tables? If so, how is it best to achieve this?

I know KML can be imported and I'm sure I've seen GeoJSON also mentioned, but this was in the API documentation rather than the regular UI documentation.


texts: How to control the priority of individual feature labels in one layer?

texts:
I have a group of city in a small area with labels.
Due to label engine it shows them in a 'random order' (and certainly optimizing) according the space the have.

There is a way, without duplicate the layer, to prioritize just ONE city name to be show whatever the conditions are, the others showed cities order around doesn't matter.


texts: Return all IDs of point features within a polygon to one field

texts:
On QGIS 2.14.14
I am trying to compile Old Land Parcel ID codes on to New Land Parcel geometries. The Old Land Parcel IDs are recorded as point data and the New Land Parcel geometries are polygons. Most of the Old Land Parcel IDs are the same as the New Land Parcel IDs but some of the New Land Parcels have multiple Old Land Parcel IDs codes associated with them as the Old Land Parcels have been combined to make the New Land Parcels.

I would like to have a field on the New Land Parcel layer showing all the Old Field Parcel IDs associated with the New Field Parcel as a string. Or more simply put - return all IDs of point features within a polygon to one field.

I am currently using AutoFields to create a live updatable field as new land parcels will be added to the project as new areas are addressed. Within this, I have used the expression

geomcontains( 'Old Land Parcels', 'Old_Parcel')

The problem is this will only return the first Old Land Parcel Id it comes to. Is there a way to repeat this so all contained features are recorded in the one field as a string?

Point Shapefile Old Land Parcel field Old_Parcel
Polygon Shapefile New Land Parcel field PARCEL_ID intended new field Old_Parcels(this has now changed to OldParcels due to character limit)


texts: Join completed, can't set field from the table add for my graduated map

texts:
I've two thing that I joined, a shapefile of PollingStation and a csv of the results of the election. When I join them by the right field, I get a new table with all the information I need.

My problem is when i go to the layer properties and I choose a graduated map I can't choose field from the dataset I add.

I've seen that my data are as string for the dataset I add, and as double for the original.

So i extracted again my data from R, but this time I first set it as characterto be sure. But when I joined it's always string.

I tried to refactor with the refactoring field tool, so I set all as double,

the tool says that it worked

but when I open the properties my data are always as string.

So Basically what I need is a way to join my data and make the field from the dataset I add usable for a graduated map.


texts: Filtering the number of pictures per month of a time series in Google Earth Engine

texts:
I'm calculating the NDVI of a time series (1 year) of Sentinel-2, but I only want two images per month to calculate that instead of the whole collection of images with less than 75% cloud.

//NDVI for Sentinel 2

//Call shapefile

Map.addLayer(table, {}, 'aneis', true)

// Create image collection of S-2 imagery for the perdiod studied
var S2 = ee.ImageCollection('COPERNICUS/S2')

//filter start and end date
.filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 75))
.filterDate('2018-08-08', '2019-08-08')

//filter according to shapefile
.filterBounds(table);

// this generates a JSON list of the images (and their metadata) which the filters found in the right-hand window.
print(S2)

// Function to calculate and add an NDVI band
var addNDVI = function(image) {
return image.addBands(image.normalizedDifference(['B8', 'B4']));
};

// Add NDVI band to image collection
var S2 = S2.map(addNDVI);
// Extract NDVI band and create NDVI median composite image
var NDVI = S2.select(['nd']);
var NDVImed = NDVI.median(); //I just changed the name of this variable ;)

// Create palettes for display of NDVI
var ndvi_pal = ['#d73027', '#f46d43', '#fdae61', '#fee08b', '#d9ef8b',
'#a6d96a'];

// Create a time series chart.
var plotNDVI = ui.Chart.image.seriesByRegion(S2, table,ee.Reducer.mean(),
'nd',500,'system:time_start', 'system:index')
              .setChartType('LineChart').setOptions({
                title: 'NDVI short-term time series',
               hAxis: {title: 'Date'},
                vAxis: {title: 'NDVI'}
});

// Display.
print(plotNDVI);

// Display NDVI results on map
Map.addLayer(NDVImed.clip(table), {min:-0.5, max:0.9, palette: ndvi_pal}, 'NDVI');


texts: Creating multiple fields with ranks using ArcPy?

texts:
While figuring out what I have to do, practically creating a field for each arrondise, and in each, rank each FID (which are part of the respective arrondise) by the size of its shape area.

Hopefully, it is a bit clear now.

EDIT
Currently, trying to rank each row which is part of a certain arrondise (ex. 1, 2) in terms of its shape area.

import arcpy

arcpy.env.workspace = "D:/M1 Geomatique/Programmation II/Dossier"

fc = "zones_rattachement.shp"

try:
    fieldRoot = "RANG_R"
    for counter in range(1,21):
        arcpy.AddField_management(fc, fieldRoot + str(counter),'SHORT')

    size_rank = 1
    numlist = list(range(1,21))
    for num in numlist:
        arcpy.SelectLayerByAttribute_management(fc, "NEW_SELECTION", "arrondisse = '%c'")
        rows = arcpy.UpdateCursor(fc, sort_fields="shape_area D")
        for row in rows:
            row.setValue("RANG_R1", size_rank)
            size_rank += 1
            rows.updateRow(row)

except:
    arcpy.GetMessages()


texts: Adding Description, Use Constraints, and Credits to Map Service when publishing with ArcPy?

texts:
I am writing a Python script to automate publishing Map Services for use in ArcGIS Online. My workflow follows the sample code in Example 1 from this ArcGIS Help page on the CreateMapSDDraft tool: Create a Service Definition Draft (.sddraft) file with CreateMapSDDraft, stage the service definition with StageService, then UploadServiceDefinition.

However, there are parameters available in the "point-and-click" process of publishing a service in ArcMap that I'm unable to add with the syntax of the CreateMapSDDraft tool. In the Item Description window of the ArcMap Service Editor, I want to add a Summary, list of Tags, Description, Access Constraints, and Credits before publishing a service.

The CreateMapSDDraft arcpy syntax only allows me to add a Summary and Tags to my item.
i.e., CreateMapSDDraft (map_document, out_sddraft, service_name, {server_type}, {connection_file_path}, {copy_data_to_server}, {folder_name}, {summary}, {tags}).
I cannot find a place within the help files of the other geoprocessing tools used in my script to add/alter the Description, Access Constraints, or Credits of my service.

What can I add to my script to allow me to define these attributes in the service that is ultimately uploaded?


texts: Force Legend to Always display in ArcGIS Online Web App

texts: Is it possible to force the legend in ArcGIS Online Web App (Web App Builder, not developer edition) to always display the legend. Perhaps in the JSON file somewhere?


texts: UpdateCursor cannot acquire lock if attribute table is opened in ArcGIS Pro

texts:
I wrote a script tool for someone to allow them to clean all whitespace errors (leading/trailing spaces, 2+ consecutive white spaces, convert empty strings to nulls), and aside from one issue it works fantastically. The issue is that, when you run the tool against a layer in ArcGIS Pro and you have that layers attribute table open, the tool cannot acquire a lock. If you close the attribute table and run the tool against the same layer in the table of contents, it runs fine.
I'm wondering if I did something wrong, or if opening the table view actually is causing this locking error and there's no way around it?
This script has much room for improvement and was created out of necessity in about an hour, feel free to point out other issues as you see fit but really my main goal is to figure out why this is causing a locking issue.
"""
Tool:               <Whitespace Cleaner>
Source Name:        <WhitespaceCleaner>
Version:            <2.7>
Author:             <John Ehlen>
Usage:              <Toolbox interface>
Required Arguments: <Input Feature>
                    <Field(s)>

Description:        <Allows user to select a feature class, and one or more target fields to clean duplicate whitespace, line breaks, and new line characters>
"""
import arcpy, os, sys, re

def ScriptTool(parameter0, parameter1):
    input_Feature = parameter0
    dataType = arcpy.Describe(input_Feature).dataType
    input_Fields = parameter1
    if dataType == 'FeatureLayer':
        arcpy.AddMessage('Input Feature: {0}'.format(input_Feature.name))
    arcpy.AddMessage('Fields to clean: {0}'.format(input_Fields))

    #Convert Params to Field Names for Cursor
    fieldNames =  input_Fields.split(";")

    #Identify workspace
    if dataType == 'FeatureLayer':
        workspace = input_Feature.connectionProperties["connection_info"]["database"]
        arcpy.AddMessage('Workspace: {0}'.format(workspace))

    with arcpy.da.UpdateCursor(input_Feature, fieldNames) as cursor:
        for row in cursor:
            attribute_list = list(row)
            new_l = [None if not item or item.isspace() else re.sub('\s{2,}', ' ', item).strip() for item in attribute_list]
            row = tuple(new_l)
            cursor.updateRow(row)
    return

if __name__ == '__main__':
    # ScriptTool parameters
    parameter0 = arcpy.GetParameter(0)
    parameter1 = arcpy.GetParameterAsText(1)

    ScriptTool(parameter0, parameter1)


texts: Python GDAL/OGR apply WGS84 CRS to Polygon

texts:
I am trying to work out how to apply WGS84 to a POLYGON I am creating from a list of latitude and longitude coordinates.  Without the correct coordinate system, I cannot calculate the area, my intersection is also failing.

# Create Polygon from formatted string of coordintes
wkt1 = "POLYGON (({0}))".format(coords_fin)
poly1 = ogr.CreateGeometryFromWkt(wkt1)

I know I can define the WGS84 geometry and have successfully done this for other tasks:

# Set up the coordinate reference system, WGS84
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)

What I can't find is how to say my polygon is in this CRS.  All the examples I have managed to find are based on exporting results as shapefiles, etc, so are applying the CRS to a layer.

I simply want to define two polygons in WGS84 (sanity check the area), then return a False if they don't overlap using:

intersection = poly1.Intersection(poly2)

I have tried SetGeometry but it doesn't appear to be an attribute of the polygon.


texts:
I have this code that works perfectly for small rasters.
var image = ee.Image("users/parivash89/Cuneo_Reclassify")
  .toInt();

Map.addLayer(image);
Map.centerObject(image);

var classes = image.reduceToVectors({
  reducer: ee.Reducer.countEvery(),
  scale: 10,
  maxPixels: 1e10
});

var result = ee.FeatureCollection(classes);

Map.addLayer(result);

// Export the FeatureCollection to a SHP file.
Export.table.toDrive({
  collection: result,
  folder: 'GEE_Polygonize',
  description:'Cuneo',
  fileFormat: 'SHP'
});

But gives me the following error while the raster input is large:
Layer 2: Tile error: Geometry has too many edges (3890235 > 2000000)

So I added this line :
// Simplify the feature
var featureSimple = ee.Feature(result).simplify({maxError: 10});

But it did not get simplified and I stil get the error. What should I change in my code? I need the resolution of 10m in my output so I cannot lower the resolution.


texts: Intervisibility tool does not consider the obstruction in ArcGIS

texts:
I have a 3D point outside a building as shown in the figure. I utilized the Construct Sight Lines tool to generate lines between the point and points which represent the vertices of the building.

Now, on running the Intervisibility tool with the building (multipatch) as an obstruction, all the lines are assigned a value of 1 as shown in the figure, which means the end and start points are visible.

The obstruction is not considered.


texts: While saving shapefile to disk, is there a way to set the geometry type and SHPT

texts:
QGIS 3.0.1

I wrote a code to create a vector layer and save it to disk as shapefile .
Is there a way I can save its geometry property as 'Point' and SHPT as 'Point'

from qgis.core import *
import qgis.utils

QgsApplication.setPrefixPath("C:/Program Files/QGIS 3.0", True)
qgs = QgsApplication([], False) ###
qgs.initQgis() ###

_layer1 = QgsVectorLayer("C:/Arundhati/Automation Test Script/Radio Signal/MasterLine.shp", "MASTERLINE", "ogr")

if not _layer1.isValid():
  print("Layer failed to load!")

else:
    print("Layer was loaded successfully!")

_writer1 = QgsVectorFileWriter.writeAsVectorFormat(_layer1,"C:/Arundhati/Automation Test Script/Radio Signal/MasterLine12.shp",'utf-8',QgsCoordinateReferenceSystem(),'ESRI Shapefile')

qgs.exitQgis()

please suggest.


texts:
I have a KML file for a NASA image. I want to create a world file from it so I can use the image in say QGIS. Is it possible to extract the key info (bounding coordinates) from the KML? It seems to be made up of lots of links to smaller tiles of the image

There is already another post on here (I can't post any more links, but search for "export-georeferenced-image-overlay-from-kml" which points to WorldFileTool and a KMZ to worldfile converter

The first (WorldFileTool) seems to never link to the all important zip file

The second throws an error when I run it.

Hence I am a bit stumped. I can georeference it by hand if needs be, but I had hoped the KML might give me the easy answer.

Any suggestions would be great.


texts:
I'm looking for a solution to snap my GPX tracks to OSM roads. Mapbox Map Matching API actually does this, but sometimes snaps to wrong roads. In this wrong case I want manually edit snapped track by dragging its vertices. When dragging vertices they should be snapped to OSM roads as well.

what could be good solution to implement this kind of task?

I do not want to download and store OSM data locally. If it is possible also not to use routing engines.


texts:
I am using QGIS 2.8 version
I am invoking QGIS from a third-party application, using the qgis-bin.exe
Along with the exe i am passing the raster image as an argument that QGIS will display on loading . In my plugin I want to know whether QGIS was invoked from a third-party s/w or it was launched in the standalone mode.
I was passing a flag as command line argument along with the raster and was checking it OnInitializationCompleted Event:

qgis-bin.exe [Raster Image] --flag

But QGIS always try to add all the arguments as layer, so I am getting a warning every time when QGIS is started from that application.

I have read some posts, they talk about calling a script on loading via --code option.
Is there any other way.


texts: Benefit of using GlobalIDs in ArcSDE geodatabase?

texts: This may be somewhere, but I can't find a clear answer on exactly the benefit of GlobalIDs are in a database.  What are the benefits of using these in an organization that has about 6 people editing the same data.


texts: Is it possible to load a WMS aerial photo from a external map into QGIS?

texts: Is it possible to load a WMS based aerial photo (like https://geoinfo.frankfurt.de/luftbild-2018) into QGIS and use it as a background-layer for internal use?


texts: ST_Intersects with degenerate LINESTRING


texts:
My script is supposed to pull weekly flu data from a department of health website, parse it for that week, and then push the data up to our ArcSDE server so that it can be displayed on our Web Mapping application.  The data is in a .csv format and the script parses the data just fine, but when it comes to doing a join with a feature class with county outlines I get a

  Traceback (most recent call last):
   File "C:\Users\username\Desktop\Projects\flu_mapping\script\flu_data0.6.0.py", >line 55, in
     arcpy.JoinField_management ("memorycounty", "ID", outputcsv, "Statecode",)
   File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\management.py", line 5383, >in JoinField
     raise e
  RuntimeError: Object: Error in executing tool

This wasn't an issues until I started using the open() command to create a new .csv file to parse the data to, instead of manually making one.

Here is my code

import os
import arcpy
from arcpy import env
import math
import numpy
import csv
import urllib2
import tempfile

#working directory

fludirectory = 'C:\\Users\\username\\Desktop\\Projects\\flu_mapping\\GIS_data\\'

#pull flu data from ND Department of health
fludata = urllib2.urlopen ('http://www.ndflu.com/data/weekenddata.txt')

#arcpy workspace
arcpy.env.workspace = "in_memory"

#load data so that numpy can use it
rows = numpy.loadtxt(fludata, int, delimiter=",",skiprows=4, usecols =(2,3,4,5,6,7,8,9))

#parse the data and write it to the CSV file
createcsv = fludirectory + "\\" + 'output_data.csv'
outputcsv = open (createcsv, 'wb')

week = max(rows [:,0])

with outputcsv as csvfile:
    flucsv = csv.writer(csvfile, delimiter=",",
                            quotechar='|', quoting=csv.QUOTE_MINIMAL)
    flucsv.writerow (["Statecode","flu_num"])
    for row in rows:
        if row[0] == week:

            countid = row[1]
            totalflu = sum(row[2:])
            csvoutput = '%s, %s' %(countid, totalflu)
            print csvoutput
            flucsv.writerow ([countid, totalflu])

outputcsv.close

#make a copy of the county dataset

#county dataset
countydata = fludirectory + "\\" + 'county_ID_code.shp'
#make a temporary feature layer
arcpy.MakeFeatureLayer_management (countydata, "memorycounty")

#join output csvoutput to the in_memory shapefile
arcpy.JoinField_management ("memorycounty", "ID", outputcsv, "Statecode",)
#test the script
arcpy.FeatureClassToFeatureClass_conversion ("memorycounty", fludirectory, 'testoutput.shp')

Also at the end I am using Feature Class to Feature Class to see if everything worked (obviously it's not) and that will be replaced by code that pushes the feature class to the ArcSDE server.


texts: "libSFCGAL.so.1: undefined symbol" upgrading to Postgis 2.2.2 (and 2.2.1)

texts:
I recently got PostGIS upgraded to version 2.2.2. I would like to upgrade a legacy database to this new version, running this query:

alter extension postgis update to "2.2.2"

The result is this error:

ERROR: could not load library "/usr/lib/postgresql/9.3/lib/postgis-2.2.so": /usr/lib/x86_64-linux-gnu/libSFCGAL.so.1: undefined symbol: _ZN5osgDB13writeNodeFileERKN3osg4NodeERKNSt7__cxx1112basic_stringIcSt11char_traitsIcESaIcEEEPKNS_7OptionsE
    SQL state: 58P01

I have the SFCGAL library installed with headers and debug symbols:

$ dpkg -l | grep sfcgal
ii  libsfcgal-dev                                               1.2.2-1                                                     amd64        Library for ISO 19107:2013 and OGC SFA 1.2 for 3D operations (dev part)
ii  libsfcgal1                                                  1.2.2-1                                                     amd64        Library for ISO 19107:2013 and OGC SFA 1.2 for 3D operations (runtime part)
ii  sfcgal-bin                                                  1.2.2-1                                                     amd64        Simple viewer for SFCGAL library
ii  sfcgal-dbg                                                  1.2.2-1                                                     amd64        Debug symbols for SFCGAL

How can I overcome this?

Update

I was not able to find a package source for a more recent version of the SFCGAL library. Therefore I decided to re-install the software stack from a single repository, the UbuntuGIS PPA. Now I have the following set up:

postgresql-9.5
postgresql-9.5-postgis-2.2 (2.2.1+dfsg-3~xenial0)
postgis (2.2.1+dfsg-3~xenial0)
libsfcgal1 (1.2.2-1)

But with this new setup I get the same error:

# CREATE EXTENSION postgis;
ERROR:  could not load library "/usr/lib/postgresql/9.5/lib/postgis-2.2.so": /usr/lib/x86_64-linux-gnu/libSFCGAL.so.1: undefined symbol: _ZN5osgDB13writeNodeFileERKN3osg4NodeERKNSt7__cxx1112basic_stringIcSt11char_traitsIcESaIcEEEPKNS_7OptionsE


texts: ArcGIS 10 export metadata to XML programmatically

texts:
I have some code that exports metadata from an iMetadata instance to a standalone XML file. This code uses the ExportXML type and the iMetadataExport interface (see below).

Private Sub ExportMetadatatoXML(ByVal Metadata As IMetadata, ByVal strXMLdocName As String, ByVal strOutputXmlFolder As String)
    Try

        Dim ExportXML As ExportXML
        Dim xml_Export As IMetadataExport
        ExportXML = New ExportXML
        xml_Export = ExportXML
        If System.IO.Directory.Exists(strOutputXmlFolder) Then
            xml_Export.Export(Metadata, strOutputXmlFolder & "\" & strXMLdocName & ".xml")
            Logging.WriteToLog("Metadata Exported to " & strOutputXmlFolder & "\" & strXMLdocName & ".xml", Serverity.SeverityLevel.Information)
        Else
            Logging.WriteToLog("Output XML Directory " & strOutputXmlFolder & " cannot be found", Serverity.SeverityLevel.Fatal)
            Result = Result - 1
        End If

        Logging.WriteToLog("Metadata was not exported to " & strOutputXmlFolder & "\" & strXMLdocName & ".xml", Serverity.SeverityLevel.Information)

    Catch ex As Exception
        Dim strRoutine As String = "Error thrown in: " & System.Reflection.MethodBase.GetCurrentMethod.Name
        Logging.WriteToLog(strRoutine & " - " & ex.Message, Serverity.SeverityLevel.Fatal)
        Result = Result - 1
    End Try

End Sub

This code worked a treat for 9.x versions of ArcGIS, but with the new metadata changes for 10+, it no longer works. The ExportXML type seems to no longer exist and I can't seem to figure out a replacement. Does anyone know how I can achieve this with 10+? I've noticed there's a metadata toolbox but would rather avoid using geoprocessing. Thanks for your time.

Edit: I have since tried using geoprocessor to export to XML, but i'm getting a hresult error on gp.execute. Can anyone suggest what the problem might be?

        Dim GP As ESRI.ArcGIS.Geoprocessor.Geoprocessor = New ESRI.ArcGIS.Geoprocessor.Geoprocessor
        Dim parameters As IVariantArray = New VarArray
        parameters.Add("C:\Temp\Temp.gdb\BLOCKS")
        parameters.Add("C:\Program Files\ArcGIS\Desktop10.1\Metadata\Translator\ARCGIS2FGDC.xml")
        parameters.Add("C:\Temp\Metadata\blocks.xml")
        GP.Execute("ExportMetadata_conversion", parameters, Nothing)


texts: Mapping bigger grid attributes to overlapping smaller overlapping grids

texts:
We have following grid structure.

There are certain attributes attached to the Bigger grid (80m). However we want to map all those attributes to smaller grids as it also contains some other attributes and is what we want to use. One approach is to use weights (by overlapping area percentage), but there are number of drawbacks. So is there any other approach we can use for this.


texts:
I'm trying to use a csv table with angle and azimuth values to compute a hillshade on the same area at different dates.

I'm using for and each date to extract the correct values, but it doesn't seem to be working out. It produces rasters with the exact same hillshade (presumably based on the values from the first line).

Date   Angle   Azimuth
Jan01  69.77   225.41
Jan06  70.50   225.56
Jan11  71.32   226.24

for(d in insol$Date){
hillShade(slope, aspect, angle=insol$Angle, direction=insol$Azimuth, filename=paste(d,".tif"))
}

Is there a way to tell R to use the values from each row at a time? The output states (sorry, in Portuguese):

> warnings()
Mensagens de aviso:
1: In cos(slp) * cos(zenith) :
  comprimento do objeto maior não é múltiplo do comprimento do objeto menor
2: In sin(slp) * sin(zenith) :
  comprimento do objeto maior não é múltiplo do comprimento do objeto menor
3: In cos(slp) * cos(zenith) :
  comprimento do objeto maior não é múltiplo do comprimento do objeto menor
4: In sin(slp) * sin(zenith) :


texts: How do I read the attribute data associated with my SHP?

texts:
I'm using DotSpatial to aid in processing some geospatial data. I'm using the code below to determine if a point (lat/lon location) is within a polygon contained in the SHP. If the point is in a polygon, I want some of the attributes stored in the DBF associated with the polygon. How do I access the attribute data?

// Reproject the strange shapefile so that it is in latitude/longitude coordinates
wierdShapefile.Reproject(DotSpatial.Projections.KnownCoordinateSystems.Geographic.World.WGS1984);

// Define the WGS84 Lat Lon point to test
Coordinate test = new Coordinate(longitude, latitude);

foreach (Feature f in wierdShapefile.Features)
{
    Polygon pg = f.BasicGeometry as Polygon;
    if (pg != null)
    {
        if (pg.Contains(new Point(test)))
        {
            ******************************************
            *** Need to access attribute data here ***
            ******************************************

            // If the point is inside one of the polygon features
            return true;
        }
    }
    else
    {
        // If you have a multi-part polygon then this should also handle holes I think
        MultiPolygon polygons = f.BasicGeometry as MultiPolygon;
        if (polygons.Contains(new Point(test)))
        {
            return true;
        }
    }
}

return false;
}


texts: Overlaying re-projected points and polygons in R?

texts:
I am trying to overlay point data with polygons in R in order to find the polygons which intersect with the points provided.

The assciated prj file from the polygons shows their coordinate reference system: "DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137.0,298.257222101]]".

The file is available here http://mapas.ine.pt/download/index2011.phtml.
Then I have sample points. The point coordinates come from the Google API and are transformed to WGS 1984.

The point data set is here https://www.dropbox.com/s/1r1rzmhjhqcsarv/Points.csv?dl=0.

But when I try to overlay both, an dataframe is returned filled with NAs.

When I plotted both files in QGIS for visual checking, I found that all most all points are aggregated in one polygon. However they should spread over Portugal,(even beyond the polygon outcrop which is Lisbon)

There should be a corruption of the coordinate system at one point, but I cannot find it!

I set a common coordinate reference system for polygons and points (WGS 1984) and repaired the polgyon with cleangeo. Below the commented code I used.

(A)  Load points and polygons

(A1) Define common CRS

CRS.new <- CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")

(A2) Read points from csv and build SpatialPointDataFrame

Points <- read.csv("Points.csv", header = TRUE, sep = ";")

Points_coords <- Points[,c(3,4)]
Population <- SpatialPointsDataFrame(coords = Points_coords, data = Points,proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))

(A3) Read polygons from BGRI dataset

Portugal_subsections <- readShapeSpatial("BGRI11_LISBOA.shp")
Portugal_subsections@proj4string <- CRS.new

(A3) Fix polygons with cleangeo

report <- clgeo_CollectionReport(Portugal_subsections) ## Check topology error
clgeo_SummaryReport(report)

Portugal_subsections_clean <- clgeo_Clean(Portugal_subsections) ## Clean topology errors

report2 <- clgeo_CollectionReport(Portugal_subsections_clean)  ## Check again
clgeo_SummaryReport(report2)

(B) Saving files

(B1) Portugals subsections

Portugal_subsections_clean@proj4string  <- CRS.new
writeOGR(obj=Portugal_subsections_clean, dsn = ".", Portugal_subsections_clean", driver="ESRI Shapefile")

(B2) Saving Population as point shape

colnames(Population@data) <- c("Index","Zip", "lon", "lat", "Population")
writeOGR(obj=Population, dsn = ".", "Population", driver="ESRI Shapefile")

(C) Overlay Population with BGRI subsections

Portugal_subsections_clean@proj4string <- CRS.new  ## Assure same CRS (WGS 84)
Population@proj4string <- CRS.new  ## Assure same CRS (WGS 84)

a <- over(Population, Portugal_subsections_clean)


texts: Strange behavior of 'qgis:joinattributesbylocation' in PyQGIS

texts:
Now, I'm trying the qgis:joinattributesbylocation in PyQGIS and I'm just facing the problem. Like the figure, I have two polygons and one is the plant community polygon (Figure1) and the other is the polygon which display the area range (Figure2).

After splitting the plant community polygon with the lines of area range polygon,  I conducted the qgis:joinattributesbylocation and wanted to assign the area range information (like 45, 46, 47 in figure).

Area range information were provided to some parts of the polygon  without any problems (pink color polygons in figure1), however, some of the polygons were failed to assign the area range information.

Why these problem was happened and any solutions?


texts: Transforming UTM layer to WGS using QGIS?

texts:
I am trying transform a UTM layer to a WGS.

I have already:

Saved it as a new layer changing the SRC
Tried the tool "define current projection"

I am using QGIS 2.14 and the original layers are in WGS84/UTM 17S EPSG:32717, and my objective is to use WGS84.


texts:
I'm trying to pass in a geometry(polygon) collected as an ESRI.ArcGIS.Client.Geometry.Geometry type on a map and passing it as an IN parameter into an oracle sql query. I'm also trying to use st_intersects to find features that intersect the passed-in geometry.

I can't seem to get my syntax correct. Here is my sql query:

WITH cons
     AS (SELECT DISTINCT culvertgid , condition AS condition_code,
                         value     AS condition_literal
         FROM   culvert_inspect
                LEFT OUTER JOIN (SELECT Extractvalue (CodedValues.column_value,
                                        'CodedValue/Code') AS Code,
                                        Extractvalue (CodedValues.column_value,
                                        'CodedValue/Name') AS VALUE
                                 FROM   sde.gdb_items_vw items
                                        inner join sde.gdb_itemtypes itemtypes
                                                ON items.TYPE = itemtypes.uuid,
                                        TABLE (
        Xmlsequence (Xmltype (DEFINITION).EXTRACT (
                     '/GPCodedValueDomain2/CodedValues/CodedValue')))
                        CodedValues
                 WHERE  itemtypes.name = 'Coded Value Domain'
                        AND items.name = 'OVERALL_CONDITION')
                coded_values
             ON culvert_inspect.condition = coded_values.code)
SELECT DISTINCT ci.condition_code, condition_literal
FROM   culverts c
       INNER JOIN  cons ci
               ON ( c.globalid_1 = ci.culvertgid )
WHERE sde.st_intersects(c.shape, sde.ST_Polygon (-10685820.278500814 5060357.6856116494, -10685820.278500814 5101327.9327725023, -10645461.527566245 5101327.9327725023, -10645461.527566245 5060357.6856116494, -10685820.278500814 5060357.6856116494)) = 1;
ORDER  BY condition_code

The code before the WHERE statement gets the domain descriptions from a table called culvert_inspect. What I'm trying to do is return only the domain values that are in the features that are within the geometry that is getting passed in. Everything works up until the WHERE statement. I've tested this same query and instead of passing in geometry, I pass in a value from a field and return expected results.


texts: How load data into file geodatabase tables with relationship classes to non simple data?

texts:
Further to my earlier question (Loading data into file geodatabase without rebuilding relationship classes manually?), is there any way to load data into GDB tables with relationship classes to composite data as geometric network without making them simple first ?

I have a file geodatabase with few M-N composite relationships classes that contain geometric network. Every few weeks or month I need to load new data into GDB tables; for this I need to delete the relationship classes delete the geometric network, and than load the data using The Simple Data Loader wizard and then rebuild the relationship classes and rebuild the geometric network. when I'm trying this without deleting first I keep getting error message:

  Cannot load into non simple data


texts: Using Search Cursor to loop through a table and use as variable

texts:
I've a fairly extensive script which requires the user to manually define a variable which populates rest of script such as file paths, file names and outputs etc. Now I'm trying get the scrip to loop through a table of file prefixes (I've these available as polygons geodatabase as well). I thought best way to accomplish this would be using Search Cursors on a table
where I've file prefix saved in each row, so if I were to set up a loop it would grab the row value, substitute it in place of the variable (coName) and then goes through rest of the steps in the script. So I've tried the following, but I'm having a difficult time trying to figure out:
(1) how to tell SC to obtain a value and then insert it as the variable
(2) in other variations it pulls the loops through the table and then ends up using the last row item in the table.

import arcpy,os
from arcpy import env
from arcpy.sa import *

# Table with prefix file names such as ABCD_10001, ABDC_10002 etc.,
fTable = "C:/test_table.txt"
field = "FileName"

cursor = arcpy.da.SearchCursor(fTable,field)

# Iterate through each rows
row = cursor.next()
while row:
    print (row.getValue(field))
    cursor.next()

    coName = row
    wFolder = "C:/Imagery/"
    coGDB = os.path.join(str(coName) + ".gdb/")

    env.workspace = os.path.join(str(wFolder), str(coGDB))
    print (arcpy.Exists(env.workspace))

    DEV = os.path.join(str(wFolder), str(coGDB), str(coName) + "_DEV")
    print (arcpy.Exists(DEV))

    outRas = Raster(DEV) * 2
    outRas.Save(str(coName) + "_x2")

UPDATE

I'm using ArcGIS 10.4. The last two lines are a psuedo map algebra equation.


texts:
I am using OL3 and I try to inspect whether two lines are intersecting or not using turf. I used GeoJSON in OL3 and fed the lines into the intersect() in turf but I get Unknown GeoJSON type.
here is my code:

lineStringdraw.on('drawend', function(e2) {
var featureEnd = e2.feature;
checkIntersection(featureEnd);
});

function checkIntersection(currentFeature){
var features = vectorLayer.getSource().getFeatures();
var geoJSON = new ol.format.GeoJSON();
geoJSON = geoJSON.writeFeature(currentFeature);
for (i=0; (features.length)>i;i++){
var feature = features[i];
geoJSON_i = geoJSON.writeFeature(feature);
var intersectionStatus = turf.intersect(geoJSON,geoJSON_i);
$('#infoList').append(intersectionStatus);
}
}

here is piece of the OL3 GeoJSON:

{"type":"Feature","geometry":{"type":"LineString","coordinates":[[-8236104.129670544,4977029.890789877],[-8236634.411554272,4974407.145256842]]},"properties":null}

Does anyone know why it is not working or what I am missing?


texts:
I have stand alone GPS points and RTK points collected at the same time and I want to compute the accuracy (1 sigma) of the GPS points. Using the haversine formula, I computed the distance between the GPS points and corresponding RTK points

De = harv(lat_rtk, log_rtk, Lat_gps, Log_gps)

then I calculated the standard deviation of

sigma_De = std(De)

I'm using the Python.

How to find the Sigma_north and Sigma_east. The mean of error is not zero should the mean of error added to sigma_DE to get the true accuracy.

I tried to use mathematics format notation $x$ but it did not work.


texts:
Did ESRI drop their high res. basemaps. I made a map in November and today I needed to make a small adjustment to it. The high resolution basemap is gone and my map looks bad. I have attached the 2 maps to show the vast difference. Has ESRI done something or is there a setting that I am missing?
I am using ArcGIS 10.4.1 for desktop.


texts: Creating training set of strike and dip data (structural geology) using ArcGIS Pro

texts:
I have been using ArcGIS Pro for the first time, and am attempting to create a training set of strike and dip data (structural geology.) I have created the polygons, after creating a new classification schema, but the files are not available for import when I got to the classification wizard. I am confused as to where I may have made a misstep, or where the program is getting confused.

For what it's worth, I have followed the Esri walkthroughs on a few deep learning/image classification projects and I think I was authentic to their directions.


texts:
I am creating a raster using GDAL and the IMAGINE file doesn't work in ERDAS because it doesn't recognize the projection.

Here is the code to generate a sample:

from osgeo import gdal

projectionWkt = """PROJCS["Albers Conical Equal Area",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AUTHORITY["EPSG","4269"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",29.5],PARAMETER["standard_parallel_2",45.5],PARAMETER["latitude_of_center",23],PARAMETER["longitude_of_center",-96],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["meters",1]]"""

outputRaster = gdal.GetDriverByName("HFA").Create("test.img",10,10,1,gdal.GDT_FLoat32)
outputRaster.setProjection(projectionWkt)

This file does work in ArcMap, but can't be used in ERDAS.

Am I doing something wrong with GDAL?


texts:
I am trying to use RemoveLayer within the mapping module, but seem to be getting mixed results. I do the following to add a layer to the map

mxd = arcpy.mapping.MapDocument("CURRENT") # define map document
lyr = arcpy.mapping.Layer(r"C:\test.lyr") # define layer object
df = arcpy.mapping.ListDataFrames(mxd) # list data frames
arcpy.mapping.AddLayer(df[0],lyr) #  add layer

Now that the layer is added, I am trying to use arcpy.mapping.RemoveLayer to remove the layer from the map. I do the following.

arcpy.mapping.RemoveLayer(df[0],lyr)

But this doesn't work. The layers still exist in the map. But if I try the following code where I list the layers and then loop through and use RemoveLayer, it works now. Why?

lyrlist = arcpy.mapping.ListLayers(mxd)
for i in lyrlist:
    arcpy.mapping.RemoveLayer(df[0],i)

Is this because Layer objects don't work with RemoveLayer for some reason?


texts:
Is there a tool to do conflation in QGIS?

ArcMap has a tool but I need to use QGIS. http://pro.arcgis.com/en/pro-app/tool-reference/editing/an-overview-of-the-conflation-toolset.htm


texts: Making Upload button for shapefile in QGIS plugin?

texts:
I am new to Plugin Builder and creating my first plugin called NearestNeighbour. The below is GUI of my plugin made using Qt Creator

I read in this link that it is important to rename Dialog Box, Toolbox because it  will be used in our code.
Label are kept as it is with object name label and label_2 respectively. The LineEdit  objectName for Upload Node Layer is called lineEdit. On a similar basis the second LineEdit  objectName for Upload Building Layer is called lineEdit_2. In a similar way the Push Buttonfor Upload Node Layer is called pushButton and the second Tool Button for Upload Building Layer is called pushButton_2. With reference to the link mentioned above, I have done the Following change in my Nearest_Neighbou.py

1) Add the below line where all imports are done

 from qgis.core import QgsProject

2) Add QFileDialog to QtWidget to the list of imports at the top of the file.

from PyQt5.QtWidgets import QAction, QFileDialog

3) Add a new method called Upload_Point_Layer with the following code just above def(run)self:

def Upload_Point_Layer(self):
  filename, _filter = QFileDialog.getSaveFileName(
  self.dlg, "Select   output file ","", '*.shp')
  self.dlg.lineEdit.setText(filename)

4) In the run method, added the following line in the block where the dialog is initialized.(This code will connect the Upload_Point_Layer method to the clicked signal of the push button widget)

 self.dlg.pushButton.clicked.connect(self.Upload_Point_Layer)

I want to know what I am doing wrong or are there any other files that also need to be edited ? I am using QGIS3.6 on 32 Bit Windows operating System and using Pycharm as Python IDE

But after saving the changes, when I tried to reload the plugin I have got the below error

Couldn't load plugin 'Nearest_Neighbour' due to an error when calling its classFactory() method

AttributeError: 'NearestNeighbour' object has no attribute 'dlg'
Traceback (most recent call last):
File "C:/OSGeo4W/apps/qgis/./python\qgis\utils.py", line 335, in startPlugin
plugins[packageName] = package.classFactory(iface)
File "C:/Users/Prachi/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\Nearest_Neighbour\__init__.py", line 36, in classFactory
return NearestNeighbour(iface)
File "C:/Users/Prachi/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\Nearest_Neighbour\Nearest_Neighbour.py", line 73, in __init__
self.dlg.lineEdit.clear()
AttributeError: 'NearestNeighbour' object has no attribute 'dlg'


texts: Importing shapefile issue PostGIS, no import made or recognized reference

texts:
I am doing an exercise where I am about to load shapefiles into a PostGIS-database. I am doing it from an already created database and schema, this will mean that the only new item is the table that is being loaded into the schema and database.

I have been following Adding shapefiles to PostGIS database and if I left out everything on the left side of the "|", it would simulate the import but not make the table itself in the database. So I could confirm it worked.

I am guessing that shp2pgsql doesn't evaluate the shapefile table and generate a copy of the table in psql but it merely imports the shapefile. This resonated because I have been using the PostGIS Shapefile Import/Export Manager.

It seems like the only way I can import a shapefile into a PostGIS DB. But it isn't very automated, which is the case of the task I am doing.

Does anyone have any comments on it?

This is my string that works:

C:\Program Files\PostgreSQL\9.5\bin>shp2pgsql -s SRID PATH_TO_FOLDER\SHAPEFILE.shp SCHEME.TABLE

But it needs the table to be created before the import, it doesn't create the table on import. Furthermore. It doesn't import. This string only simulates the import.

I then tried to do as the cheatsheet made an example about. It includes the "|"

But it only gives me this output:

  ERROR: relation "schema.table" does not exist.

Which I read into that fact that shp2pgsql only imports the shapefile if there is a table already.

I then tried to import it with PostGIS Shapefile Import/Export Manager as before mentioned. It imports it (In absolute no time) and creates a table with no more information than otherwise needed. But now the table is there in the database.

So I try again with the table created already in PostGIS, importing the shapefile from the string mentioned before. This time with the credentials left of the "|".

so it looks like this:

 C:\Program Files\PostgreSQL\9.5\bin>shp2pgsql -s SRID Path_To_folder\SHAPEFILE.shp schema.table | psql -h localhost -d postgres -U postgres

And it outputs this:

  ERROR:  current transaction is aborted, commands ignored until end of transaction block ROLLBACK ANALYZE

Is there another way I could do this. Already looked into GDAL's SQL statement but no luck either. Maybe there is a methodology that I could try differently in order to narrow down possibilities of errors occurring?


texts:
EDIT: I have cross-posted this to Drupal Answers. I hope this isn't too much against the rules --- I think the answer will be something to do with both Drupal and OpenLayers.
https://drupal.stackexchange.com/questions/50733/openlayers-toggle-visibility-of-layer

I'm using OpenLayers in Drupal and am stuck while trying to make a layer toggle on and off (as it does from the right-side layer switcher).

See this map.
http://dev.lambeth.coop/map

Click on food growing projects.

EDIT: I've changed the code to show the Drupal-specific elements.

(function ($) {
  Drupal.behaviors.openlayers_switcher = {
    attach: function(context, settings) {
      $(document).ready(function() {
        $('[id^=openlayers-switcher-block--]').click(function(e) {
          $(this).toggleClass('active');
          // 27 is the length of 'openlayers-switcher-block--'
          // layername is the string following that.
          var layername = e.target.id.substring(27);
          var ol = $('.openlayers-map').data('openlayers');
          var layers = ol.openlayers.getLayersByName(layername);

          // this doesn't work
          // if (layers[0].visibility) {
            // layers[0].setVisibility(false);
          // }
          // else {
            // layers[0].setVisibility(true);
          // }
        });
      });
    }
  }
}(jQuery));

When I try to toggle visibility - the visibility doesn't change.

EDIT: So it appears that the entire click function is being called twice.

I've tried putting these lines (in various combinations) at the beginning of the function, as these have been suggested in google searches for 'click events firing twice' etc.

$('[id^=openlayers-switcher-block--]').click(function(e) {
  e.preventDefault();
  e.stopPropagation();

... but these don't help - I still get the function called twice.

In any case, I'm new to openlayers and would be grateful for a way forward!

Thanks!


texts: QGIS 3.18 Slight misalignment of two vector layers

texts:
Per the attached image, my vector layer of political boundaries (in Nepal) is slightly offset to the east, by 300 m. This offset doesn't change from one end of the country to the other, so it doesn't seem to be a projection issue. In reviewing posts on similar topics, they don't seem to address this simple case. There doesn't appear to be any north-south displacement. We want these layers to line up with OSM, which we are using as an underlay.
When we receive a vector layers for, say, administrative boundaries, some are offset 300m to the East, others 100m to the South. So, as Vince suggests, looks like we have a datum transformation issue. Or, rather, several. Now, reading up on the subject, it looks like there are various parameter strings that one enters in the appropriate QGIS dialogue Datum transformation in QGIS 3
My basic understanding is that these parameters tweak a given projection, possibly for a large geographic region, to a particular location.
My larger question is that we want to do this in a robust manner, and it looks like each new set of boundary layers will be slightly offset from the OSM underlay, so we'll need to do this repeatedly.


texts: Highlight all countries a person has visited

texts:
I am interested if it is possible to select a Person in the Attribute table and directly see
which coutries are affected by highlighting them.
Example:
(Here Person 07 is selected and the countries Czech Republic and Estonia should be highlighted in
the Map View.)

I am working with ArcGIS 10.7.
I assume I have to do it via Python/SQL!?


texts: Object does not exist in the namespace "http://schemas.esri.com/arcgis/client/2009"

texts:
I am trying to move all of my GIS code over the Visual Studio 2012. I have added the following references:

C#:

using ESRI.ArcGIS.Client;
using ESRI.ArcGIS.Client.Geometry;
using ESRI.ArcGIS.Client.Tasks;

XAML:

xmlns:esri="http://schemas.esri.com/arcgis/client/2009"

The problem is that many of the objects are not being recognized because they do not exist in the namespace...
Some of these include the InfoWindow and Legend.
However it seems to be OK with objects such as the PictureMarkerSymbol and the Map and Layers.
I don't know if something has changed since I last looked at it, if so I think ESRI need to update their SDK samples. Any help would be appreciated. Thank you in advance.


texts: Date format in GeoServer WMS requests strange

texts:
I got some date values in layer properties.

When i use GetFeatureInfo request to WMS, i got date on day back, for example in DB or in WFS requests or in layer preview it is '2015-04-10', but i got '2015-04-09Z' in my WMS request.

How can i fix or change this output on server side?

Geoserver + PostGis


texts: Splitting DEM based on elevation using ArcMap?

texts:
I am trying to split a DEM TIFF file into three separate shapefiles (upland, middle and lowland) on ArcMap 10.6. I have used natural breaks to classify the DEM into three elevation categories.
How do I use this and split it further into different rasters which I can then create shapefiles from?


texts:
How can I fix an MXD that was created with ArcMap 10.5 and then a copy saved as 10.4?

Background: I worked for 2 hours at school creating MXDs I need for a project. After each MXD was saved using ArcMap 10.5, I saved a copy for use with 10.4 at home. I can open the 10.4 copy at home, but none of the layers appear in the TOC. I do not mean that the link to the data is broken which I expect. I mean it shows no layers at all in the data frame or TOC. Only 1 of the MXDs shows the basemap, but the 40 raster files are not there. I do notice the data frame still has the proper projection. Is there a fix or do I have to repeat 2 hours of work?


texts: Getting subbasin area of each specified river segment？

texts:
I have a DEM file of one district, and I have used those following tools in spatial analyst tools in ArcGIS 10.1 to get files related to the original DEM file：

Fill
Flow direction
Flow accumulation

Besides, I already have the point shapefile that these points are vertices in the river the one belongs to the same district as DEM file, so I use the point shapefile as the pour points in the river.

I use the watershed tool, setting my flow direction file and point shapefile as input files, try to get the subbasin area in this river.

I just got this：

But here's the schematic photo that I want：I want to get every subbasin area and calculate them, assigning them to the attribute of the designated segment

In what way should I use tools in ArcGIS to achieve my goal, thanks for giving me any suggestion!

Although I have a great answer to this question, I got another problem afterwards：Getting additional subbasin areas designated by points in river network？


texts: QGIS batch processing does not allow the output from an algorithm as the input to a subsequent algorithm

texts:
Using QGIS 2.18.14, Windows 10, I need help with Batch Processing.  I have created this simple model with just two algorithms, a reproject followed by polygon centroids:

When I try to run this model as a batch process, I follow the QGIS instructions for filling out the Parameter Table.  The autofill dropdown for the first algorithm works correctly, allowing me to select the input layer.
However, there is a problem with the second algorithm's dropdown; just like the first algorithm it only includes the input layer as an option.  It does not allow me to select the output from first algorithm, which is necessary to complete the model.

How do I enable the parameter table to offer intermediate output as a selection option?  Without this ability, only one algorithm is possible in batch processing.
EDIT: For clarity, my problem is not with the model itself, but rather the batch processor: it does not provide the ability for an algorithm to accept, as input, the output of a previous algorithm.
inability of the batch processor to accept the output of an algorithm


texts: User defined folder location for QGIS3 temp files in Pycharm

texts:
I set up user defined paths for temporary files in the QGIS version that I'm using (settings -> options -> systems). When I run code within QGIS3 the temporary files are saved in the indicated folder destination as expected. When I run PyQGIS code in Pycharm however, temporary files are simply stored in the default local temporary file folder (C:\Users\XXX\AppData\Local). I configured the qgis.bat of the QGIS version that I want to use to launch Pycharm and created the QGIS environment at the beginning of my script.
Is there anything that I'm missing in either the batch file or the QGIS environment setup? I'm using QGIS 3.12 and PyCharm Community Edition 2020.1.1
This is my batch file:
@echo off
call "%~dp0\o4w_env.bat"
call "%OSGEO4W_ROOT%\apps\grass\grass78\etc\env.bat"
call qt5_env.bat
call py3_env.bat
@echo off
path %OSGEO4W_ROOT%\apps\qgis\bin;%OSGEO4W_ROOT%\apps\grass\grass78\lib;%OSGEO4W_ROOT%\apps\grass\grass78\bin;%PATH%
set QGIS_PREFIX_PATH=%OSGEO4W_ROOT:\=/%/apps/qgis
set GDAL_FILENAME_IS_UTF8=YES
rem Set VSI cache to be used as buffer, see #6448
set VSI_CACHE=TRUE
set VSI_CACHE_SIZE=1000000
set QT_PLUGIN_PATH=%OSGEO4W_ROOT%\apps\qgis\qtplugins;%OSGEO4W_ROOT%\apps\qt5\plugins
start "PyCharm aware of Quantum GIS" /B "C:\Program Files\JetBrains\PyCharm Community Edition
2020.1.1\bin\pycharm64.exe"

And this is my QGIS environment setup in Pycharm:
from qgis.core import QgsApplication, QgsProcessingFeedback
from qgis.analysis import QgsNativeAlgorithms

QgsApplication.setPrefixPath(r'C:\QGIS_Installations\QGIS_Aktuell\apps\qgis', True)
qgs = QgsApplication([], False)
qgs.initQgis()

sys.path.append(r'C:\QGIS_Installations\QGIS_Aktuell\apps\qgis\python\plugins')

import processing
from processing.core.Processing import Processing
Processing.initialize()

QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())
feedback = QgsProcessingFeedback()


texts: ArcGIS Pro Python Validate and Alter Field Names

texts:
I want to check if the field names in shapefile start with a number (or other forbidden character) or if they contain spaces. If so I want to update the names into some which are correct.

I tried:

Fields = ['abc', '0-9 abc' ]
list_len = range(len(Fields))
        for i in list_len:
            arcpy.AlterField_management(FromPolygon, Fields[i], arcpy.ValidateFieldName(Fields[i]))

But I get an error message:

  Failed to execute. Parameters are not valid.
  ERROR 000664: Invalid input: The type of dataset is not supported.
  Failed to execute (AlterField).

Fields is a list of strings.


texts: Import OS Mastermap to PostgreSQL database

texts:
I have ordnance survey master map a coverage for the whole of the UK. These are currently received on a disk every 3 months and arrive in .gz format.

Does anyone have any recommendations on how to import this into a PostgreSQL database in bulk so that multiple users can access this?


texts: Can someone give me an example each where these two tools would be appropriate to use. I kind of understood what each of them does but can't seem to figure out examples where one would be more appropriate than the other. Both seem to show density just fine.


texts: Getting Leaflet and image popup from QGIS export working?

texts:
I am a novice in using leaflet, I exported my map layers from QGIS to a working Leaflet map. I am trying to modify the popup to use the attribute so as to get an image popup. i am however stuck as i have read documentation but still unsure. Any leads will help. Currently my popup does a get feature info to show the following.

function pop_dwelling(feature, layer) {
                                        var popupContent = '

house_numb	' + Autolinker.link(String(feature.properties.house_numb)) + '

photo	' + Autolinker.link(String(feature.properties.photo)) + '

<th

';
                    layer.bindPopup(popupContent);

                }

where the photo is a column in the vector that has the path to the image


texts: KeyError when trying to set attribute in QGIS?

texts:
I have two input vector layers and I want to create an Output vector layer by doing some calculations. But somehow, the setAttribute method does not work.
In QGIS 2.18, I was able to set an attribute like this:

ID = 'id'

inEdges = self.parameterAsSource(parameters, self.INPUT_VECTOR_LAYER_EDGES,
                                          context)

# get field indices for inputlayer
fieldsEdges= inEdges.fields()

idxEdgeId = fieldsEdges.indexFromName(ID)

edgeIdTrimmed = some calculation ...
(...)
newEdge = QgsFeature()
newEdge.setAttribute(idxEdgeId, edgeIdTrimmed)   # <-- the error is here

But now, using QGIS 3, I get the following error message:

  Traceback (most recent call last):
  File "", line 125, in processAlgorithm
  KeyError: '0'

Any ideas on how to set an attribute value in QGIS 3?

Is there any other method to do this now?

And what does that KeyError mean exactly?


texts: Hillshade for 30m NED at continental scale

texts:
I'm new to terrain analysis. I have a 30m NED DEM (.tif) for the conterminous US. I would like to generate a nice looking hillshade for the US, but it looks bad - as in pixelated and rough (See image). I have tried different classification schemes and number of classes, which helps some, but not enough. I have tried majority filter, but the output looks exactly the same. Is there some sort of smoothing function that should be applied for this scale? Any suggestions much appreciated.

Edit: Using ArcMap 10.1, projected w/ NAD83 Albers
Edit 2: The hillshade could use more "depth", but looks good for now (image 2).


texts: Interpolating points along line in Python for ArcGIS Desktop?

texts:
I am building a Python geoprocessing script for ArcGIS 9.3. In the script, I have a simple LINE feature class with one row, one line.

How can I use a linear referencing system to interpolate a point along the path?

I can use a cursor to access the feature geometry, and get the vertex geometries. But I can't find anything that could help me interpolate points along it.

Note: this is ridiculously simple with Shapely:

import shapely.wkt
line = shapely.wkt.loads('LINESTRING (380 60, 360 130, 230 350, 140 410)')
line.interpolate(0.0).xy # (array('d', [380.0]), array('d', [60.0]))
line.interpolate(100.0).xy # (array('d', [346.16312174913622]), array('d', [153.41625550146173]))
# ... etc

Is there any equivalent in ArcGIS?

I have all the common ArcGIS extensions to my disposal.

Or should I just bring the geometry over to Shapely to do the work?

The geometry processing needs to eventually go back to ArcGIS.


texts: How can I measure linestring on postgis?

texts:
I have a linestring and I want to measure this line.
I try it on google maps. the result is that 0,65 km
https://www.google.com/maps/dir/40.144048,26.45235/40.138443+26.450246/@40.1412448,26.4467258,16z/data=!3m1!4b1!4m6!4m5!1m0!1m3!2m2!1d26.450246!2d40.138443

Also I try it bing maps. the result isthat 0,60 km.
http://www.bing.com/maps/?FORM=Z9LH3#Y3A9NDAuMTQxMjM3fjI2LjQ1MTMyMSZsdmw9MTYmc3R5PXImcnRwPXBvcy40MC4xNDQwNDhfMjYuNDUyMzUwX0QtNTUwJTIwJTI2JTIwJUMzJTg3YW5ha2thbGUlMjBLYXYlQzUlOUZhJUM0JTlGJUM0JUIxJTJDJTIwMTcwMDAlMjAlQzMlODdhbmFra2FsZSUyQyUyMFQlQzMlQkNya2l5ZSUyMHlhayVDNCVCMW5sYXIlQzQlQjFuZGFfX19hX35wb3MuNDAuMTM4NDQzXzI2LjQ1MDI0Nl9ELTU1MCUyQyUyMDE3MDAwJTIwJUMzJTg3YW5ha2thbGUlMkMlMjBUJUMzJUJDcmtpeWUlMjB5YWslQzQlQjFubGFyJUM0JUIxbmRhX19fYV8mbW9kZT1EJnJ0b3A9MH4wfjB+

and I measure it with postgis sql;

SELECT ST_Length(ST_GeomFromText('LINESTRING(40.144048 26.45235,40.138443 26.450246)'));

the result is :0.00598688909200328  km

SELECT ST_Length_Spheroid('LINESTRING(40.144048 26.45235,40.138443 26.450246)',
                          CAST('SPHEROID["WGS 84",6378137,298.257223563]' As spheroid) );

the result is :605.657141879302 km

I wonder that which one is correct lenght? please say some advice for me ?


texts: Looping through CSV files to write Lists

texts:
I am very new to Python and am attempting to use it for some relative humidity data. I am looking at 60 months of data. The files include a Lat and Long file (CSV in a Matrix) and Relative Humidity data (60 CSV Matrix files). What I'm trying to do is correspond the Lat and Long to the value in the Relative Humidity files by turning them into 60 CSV XYZ (latitude, longitude. RH) files - essentially merging the data and converting from matrix to list form. My problem is I don't know how to loop through the RH files.

# Create XYZ CSV file from Latitude, Longitude, and relative_humidity files
import csv

path = "/Users/msayre1/Desktop/CSV_Input/"
longitude = []
latitude = []
humidity = []

with open(path + 'LATITUDE.csv', 'rb') as f:
    reader = csv.reader(f)
    latitude = list(reader)

with open(path + 'LONGITUDE.csv', 'rb') as f:
    reader = csv.reader(f)
    longitude = list(reader)

with open(path + 'relative_humidityApr09.csv', 'rb') as f:
    reader = csv.reader(f)
    humidity = list(reader)

with open(path + 'RH_200904.csv', 'wb') as newFile:
    newWriter = csv.writer(newFile, delimiter = ',')
    for row, value in enumerate(latitude):
        for col, lati in enumerate(value):
            print lati + "," + longitude[row][col] + "," + humidity[row][col]
            newWriter.writerow([lati, longitude[row][col], humidity[row][col]])


texts: Installing scikit-learn python package for QGIS on Win

texts:
I want to install scikit-learn package to be executed in Python-QGIS. I've tried to use OSGeo4W Shell to install with pip. But I got this error.
I've tried to run osgeo4w-setup-x86_64.exe and install all Python packages but it doesn't work.


texts: K-Means clustering in QGIS (Attribute based clustering)

texts:
In the 2.18.4 and 2.18.5 version of QGIS K-means clustering was moved to a plugin called 'Attribute based clustering'. However, I am working with a virtual raster dataset (vrt.) and the current version of he plugin works with vectors.

Is there another method i can use to  perform an unsupervised K-means clustering in QGIS?


texts: Using a parametric SQL View in Geoserver

texts:
In Geoserver i create parametric view:

SELECT *
FROM MyTable
LEFT JOIN Track On MyTable.SegmentID=Track.SegmentID
WHERE (DATEPART(yy, MyTable.Date)=%year%
AND DATEPART(mm, MyTable.Date)=%month%
AND DATEPART(dd, MyTable.Date)=%day%)

i set default values 2009-09-18

Now i go in layer preview pick my layer and see layer with features with Date=2009-09-18.
Now i want to see features with another date. In address line i add:

viewparams=year:2009

But still see same features. I try another dates but not help.
Whats wrong?


texts: Import waypoints from Garmin GPSmap into R

texts:
I would like to import waypoints from my Garmin GPSmap 62 into R. I have about 270 waypoints saved on the GPS unit. I can use Garmin's BaseCamp program to export them, but I would like to do this with R.

Once I plug the Garmin unit into my computer (Mac OS X) via USB, the file system I see is the following:

~ $ cd /Volumes/GARMIN/
/Volumes/GARMIN $ ls
Documents   Garmin      autorun.inf
/Volumes/GARMIN $ cd Garmin/
/Volumes/GARMIN/Garmin $ ls
BirdsEye             Garmintriangletm.ico gmapbmap.sum
CustomMaps           JPEG                 gmaptz.img
CustomSymbols        Profiles             startup.txt
ExtData              SQL                  system.xml
Filters              Text                 wytopo11.img
GPX                  fs_image.ver
GarminDevice.xml     gmapbmap.img
/Volumes/GARMIN/Garmin $ cd GPX
/Volumes/GARMIN/Garmin/GPX $ ls
Archive                 Waypoints_12-MAY-13.gpx Waypoints_19-MAR-16.gpx
Current                 Waypoints_13-APR-15.gpx Waypoints_20-OCT-16.gpx
Nav                     Waypoints_13-JUN-13.gpx Waypoints_21-MAR-16.gpx
Waypoints_01-MAY-14.gpx Waypoints_13-MAY-13.gpx Waypoints_21-MAY-14.gpx
Waypoints_02-APR-15.gpx Waypoints_13-MAY-14.gpx Waypoints_21-OCT-16.gpx
Waypoints_02-MAY-15.gpx Waypoints_14-JUN-13.gpx Waypoints_22-MAR-16.gpx
Waypoints_03-APR-15.gpx Waypoints_14-MAR-16.gpx Waypoints_22-MAY-14.gpx
Waypoints_03-JUN-13.gpx Waypoints_14-MAY-13.gpx Waypoints_22-OCT-16.gpx
Waypoints_03-MAY-15.gpx Waypoints_14-MAY-14.gpx Waypoints_23-MAR-16.gpx
Waypoints_04-MAY-14.gpx Waypoints_15-APR-15.gpx Waypoints_23-MAY-14.gpx
Waypoints_05-APR-15.gpx Waypoints_15-JUN-13.gpx Waypoints_24-MAR-16.gpx
Waypoints_05-JUN-13.gpx Waypoints_15-MAR-16.gpx Waypoints_24-OCT-16.gpx
Waypoints_05-MAY-14.gpx Waypoints_15-MAY-14.gpx Waypoints_25-APR-15.gpx
Waypoints_06-APR-15.gpx Waypoints_16-MAR-16.gpx Waypoints_25-OCT-16.gpx
Waypoints_06-MAY-14.gpx Waypoints_16-MAY-13.gpx Waypoints_26-APR-15.gpx
Waypoints_07-JUN-13.gpx Waypoints_16-MAY-14.gpx Waypoints_26-MAR-16.gpx
Waypoints_08-APR-15.gpx Waypoints_17-JUN-13.gpx Waypoints_26-MAY-13.gpx
Waypoints_08-JUN-13.gpx Waypoints_17-MAR-16.gpx Waypoints_26-OCT-16.gpx
Waypoints_10-APR-15.gpx Waypoints_18-JUN-13.gpx Waypoints_27-APR-15.gpx
Waypoints_10-JUN-13.gpx Waypoints_18-MAR-16.gpx Waypoints_28-APR-14.gpx
Waypoints_10-MAR-16.gpx Waypoints_18-MAY-14.gpx Waypoints_29-APR-14.gpx
Waypoints_11-JUN-13.gpx Waypoints_18-OCT-16.gpx Waypoints_29-APR-15.gpx
Waypoints_11-MAY-13.gpx Waypoints_19-APR-15.gpx Waypoints_30-APR-14.gpx
Waypoints_12-MAR-16.gpx Waypoints_19-JUN-13.gpx Waypoints_30-APR-15.gpx

Within the GPX directory, there is a Current directory that contains Current.gpx.

I'm not sure if that is the file that contains my waypoints, but this is what I started with.

I tried to load this into R with rgdal:

> readOGR('/Volumes/GARMIN/Garmin/GPX/Current/Current.gpx', layer="track_points")
OGR data source with driver: GPX
Source: "/Volumes/GARMIN/Garmin/GPX/Current/Current.gpx", layer: "track_points"
with 9759 features
It has 26 fields

There are way more entries than waypoints that I expect to have. Also, none of the columns in the data slot of this SpatialPolygonsDataFrame have the names of the waypoints.

So, am I looking at the wrong file? Is there a better way to load the waypoints into R?


texts: Georeferencing PNG against Google maps with OpenLayers plugin

texts:
I am trying to georeference blueprints to the footprint of a building. I'm using the Google Satellite image from the OpenLayers plugin to do this. The blueprints are rotated about 30 degrees. When I run the georeferencer after placing GCP, the blueprints load in vertically, but the rotation doesn't occur...

Why is this?

I figured it would line up the GCPs properly...


texts: Using a PostGIS raster catalog in GeoServer?

texts:
I have a PostGIS table cataloging 10,000+ GeoTIFFs. Each table row has a name, file path, and a polygon as the bounding box. The rasters are spread out all over the globe, they are not adjacent.

I would like to be able to search them, and add them to my map one at a time.

I thought about using imageMosaic, and including the name attribute in "AdditionalDomainAttributes" in my indexer.properties file like I saw in this question, but I don't need them as a mosaic as I never need to serve more than one or two at a time, and imageMosaic doesn't seem allow rasters to be added that  not adjacent or overlapping to existing ones. (Also this database is added to frequently.)

Is there a different plugin I should use, or a way to handle this effectively? The best I have come up with so far is just using image mosaic, and creating a global "background" raster so that they are all overlapping, and then just only allowing my client to request one "name" at a time when requesting tiles for that layer.

I'm using GeoServer 2.15.0