texts: Removing existing backup layers and renaming current live layer in ArcSDE

texts:
Looking at the attached screenshot of Catalog, there are the file "CRT_CANAL" and many backups. The aim of this code is deleting old backups in the ArcSDE database,
renaming the existing "CRT_CANAL" as backup + the current date (es CRT_CANAL_bck20151222) and copying the new layer in the database that will get the original name "CRT_CANAL" without date.

Running my code I got a generic error in the deleting layers section. So that the machine overcomes that part of the code and goes to the next stage renaming current live layer.

This is the message I get:

"delete old backup layer ERROR: "

without any other specification of the kind of error

This is the code:

                if featCount > 0:
                    sde_Connection = "C:\Users\<user>\AppData\Roaming\ESRI\Desktop10.1\ArcCatalog\GISAdmin DC to GISTest.sde\GISTest.GISADMIN."
                    layer_sde_output = sde_Connection + wfsOutputName
                    layer_bck_delete = sde_Connection + wfsOutputName + "_bck*"
                    layer_bck = sde_Connection + wfsOutputName + "_bck" + strftime("%Y%m%d")
                    data_type = "FeatureClass"

                    # Check if backup layer exists in SDE and remove it
                    if arcpy.ValidateTableName(layer_bck_delete):
                        try:
                            fcList = arcpy.ListFeatureClasses(layer_bck_delete)
                            for fc in fcList:
                                data_type = "FeatureClass"
                                arcpy.Delete_management(fc, data_type)
                                if debugMode:
                                    msg = msg + "delete old backup layer success" + "\n"
                        except:
                            msg = msg + "delete old backup layer ERROR: " +  arcpy.GetMessages(2) + "\n"
                            process = 0

                    if arcpy.Exists(layer_sde_output):
                        if process == 1:
                            # Rename current live layer
                            try:
                                arcpy.Rename_management(layer_sde_output, layer_bck, data_type)
                                if debugMode:
                                    msg = msg + "create backup layer success" + "\n"
                            except:
                                msg = msg + "create backup layer ERROR: " +  arcpy.GetMessages(2) + "\n"
                                process = 0

                    if process == 1:
                        # Import into SDE
                        try:
                            sdeLayer = sdeConnection + sdeOwner + wfsOutputName
                            arcpy.CopyFeatures_management(gdbLayer, sdeLayer)
                            if debugMode:
                                msg = msg + "SDE import success" + "\n"
                        except:
                            msg = msg + "SDE import ERROR: " +  arcpy.GetMessages(2) + "\n"
                            process = 0

Anyone can help me in correcting this code?


texts:
To summarize, trying to read a raster keeps giving me an UnsatisfiedLinkError that makes sense, but looking into the compiled GDAL jar, I don't know how to fix.

gdal.AllRegister();
TheFile = gdal.Open(filepath);

double[] attributes = TheFile.GetGeoTransform();
extent = new int[] {TheFile.GetRasterXSize(), TheFile.GetRasterYSize()};
pixel_size = new double[] {attributes[1], attributes[5]};
corner = new double[] {attributes[0], attributes[3]};

Band TheBand = TheFile.GetRasterBand(1);
TheData = new double[extent[0]][extent[1]];
for(int y = 0; y < extent[1]; y++) {
    TheBand.ReadRaster(0, y, extent[0], 1, TheData[y]);
}

So the preceding code should just open and load the raster (specifically a tiff), line by line, into an array. This however fails at ReadRaster() with the following error message:

java.lang.UnsatisfiedLinkError: org.gdal.gdalconst.gdalconstJNI.DMD_LONGNAME_get()Ljava/lang/String;
    at org.gdal.gdalconst.gdalconstJNI.DMD_LONGNAME_get(Native Method)
    at org.gdal.gdalconst.gdalconstConstants.<clinit>(gdalconstConstants.java:75)

Now I notice that the problem is the method called DMD_LONGNAME_get() can't be found. Looking into the gdalconstJNI.java in the GDAL jar file I see there is a similar but slightly differently named method called GDAL_DMD_LONGNAME_get(). So I understand the unsatisfied link error, but why is it named differently and how do I fix it?

I used the GDAL wrapper from here, in particular the latest MSVC2010 Win64 stable release (as of now: release-1600-x64-gdal-1-10) and followed all the proper installation instructions, updating my 2010 Visual C++ Redistributable x64 (and x86 for good measure), running the bat file, creating all the environment paths, adding the -Djava.library.path to the project's VM options, making sure everything's compatible with 64-bit Win7, so I don't think there's a problem there.

Is this just a GDAL badly written to Java issue? (I mean it seems that way with the mismatched naming issue). Is there a previous release version I should go with? If so, which version? as I don't want to have to try it one by one.

UPDATE Just tried with the oldest GDAL version in the previously mentioned site with compiled GDAL wrappers (gdal-1-7-3), in which the DMD_LONGNAME_get() exists. However now I get the same error but instead the method missing is now GARIO_PENDING_get()I which I don't find any version of it in gdalconstJNI.java. Unless I get a response I suppose I'll just have to keep trying different versions and see if there is a sweet spot, but this seems clumsy/time-consuming so hopefully someone has run into this before.


texts: QGIS version 2.18.14 exporting georeferenced tiff file to kml

texts:
It seems everybody can use right-click on layers, and "save as" to export to .kml. I only get the opportunity to save as "GTiff".

What am I missing out here?


texts: Changing Spatial Reference System (SRS) of existing GeoTiff file?

texts:
How do I change the Spatial Reference System (SRS) of an existing GeoTiff file?

Or more fully, what can I use to convert the meta data for the file to a different SRS?


texts:
I tried to do calculations for 2 dataset of EVI on different year. I got an error message evi2020.subtract is not a function.
The code is below
    var countries = Indonesia

var Nama_Provinsi = ['Kalimantan Selatan']
var geometry = countries.filter(ee.Filter.inList('NAME_1', Nama_Provinsi));
var area1 = ee.FeatureCollection(geometry);

var dataset = ee.ImageCollection('MODIS/006/MOD13Q1')
                  .filter(ee.Filter.date('2020-01-01', '2020-12-31'));
var evi2020 = dataset.select('EVI');
var eviVis = {
  min: 0.0,
  max: 8000.0,
  palette: [
    'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
    '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
    '012E01', '011D01', '011301'
  ],
};
Map.setCenter(116, -3, 8);

Map.addLayer(evi2020, eviVis, 'EVI20');

var dataset = ee.ImageCollection('MODIS/006/MOD13Q1')
                  .filter(ee.Filter.date('2010-01-01', '2010-12-31'));
var evi2010 = dataset.select('EVI');
var eviVis = {
  min: 0.0,
  max: 8000.0,
  palette: [
    'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
    '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
    '012E01', '011D01', '011301'
  ],
};
Map.setCenter(116, -3, 8);

Map.addLayer(evi2010, eviVis, 'EVI10');

var selisihevi = evi2020.subtract(evi2010);
Map.addLayer(selisihevi.clip(area1), {min: -0.5, max: 0.5, palette:'blue, green, yellow, red'}, "Selisih EVI");


texts: How do I change the style of all markers in OpenLayers3?

texts:
This is probably a newbie question, but am at my wits end.

I have a map that has several markers (see attached code).  What I am trying to do is alter ALL the markers depending on a map event (e.g., resizing marker for a zoom event). I am trying to do this without embedding the entire for loop creating markers within the event handler.  Any idea how I can change all markers together?

JS Fiddle

<!DOCTYPE html>

<link rel="stylesheet"
href="//cdnjs.cloudflare.com/ajax/libs/ol3/3.5.0/ol.min.css" type="text/css">


texts: How to replicate ArcGIS geometric intervals in QGIS

texts:
Is it possible to replicate ArcGIS's Geometric Intervals in QGIS 3 layers symbology?
I am using QGIS 3.10.9 LTS and would like to assign a geometric interval style to a simple shapefile of UK counties which contains a numeric field. In QGIS-layers-symbology one can select 'graduated' as the style and select the numeric field as 'value', but none of the six modes QGIS provides is appropriate to my analysis: What I need is a Geometric Interval classification like the one available in ArcGIS. I cannot see a Plugin that attempts to address this and the one online article I can find is from 2014 and uses an experimental Python script on a now unavailable website. Is there a way of solving this task? My knowledge of Python is very limited.


texts: How to add up election zone voting upto county level?

texts:
Sorry if it is a trivial question.

I am trying to add voting data from election zone up to county-level using spatial join tool, but I cannot figure out which relationship I should use. When I use Intersect, it seems data are being repeated (original sum:13198656 After spatial join:13921985). I also tried Contains, it seems data are being left then (original sum: After spatial join:12262146). I guess some election zones are intersecting between two counties. What should I do about it?


texts: WriteArray method not updating raster pixel value

texts:
I have created a raster.tiff file using GDAL driver with data type as Byte and have also added a palette of unique values to it (0,1,2).

When i try to write a pixel value for a particular Raster Band the value is not getting updated.

I have tried using FlushCache() and del variables, but nothing seems to work.

rasterDS=gdal.Open("C:/Anuj/fyp/Delhi_Sample/Clutter/raster.tiff")
rasterBand=rasterd.GetRasterBand(1)
rasterArray=rasterb.ReadAsArray()
rasterArray[0][1]
1
rasterArray[0][1]=2
rasterBand.WriteArray(rasterArray)
0
rasterBand.FlushCache()
rasterDS.FlushCache()
rasterBand=None
rasterDS=None
del rasterBand
del rasterDS
rasterDS=gdal.Open("C:/Anuj/fyp/Delhi_Sample/Clutter/raster.tiff")
rasterBand=rasterd.GetRasterBand(1)
rasterArray=rasterb.ReadAsArray()
rasterArray[0][1]
1

The pixel value is not getting updated after writing the array to the raster band.


texts: Connecting Problem: ODBC versus ESRI Personal Geodatabase (*.mdb)

texts:
I try over 2 days to connect an ESRI Personal Geodatabase with QGIS 2.8.1

Here is my testing report.

Testing environment: Windows 7 64-bit; MS Office 32-bit; QGIS 32-Bit

1. Step: Creating an ODBC-System-DSN
We have a 32-bit Office installation. We find the 32-bit ODBC-Administration-Tool under:

C:\Windows\SysWOW64\odbcad32.exe

The mdb-file is an ESRI Personal Geodatabase.

2. Step: Testing the ODBC-Connecting work
With the database4.exe, which is free, I tested the new ODBC-Connecting does work.
(http://fishcodelib.com/Database.htm)

Yes, very good!

And here you can see all tables in the mdb-File.

3. Step: Open with ODBC-Connection in QGIS

Yes, very good. The ODBC is working!

Oh no, the geometry type is Unknown!!

4. Step: Open with ESRI-Personal-GeoDatabase-Connection in QGIS

Oh no, the ESRI-Personal-GeoDatabase-Connection doesn't work! Why?

Reference 1:
Should I try the QGIS 64-bit version? Like this post: https://gis.stackexchange.com/a/139716/53879

My opinion is, that in this post with drag the mdb file into your workspace is NOT an ODBC-Connection. The ESRI Personal Geodatabase will add with "Add Vector-Layer.." in File-Modus. Just try, drag a mdb-file into the workspace. It will open...
But I want to connect the ESRI Personal Geodatabase in Database-Modus.

Reference 2:
In this post: http://www.northrivergeographic.com/qgis-accessing-personal-geodatabase
the workaround work. I think here is QGIS 64-bit installed.

QUESTIONS:
-Why work the ODBC-Test and the ESRI-Personal-GeoDatabase-Test not?

-Should I try my workaround in QGIS 64-bit installation? Have anybody some experience?

PROBLEM:
I want to connect an existing ESRI-Personal-GeoDatabase (*.mdb) over ODBC and edit the data. In 32-bit installation of QGIS is this not possible??

INFORMATION:
I tryed over 4 days to connect an ESRI-Personal-GeoDatabase (*.mdb) over the Database-Radio-Button.
UNSUCCESSFUL!
I tested extra on a 64Bit PC: 64bit Windows 7 with SP1; 64Bit QGIS; 64Bit Office 2010.
I received the same error.

And it was an Original ESRI-Personal-GeoDatabase (*.mdb). I opened with ESRI-Software.


texts:
I am a new user of Python on ArcGis10 i am trying to apply symbology settings to multiple files but to no success. I have tried my own ways copied syntaxes from various places still getting error. The following is the one that make sense to me but What am i doing wrong??

import arcpy
... from arcpy import env
... env.workspace = r"C:\Experiments"
... inputLayer = ["b","c","d","e","f","g","h","i","j","k","l"]
...
... symbologyLayer = "a"
... arcpy.ApplySymbologyFromLayer_management (inputLayer, symbologyLayer)
...
***Runtime error <type 'exceptions.RuntimeError'>: Object: Error in executing tool***


texts: How to find intersection polygons in single layer on QGIS version 3.4.2?

texts: I have polygons intersection in single layer. I want polygon in yellow frame and use polygon intersection to calculate area , how to do this?


texts: How to add marker in Leaflet base on Google Maps layer?

texts:
I have tried many way to add marker base on Google Maps but it doesn't work at all, but if I changed my base map to such as open street map which provide by Ersi or Leaflet, it does work.

Why is it that the only Google Maps layer doesn't work?

Here is my code:

<!DOCTYPE html>

    <title>Leaflet</title>
    <link rel="stylesheet" href="https://unpkg.com/leaflet@1.0.3/dist/leaflet.css" />

    <link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-0.3.1/leaflet.css" />


texts: OpenLayers getfeatureinfo shows all features, even with CQL filter

texts:
I have a map with several CQL filters that a user can use in my website. When the user applies the CQL filter, I want them to be able to use getfeatureinfo which only shows features from the CQL filter. However, after applying a CQL-filter, it still gives featureinfo of all features, even the ones which are not visible. I looked for solutions here: OpenLayers CQL filter gets info about NOT VISIBLE features and here: How to add CQL filter into GetFeatureInfo request in OpenLayers, but I couldn't solve it myself.

My code looks something like this:

Ext.onReady(function() {

var control_zoom_max = new OpenLayers.Control.ZoomToMaxExtent({title:'Zoom naar het volledige kaartbeeld'});
    var control_zoom_in  = new OpenLayers.Control.ZoomIn({title:'Zoom gecentreerd in op het kaartbeeld'});
    var control_zoom_out = new OpenLayers.Control.ZoomOut({title:'Zoom gecentreerd uit op het kaartbeeld'});
    var control_zoom_box = new OpenLayers.Control.ZoomBox({title:'Zoom in met een te slepen rechthoek'});
    var control_nav_hist = new OpenLayers.Control.NavigationHistory();
    var control_drag_pan = new OpenLayers.Control.DragPan({title:'Verplaats het kaartbeeld met de linkermuisknop'});
    var control_identify = new OpenLayers.Control.WMSGetFeatureInfo({
        title : 'Klik op een punt of lijn voor informatie',
        layers : [ AmiceHerkomst, AmiceVerwerking],
        maxFeatures: 1000,
        queryVisible : true,
        infoFormat : 'application/vnd.ogc.gml',
        vendorParams: {
            buffer: 5  // geoserver buffer in pixels
        },
        eventListeners : {
            getfeatureinfo : function(event){
                showInfo(event);
            }
        }
    });
});

function updateFilter() {
    var plaatshtml = document.getElementById('filter').value.toLowerCase();
    var volumehtml = document.getElementById('volume').value;
    var volumeType = document.getElementById('filterType').value;
    var jaar = document.getElementById('Jaar').value;
    var maand = document.getElementById('Maand').value;
    var bedrijfhtml = document.getElementById('filternaam').value.toLowerCase();

    var plaatsnaam = "strToLowerCase(plaats) LIKE '%"+ plaatshtml +"%'";
    var volume = "totgw_1 "+ volumeType + " " + volumehtml;
    var datum = jaar + maand;
    var bedrijfsnaam = "strToLowerCase(straatnm) LIKE '%"+ bedrijfhtml +"%'";

    var filterParams = {
                    cql_filter: null
                };

    if (OpenLayers.String.trim(plaatshtml) != "") {
                    filterParams["cql_filter"] = plaatsnaam;
                }
    if (OpenLayers.String.trim(volumehtml) != ""){
                    filterParams["cql_filter"] += "AND " + volume;
                }
    if (OpenLayers.String.trim(bedrijfhtml) != ""){
                    if (filterParams ["cql_filter"] != null){
                        filterParams["cql_filter"] += "AND " + bedrijfsnaam;
                    }
                    else {
                        filterParams["cql_filter"] = bedrijfsnaam;
                    }
                }
    if (OpenLayers.String.trim(jaar) != ""){
                    filterParams["cql_filter"] += "AND " + datum;
                }

    mergeNewParams(filterParams);
}

How do I add my CQL-filter parameters into my getfeatureinfo request?

EDIT
this is what the getfeatureinfo request does looking via firebug:


texts: MakeFeatureLayer_management parsing error (for loop)

texts:
I want to make a Feature Layer for each month in the table "Sightings1995-2014".

Why doesn't it work?

      for month in range (01,12):                                                            arcpy.MakeFeatureLayer_management("Sightings1995-2014", "test%s" %month,"DATE >= date'01.%s.1995' AND DATE < date'01.%s.2014'" % (month, month+1)

Parsing error SyntaxError: invalid syntax (line 1)

for the output in each year this code worked

      for year in range(1995,2014):
      arcpy.MakeFeatureLayer_management
      ("Sightings1995-2014", "test%s" % year, "DATE >= date'01.01.%s' AND DATE < date'01.01.%s'" % (year, year+1))


texts:
I recently downloaded some soils data using a WMS server into ArcMap.  Here is a link to these data http://maps.isric.org/. The data gets imported into ArcGIS as a .lyr file composed of about 20 layers (for example the clay dataset in the link) . I want to convert and save each of these layer files as a raster file.

How do you do this using ArcGIS Desktop?

I eventually want to resample a few of these layer files and save them as .asc (ASCII) files.


texts: Shapefiles do not store times so how do I create time intervals for use in Torque?

texts:
I created a bunch of time intervals in a file geodatabase with the intention of using them in Torque to create an animation based on those time intervals. However I can't upload a .gdb to my visualisation in CartoDB right? So I exported my .gdb features to shapefiles but, ahhh, shapefiles do not store times! I tried converting my time intervals to integers thinking I was clever but then my animation has annoying pauses when one hour changes to the next such as from 10.59 to 11 and since that is the actual time that I would want to show on the time slider I really can't change it.

For now my animation Time Column is set to the UID which is in sequential order according to the times so it looks alright but its not what I was hoping for. There is no way to see the progression of time on my map. Lastly I am using the Torque CAT option because I have several events happening that start at different times albeit on the same day and sometimes these events overlap as well - so you see I really need to be able to use my time intervals to achieve the desired effect.


texts: How can I update geojson properties using popup leaflet?

texts:
I have GeoJSON file named json_point2. This file holds point data as follows:

var json_point2={
"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
"features": [
{ "type": "Feature", "properties": { "Class": null }, "geometry": { "type": "Point", "coordinates": [ 40.193168611061743, 40.904508523555492 ] } },
{ "type": "Feature", "properties": { "Class": null }, "geometry": { "type": "Point", "coordinates": [ 40.191358449981372, 40.902051540140285 ] } },

When I click this point on my web apps, I want to enter input value on popup and save this value geojson Class properties.

I can click on the point with the following code and enter the attribute. But he does not record it in his database.

My code:

map.createPane('pane_point2');
map.getPane('pane_point2').style.zIndex = 602;
map.getPane('pane_point2').style['mix-blend-mode'] = 'normal';
var layer_point2 = new L.geoJson(json_point2, {
    attribution: '',
    pane: 'pane_point2',
    onEachFeature: function (feature, layer) {
        var input = L.DomUtil.create('input', 'my-input');
        input.value = feature.properties.Class;
        L.DomEvent.addListener(input, 'change', function () {
            feature.properties.Class = input.value;
        });
        layer.bindPopup(input);

    }
}).addTo(map);


texts: Selecting pixels in raster using polygon selection tool in QGIS?

texts:
I want to reclassify landcover image manually with respect to online imagery. I am looking for a software/plugin/script to manually select raster pixels using polygon selection tool and could be able to enter new value.
Can you Please explain the process in QGIS?


texts: Replace the empty cells of one raster with the values of a second raster

texts:
I have 2 rasters with the same extent. In the first raster, there are a lot of empty cells that I want to replace by the values of the second raster.
I have tried Raster calculation with the following expression...

Con(IsNull("DTM_Clip_Clip1"),"DSM_Clip","DTM_Clip_Clip1")

I don't get an error message, but no output is created?

I also tried Mosaic to new raster, but also here there is no output and no error message?

I use ArcGIS 10.6 Desktop.


texts:
I have a KML-file with a path in Google Earth that is slightly displaced.  And I don't want to move each coordinate manually.

How can I move all of the coordinates together?

I think this site is using the same kml as I do. It depicts a city wall that is off. It should be moved like 100meters to the south so it don't crosses the major road that is north of it (as it does now). Or if you want to be more precise, each nod is of by longitude 0,000182 and lattitude 0,001217 if my calculation is right.


texts:
Is there a simple way to create a "CAD template" for map outputs from ArcGIS as in the image below:

I'd like to create a standard template and then populate some simple 'fields' with map information (title, scale, date).


texts: Plotting proportional symbol map of horse mackerel species using R

texts:
library(sf)
library(rgdal)
library(cartography)
hm_SHP <- st_read(file.choose())

Reading layer Oman -coastline' from data source /Users/salmanmaf/Documents/OMAN map_shapefiles/OTrIs/Oman
-coastline.shp' using driver `ESRI Shapefile' Simple feature collection with 2279 features and 2 fields geometry type:  POLYGON
dimension:      XY bbox:           xmin: 46.2 ymin: 14.5 xmax: 65.6
ymax: 29.1 CRS:            NA

hm_sp <- read.csv(file.choose())

plot(st_geometry(hm_SHP))

plot(hm_SHP)

# Total length plot on proportional symbols
propSymbolsLayer(x = hm_SHP, var = "TL",
                  symbols = "circle", col =  "white",
                  legend.pos = "right", border = "grey",
                  legend.title.txt = "TL",
                  legend.style = "c")

Error in abs(x[[var]]) : non-numeric argument to mathematical function

I got the above error but what is the problem?


texts: Eliminating thin polygon connection in QGIS

texts:
I'm facing a problem concerning the fact that 2 polygons are connected by a thin connection.

Is there an automatic method/process to deal with this situation, separating the 2 polygons and maintaining the same information for both?
I've already tried the "Correct Geometries" process.
I know that on QGIS 3.0 or later, there is the verify geometries toolbox, is there any specific configurations inside there that is able to deal with this?


texts: Which GeoTIFF tags tell me what I want to know?

texts:
I'm reading some GeoTIFF files into my app using libgeotiff. After the model of the listgeo program that comes with libgeotiff, I want to extract the extents (corners) of the GeoTIFF so I can position it.

I only really want to deal with coords in WGS84 and UTM (ultimately, I'll convert WGS84 to UTM). However, I don't know which bit of the GeoTIFF info tells me what I want to know about the 'corner' coords. Obviously, I can use the magnitude of the numbers to guess if the resulting coords are lats and longs, but surely there must be an official way?

Something like (pseudo-code):

If (geotifftag == UTM || geotifftag == WGS84)
    // Proceed
Else
    Msg("Unrecognised data format")
End If

Could anybody help me out with what the GeoTIFF tags actually tell me?


texts: Preserving polygon boundaries when vectorizing raster according to polygon vector?

texts:
I have a raster which I clip by some polygons and then convert to polygon. I would like to obtain the same number of polygons in both my vector files.

The problem is that some areas in the vector clipper are separated by a line (i.e. adjacent polygons) are now together (as I would expect). When I perform the vectorization of my cropped raster, I get less polygons than what I want (for the reason I just mentioned).

Is there a way to preserve polygon boundaries when performing raster to vector transformation according to a given polygon vector (that is, my clipper file)?

I use ArcGIS, QGIS, GRASSgis, python, MATLAB/octave.


texts: Raster ST_UNION Losing Cells with Zero as Value

texts:
I am encountering a bit of frustration in trying to get 3 rasters to merge through the use of ST_UNION(rast, 'MAX'). It appears as though it is removing anything with a value of zero from the final result, which is behavior that breaks things downstream.
SQL to check before:
SELECT ST_SummaryStats(rast) FROM table WHERE blah='foo'

Results:
"(1048576,232,0.00022125244140625,0.014872911240353,0,1)"
"(1048576,26716,0.0254783630371094,0.157572891241033,0,1)"
"(1048576,48659,0.0464048385620117,0.210360237497599,0,1)"

Union Attempt:
SELECT ST_SummaryStats(ST_UNION(rast, 'MAX')) FROM table WHERE blah='foo'

Results:
"(75607,75607,1,0,1,1)"

So in this union it looks like it was only keeping items with a value of 1, which is not the desired behavior for what is effectively being used as a Boolean layer.


texts: Get FeatureDataSet of a invalid FeatureClass

texts:
I have a project that repairs broken datasources. One goal of the project is to add the datasource of a featuredataset to repair any broken feature classes. However, so far the only way I've seen to get the FeatureDataSet is through the FeatureClass. However, since these layers are invalid, there is no featureclass.

Normally the parent featuredataset of a featureclass is held in the FeatureClass. So is it accessible anywhere else?

I have already checked the properties in the IProperties interface, but no such luck. I have also tried casting the IFeatureDataset from IDataset, but obviously, that wont' work since I'm using a featureclass to get to the featuredataset that it's from and it's not a featuredataset to begin with.


texts: Using List field function for Python script tool?

texts:
I am trying to write a simple script for a toolbox, similar to the Addfield_datamanagement toolbox. The script adds a new field to a table of a feature class or a shapefile with attribute tables.
So far so good

# Import system modules
import arcpy
#Workspace
arcpy.env.workspace="C:\\Users\\xxx.gdb"
# Set feature class
inputFC = arcpy.GetParameterAsText(0)
# Get List of Fieldnames
inputString =arcpy.GetParameterAsText(1)
inputString =arcpy.ValidateFieldName(inputString,arcpy.env.workspace)
fieldList =inputString.split(";")
fieldType = arcpy.GetParameterAsText(2)
for fieldName in fieldList:
          arcpy.AddField_management(inputFC , fieldName, fieldType)
arcpy.AddMessage("Field generated:" + fieldName)
arcpy.AddMessage ("Script executed")

So this seems to work, but now i want to check first if the fieldname already exists before it´s created and print the existing name!
I thought about the list fields together with the built-in len() function:

if len(arcpy.ListFields(fieldName(?) or inputString,?)==1:
             arcpy.AddMessage("Field already exists"+fieldName)

Not sure about the synthax of the ListFields command and how to integrate it in the code above!


texts:
I'm having an attribute table like this of a shapefile on QGIS:

Movie / Actor

A   /    1

A    /   2

A    /   3

B    /   4

I want to get the name of a movie and all actors in that movie, and I want the result to be in a format like this:

Movie /  ActorList

A   /    1, 2, 3

B    /   4

How can I do it?


texts: Mosaic data looks like checkerboard using ArcMap

texts:
I have copied mosaic data from the external disk to the local disk. When I try to load the mosaic data into a new ArcMap document, it is displayed as a checkerboard. But if I first load the mosaic data from the external disk and then the mosaic from the local disk, the mosaic from the local disk is also displayed correctly.
What could be the reason for this and how can I fix it?
When i use ArcCatalog , and then use the repair path (or i have tried tool repair mosaic dataset in ArcMap) the old path is the one that is on internal drive


texts: Maybe I'm just missing something here, but i can't figure out how to set a WMTS Layer as an Input Parameter in QGIS Modeler. I just wanted batch the process of exporting a GeoTiff from a WMTS Layer. Now it takes like 4 steps or so to do it. Since I use that a lot, I wanted to make it a little Model. Anyway I am not to good with Python, so I tried Modeler.


texts: Unexpected EOF while looking for matching

texts:
I was wondering if someone could help me with this error. I am new to writing shell script so here is my problem. When I run the following script it gives the error,
line 31: unexpected EOF while looking for matching
line 33: syntax error: unexpected end of file

g.remove rast=tmp
i=1
while [ $i !=366 ]
do
     j=1951
     r.mapcalc "tmptot = 0.0"
     r.mapcalc "tmpsum = 0.0"
     r.mapcalc "tmptot2 = 0.0"

     while [$j !=2008 ]
     do
         g.copy rast=prec$j.$i,tmp
         r.mapcalc "tmp2 = if(isnull(tmp),0.0,tmp) + tmptot"
         r.mapcalc "tmptot = tm2"
         r.mapcalc "tm3 = if(isnull(tmp),0.0,tmp*tmp) + tmptot2"
         r.mapcalc "tmptot2 = tm3"
         r.mapcalc "tmp4 = if(isnull(tmp),0.0,1) + tmpsum"
         r.mapcalc "tmpsum = tmp4"
         g.remove rast=tmp,tmp2,tmp3,tmp4
         j=`expr $j + 1~
     done

    r.mapcalc "tmp5 = tmptot / tmpsum"
    r.mapcalc "tmp6 = (tmptot * tmptot) / tmpsum"
    r.mapcalc "tmpvar = sqrt((tmptot2 - tmp6) / (tmpsum - 1))"
    g.rename rast=tmp5,av$i
    g.rename rast=tmpvar,sd$i
    g.rename rast=tmpsum,sum$i
    g.remove rast=tmp6

   i=`expr $i + 1`
done


texts:
I want to work with PyQGIS Python 3 using QGIS 3.2 standalone install. I created a simple batch file to work with PyQGIS but I have some errors.

Here's the batch file:

@echo off
SET OSGEO4W_ROOT=C:\Program Files\QGIS 3.2

call "%OSGEO4W_ROOT%"\bin\o4w_env
call "%OSGEO4W_ROOT%"\apps\grass\grass-7.4.1\etc\env.bat
@echo off
path %PATH%;%OSGEO4W_ROOT%\apps\qgis\bin
path %PATH%;%OSGEO4W_ROOT%\apps\grass\grass-7.4.1\lib
path %PATH%;%OSGEO4W_ROOT%\apps\Qt5\bin
path %PATH%;%OSGEO4W_ROOT%\apps\Python36\Scripts

set PYTHONPATH=%PYTHONPATH%;%OSGEO4W_ROOT%\apps\qgis\python
set PYTHONHOME=%OSGEO4W_ROOT%\apps\Python36

SET PYCHARM="C:\Program Files\JetBrains\PyCharm 2018.1.4\bin\pycharm.exe
start "PyCharm aware of QGIS" /B %PYCHARM% %*

Error is:

  File "C:\PROGRA~1\QGIS3~1.2\apps\Python36\lib\site.py", line 177
    file=sys.stderr)
        ^
SyntaxError: invalid syntax


texts: Write World File option blocked when exporting Layout View in ArcGIS Desktop?

texts:
When i try to export map to jpg file in ArcGIS DEsktop 10.3, the "write world file" is blocked.

How can i solve this problem?


texts: Export Shapefile with UTF-8 encoding for non-english characters in GeoPandas?

texts:
This part of code puts the EΞ characters in a column and when exported it shows ?? in their place.

dff['KAEK']=dff['CAD'].str[:5].add('ΕΞ').add(dff['CAD_AD_F'].str[5:]).add(dff['rank'])
dff.to_file(out)

Where can i set the utf-8?


texts: Transforming Lambert to WGS84 in SQL Server

texts:
Need to create a SQL view of calculated Lat/Lon values from an existing SQL view of County Coords.
Loaded up SQL Server Spatial Tools and registerd CLR types. Lambert generated a resource error.
-- Project point and linestring using Lambert Conformal Conic projection

declare @lambert Projection

set @lambert = Projection::LambertConformalConic(0, 0, 0, 0)

select @lambert.Project('POINT (45 30)').ToString()

select @lambert.Unproject(@lambert.Project('LINESTRING (10 0, 10 10)')).ToString()

select @lambert.ToString()

--
XY going in=(521387.001,204215.954) Need LatLon coming out.
Do you have any ideas about how to do this?


texts: How to install PostGIS 8.4 in Windows 7 64 bit parallel to an existing PostGIS 9.2?

texts:
PostGIS 9.2 is installed along with OpenGeoSuite using one click setup. But 9.2 doesn't support pgRouting. So I need to install 8.4. I don't how to include binaries and all.

So someone please tell me how to install it in Windows 7 64 bit in detail including what port etc.


texts: Why is toolbox of ArcGIS Editor for OpenStreetMap not in ArcToolbox?

texts:
I'm using ArcGIS 10.0 and I tried to install ArcGIS Editor for OSM (version 2.0) on ArcGIS

But I couldn't make my ArcToolbox list the tools that ArcGIS Editor for OSM should've installed.

Any ideas about that?


texts: ERROR 000582 from clip processes in ArcPy

texts:
Durring execution of clip in arcpy, an error has occurred. It is showing ExecuteError: ERROR 000582: Error occurred during execution.
import arcpy

arcpy.env.workspace = (r"L:\\Arcpy\\Data\\vector_data\\selected_district")

rasterlayer = (r"L:\\Arcpy\\Data\\Mosaic")

output = (r"L:\\Arcpy\\Data\\clip")

vector_layer = arcpy.ListFeatureClasses()

for vectordata in vector_layer:
        # the first part of the file name before extension (if present) then add .tif
        output = vectordata.split('_')[0] + '.tif'
        arcpy.Clip_management(rasterlayer, "", output, vectordata, "0", "ClippingGeometry", "NO_MAINTAIN_EXTENT")


texts: Dissolve multiple rows straight into one layer and add OIDs into table

texts:
I am currently running this script to direct multiple dissolves of a polygon layer

FC = "BLPUsCavPlace"
TABLE = "BLPUsCavPlaceContigBLPUs"
FIELDS = ["blpuid", "contigid"]
with arcpy.da.SearchCursor(TABLE, FIELDS) as cursor:
    for row in cursor:
        OID1 = str(int(row[0]))
        OID2 = str(int(row[1]))
        SQL = "OBJECTID = %s OR OBJECTID = %s"%(OID1, OID2)
        arcpy.MakeFeatureLayer_management(FC, "TEMP_LYR", SQL)
        OUTPUT_NAME = "Dissolve_%s_%s"%(OID1, OID2)
        arcpy.Dissolve_management("TEMP_LYR", OUTPUT_NAME)
        arcpy.Delete_management("TEMP_LYR")

However, this knocks out individual dissolved layers. I can always append them afterwards, but I would like to keep it clean and write them straight into another polygon layer called BLPUsCavPlaceDissolved. I would also like each row to contain the two OIDs that were dissolved, so there would be two extra columns called OID1 and OID2 that record them in the layer. Is this possible? Or will dissolve only ever create individual layers?

Thanks


texts:
I compute shortest paths between ports. However, some of the resulting paths are simply incorrect.
Consider the following example where coast is multi-polygon object of the planet's landmass (e.g. based on GSHHG).
# Load packages
packs <- list("tidyverse", "raster", "sf", "tmaptools", "gdistance", "fasterize")
lapply(packs, require, character.only = T)

# Geocode ports
brisbane <- geocode_OSM("Brisbane, Australia", as.sf = T, geometry = "point")
liverpool <- geocode_OSM("Liverpool", as.sf = T, geometry = "point")

# Derive transition layer
tr_paths <- raster(crs = "+proj=longlat +datum=WGS84 +no_defs", vals = 1, resolution = c(1/8, 1/8), ext = extent(c(-180, 180, -90, 90)))
tr_paths <- fasterize(coast, tr_paths) %>%
   mask(tr_paths, ., inverse = T, updatevalue = 1000000) %>%
   transition(., function(x) 1/mean(x), directions = 8, symm = T) %>%
   geoCorrection(., type = "c", scl = F)

# Derive shortest path
path <- shortestPath(tr_paths, as(brisbane, "Spatial"), as(liverpool, "Spatial"), output= "SpatialLines") %>%
   st_as_sf()

This is the resulting path:

The problem appears to mostly affect particularly long paths. Shorter connections, e.g. between Recife in Brazil and Liverpool in the UK, are correctly computed.
The gdistance::shortestPath function that I use employs Dijkstra's algorithm implemented via the igraph package. Any ideas why that function produces such odd results? Does igraph change its optimization above a certain path length?


texts: Recently-opened documents option in QGIS

texts:
Is there a recent documents option somewhere in QGIS similar to other software? Or do people always use add layers / drag and drop from the Browser, or ensure that  all is saved as a project each time? Looks like the Project option is basically the same as the .WOR file in MapInfo which I'm very familiar with. Very handy but can be a bit of a pain if files start moving round...

I've just had QGIS hanging up on me a few times when experimenting and requiring restart of program. A recently opened option would be nice. Not a priority and I know the spinning wheel usually means lose of record of last file opened, not saved, etc.


texts:
We exported DEM, point, line and polygon data as .mdb or .gdb formats from ArcGIS. When we import these data to MapInfo, present database and when choosing the data from feature class there is only the point, line and polygon data shown, but there is no DEM data. I could try to import the DEM in raster format (.tiff or .bmp) but I believe that the image will get spoiled.

Is there any other way to import a DEM on MAPublisher?


texts: Modifying Time attribute using ModelBuilder?

texts:
I need to remove the DD-MM-YY from the "Time" attribute and only keep the HH:MM:SS using the modelbuilder. The model is displayed below and I use ArcMap 10.2.

I have tried the "Calculate Field" and the "Right function" but it only works on string.

Is this possible?


texts: Determining Min/Max values from feature class using ArcPy?

texts:
I am very rusty/inexperienced in Python.

I would like to create a more compact script that prints min and max values of each numeric field from a feature class. Right now I have the code below repeated for each field that contains numeric values.

How would I code this to identify if a field has numeric values and then print the min and max values to the python shell or even better to a .csv/text file?

import arcpy
data = "C:/Workfolder/Input.gdb/FeatureClass"

listname = []
cursor = arcpy.da.SearchCursor(data, "field1")

for row in cursor:
    listname.append(int(row[0]))

print 'Field1_min, {0}'.format(min(listname))
print 'Field2_max, {0}'.format(max(listname))


texts:
Is there a way to speed up rendering data into ArcMap ?

When I converted the raster to polygon and some data would render quicker and some other renders slow.

I am wondering if I could add the spatial index to the data that renders slow.

Would that work or is there anything that I could make the data render quicker?


texts:
Is there any way to import the OSM data into pgRouting with the toll data?

I already have a database with OSM data with topology, and working with pgRouting, but I need to add some nodes to it.

I can see that in the barrier=toll_booth tag I can have the information that I need, and I can download this points from OSM separately.

The key is:

Can I add this nodes to my current network,  in order to get it if my route crosses it?

Or if not,

Can I re-download all the data with osm2pgrouting or osm2po including this nodes too?


texts: Why are the holes in my polygons not generalized?

texts:
I have one polygon layer with many holes. I would like to generalize it. I used v.generalize of grass. With it, I used chaiken algorithm. But, only the borders of polygons were generalized, the holes in them haven't smoothed their edges.

What can I do to smooth all edges (adding hole edges) of polygon?


texts: Using CASE Statement with two strings in QGIS

texts:
I am trying to fill out a field based on information in a different field in the same layer using a CASE Statement in the Field Calculator. My table structure looks like this:

SPECIES
SHORT

Toad
NULL

NULL
B

Fish
NULL

Plant
NULL

My goal looks like this:

SPECIES
SHORT

Toad
T

NULL
B

Fish
F

Plant
P

The Expression I am using and does not work looks like this:
CASE
   WHEN "Species" LIKE 'Toad' THEN 'T'
   WHEN "Species" LIKE 'Fish' THEN 'F'
   WHEN "Species" LIKE 'Plant' THEN 'P'
END

The result is that everything in "SHORT" ends up being NULL.
The Field "SPECIES" is a VALUE MAP widget type. "SHORT" is a TEXT EDIT.
I am using QGIS 3.1


texts: Converting Python GDALDataset to C++ GDALDataset

texts:
I'm currently writing a C++ extension to python to speed up a raster data viewer which is running extremely slowly with our large datasets. I've managed to get a small speedup by writing a C++ extension, but after some investigation the slowest part in the extension is a python callback to a GDALRasterBand. Alot of this time is due to the overhead of doing a callback to python from a C-Extension,
So I was wondering if anyone knew of a way to easily convert a GDALDataset from python to a C++ GDALDataset.
From my research into it I believe they use the same backend structure, but the python one has the additional python headers, so I don't think using the pointer reference to the python version will work.


texts: Adding elevation data from contour lines to point layer using QGIS?

texts:
I want to calculate the number of vacant lots within 300m from a given building. My source file contains point layers of buildings and vacant lots.

Here's a short description of my workflow. I've successfully completed steps 1 and 2. How do I complete step 3 ?

Create heatmap of vacant lots @300m radius.
Create a layer of contour lines based on the heatmap.
On "buildings" layer, add attribute indicating which contour line each building is located in.


texts: How to get gdaltranslate to create world file

texts:
Am trying to recreate "gdal... -co worldfile=yes"

However whenever I do I get error: "Warning 6: driver GTiff does not support creation option worldfile"

Though I think it may be working anyway. Any ways to prevent error message, or workarounds?


texts: How to get the result of WMS GetFeatureInfo request result into an array?

texts: I am using OpenLayers 2.1.3 and using OpenLayers I obtain the GetFeatureInfo result but it is a web page, I want to get the output fields into a variable so I can display them in grid.


texts: How to use VBA in the field calculator, to copy the last 7 characters to a new field?

texts: Im using ArcMap 9.3.  Im looking for a VBA code to use in the field calculator to make a new field of the last 7 digits of a "pin" field.  For example.  my pin is 00-00-999-999 and i want it to display as 999-999 in my new field.  Anybody have code for this?


texts:
I'm trying to add a point to a map that already has points (table imported from a public source).  However I do not have the 'button' on the right hand side in my 'map view' to be able to add in my new point.  Any reason why this button might be absent?

I'm using CartoDB.


texts: Qgis Web Client change the default EPSG 3857 to 29194 SAD 69 UTM ZONE 24

texts:
The Qgis standard Web client's initial installation is EPSG default +3857 like to know how to change EPSG 29194 SAD 69 UTM ZONE 24 because all my projects are in this coordinate system.

What is the procedure and changes should I make?


texts: How do I edit shapefile attribute tables so they all have the same information?

texts:
I am trying to merge several pipeline shapefiles using GEOMERGE but it cannot complete the function because they have different column definitions.

So I need to edit the attribute tables so that they all have the same data/layout.

For e.g one pipeline has the following columns in the attribute table: ID   REF   TYPE

and in another pipeline attribute table it has:  Leg-No   From   To   Length_Km

What is the best way for me to overcome this?

If it is relevant some of the shapefiles have been converted from KML to shapefiles using GQIS.


texts: Portraying quarter circle symbology with multiple attributes in ArcGIS Desktop?

texts:
I have four columns, (A, B, C, D) and each can be ranked 1, 2, or 3 (point data).

I would like to portray each category at once based on colors assigned to the rankings (1 = green, 2 = yellow, 3 = red).

I am imagining this would look like a pie chart with four quarters, and each quarter would be colored based upon the ranking.

I'm using ArcGIS Desktop.


texts: Volume using VB Script of the field calculator

texts:
I am trying to compute the volume of the area given certain condition. I wrote the following VBScript code, however I am getting error that it has failed to do so.

Can anyone look into it and suggest a way forward.

if [Area_Ha]<=70 Then
 [Volume__De]=2.85* [Area_Ha]^1.22
elseif  [Area_Ha] >70 Then
 [Volume__De] =7.1* [Area_Ha]+9.97
endif

The picture of the attribute table is attached.

Attached is the error message.

Matte,
The modified code didn't work


texts: Finding distances associated with change in latitude and longitude at given location

texts:
I am writing a program that uses GPS lat and long coordinates and a time change to determine speed with the output in units of m/s (precise to 0.1 m if feasible, but 1 m is sufficient).  When researching the distances associated with a change in latitude and longitude, however, I saw just how variable it was.  When continuing my research to find the distances associated with the general location I'm interested in, I could not find such a resource.

I know what area this program will be used in (Boston, Massachusetts) and the distance travelled during each time step is short enough that I can assume it moved across a plane.  As such, I am not looking for a formula to use in the general case over long distances but a look-up table or data base that has the relevant distances I need.

Where should I look online (or elsewhere) for this information?


texts: Dividing polygons into equal parts by the number of points inside

texts:
I would like to divide many polygons into 2, 3 or 4 parts.

I have vectorized plants. Sometimes 2, 3 or more plants hide under one polygon. Depending on the size I just have to divide the polygon into 2, 3 or 4 equal parts. Target is to set polygon centroids and count point (plants).

  Enhanced from a comment

Actually I only need points in the polygon depending on the size of the polygon.

Polygon size < 0.004m² => Single plant => One point
0.004m² < Polygon size < 0.008m² => Two plants => Two points
0.008m² < Polygon size < 0.012m² => Three plants => Three points

It does not matter what position has the points in the polygon. They should only be within the polygon. In the end, I want to count the points in my field (that's easy)


texts: Cloud-Optimized GeoTIFF overlay on Leaflet?

texts:
I want to load external Cloud-Optimized GeoTIFF files onto a Leaflet map as overlays, and manipulate the colors based on user input. I've found examples online that appear to to this exact thing, but with OpenLayers: https://geotiffjs.github.io/cog-explorer/
I've also found a few Leaflet plugins, of various age and stability, but they all appear to have caveats:

My COG files are big, so it would need to load parts of the tiff file on demand depending on map position/zoom, (rather than the whole file at once as you would with a regular raster overlay). The OpenLayers demo above does this.
Ideally, we could do this with EPSG 3857 projection

Is Leaflet capable of doing this?


texts: How many GIS projects deal with just basic "putting markers on a map"?

texts:
Is there an academic paper or conference contribution, where someone estimated how many GIS projects actually deal with only a very basic usage of GIS data? Something equivalent to the infamous "80% of all data is spatial" quote. such as

"X% of all use-cases of GIS data are concerned with plotting simple point data on maps"

There was a short rant about this question in one of the weekly podcasts from Directions magazine, but I don't remember which podcast, and it was only a qualitative statement. If I remember correctly, someone said "here we are in 2012 and most of us in the GIS community are still mainly putting data on a map.'

I am interested in getting a quantitative answer to this. If this is possible. The answer does not have to be academic. A back-of-the-envelope calculation will do as well.

Edit (w.r.t. comment): I don't intend to trivialize typical GIS programming tasks here, but I cannot really draw a line where "serious" work starts. I 'd like to emphasize that putting even just a single point on a map accurately can indeed be quite complex. (Understand coordinates and determine coordinates of a location, choose a map projection, perform a coordinate transform, then add many more points on map when first is right, classify etc). But most projects even don't work with "line"  data type, or let alone polygon data types.

That is not just typical for the GIS community. There are equivalent studies from the relational-database community. Market researchers have found out that surprisingly high number of database-backed applications contain only a single table.


texts:
I want to use GeoTools to convert collection of SimpleFeature objects to GML 3.2.1 geometry.

However, I'm dealing with this:

Creation of XSD works with both GML 2 and 3.
Creation of GML works only with GML 2.
With version 3, it creates empty string.

This is the piece of code:

    import org.geotools.GML;
    import org.geotools.data.DataUtilities;
    import org.geotools.data.collection.ListFeatureCollection;
    import org.geotools.feature.simple.SimpleFeatureBuilder;
    import org.geotools.geometry.jts.WKTReader2;
    import org.opengis.feature.simple.SimpleFeature;
    import org.opengis.feature.simple.SimpleFeatureType;

...

    SimpleFeatureType featureType = DataUtilities.createType("location", "geom:Point,name:String");

    File xsdFile = new File("C:/Dev/VTS/GML_samples/location.xsd");
    xsdFile = xsdFile.getCanonicalFile();
    xsdFile.createNewFile();
    URL xsdURL = xsdFile.toURI().toURL();
    URL baseURL = xsdFile.getParentFile().toURI().toURL();

    // Creation of XSD file using GML 3 works:
    GML xsdEncoder = new GML(GML.Version.GML3);
    FileOutputStream xsd = new FileOutputStream(xsdFile);
    xsdEncoder.setBaseURL(baseURL);
    xsdEncoder.setNamespace("location", xsdURL.toExternalForm());
    xsdEncoder.encode(xsd, featureType);
    xsd.close();

    // Simple Feature collection:
    WKTReader2 wkt = new WKTReader2();
    List<SimpleFeature> collection = new LinkedList<SimpleFeature>();
    collection.add(SimpleFeatureBuilder.build(featureType, new Object[] { wkt.read("POINT (1 2)"),"name1" }, null));
    collection.add(SimpleFeatureBuilder.build(featureType, new Object[] { wkt.read("POINT (4 4)"),"name2" }, null));

    // Creation of GML file using GML 3 doesn't work. It does work with GML 2 though:
    // GML gmlEncoder = new GML(GML.Version.GML2); // this works
    GML gmlEncoder = new GML(GML.Version.GML3); // this doesn't work
    ByteArrayOutputStream xml = new ByteArrayOutputStream();
    gmlEncoder.setLegacy(true);
    gmlEncoder.setBaseURL(baseURL);
    gmlEncoder.setNamespace("location", "location.xsd");
    gmlEncoder.encode(xml,  new ListFeatureCollection(featureType, collection));
    xml.close();

    System.out.println(xml.toString());


texts: Possible for there to be same street, city, and country in two locations

texts:
I am working on a project, and need to know if it is at all likely in most of the world that there can be two streets in a city with the same name, in the same country, but in two different locations.

Eg. Example Street, Chicago, USA & Example Street, Chicago, USA. These would be different streets. It doesn't matter to me if they have different address numbers on them.

Is this likely in the USA or Europe?


texts: Regression between reflectance at pixel value vs input daily data over image collection

texts:
I have been trying to make a raster of regression values between reflectance value vs input data at each pixel position in a collection (So, for example the R^2 relationship between daily average precipitation and reflectence value at every pixel position within a given extent). I have tried doing this by filtering MODIS by a list of dates and applying a constant band to each image with a respective input daily data value from CSV data, but I am having trouble extracting a specific CSV value for an individual date and applying it to an individual, respective image. The following code has returned the error "single positional indexer is out-of-bounds" when using the python Google Earth Engine API :
    import numpy as np
    import pandas as pd
    import ee
    ee.Initialize()
    df=pd.read_csv('.csv')
    ROI=ee.Geometry.Polygon([coords])

    dates=df['Date'].tolist()
    df['DateMillis'] = df['datetime'].astype(np.int64) / int(1e6)
    data=df['Data'].tolist()
    date = ee.List(dates).map(lambda date: ee.Date(date).millis());
    d=ee.List(data)

    def maskModisClouds(image):
        return image.updateMask(image.select('state_1km').bitwiseAnd(1 << 10).eq(0))

    def timeband(image):
         return image.addBands(image.metadata('system:time_start'))

    def databand(image):
        return image.updateMask(ee.Image(((df.loc[df['DateMillis']==image.get('system:time_start'),'Data']).iloc[0])).addBands(image))

    col = ee.ImageCollection('MODIS/006/MOD09GQ').filter(ee.Filter.inList("system:time_start", date));
    Image = col.filterBounds(ROI);
    image=Image.map(maskModisClouds).map(timeband).map(Databand)


texts: Accessing raster dataset's attribute table?

texts: I need to loop through the rows in the raster attribute table and add a date for each row.


texts: Weird output when projecting/interpolating points on lines

texts:
I am trying to project a points' shapefile on a multilinestrings' shapefile. (Project each point to the closest line)

Here's my code :

import geopandas as gpd
import time

start_time = time.time()

points4326 = gpd.GeoDataFrame.from_file('./Outputs/points.shp')
points = points4326.to_crs({'init': 'epsg:2269'})
lines = gpd.GeoDataFrame.from_file('./Shapefiles/lines.shp')
projected_points = gpd.GeoDataFrame(columns=['ID','Arrow','Comments','ZoneName','geometry'], crs={'init' :'epsg:2269'})

for index, row in points.iterrows():

    pt = gpd.GeoSeries(row['geometry'])
    # Project
    projOnLines = lines.project(pt)
    # print projOnLines

    #  Interpolate
    ptInterpolate = lines.interpolate(projOnLines)
    # print ptInterpolate
    projected_points.set_value(index, 'geometry', ptInterpolate.iloc[0])
    # print projected_points

    # Set the value of each column of projected_points with the value of the corresponding row in the input points' shapefile
    projected_points.set_value(index, 'ID', row['OBJECTID_l'])
    projected_points.set_value(index,'Arrow',row['Arrow'])
    projected_points.set_value(index,'Comments',row['Comments'])
    projected_points.set_value(index,'ZoneName',row['ZoneName'])

print projected_points.head(5)
projected_points.to_file('./Outputs/pointsOnLines.shp')

print "time elapsed: {:.2f}s".format(time.time() - start_time)

If I just put :

projOnLines = lines.project(pt)
ptInterpolate = lines.interpolate(projOnLines)

I get this error :

ctypes.ArgumentError: argument 3: <type 'exceptions.TypeError'>: wrong type

The problem comes from the fact that projOnLines's result is a Geoseries, and interpolate expects one float. Here's what print projOnLines shows :

0        199.894
1          False
2          False
3          False
4          False
      ...
69955      False
69956      False
69957      False
dtype: object

Now If I take only that first value and pass it to interpolate like this :

ptInterpolate = lines.interpolate(projOnLines[0])

The script takes a long time to execute and when I finally get to display the result on QGIS I see all points are overlapped on one same location! (the attribute table has all features but they're somehow gathered in one place although I'm looping through all the points in the shapefile)

Any ideas?


texts: Cluster spatial - cannot allocate vector

texts:
I'm trying use this code for create cluster spatial in R, I get from @ssanch:

# example data from the thread

x <- c(-22.83982, -23.49961, -27.08745, -22.98890, -19.92278)
y <- c(-47.20671, -46.85583, -52.62162, -43.22838, -43.95431)

# convert data to a SpatialPointsDataFrame object
xy <- SpatialPointsDataFrame(
  matrix(c(x,y), ncol=2), data.frame(ID=seq(1:length(x))),
  proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))

# use the distm function to generate a geodesic distance matrix in meters
mdist <- distm(xy)

# cluster all points using a hierarchical clustering approach
hc <- hclust(as.dist(mdist), method="complete")

# define the distance threshold, in this case 40 m
d=275000

# define clusters based on a tree "height" cutoff "d" and add them to the SpDataFrame
xy$clust <- cutree(hc, h=d)

But I change the latitude and longitude for my data with 600k registers and the R can't processing with this errors:

  Error: cannot allocate vector of size 8.0 Gb

Is there any solution?


texts:
I have a population density raster data, but I'm trying to calculate the density per pixel? how do i go about this?

working with ArcGIS


texts: How to obtain road classification type at a point via nominatim?

texts:
This is my first time posting here so please forgive me if i ask a trivial question. Is there a way of obtaining the type of road an object is settled in from a request via nominatim? (i have the geocoded sample)

For example:
http://nominatim.openstreetmap.org/reverse?format=xml&lat=52.5487429714954&lon=-1.81602098644987&zoom=18&addressdetails=1

does provide me with loads of info. yet, there is no info on the road the address is settled in. would like to obtain something like: rural road or even just grade*N* (where N stands for the level of hierarchy.) i hope what i'm looking for is understandable :) thanks in advance!


texts:
This is my code:
GeoRef = georasterref('Rastersize',[1440,721],'Latlim',[-90,90],'Lonlim',[-180,180]);

data=ncread('G:\ERA_humid\humid.nc','r');

for year=1979:2020
        data1=data(:,:,1+12*(year-1979):12*(year-1978));

        data3=sum(data1,3)/12;

        data4=rot90(data3);

        filename=strcat('G:\ERA_humid\year_mean_humid\',int2str(year),'year_mean_humid.tif');

        geotiffwrite(filename,data4,GeoRef)

    for mon=1:12

        data2=data1(:,:,mon);
        data4=rot90(data2);
        filename=strcat('G:\ERA_humid\month_humid\',int2str(year),'_',int2str(mon),'month_humid.tif');
        geotiffwrite(filename,data4,GeoRef)
    end
end

The error is:

Untitled23 (line 9)
geotiffwrite(filename,data4,GeoRef)


texts:
I m triying to add a secured feature service to CMV using a proxy page.

My viewer.js file is the following:

esriConfig.defaults.io.proxyUrl = 'myProxyURL';
esriConfig.defaults.io.alwaysUseProxy = false;

operationalLayers: [
{
            type: 'feature',
            url: 'myURL',
            title: i18n.feature1,
            options: {
                id: 'feature_layer_1',
                opacity: 1.0,
                visible: true,
                outFields: ['*'],
                mode: 0
            }
},

And my proxy.config file:

<?xml version="1.0" encoding="utf-8" ?>
<ProxyConfig allowedReferers="*"
         logFile="proxy_log.log"
         mustMatch="false">
    <serverUrls>
        <serverUrl
            url="myURL"
            username="myusername"
            password="mypassword"
            rateLimit="120"
            rateLimitPeriod="60"
            matchAll="true"/>
    </serverUrls>
</ProxyConfig>

But when the CMV is starting, it prompts me to put username and password... What am I doing wrong?


texts:
I am new with leaflet.

I am trying to implement easyButton plug-in: https://github.com/CliffCloud/Leaflet.EasyButton/tree/v1.x

I am using version 1.3.2 because on my project I use leaflet 0.7.7.

I use these instructions to download:

cd your/project/javascript-files/
curl -O https://raw.githubusercontent.com/CliffCloud/Leaflet.EasyButton/master/src/easy-button.js
# saved at your/project/javascript-files/easy-button.js

cd your/project/css-files/
curl -O https://raw.githubusercontent.com/CliffCloud/Leaflet.EasyButton/master/src/easy-button.css
# saved at your/project/css-files/easy-button.css

Calling the library on my code:

<script src="assets/js/easy-button.js"></script>
<link rel="stylesheet" href="assets/css/easy-button.css">

I run this code:

var buttoms = L.easyButton('fa-star', function(btn, map){
    var antarctica = [-77,70];
    map.setView(antarctica);
});
L.easyBar(buttoms).addTo(map);

And I got this error:

ReferenceError: L is not defined       [Learn More] easy-button.js:6:1
TypeError: L.easyButton is not a function

I cannot understand why I got this error and how to solve it..


texts:
I am working on a web map using ArcGIS JavaScript API. The map should be able to query a layer and show the result. Since this is my first time using query task of the API, I got some issue. Below is the syntax of the query part. When I clicked the button in webpage, It always give the error:

init.js:182 Error(…)y
@ init.js:182(anonymous function)
@ init.js:181c
@ init.js:66d
@ init.js:66resolve.callback
@ init.js:67c
@ init.js:66d
@ init.js:66resolve.callback
@ init.js:67c
@ init.js:66d
@ init.js:66resolve.callback
@ init.js:67c
@ init.js:66d
@ init.js:66resolve.callback
@ init.js:67(anonymous function)
@ init.js:185k
@ init.js:189q
@ init.js:189resolve
@ init.js:191a
@ init.js:190k
@ init.js:190q
@ init.js:189resolve
@ init.js:191a
@ init.js:190k
@ init.js:190q
@ init.js:189resolve
@ init.js:191a
@ init.js:190k
@ init.js:190q
@ init.js:189resolve
@ init.js:191n
@ init.js:150f
@ init.js:154

As you can see, I put some console in the syntax and here is the result:

executequery started
where statement completed in executequery function
myquerytask completed in executequery function

It looks like the queryTask has already been excuted but the showResults, which is the callback function, doesn't.
Could anybody help me with the entry level problem?

var whereExpression = "";
myQueryTask = new QueryTask("http://arcgis-cabhpgis-2026970499.us-west-1.elb.amazonaws.com/arcgis/rest/services/HCIC/160524_interactive_map_jzz_v02/MapServer/0");
myQuery = new Query();
myQuery.returnGeomery = true;
myQuery.outFields=["*"];
myQuery.where = whereExpression;
function executeQuery(){
    console.log("executequery started");
    whereExpression = "";
    for(i=0;i<serviceList.length;i++){
        if(i==0) {
            whereExpression += ("ServiceProvided LIKE '%" + serviceList[i] + "%'");
        }
        else{
            whereExpression += (" AND ServiceProvided LIKE '%" + serviceList[i] + "%'");
        }
    }
    console.log("where statement completed in executequery function");
    myQueryTask.execute(myQuery,showResults);
    console.log("myquerytask completed in executequery function");
}
function showResults (results) {
    console.log("showresults started");
    document.getElementById("info").innerHTML = results.features.length;
}
on(dom.byId("applyServiceSelection"),"click",executeQuery);


texts: Using Near (proximity) tool in ArcGIS for Desktop?

texts:
I've tried changing GCS and PCS multiple ways but can't seem to figure out how to get the Near (from Proximity tools) to output the Near_DIST in feet. All it outputs is what looks like decimal degrees. Here is my data:

Data frame: NAD_1983_StatePlane_Texas_North_FIPS_4201_Feet

Input features (points): GCS_North_American_1927

Near features: GCS_North_American_1983

I'm not fond of the Data frame features but it's for my geographic area and ensures that when I draw buffer circles they actually look like circles and are not skewed to appear oblong. Based on the tool, it says the output Near_DIST will return units of measure based on the Input features. I'm not sure what I should change to get it to return the distance in feet.


texts: QNEAT3 Iso-Area as polygons algorithm "float division by zero" error

texts:
I tried to generate isochrones using Iso-Area as polygons algorithm in QNEAT3 plugin. I fed the algorithm by the speed field of the road network attribute table, I set the default speed to 0 and topology tolerance to 10, I have also selected the fastest path for generating isochrones. After running the algorithm I got the following error, ZeroDivisionError: float division by zero
It seems that despite the fact that I fed the algorithm by the speed field of the network layer attribute table, it still uses the value of default speed for generating isochrones.
How can I solve the problem?
The error message is:
Traceback (most recent call last):
File "C:/.../.../AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\QNEAT3\algs\IsoAreaAsPolygonsFromLayer.py", line 232, in processAlgorithm
list_apoints = [Qneat3AnalysisPoint("from", feature, id_field, net, net.list_tiedPoints[i], entry_cost_calc_method, feedback) for i, feature in enumerate(getFeaturesFromQgsIterable(startPoints))]
File "C:/.../.../AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\QNEAT3\algs\IsoAreaAsPolygonsFromLayer.py", line 232, in <listcomp>
list_apoints = [Qneat3AnalysisPoint("from", feature, id_field, net, net.list_tiedPoints[i], entry_cost_calc_method, feedback) for i, feature in enumerate(getFeaturesFromQgsIterable(startPoints))]
File "C:/.../.../AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\QNEAT3\Qneat3Framework.py", line 542, in __init__
self.entry_cost = self.calcEntryCostEllipsoidal(feedback)
File "C:/.../.../AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\QNEAT3\Qneat3Framework.py", line 559, in calcEntryCostEllipsoidal
return dist/(self.entry_speed*(unit_to_meter_factor * 1000.0 / 3600.0)) #length/(m/s) todo: Make dynamic
ZeroDivisionError: float division by zero


texts:
When I activate map click event, I want to suppress map move event

Do you have any idea how can I do that?


texts: Create atlas series with yearly maps from 1920 to 2018

texts:
I have a lot of maps I need to create. One for each year that shows what have happened in the map area. First map are supposed to shot 1920. Then the next one are supposed to shot 1920+1921. And the next after that are supposed to show 1920+1921+1922.

I'm wondering if it's possible to create an atlas that can do this ?

I have to create one each year until 2018.


texts: How to downgrade from Processing 2.6 to 2.2?

texts:
I have QGIS 2.8.2 installed and Processing Toolbox 2.6.

I want to test a thing and I need to install processing Toolbox 2.2 in QGIS 2.8.2. Is it possible? How can I do that?


texts: Creating Slope from DEM in ArcGIS Desktop so that information is in feet per mile?

texts:
I'm trying to get the Time concentration for the ModClark method in Hec-HMS, and one of the variables is Slope in feet per mile.

I have a DEM for the sub basin that I've download, the elevation was in metres but I've used ArcGIS Spatial Analyst - Math to convert to feet.

How could I get the slope value (preferably an average for the basin) that is in feet per mile?


texts: How do I know what caused the extent of my axMapControl to be updated?

texts:
I have an axMapControl in a WinForms application. According to the ArcGIS documentation,

  [The]OnExtentUpdated [event] is triggered whenever the IMapControl2::Extent
  property is changed. Setting the extent explicitly, zooming, panning
  or using the IMapControl2::CenterAt method will change the extent.

But how can I know (programmatically) inside the event handler that it was zooming or panning that caused the extent to change? I am looking for some C# code which would enable me to know that the extent change was caused by one of those two things.


texts: Create Shapefile of all overlapping features from single Shapefile

texts:
I need to create a Shapefile out of another Spahefile that consists all overlaps of features within the initial Shapefile.

I'm sure there is a simple solution using ArcMap or QGIS, but at the moment I am just lost and somehow I struggle to find something via search engines although I am aware that this question must have been asked a thousand times.
The screenshot shows the desired result. My new Shapefile should only consist the darker parts independent on the number of overlapping features.


texts: Replace values in one DEM with values from another DEM

texts:
I have two DEMS. One is a 10m DEM (let's call it DEM-A) covering a large area (projection EPSG 4326). The other is a higher resolution DEM (DEM-B) that was constructed from a structure-from-motion point cloud, but has since been down-sampled to 10m and projected into the same coordinate system (EPSG 4326).
I need to replace the values in DEM-A with the values in DEM-B where there are such, and keep all values of DEM-A where there there is nothing (i.e. nodata in DEM-B or no coverage of DEM-B).

In the image below, you can see parts of DEM-A on the right, and parts of DEM-B on the left. Black cells are nodata (-9999). You can see that the grids don't quite match, so DEM-B will need to be resampled to fit the grid in DEM-A.

I feel like this should be a straigthforward conditional statement, but I can't get it to work. Happy to get suggestions for QGIS, gdal, or anything open source.


texts: Create raster in R - distance to line feature

texts:
I have some line features (eg. 1000 meter isobath), and I would like to make a raster file for my study area where the value of each cell is the distance to the line feature.
Any idea what function I would use in R to accomplish this?


texts: Update geometry(polygon) of shapefile based on a point layer

texts:
I have two shapefiles.

Polygon shapefile
Point file contains some points.

Now I want to delete all the vertex of polygon shapefile which are available in point file.

I am using ArcGIS (ArcEditor) platform.

Any ideas?

This is result of procedure you have mentioned via GeoWizard.

Till now I have not deleted any point, just converted polygons to point and then point to polygons.


texts: Error using 'Clip raster by mask' in QGIS v3.8.2

texts:
I am attempting to clip a raster to a shapefile boundary (country boundary, 1 polygon) using the 'clip raster by mask' tool in QGIS v3.8.2.
Both layers are in WGS84, so projection mismatch shouldn't be the issue.

I'm getting the following error message that I haven't been able to sort out:

  2019-09-09T16:23:13 INFO GDAL execution console output
  Creating output file that is 435P x 193L.
  Processing R:\exampledata\wc20_25m_tmin_01.tif [1/1] : 0ERROR 1: IllegalArgumentException: point array must contain 0 or >1 elements

  ERROR 1: Cutline polygon is invalid.

  Using internal nodata values (e.g. -3.4e+38) for image R:\exampledata\wc20_25m_tmin_01.tif.

  Copying nodata values from source R:\exampledata\wc20_25m_tmin_01.tif to destination
  C:/Temp/processing_7b281fab84954f6ca3dd66d12d696f5c/0b2ef639021f42f0b744260b98b85235/OUTPUT.tif.

To correct this I've tried:

1) checking the geometry of the polygon file, no errors were found. Geometry is valid

2) correcting the geometry anyway. No change.

3) creating a zero distance buffer, didn't fix it.

4) tried a different vector shapefile, same issue.

5) tried a different raster with both shapefiles, same issue.

6) tried simplifying the geometry, no change to the issue.

7) discovered the single feature is in fact a 'multi polygon'. tried deleteing all the small features and saving the main area as a single polygon. This did not correct the issue.

8) Using the SAGA 'clip raster with polygon' tool. This produces a usable output with a warning:

  Loading resulting layers
  The following layers were not correctly generated. C:/Users/wkessler/AppData/Local/Temp/processing_1b387f0cf3564641840068eedfe5ec05/cf6aaedbd8e44f09b1eeb194d5d23a36/OUTPUT.sdatYou can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm.

Now, I have tried the same operation in both QGIS 2.18 (Clipper tool) and ArcGIS 10.6 (Clip) and have no issues. I have also managed to get the SAGA 'clip raster with polygon' function to work on this data in QGIS 3.8.2, although it does give a warning that the result was generated incorrectly. Can anyone tell me what's going on?

EDIT: Since posting this question I have also sent the data to a friend who tried the operations using QGIS 3.8.2, 3.4, and 2.18. They found the same thing; that it works in 2.18 but not in 3.4 or 3.8.2.

EDIT 2: tried additional solutions suggested by comments. None have solved the issue.


texts: Creating feature class in memory and adding to map using ArcObjects

texts:
I have the code below, which is supposed to create a new featureclass in memory then add the FC to the map, but I get a Automation error when I use the createfeatureclass method

  Set pFC = pFeatWS.CreateFeatureClass("NewBus", pFields, Nothing,
  Nothing, esriFTSimple, "Shape", " ")

Here' s the full code

Public Sub InMemoryWorkspace()

'Create Workspace factory
Dim pWorkspaceFactory As IWorkspaceFactory2
Set pWorkspaceFactory = New InMemoryWorkspaceFactory

'Create workspace
Dim pWorkspaceName As IWorkspaceName2
Set pWorkspaceName = pWorkspaceFactory.Create("", "MyInMemoryworkspace", Nothing, 0)

'Get a handle on Name
Dim PName As IName
Set PName = pWorkspaceName

'Get handle on Featureworkspace
Dim pFeatWS As IFeatureWorkspace
Set pFeatWS = PName.Open

'Create field object and set number of fields
Dim pFields As IFields
Set pFields = New Fields
Dim pFieldsEdit As IFieldsEdit
Set pFieldsEdit = pFields
pFieldsEdit.FieldCount = 3

'Create objectID field
Dim pField As IField
Set pField = New Field
Dim pFieldEdit As IFieldEdit
Set pFieldEdit = pField
With pFieldEdit
    .Name = "OBJECTID"
    .Type = esriFieldTypeOID
End With
Set pFieldsEdit.Field(0) = pField

'Create postcode field

Set pField = New Field

Set pFieldEdit = pField
With pFieldEdit
    .Name = "POSTCODE"
    .Type = esriFieldTypeString
End With
Set pFieldsEdit.Field(1) = pField

'Create Shape field
Set pField = New Field
Set pFieldEdit = pField
With pFieldEdit
    .Name = "SHAPE"
    .Type = esriFieldTypeGeometry
End With
Set pFieldsEdit.Field(2) = pField

Dim pGeomDef As IGeometryDef
Dim pGeomDefEdit As IGeometryDefEdit
Set pGeomDef = New GeometryDef
Set pGeomDefEdit = pGeomDef
With pGeomDefEdit
  .GeometryType = esriGeometryPolygon
  Set .SpatialReference = New UnknownCoordinateSystem
End With
Set pFieldEdit.GeometryDef = pGeomDef
pFieldsEdit.AddField pField

' Create featureclass
 Dim pFC As IFeatureClass
 Set pFC = pFeatWS.CreateFeatureClass("NewBus", pFields, Nothing, Nothing, esriFTSimple, "Shape", " ")

 Dim pmemLayer As IFeatureLayer
 Set pmemLayer.featureClass = pFC

 Dim pmxdoc As IMxDocument
 Set pmxdoc = ThisDocument
 Dim pmap As IMap
 Set pmap = pmxdoc.FocusMap

 pmap.AddLayer pFC

Can someone explain where I am going wrong?


texts:
I am trying to convert a MapInfo .TAB file to a Shapefile while setting the projection of the target to NAD83 UTM Zone 19. The latest version of OGR seems to not have this projection available in the drop down for the target file. There are other options for NAD83 but when I check the output files in ArcCatalog they seem to be projecting in a geographic coordinate system.

Does anyone know if NAD83 UTM Sone 19 is available?


texts: How to quickly get coverage of multiple services at once using pgrouting driving distance

texts:
I'm trying to get some useful information out of openstreetmap data to choose where to rent an apartment in my city.
I used osm2pgrouting and osm2pgsql to get OSM data into my postgis database.
So far I've managed to create a view (nearest_poi) with the id of the nearest neighbor of my POIs (some type of shops, bus stations and parks from OSM points and polygons), along with the category of the poi and the maximum distance I would like to be from that category (near_vertex_id, category, distance).

To be able to get the coverage in a single view I used pgr_drivingdistance from all the points of nearest_poi using the absolute maximum distance (max_distance is a table that stores the couple category-distance) and then I get rid of the points with distance greater than the maximum distance for that category:

CREATE OR REPLACE VIEW coverage_pts AS
SELECT id, the_geom, min(distance) as distance, category, max_distance
FROM (
  SELECT v.id, v.the_geom, di.agg_cost as distance, n.category, n.distance as max_distance
  FROM pgr_drivingdistance(
    'SELECT gid as id, source, target, length_m as cost FROM ways',
    array(SELECT DISTINCT near_vertex FROM nearest_poi),
    (SELECT max(distance) FROM max_distances),
    false, false
  ) as di
  JOIN ways_vertices_pgr v ON di.node = v.id
  JOIN nearest_poi n ON di.from_v = n.near_vertex
  WHERE di.agg_cost<=n.distance
) coverage
GROUP by id, the_geom, category, max_distance

Finally I get the points covered by all the services using:

SELECT * FROM (
  SELECT  id, the_geom, count(category) AS tot_services
  FROM  coverage_pts
  GROUP by id, the_geom
) AS hottest_point
WHERE tot_services = 20 --I have 20 categories

All is working fine except the last query takes about 10 seconds for 427 POIs on my 2008 MacBook Pro. When I try to visualize the points in QGIS, the refresh time when panning/zooming is very high.

Is there something I can do to speed up the operations? I also would like to keep it as "compact" as possible (avoid to break data into multiple views by category, the only alternative I thought of).


texts: Safely replacing ArcSDE feature classes using transactions in ArcObjects?

texts:
I am using C#, ArcGIS Server 10, MS-SQL Server 2008, and I want to replace an ArcGIS SDE feature class inside of a transaction. After the deletion of the feature class, I have to create a new one with the same name (as a copy of another feature class).

After some research I have implemented the following:

    public void deleteFeatureClassInTransactionTest(string featureClassName) {
        using (WebObject webObj = new WebObject())
        {
            // Connect to server
            AGSServerConnection agsServerConnection = ConnectServer();
            agsServerConnection.Connect();

            // Initialize workspace
            IWorkspace workspace = GetSdeWorkspace();
            webObj.ManageLifetime(workspace);
            // Initialize edit workspace
            IWorkspaceEdit workspaceEdit = (IWorkspaceEdit)workspace;
            IMultiuserWorkspaceEdit muWorkspaceEdit = (IMultiuserWorkspaceEdit)workspace;
            // Initialize feature workspace
            IFeatureWorkspace featureWorkspace = (IFeatureWorkspace)workspace;
            // Open feature class
            IFeatureClass featureClass = featureWorkspace.OpenFeatureClass(featureClassName);
            webObj.ManageLifetime(featureClass);

            // Start transaction
            muWorkspaceEdit.StartMultiuserEditing(esriMultiuserEditSessionMode.esriMESMNonVersioned);
            workspaceEdit.StartEditOperation();
            try
            {
                // Delete feature class
                ESRI.ArcGIS.Geodatabase.IDataset pdataset;
                pdataset = (ESRI.ArcGIS.Geodatabase.IDataset)featureClass;
                pdataset.Delete();

                // ... CREATE THE NEW FEATURE CLASS ...
                // Simulate exception
                throw new Exception("DEBUG");

                // Commit
                workspaceEdit.StopEditOperation();
                workspaceEdit.StopEditing(true);
            }
            catch (System.Exception ex)
            {
                // Rollback
                workspaceEdit.AbortEditOperation();
                workspaceEdit.StopEditing(false);
                throw ex;
            }
            finally ..

However, the rollback does not work. The table gets deleted immediately.

How can I make this whole process atomic?


texts:
I don't see the reason why multiversioned views exist.

It's not possible to use joins like in normal views. And editing versions is possible by directly using the feature class itself. So when am I supposed to use them?

So if I understand:

When I query the feature class, I query the compressed DB. And the only way to access data in the versions, is using mv views.
Am I wrong?


texts: Modifying projection in Arc-changing standard parallel in N. Pole Lambert Azimuthal

texts:
I'm working in a few targeted areas, and I wanted to project them all similarly. ROIs are on the order of 200 km x 200 km. For various reasons I'd like to use an equal area projection.

I like the Lambert Azimuthal Equal Area projection, but the only preset Lambert Azimuthal Equal Area projection in Arc is the North Pole Lambert Azimuthal Equal Area.

Are there any drawbacks to simply modifying the standard parallel and meridian of this projection for each target area?


texts: Automatically selecting main rivers of very branched river network using Arc Hydro?

texts:
I'm using ArcGIS 10.1 for Desktop and Arc Hydro.

Is there any way to automatically select the main rivers of a very branched river network?


texts:
I have a shapelayer with lines (streets), each line has a value for the traffic on it. The value is shown by the width of each line. Now I want to smooth the transition between one segment that has e.g. the value 500 which goes into a segment with value 100. Any ideas?
I want to achieve more or less the same result, as in this post, only automatically for a larger dataset …
Blend differing line thicknesses in QGIS?

Maybe a little graphic helps: Above what I got. Bottom, what I want to have. Each line segment has only 1 value, and I want to smooth the segments into each other.


texts:
I have loaded GPS points into a QGIS vector layer.
I also imported a vector line from a dxf file which doesn't have georeferencing.

I would like to fit in line geometry between the GPS points - how can I do this?


texts: Calculating change between two sets of district lines?

texts:
I'm using the Census Shape files for state legislative district lines for 2008 and 2012. I want to create a measure that indicates how much these lines have changed between these two time periods (because of redistricting). At first, I thought the best approach was to calculate the percentage of overlap between the two different layers by each district. But, I'm not sure if that's the appropriate way to do it.

Is there a better measure that tells me how similar (or different), a 2012 district is to its 2008 counterpart?


texts: Using Python to label points with format dependent on date of measurement in ArcMap

texts:
I have a database of wells and water level measurements which I would like to label in a stacked manner with the ID on top and the measurement below it. If the measurement is from 2020, it should appear normally but if it is from before 2020 I want to put the value in parentheses. It seems what I have written so far ignores the else part of the function (everything comes out normal). If I switch the >/< symbol, everything comes out in parentheses.

Does anyone know what I am doing wrong in ArcMap here?

def FindLabel ( [IDMINE], [Date], [Q1_2020_WL] ):
    if [Date] > 12/30/2019:
        return [IDMINE] + '\n' + [Q1_2020_WL]
    else:
        return [IDMINE] + '\n' + "(" + [Q1_2020_WL] + ")"


texts: I am making a map where I want to label two points with only one label. How can I do that? I've attached the image of my map.


texts: Getting NameError when using Cursors to read vertex coordinates in ArcPy?

texts:
I am new to python and having trouble Looping each feature in the new polygon shapefile, then use geometry cursors to read vertex coordinates of this lake polygon, identify the city that is the closest to the lake’s vertex, and add the city information into the fields of this lake feature.

for lake in arcpy.da.SearchCursor(In_Lakes,["OID@","SHAPE@"]):
    for feature in lake[1]:
        for vertex in feature:
        print vertex, vertex
Near_analysis(lake, In_Cities)

The error I get is:

  NameError: name 'Near_analysis' is not defined


texts: GeoServer using GeoWebcache and SQL Server whats the process?

texts:
I'm using GeoWebCache and GeoServer to tile cache a layer of data from a SQL Server database. If I seed the cache and I send a tiled WMS request what occurs?

Does GeoServer make a call to SQL Server or are the tiles taken straight from the GeoWebCache?
If the underlying SQL server data changes how does GeoServer know to regenerate the relevant tile?

Any help apprciated.


texts: Can the JSON output for an ArcGIS Rest Server be opened in a similar way to a WFSServer?

texts:
I would like to construct a python script to obtain JSON vector data from the following ArcGIS Rest Server:

http://maps.six.nsw.gov.au/arcgis/rest/services/sixmaps/LPIMap/MapServer/69

For a WFSServer output for an ArcGIS Rest Server hosted by a different organisation, I am able to open the layer and view information about the data using the following Python code:

import ogr

url = 'WFS:http://idpgis.ncep.noaa.gov/arcgis/services/NWS_Observations/ahps_riv_gauges/MapServer/WFSServer'
ds = ogr.Open(url)
lyr = ds.GetLayer(1)
name = lyr.GetName()
print (name)

Is it possible to "open" the url of the JSON output of an ArcGIS Rest Server in the same way?  For the NSW six maps ArcGIS Rest Server the following python code returns a "OSError: Could not open " error:

import ogr

url = 'http://maps.six.nsw.gov.au/arcgis/rest/services/sixmaps/LPIMap/MapServer/69/JSON'

ds = ogr.Open(url)

The following variations of the url also return a "OSError: Could not open " error:

url = 'JSON:http://maps.six.nsw.gov.au/arcgis/rest/services/sixmaps/LPIMap/MapServer/69/JSONServer'
url = 'http://maps.six.nsw.gov.au/arcgis/rest/services/sixmaps/LPIMap/MapServer/JSONServer'
url = 'JSON:http://maps.six.nsw.gov.au/arcgis/rest/services/sixmaps/LPIMap/MapServer/JSONServer'


texts: Get Path of my toolboxes ESRI addins Python 2.7

texts:
I need my toolboxes path. These are in the installation folder (toolboxes) but when using the function:

dir_tools = os.path.join(os.path.dirname(__file__),r'toolboxs\parametros.tbx')

  result:
  C:\Users\xxx\AppData\Local\ESRI\Desktop10.4\AssemblyCache{8F89910B-C991-4EAA-B2A5-294D1F72C1A4}\toolboxs\parametros.tbx

Return path of system, not path installation folder.

How do I get the path of the installation folder?


texts: Why OSM doesn't provide administrative division for national parks?

texts:
I noticed that for many national parks OSM doesn't provide administrative division (state, province) where they belong to (only country). Some examples:

https://nominatim.openstreetmap.org/details.php?place_id=2597516462
https://nominatim.openstreetmap.org/details.php?place_id=145053908

Meanwhile for some others it is specified:

https://nominatim.openstreetmap.org/details.php?place_id=145091198

Any idea why? And how can I figure out where they belong to using OSM data?


texts: Including chart based on fields in CartoDB infowindow?

texts: I would like to make graphs, pie charts, etc based on fields in my table within my infowindows. Is there a way to accomplish this?


texts: How to use image_url from dataset in Carto Builder HOVER pop-up window?

texts:
I am creating a map in Carto Builder with several points. I want to set the pop-up (hover) to display information from my dataset, including an image which is dictated by the corresponding image_url column in my dataset.

Below is the new HTML for the hover popup window. How do I get the image to actually appear? No matter what I do or what small syntax changes I make I just cannot get the image to show up!

      >

name

{{name}}

university

{{university}}

class_year

{{class_year}}


texts: issues with updating data source, script runs but nothing happens, arcpy

texts:
I have seen many questions posted that are similar to mine but still havent found the answer so figured I would post my script. I want to change all the data sources found in all mxds in a location to point to one geodatabase. The mxds are currently referencing feature classes, shapefiles and rasters. I am wondering if there is an issue with the shapefiles that have now been moved to the gdb and no longer have the file extension.. even if so why wouldn't the feature classes sources change? Here is what I have, it runs, but there are no changes to the data sources that I can see.

UPDATE: I've now modified it to use lyr.ReplaceDataSource, which works it changing the data source for the first shapefile it encounters in the mxds, but it doesn't change the source for the other feature classes and rasters.

UPDATE 2: Im still struggling to find the answer here. I now am trying to use findandreplaceworkspacepaths on each mxd, and again the script runs, but the data sources don't change. I am printng out the old and new paths in the script and they are printing correctly.. so what is happening?? updated code below:

import arcpy, os
mxd_path = r'C:\GIS\Assimilate_Mxd_Data\Assimilate_Mxd_Data_TESTING\Input Mxd'
new_path = r'C:\GIS\Assimilate_Mxd_Data\Assimilate_Mxd_Data_TESTING\Output Data\test_output.gdb'
mxd_list = []

for dirpath, dirnames, filenames in os.walk(mxd_path):
    for filename in filenames:
        fullPath = os.path.join(dirpath, filename)
        basename, extension = os.path.splitext(filename)
        if extension == ".mxd":
            each_mxd = fullPath
            mxd_list.append(each_mxd)
for m in mxd_list:
    mxd = arcpy.mapping.MapDocument(m)
    for lyr in arcpy.mapping.ListLayers(mxd):
        if lyr.supports("DATASOURCE"):
            old_path = lyr.dataSource
            mxd.findAndReplaceWorkspacePaths(old_path,new_path,validate=True)
            print old_path,new_path
mxd.save()
del mxd


texts: Leaflet API - select one overlay at a time like base layers

texts:
I have a set of basemaps and overlays on my simple web mapping application with Leaflet. I would like to be able to select one overlay at a time like base layers so that multiple selection of overlays is not allowed.

I cannot make those overlays into base layers as at least one base layer should be visible on the map at any time.

Currently I have the following code snippet.

var baseLayers = {baselayer1, baselayer2, baselayer3};
var overlayMaps = {overlay1, overlay2, overlay3};
L.control.layers(baseLayers, overlayMaps).addTo(map);


texts: Pass additional options to GDAL in QGIS v3

texts:
I'm using gdalwarp inside a QGis v3 model (using qgismodeler) and I'd like to pass the -tap option to the command.

I can't find where to place this. Is it possible to add additional options to gdal command


texts: How to change the symbol and/or geometry in ArcGIS API for JavaScript 4?

texts:
The below code works to change the Graphic Symbol and Geometry in ArcGIS 3.18 JavaScript API.

        /**********************
        * Create a point graphic
        **********************/
        var pointX = new Point({
            longitude: location.long,
            latitude: location.lat
        });

        var templateX = new InfoTemplate({
        title: "${Name}",
        content: "Id: ${Id}"
        });

        var pointGraphicX = new Graphic(
            pointX,
            markerSymbolGreen, {
                Id: deviceId,
                Name: name,
            }, templateX);

        assets.graphic = pointGraphicX;
        map.graphics.add( pointGraphicX);

        ....
        ....
        if (assets.graphic != null){
            if(data.sensors.temp > 25) {
                assets.graphic.setSymbol(markerSymbolRed);
            } else {
                assets.graphic.setSymbol(markerSymbolGreen);
            }
        }

        ....
        ....
        var newPoint = new Point({
            longitude: newlocation.long,
            latitude: newlocation.lat
        });
        assets.graphic.setGeometry(newPoint);

But setSymbol() or setGeometry() is not a method in the ArcGIS 4.2 JavaScript API Class.

According to https://developers.arcgis.com/javascript/latest/guide/working-with-props/index.html, I've tried

assets.graphic.set(markerSymbolRed)
assets.graphic.set(newPoint)
assets.graphic.symbol = markerSymbolRed

But it seems it's not working.

How do I change the Symbol  and Geometry of existing graphic Symbol in ArcGIS 4.2 JavaScript?


texts: GPS points of a vehicle snapped to a specific polyline path

texts:
So I have been capturing bus GPS locations for a while now, and I am became interested in snapping them to their corresponding polyline (route). While at first this seemed like an easy enough task I began to realise how complex this task is, particularly in situations where the route doubles back on itself and a simple snap to closest isn't always accurate.

Now I assume that I am not the first person in the world of GIS to want to do something like this, I was wondering if you have come across any algorithms, scripts, functions, programs or whatever that could help me with this task. I am currently using PostGIS and QGiS but am open to any suggestions.


texts: How to define inputs that are the result of an iteration? Python ArcGIS

texts:
The script (helped by @BERA and badly added to by me *completely new to coding) iterates through the fields in a polygon attribute table to convert to a temporary raster (for different scenarios in the cost analysis) and then do a cost analysis on each raster.  I need the resulting accumulated cost surface to be the input of the next part of the analysis – namely set null values and the resulting raster to be the input for summarize zones  - and repeated for each cost surface.

How do I define the inputs and outputs please (of null value and summarise zones)?

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

 # Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Set overwrite option
arcpy.env.overwriteOutput = True

# Local variables:
habitat = "D:\gis_data\habitat.shp"
fodyrange_current = "D:\gis_data\fodyrange_current.shp"
AQ = "D:\gis_data\AQ.shp"
R_fody_exp1= ' D:\gis_data\model_outputs\raster_f_1'
backlink_fmaq = ""
nullc_f1 = "D:\gis_data\model_outputs\nullc_f1"
zone_fody_exp1 = "D:\gis_data\model_outputs\zone_fody_exp1"

fieldlist=['fody_exp1','fody_exp2','fody_exp3','fody_exp4'] #extend/edit to include all your fields

for item in fieldlist: #For each field specified in fieldlist perform all below:
arcpy.FeatureToRaster_conversion(in_features=habitat, field=item, out_raster=r'in_memory\tempraster')
#Dont know what you want to do next, cost distance?:
outCostDistance = CostDistance(AQ, r'in_memory\tempraster')
outCostDistance.save(os.path.join(r'D:\gis_data\model_outputs','R_'+item)) #The outputs will be named after field, for example R_fody_exp1
#and so on...

# I don’t know how to set the input in the following so that it iterates through all the outputs of the cost analysis

# Process: Set Null
arcpy.gp.SetNull_sa(R_fody_exp1, R_fody_exp1, nullc_f1, "\"value\" = 0")

# Process: Zonal Statistics as Table
arcpy.gp.ZonalStatisticsAsTable_sa(AQ, "id", nullc_f1, zone_f_1, "DATA", "ALL")

This is the ModelBuilder model exported as python (sorry about long pathnames)- but does not include any iterations.  It is necessary to iterate through the fields in an an attribute table to get the values to be used as cell values when converting to raster in the first step in this process and then use the outputs in the rest of the process.

# Import arcpy module
import arcpy

# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")

# Set Geoprocessing environments
arcpy.env.scratchWorkspace = "C:\\Users\\willow\\Documents\\ArcGIS\\Default.gdb"
arcpy.env.snapRaster = "scores_fody_m"
arcpy.env.extent = "534630.586094391 7813776.96082858 552599.503136767 7825252.24043838"
arcpy.env.workspace = "C:\\Users\\willow\\Documents\\ArcGIS\\Default.gdb"

# Local variables:
habitat = "habitat"
fodyrange_current = "fodyrange_current"
AQ = "AQ"
hab_sn_1 =
"D:\\gis_data\\rodrigues\\Rod_analysis\\iterate_model\\hab_ras\\hab_sn1"
dist_sn_1 =

"D:\gis_data\rodrigues\Rod_analysis\iterate_model\dis_raster\dist_sn_1"
    backlink_fmaq = ""
    nullcost_sn_1 =

"D:\gis_data\rodrigues\Rod_analysis\iterate_model\setnull_cost_r\nullcost_sn_1"
    zone_sn_1 =

"D:\gis_data\rodrigues\Rod_analysis\iterate_model\zonal_tables\zone_sn_1"

# Process: Polygon to Raster
arcpy.PolygonToRaster_conversion(habitat, "warb_mean", hab_sn_1, "CELL_CENTER", "NONE", "20")

# Process: Cost Distance
arcpy.gp.CostDistance_sa(fodyrange_current, hab_sn_1, dist_sn_1, "", backlink_fmaq)

# Process: Set Null
arcpy.gp.SetNull_sa(dist_sn_1, dist_sn_1, nullcost_sn_1, "\"value\" = 0")

# Process: Zonal Statistics as Table
arcpy.gp.ZonalStatisticsAsTable_sa(AQ, "id", nullcost_sn_1, zone_sn_1, "DATA", "ALL")