texts: Improving processing time of ArcPy script?

texts:
Background: I work for an electrical utility, and we have a point feature class called "meters" that contains approximately 35,000 records. These records represent the actual power meter, and contain the customer's account number among other things.

Approximately 2000 of these records are missing account numbers, and I've been tasked with filling in this information which is required for another project the utility is working on.

I selected out the meters that were missing accounts, and created a feature class containing only those problem meters. I then filled in the missing information using a separate script that I wrote.

Now, I'm trying to take the completed account numbers from the subset feature class and transfer them into the "live" Meter feature class. I wrote a fairly simple script to accomplish this, however it seems to be taking a very long time to run. It ran all night last night without finishing (over 12 hours long) and stopped responding after my PC's connection to the server was booted for backup/synch. So this morning I created a local (non SDE) copy of the Meters feature class and am running the script against that to see if it makes a difference. It's been running now for three hours and I'm feeling like this is taking too long.

I'm hoping someone might be able to take a look at the code below and give me some tips on how to possibly improve the processing time.

import arcpy, time, csv
from arcpy import env
import datetime as dt

logfile = r"U:\Python Scripts\scriptLog.csv"
scriptName = "Meters_AttributeTransfer"
startTime = time.clock()
clockTime = time.time()

workspace = arcpy.GetParameterAsText(0)
arcpy.env.overwriteOutput = True

targetMeters = arcpy.GetParameterAsText(1)
goodMeters = arcpy.GetParameterAsText(2)

tMeterFields = ["LOC_NBR"]
gMeterFields = ["LOC_NBR"]

arcpy.MakeFeatureLayer_management(targetMeters, "tMeters_lyr")
arcpy.MakeFeatureLayer_management(goodMeters, "gMeters_lyr")

with arcpy.da.UpdateCursor(targetMeters, tMeterFields) as uCursor:
    for uRow in uCursor:
        arcpy.SelectLayerByLocation_management("gMeters_lyr", "INTERSECT",     targetMeters, None, "NEW_SELECTION")
        with arcpy.da.SearchCursor("gMeters_lyr", gMeterFields) as bCursor:
            for bRow in bCursor:
                if uRow[0] == bRow[0]:
                    pass
                elif uRow[0] != bRow[0]:
                    uRow[0] = bRow[0]
                    uCursor.updateRow(uRow)

#----Write Log File---------------------------------------------------------------------------------------------------------------
duration = ((time.clock() - startTime))/60
dateObj = dt.datetime.now()
dateStr = dateObj.strftime( '%Y-%m-%d' )
newRow = [scriptName, str(duration), dateStr, clockTime, "V0.1"]
with open(logfile, 'a') as edit:
    writer = csv.writer(edit)
    writer.writerow(newRow)

The idea for this script is fairly simple:

An update cursor is created against the "target meters", and for each row in the target meter feature class, the matching/intersecting point in the "good meters" feature class (which contains the missing account information) is selected and searched to see if there is a matching number in the target feature class. If there is no match, then it copies the account number from the good meter to the target meter and moves on to the next one.

There is also code included which updates a log file that I maintain, which tracks when I run scripts and how long they take, and which version of the script I'm running.

I hope someone has some insight into this.

PS: I'm using ArcGIS 10.2.1 with Python 2.7.


texts: Select neighboring cells from a raster under conditions

texts:
I generated river networks (raster format) from a DEM and I also want to select the neighboring cells of the river networks from the original DEM, if the neighboring cells have a lower or equal elevation. Is there any way to achieve this with ArcGis 10?

In this plot, a branch of river network (the blue line) crosses cells no.2, no.5 and no.8. Let's pick no.5 as example. All the cells except 5 itself are neighboring cells of no.5. I want to filter cells that has an elevation no higher than cell no.5 from the neighboring cells. Assuming that the number on the plot denotes elevation, the cells I need would be cells no.1, no.3 and no.4.


texts: Adding values from list to field in feature class using ArcPy?

texts:
I have the code below. fc is a shapefile. I add a new field called "RANK". I want to populate the RANK column with values from Final_Rank list. Fortunately, they are of the same size. Field RANK has got 11 rows and Final_Rank has 11 values. The code below creates the field RANK and then adds 11 additional rows, below the existing 11 rows in the attribute table and inserts the values to RANK beginning notfrom the start but from the 11 additional rows. I want the code to add the values to the existing RANK field without adding any rows.

print Final_Rank

arcpy.AddField_management(fc, "RANK", "Double")

with arcpy.da.InsertCursor(fc, 'RANK') as cursor:
    for i in Final_Rank:
    cursor.insertRow(i)


texts: How to install older GDAL v1.11 on Ubuntu 16.04 instead of v2.1.3

texts:
I have another question going on this forum regarding GDAL 2.1 ogrinfo not working at Trying to get GDAL 2.x to work, but ogrinfo reports ERROR 3: and ERROR 4: and "Failed to open table res_t48.PAT", but in parallel, as a fix to an apparently "broken" ogrinfo 2.1 reading coverages, I'd like to see if I can just go back to gdal v1.11 on Ubuntu 16.04.

Can anyone tell me how to install an older version of gdal on Ubuntu 16.04 when "apt-get install gdal-bin" installs v2.1.3?  What repositories do I have to point to, or what packages do I install?


texts:
I created a virtual field in my layer where I want a specific value to be calculated. In my attribute table I have several columns. The expression I am trying to create is the following. I want to iterate through all features and add all values of the column age who have a specific value, let's say the name of the features should be marc. In my first attempt I created this:

CASE WHEN  "name"  = "marc" THEN  sum(  "age" ) END

But it kinda doesn't work and I don't know how to iterate though all features in the layer and only choose those features with the name marc. Any advice?


texts: How to add a layer with city names into QGIS project?

texts:
I am new to GIS technologies and I wonder how to add a layer with city names to the QGIS project. I have two layers: (1) Base cartography that is a map without any labels, (2) Delimited text layer that shows some points on the map.

I assume that I should add another layer, but I don't know where can I find appropriate data. Which format should I search for? SHP or any other?

To be more concrete, I need this information for Spain.


texts: External GeoJson file integration and showing a popup

texts:
I have found a Github repository which contain sample data set to include a external geojson file and showing a popup when click on an feature.

BUt I dont know what I am doing wrong.

I have found similar questions regarding this issue in gis.stackexchange but I couldn't find a solution for my problem in those answers.
https://stackoverflow.com/questions/34751488/openlayers-3-popup-with-external-geojson

Load external geojson file into leaflet map

This is the code I have used(Which included in github repository). I only changed the mapbox token.

<!DOCTYPE html>

    <meta charset="utf-8">
    <meta http-equiv="X-UA-Compatible" content="IE=edge,chrome=1">
    <meta name="viewport" content="width=1024, user-scalable=no">

    <link rel="stylesheet" href="leaflet.css" />

    <title>Leaflet AJAX</title>

When I open the index.html It shows only the basemaps (No any other map layers)

Base in the code output should be like this


texts: OGC Filter set but symbolizer expected in GeoServer SLD

texts:
I have a problem with one of my SLD files and do not know what is causing the problem.
The error message is pretty clear: line 18: cvc-complex-type.2.4.a: Invalid content was found starting with element 'ogc:Filter'. One of '{"http://www.opengis.net/sld":Symbolizer}' is expected. However, when I look into the example of the cookbook for attribute based
or for zoom based polygons then I cannot figure out why I am getting the error. Have a look at the code below... I don't know why a symbolizer is expected in line 18, even though the example in the cookbook has the structure the same way as I do, or not!

Any suggestions on how to overcome the problem?

I am running Geoserver in version 2.3.2

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0"
 xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
 xmlns="http://www.opengis.net/sld"
 xmlns:ogc="http://www.opengis.net/ogc"
 xmlns:xlink="http://www.w3.org/1999/xlink"
 xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <NamedLayer>
    <Name>Baggrund NO2</Name>
    <UserStyle>
      <Title>Baggrund NO2</Title>
      <FeatureTypeStyle>
        <Rule>
          <Name>4,75 - 5,31</Name>
          <Title>4,75 - 5,31</Title>
          <MinScaleDenominator>23000</MinScaleDenominator>
          <MaxScaleDenominator>1000000000</MaxScaleDenominator>
          <ogc:Filter>
            <ogc:PropertyIsEqualTo>
              <ogc:PropertyName>Type</ogc:PropertyName>
              <ogc:Literal>1.0000</ogc:Literal>
            </ogc:PropertyIsEqualTo>
          </ogc:Filter>
          <PolygonSymbolizer>
            <Fill>
              <CssParameter name="fill">#1a9641</CssParameter>
              <CssParameter name="fill-opacity">0.40</CssParameter>
            </Fill>
          </PolygonSymbolizer>
        </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>


texts:
I'm trying to create elevation graphs from a set of line features and raster elevation file, while the Stack Profile tool from ArcMap can do that, I need to modify the min/max values of the y-axis to make them uniform. I know I can manually change the values but I need to be able to automate the process of changing those values.

I have looked all over to find if there's way to do that, and the closest answer I got was the addSeriesBarMinMax function, which does not work for my graph. This is a test run of the code:

import arcpy
from arcpy import env

env.workspace="D:"

graph="test1.grf"
output_graph="test.grf"

graph=arcpy.Graph()

graph.addSeriesBarMinMax(graph, 400, 510)
arcpy.SaveGraph_management(graph,output_graph)

And I get this error:

  arcpy.SaveGraph_management(graph,output_graph)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\management.py", line 5473, in SaveGraph
    raise e
RuntimeError: Object: Error in executing tool

Is there any way for me to modify the min/max y-values from ArcPy?


texts:
GIS: Geomedia 2015, Q-Gis 2.14

Avaible Data: Raster files (tif) with a seperate saved georeference in Geomedia mdb, Shapefiles with the outline of every tif file.

I'm using Geomedia 2015 and have more than 50k small Tiff files to export.
There is an option to do this in Geomedia, but only for one image at once.

I know the Coordinate System/Datum and also the 4 node coordinates for every tiff.

So, does anybody know a way to batch that process, either in Geomedia or with QGIS?


texts: FeatureCollection.geometry() shortcut no longer working to return a single geometry from the features of a FeatureCollection

texts:
Up until sometime in the last week or so, the following code was working to convert the features of a filtered FeatureCollection into a single multipolygon geometry:
var predicates = ee.Filter.inList('GEOID', countylist);
var unioncounties = ee.FeatureCollection(counties.filter(predicates).geometry());

Now however, the features all remain separate and I have to use the following to get a single geometry, which is much slower:
var predicates = ee.Filter.inList('GEOID', countylist);
var unioncounties = ee.FeatureCollection(counties.filter(predicates).geometry().dissolve());

Anybody else seeing this change of behavior?
UPDATED TO EXCLUDE COUNTIES WITH LINESTRINGS from Noel's answer below.  Still,  .geometry() no longer is merging the features into a single geometry like it did a month ago.  It does not matter what FC is used, the behavior has somehow changed:
Full working code example here: https://code.earthengine.google.com/6d3c7110f0994cfde2a467cd9ee3911b


texts:
I have two raster layers.

(A) is state-sized, 30m.

(B) is county-sized, 3m.

I'd like to get a new raster (C) which is the part of (A) which coincides with (B), preferably scaled to the same resolution as (B).

Is there a way to do this from the command line? From Python or C++?


texts: How to create report using arcpy in arcgis 10.2 (result the same as using view -> Reports -> Create Report)

texts:
I need to create reports in arcgis 10.2. because I am not intent to request user to go through all the steps for generating reports, so I want to create an add-in tool to do this. The result is the same as the picture, Can it be realized?
Thanks in advance.


texts: Displaying progress bar when using QGIS processing tool

texts:
I'm building a QGIS plugin and I want to improve the user interface. I'm able to use processing tools like buffer, but with heavy layers I want to display a progress bar, users can easily see the treatment progression.

In documentation, there is some words about it but when I run simples examples, QGIS systematicaly not responding and there isn't any progress bar. Note that this treatment works perfectly without this progress bar issue.

Here is my code :

import time
from qgis.PyQt.QtWidgets import QProgressBar
from qgis.PyQt.QtCore import *

# Progress bar
progressMessageBar = iface.messageBar().createMessage("Work in progress...")
progress = QProgressBar()
progress.setAlignment(Qt.AlignLeft|Qt.AlignVCenter)
progressMessageBar.layout().addWidget(progress)
iface.messageBar().pushWidget(progressMessageBar, Qgis.Info)

# Layer
vlayer = qgis.utils.iface.activeLayer()
count = vlayer.featureCount()
features = vlayer.getFeatures()

# Progression loop
for i, feature in enumerate(features):
    percent = i / float(count) * 100
    # My treatment
    processing.runAndLoadResults('native:buffer', {"INPUT": vlayer, "DISTANCE": 20, "OUTPUT": 'memory:buffer'})
    # Progression status
    iface.statusBarIface().showMessage("Processed {} %".format(int(percent)))
    progress.setValue(percent)

iface.statusBarIface().clearMessage()

Am I missing something ? Maybe something is wrong in my loop. I'm working with QGIS 3.10.1.


texts: Combining Geological Grids and extruding between them?

texts:
I have a situation where I have 3 Geological grids, each showing a different horizon + the Seabed Grid.

I need to extrude between these grids in order to make the 3D polygon for the geological unit which is between these horizons. The problem occurs however when (due to Shallow gas or lower geological units protruding upwards) there are gaps in my grid; sometimes to the extent that one of the horizons only covers a small segment of the sea bed area at the surface. The solution therefore (I think) is to combine say, H3 + H2 for example, and H2 +1, and then: Extrude between [(H3+H2) (where H2 fills in the gaps of H3)] and [(H2 +H1)(Where H2 Fills in the Gaps of H1)].

The Situation becomes very complicated when some of H3 actually rises to the Seabed and the Seabed fills in some gaps in H3. A very complex series of combinations of grids is therefore needed involving probably H3, H2, H1 and the Seabed.

The only way I have found to do this is in Global Mapper, but version I have does not have the combine grids tool, unlike the expensive newest version. The software I Have access to include Arc GIS 10.2 and freeware such as Map window, QGIS, GRASS etc. The solution may involve raster calculator conditional function, however I have no experience in how to accomplish this task using this.


texts: Unable to get correct @row_number in Attribute Table

texts: I am working with QGIS 3.10 and would like to create a field in the attribute table which is a unique ID for a feature.  I have read that using the @row_number as the expression should provide this capability, but the resulting value is always NULL (both in the Preview and once a feature is created).  If I select the field from the Attribute table and apply the @row_number in the expression box and select 'Update All' the table is correctly updated with the row numbers.  How can I get the correct row number to show up when the feature is selected and then displayed in the Attribute table?


texts: Error to save Custom coordinate system definition in QGIS Desktop 2.14.7

texts:
I'm trying to enter a coordinate reference system applicable to El Salvador in QGIS Desktop, but it doesn't save the coordinate system.
I am currently using version 2.14.7, however, I tried before using version 2.16.3., and the same error occurrs.

QGIS allows me to perform the test with the coordinate system. And then, I apply the changes. No error message was displayed, but when I back and looked for the new coordinate system, it wasn't found.


texts: Why pgRouting GUI panel may not work in QGIS?

texts:
What could be the reason behind the mal-functionality of the pgRouting Layer GUI panel?

In order to route shortest path via dijkstra function code I'm able to conduct the task using the following code:

SELECT seq, id1 AS node, id2 AS edge, cost, geom
  FROM pgr_dijkstra(
    'SELECT id, source, target, st_length(geom) as cost FROM schema.cvl',
    1, 100, false, false
  ) as di
  JOIN schema.cvl pt
  ON di.id2 = pt.id ;

but unable to do the same using the pgRouting Layer GUI panel!

  An error has occured while executing Python code: .............

I searched for this error on the net and found out that it may happen because the database password is not saved.. So, I saved the password, but still getting same error when I try to locate source_id or target_id on the map

I'm using QGIS 2.0.1 and pgRouting pg93-binaries-2.0.0


texts: Configure QGIS in a project on PyCharm on Ubuntu

texts:
I'm creating a PySide2 project and i'd like to know if it's possible to add QGIS in it.
When i puted "from qgis" it does not work.

from PySide2 import QtCore, QtGui, QtWidgets
from qgis._core import QgsVectorLayer, QgsMapLayer


texts: Adding geojson features to OpenMapTiles displayed with MapboxGLJS is not working

texts:
The OpenMapTiles example for display using Mapbox GL JS works nicely to display the map.  However, when trying to add lines and features to the map using Mapbox GL JS functions like this:

map.on('load', function () { map.addLayer({ .... }) })

the desired added lines and features do not display on the map.  The full example HTML is shown below.

Questions:

Is there something about the openmaptiles style osm-bright that makes it incompatible which the specified data as added below?
Could it be made to work by changing something in the code below?
Could it be made to work by changing something in the osm-bright style file?

The added code for geojson features are mapbox examples:

https://www.mapbox.com/mapbox-gl-js/example/geojson-markers/
https://www.mapbox.com/mapbox-gl-js/example/geojson-line/

The features should appear in and around San Francisco.

Here is the whole HTML file:

<!DOCTYPE html>

    <meta charset='utf-8' />
    <title>OpenMapTiles OSM Bright style</title>
    <meta name='viewport' content='initial-scale=1,maximum-scale=1,user-scalable=no' />

    <link href='https://api.tiles.mapbox.com/mapbox-gl-js/v0.29.0/mapbox-gl.css' rel='stylesheet' />


texts: Convert a multi-band Landsat RGB to a KML using MapTiler

texts: I was able to use MapTiler for converting a single raster (ortofoto GeoTIFF obtained with UAV) into a KML file. Now I want to convert a RGB Landsat into KML for visualisation in Google Earth. It is not clear to me how this can be done as the raster has three bands now. In addition, how I can assure that the nice visualisation of the LANDSAT in QGIS, obtained after setting manually min/max values with'Colourmap', is maintained during conversion MapTiler?


texts:
I've investigated the issue of coordinate order in GML and each time I think I've got it but I find my self confused once again.

I think some of my confusion may be to do with PostGIS.  Doing a simple query on my database I can see the coordinates are stored correctly.  The st_asgml function is outputting these coordinates in long,lat order which I think is correct for this srsName:

select title, ST_AsGML(3, locationgeometry, 15, 1), xmin(locationgeometry), ymin(locationgeometry) from tblobject where title='Chuska East Ridge';
-[ RECORD 1 ]-------------------------------------------------------------------------------
title    | Chuska East Ridge
st_asgml | <gml:Point srsName="urn:ogc:def:crs:EPSG::4326"><gml:pos srsDimension="2">-108.848780000000005 36.078707000000001</gml:pos></gml:Point>
xmin     | -108.84878
ymin     | 36.078707

However, if I take the same GML output above and run it through st_geomfromgml, the coordinates appear to have been reversed:

select xmin(st_geomfromgml('<gml:Point srsName="urn:ogc:def:crs:EPSG::4326"><gml:pos srsDimension="2">-108.848780000000005 36.078707000000001</gml:pos></gml:Point>'));
-[ RECORD 1 ]---
xmin | 36.078707

Am I missing something, or is this a bug in PostGIS?  Can someone please confirm for me once and for all whether I should have lat,lon or lon,lat in GML when using the srsName urn:ogc:def:crs:EPSG::4326?

PostgreSQL version: 9.1.14
PostGIS version   : POSTGIS="1.5.3" GEOS="3.2.2-CAPI-1.6.2" PROJ="Rel. 4.7.1, 23 September 2009" LIBXML="2.7.8" USE_STATS


texts: how to build a reverse geocoding service

texts:
We want to build a web service like google reverse geocoding.

Of course, we will only provide the service for only a specified region.

Then we want to know things we need to implement it.

By now, we have the pois,roads data of the region with ESRI Shapefile data.

Any suggestions?


texts: Adding fields to shapefile from standalone PyQGIS script?

texts:
When I run the following script via OSGEO4W64 Shell, I receive the following error when trying to run the qgis:fieldcalculator algorithm:

  'NoneType' object has no attribute 'mapCanvas'

I read that the error message relates to the script trying to access an interface although I'm not sure. The algorithm does what it's supposed to with the output layer containing a new ID column, but the layer has no features saved.

Basically, I would like to add a new ID field using the formula $rownum and remove the longitude and latitude fields which are generated by the Create grid algorithm. I would also be happy if this can be achieved without having to use the Field calculator algorithm (ie. by directly accessing the data).

import os, sys, glob

from qgis.core import *
from PyQt4.QtGui import *
app = QApplication([])
QgsApplication.setPrefixPath("C:\\OSGeo4W64\\apps\\qgis", True)
QgsApplication.initQgis()

from os.path import expanduser
home = expanduser("~")

#   Folder path of the Results for shapefiles
path_dir = home + "\Desktop\Test\\"
path_res = path_dir + "Results\\"

# Prepare processing framework
sys.path.append( home + '\.qgis2\python\plugins' )
from processing.core.Processing import Processing
Processing.initialize()
from processing.tools import *

def run():

    outputs_1=general.runalg("qgis:creategrid", 1000, 1000, 24108, 18351.157175, 258293.802316, 665638.226408, 1, 'EPSG:7405', None)
    outputs_2=general.runalg("qgis:fieldcalculator", outputs_1['SAVENAME'], 'ID', 1, 10, 0, True, '$rownum', None)
    outputs_3=processing.runalg("qgis:deletecolumn", outputs_2['OUTPUT_LAYER'], 'longitude', None)
    outputs_4=processing.runalg("qgis:deletecolumn", outputs_3['SAVENAME'], 'latitude', path_res  + "/" + "grid.shp")

#   Paths of the shapefiles in the Result folder with list comprehension
    output = [shp for shp in glob.glob(path_res + "*.shp")]

run()
QgsApplication.exitQgis()


texts: GRASS 7 tools in QGIS results in EPSG:3397 (PD/83 GK4) instead of 31468 (DHDN GK4) causing dislocation

texts:
GRASS 7 tools used from QGIS produce a result in a projection differing from the input projection which might lead to inconsistences:

Example v.net.connect, input EPSG of both layers is 31468 (DHDN GK4):

Result is in EPSG:3397 (PD/83 GK4):

Same with v.clean.

The resulting network is located ~1.5m away from the original, which is of course unacceptable (blue: original network, red: result after v.clean rsp. v.net.connect):

Why this unauthorized change of coordinate system and how to force GRASS 7 tools executed from QGIS (2.18.16) to output results in the same CRS as the inputs?

UPDATE

Trying the same in QGIS 2.99.0-280 leads to the same change in CRS, but does not cause erroneous dislocation of the data:


texts: Using Access Hyperlink Field Type in ArcGIS

texts:
I've been searching for an answer to my particular issue, but haven't turned up any answers after a couple days of searching. I am using Access 2013 and ArcInfo 10.X

My issue is:

I have an Access 2013 database containing land transaction records in a hyperlink field that is concatenated from several different fields to generate a file path. The hyperlink field displays text (i.e., a person's name) which hyperlinks to the appropriate folder. I have underlying information organized in the database where personal information broken up into several fields that auto-generates the file path to a folder that will contain PDFs.

Everything functions well in Access. However, when I establish the OLE DB connection in ArcCatalog, the field type is not recognized as a hyperlink field and in ArcMap (as expected) and is imported as a string.

My question is:

In ArcMap, can you import Access functionality, such as the hyperlink field generated from concatenated fields represented by text, and not the actual file path?

I have tried a couple workarounds, such as inputting the base URL in the MXD properties and designating a field to specify the exact folder. This works just fine, but I will eventually have several additional paths per parcel that I'll need to add for historical land transaction data. Ideally, a given tract of land will have several hyperlinked PDF's with the hyperlinks represented by the concatenated field (auto-generated from the Access form).

I'm under the impression this could probably be solved through a python script, but am at a loss of how to do it.


texts: How to make input data in MaxEnt with several occurences at the same coordinates?

texts:
I have insect traps data with several hundreds insects caught every day.

In this paper, authors say "Running the model using the number of mosquitoes per trap-night resulted in clustering of high probabilities of occurrence around individual collection sites."

How comes that they used the number of mosquitoes at a single individual trap location since so far what I understand from Maxent is that you only need three columns, occurence (species), x and y ?

How can I manage to tell MAxent to take into account the number of insects instead of their presence only ?

I am very new to MaxEnt and I really fail to resolve that.


texts: Add GeoServer WMS layer on web page using OpenLayers

texts:
I am trying to add GeoServer WMS Layer using OpenLayers on the web page.

And here is my code:

<%@ Page Language="C#" AutoEventWireup="true" CodeBehind="WMSLayer.aspx.cs" Inherits="OpenLayerMap.WebPages.WMSLayer" %>

<!DOCTYPE html>

    <link rel="stylesheet" href="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.2.1/css/ol.css" type="text/css">

    <title>OpenLayers WMS</title>

WMS Layer

I am getting the base map, but layer is not added on the basemap. Can somebody help me to solve this problem?


texts: QGIS3 Atlas: expression to pull values from a non-coverage layer at current location?

texts:
I have an atlas composition where I would like to, in addition to attributes from my coverage layer which consists of polygons that are scattered throughout the world, place other values from another layer.

I would like an expression in a text box that displays the name of the country that the current feature is in.


texts: Buffering vector polygons without projecting in R?

texts:
Nowadays, it is possible to buffer without projecting vectors in ArcMap and I want to do the same in R. So, could somebody suggest me some way outs to buffer vector polygons without projecting in R?

Function descriptions for "gbuffer" and "buffer" say that they can do that but ends up with error message that they expect metric CRS.


texts: Fitting polylines to their true position using ArcGIS Desktop?

texts:
I have a satellite image and CAD data of the same area. The polylines from CAD to not overlay correctly on the satellite image.  I want to fix polylines to their true coordinates but I could not fix them.

I'm trying spatial adjustment but they did not place truly. Is there anything I can do to correct this problem, or  could my CAD drawing be wrong?

In Spatial Adjustment, which adjustment type should I use ?

Affine
Projective
Similarity
Rubbersheet
Edge snap


texts:
I have been working with creating a story map on ArcGIS Online and I have been struggling to figure out how I can embed the webpage of my web scene to display a limited UI. When trying to embed the web scene into the story map, the display also includes my account logged in, and tools that I do not need.

I want to limit the UI so that the user of my story map gets only the map itself, and not additional tools, login options, and editing privileges.


texts: Blank values in attribute table don't work as null in QGIS?

texts:
I am using QGIS 3.4.

I have a layer generated from the QuickOSM plugin with various pieces of OSM data.

For lines 17-19 which do not have a value in the "name" attribute, I don't know why  null isn't showing up. I would like to select all the rows that have the name value empty. I tried using name is null, but that did not select anything. In settings, I have the option selected to say null in any blank values. My only guess is that maybe these values somehow are not empty, but they seem to be blank and I have tried writing something and then deleting it.

So what I would like to know is how I would select all rows of the table that don't have a value in the name column. Also, if someone could explain why NULL isn't showing up in the values that would be helpful.

I am pretty new to QGIS and GIS in general.


texts:
Possible Duplicate:
  Standard SLDs for OpenStreetMaps?

I can't find a good SLD for OSM.
I downloaded the OSM data and imported it to PostGIS. With GeoServer I use it as a data source to create my layers.

Is there no open source SLD for OSM out there?


texts: Why can't I make the blue postgis based multipolygons render like the black shapefile polygons in QGIS? Am I missing a button or render setting? I've compared properties on both in QGIS and they are identical simple fills with transparent centers and colored outlines.


texts:
I've been playing around with C++ and the API for GDAL to create some functions that I have done using Python. One of them was to create a virtual raster based on a list of files found in certain directory. I accomplished the later using the boost library. And when it comes the turn to write the virtual raster, I find that the vrt generated only contains one band.
I checked some questions regarding this problems, and the most close to my issue is found in the following link:

https://gis.stackexchange.com/questions/269976/c-gdal-api-add-more-than-one-rasters-to-vrt

However when compiling my program using
g++ -o stack_raster.exe stack_raster.cpp -lboost_system -lboost_filesystem -lgdal

I get the following message:

stack_raster.cpp:35:99: error: 'YES' was not declared in this scope ver->CreateCopy("media/B02_stack.vrt", poSrcDS,RALSE,APPEND_SUBDATASET=YES,NULL,NULL,NULL);

Before I was having no issues during the compilation by replacing APPEND_SUBDATASET=YES,NULL,NULL,NULL), for NULL. The resulted Virtual raster generated only contains one band, as expected. Could someone gives me a hint about my problem?

This is my code:

#include <iostream>
#include <boost/filesystem.hpp>
#include <boost/algorithm/string/predicate.hpp>
#include "gdal/gdal.h"
#include "gdal/cpl_conv.h" /* for CPLMalloc() */
#include "gdal/vrtdataset.h"
#include <iostream>
#include <string>
#include "gdal/cpl_string.h"
#include <stdbool.h>

using namespace std;
using namespace boost::filesystem;

int main()
{
    GDALAllRegister();

    path p("/home/roger/Documents/09_Stack_Raster/media");
    GDALDriver *poDriver = (GDALDriver*) GDALGetDriverByName( "VRT" );
    GDALDataset *poSrcDS, *poVRTDS;
    string bdir = "media/";

    for (auto i = directory_iterator(p); i != directory_iterator(); i++)
    {
        string filename = i -> path().filename().string();
        cout << filename << endl;
        if (!is_directory(i->path()) && boost::ends_with(filename, ".tif"))
        {
            cout << i->path().filename().string() << endl;
            std::string str_path = "media/" + i->path().filename().string();
            poSrcDS = (GDALDataset *) GDALOpenShared(str_path.c_str(), GA_ReadOnly);
            poVRTDS = poDriver->CreateCopy("media/B02_stack.vrt", poSrcDS,FALSE,APPEND_SUBDATASET=YES,NULL,NULL);
            GDALClose((GDALDatasetH) poSrcDS);
        }

        else
            continue;
    }

    GDALClose((GDALDatasetH) poVRTDS);
}


texts:
OgrGML and OracleGML format are not the same, I want to edit the parcel on QGIS which getting from OracleGML  format.

Is there any way to covert OgrGML to OracleGML and vise versa? (OracleGML format the geometry coordinates x and y separates by comma (,) but the OgrGML x,y coordinates separated by space.)


texts: Create dynamic layer for each unique value

texts:
I'm new to ArcPy.
I'm trying to write a script that will loop through the unique values in a field in a FC and create dynamic layers, with definition queries, for each unique value.

Each of the dynamic layers would all point to the same master feature class.
I do not want to output static feature classes. For example, I can't use the split by attributes tool for this.

I've attempted to adapt an ArcPy script from a related answer to do it:
from arcpy import *

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
fc = "WO_VW_FGDB"
field = "CLASSIFICATIONID"

env.overwriteOutput = True

#Create cursor to iterate rows
cursor = da.SearchCursor (fc, field)
for row in cursor:
    #sql statement for a single feature
    sql = '"{0}" = \'{1}\''.format (field, row[0])
    #Make layer with sql for one feature only. Use the field value as the layer name.
    MakeFeatureLayer_management (fc, row[0], sql)
    #Make mapping layer object
    lyr = mapping.Layer(row[0])
    #Add mapping layer object to map
    mapping.AddLayer (df, lyr)
del cursor
mxd.save ()

When I run the script in ArcGIS Desktop 10.7.1, it successfully creates the first layer. In fact, it actually creates duplicate layers.
But then it produces an error:
Runtime error  Traceback (most recent call last):   File "<string>", line 17, in <module>
File "c:\program files (x86)\arcgis\desktop10.7\arcpy\arcpy\management.py", line 6986, in MakeFeatureLayer
raise e ExecuteError: ERROR 000622: Failed to execute (Make Feature Layer). Parameters are not valid.
ERROR 000628: Cannot set input into parameter in_features.

What am I doing wrong?


texts: Is there a way to turn a data frame off in ArcMap 10.2?

texts:
I have an mxd with three data frames. I've written a python script that loops through a list of user-specified coordinates and scales for each data frame before exporting the mxd as a pdf.

Is there a way to "turn off" one of the data frames using python? (e.g. set data frame 1's visibility to False leaving just data frames 2 and 3 visible). The best idea I can come up with is to move data frames off screen when they're not needed and bring them back when they are.


texts: Removing black background when publishing Image Service?

texts:
I have a mosaic dataset that I want to publish as Image Service via ArcGIS Server.

I know how to set a background RGB value and display it as transparent in ArcMap.

Is there a way to do that when I publish the Image Service?


texts:
I am just trying to get the admin area out of FAO/GAUL/2015/level2 for Romania when I click in the map and use it as a label in a panel.

It works in the console but in the panel it prints out way too much. What am I doing wrong??

See code below:

// Create a FeatureCollection from the list and print it.
var fromList = ee.FeatureCollection("FAO/GAUL/2015/level2");

var filter = ee.Filter.inList('ADM0_NAME', ['Romania']);
var filteredArea = fromList.filter(filter);

//print(filteredArea);

Map.setCenter(25.1, 45.5, 6);
Map.addLayer(filteredArea, {}, 'Regions');

// Create a panel to hold the chart.
var panel = ui.Panel();
panel.style().set({
  width: '400px',
  position: 'bottom-right'
});
Map.add(panel);

// Register a function to draw a chart when a user clicks on the map.
Map.onClick(function(coords) {
  panel.clear();
  var label = ui.Label('This is a panel Example',{textAlign: 'center'});
  panel.add(label);
  var point = ee.Geometry.Point(coords.lon, coords.lat);
  var filteredPoint = filteredArea.filterBounds(point);
  print(filteredPoint.aggregate_first("ADM1_NAME"));
  var AreaName = filteredPoint.aggregate_first("ADM1_NAME");
  print(AreaName);
  var label = ui.Label(AreaName);
  panel.add(label);
});


texts:
Is there any REST API for remove wms cache of a layer in GeoServer?

I use gwc/rest/masstruncate but it returns status 500. I can clear cache of a layer manually in TileCashign in GeoServer but I need API for this.

my request is: curl -v -u user:pass -POST -H "Content-type: text/xml" -d "<truncateLayer><layerName>workspace:layer</layerName></truncateLayer>"  "http://localhost:8080/geoserver/gwc/rest/masstruncate/"

and my response is: Truncation failed


texts:
I am tracing routes from a road segments layer with the goal graphing the elevation profile of the route in JavaScript with GeoJSON. This is my "real world" hobby project to get familiar with GIS and JavaScript. Open to any suggestions.
After picking a route, from the processing toolbox I "Explode lines" then "Segmentize by maximum distance". I would like to export this with the Z coordinates, the closest thing I have found so far is "Join attributes from closest" where I select a layer created from a CSV with a grid of XYZ data. The result is close but I would like to have the Z value in the coordinates rather than the attributes.
All my data is from https://maps.canada.ca/czs/index-en.html


texts: How to make geoserver return JSON for WMS GetFeatureInfo?

texts: I wish to get the results of a WMS GetFeatureInfo request as a JSON instead of the table that comes up by default.  Out of the box, the only options for the response from a WMS GetFeatureInfo are the table, XML and text.


texts:
I have a set of lines that were created as lines between point pairs. Now I need to know, which points the lines connect. I thus want to join the Ids of both points to the line.
So far I tried to use spatial join, but one ID is always right, while the other ID is always wrong. The lines were created from the points, so they mark the vertices of the lines.
In the end I want to have a single ID, derived from both IDs. For this question I only need to know how to get 2 columns in the line table with the IDs from the vertice points.

How can I (spatial) join two points to one line?

I am using ArcGIS 10.3 on Windows 7.


texts: Image of police car -> dynamically draw police car number

texts: Let's say I have 10 police cars using the same image, and I want each one to be numbered 1 - 10, but the number 1 is dynamically drawn on top of the police car image.  What is the class / interface  I use for that with C# and ArcEngine 10.


texts: Why does the AddGeometryColumn function not exist in RSQLite?

texts:
I am trying to add the geometry to a new table with the following statement:

library(RSQLite)
# db.sqlite was created with spatialite-gui
con <- dbConnect(RSQLite::SQLite(), dbname  = "C:/GIS/db.sqlite")

rs <- dbSendQuery(con,
  "CREATE TABLE test_geom (
    id INTEGER NOT NULL
      PRIMARY KEY AUTOINCREMENT,
    name TEXT NOT NULL,
    measured_value DOUBLE NOT NULL);")
dbClearResult(rs)

rs <- dbSendQuery(con,
"SELECT AddGeometryColumn('test_geom', 'Geometry', 4326, 'POINT', 'XY');")

But it generates the following error:

Error: no such function: AddGeometryColumn

If I execute the same statement in spatialite-gui if it works.


texts:
I am trying to display the query results in an Infowindow with a click on my map. But the infowindow displaying on a second click not with a firstclick on map.

Can you run the following code and suggest where I have to modify the code?

<!DOCTYPE html>

    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
    <meta name="viewport" content="initial-scale=1, maximum-scale=1,user-scalable=no"/>
    <title>Simple Map</title>
    <link rel="stylesheet" href="https://js.arcgis.com/3.17/esri/css/esri.css">


texts: How to select one band in image using rasterio

texts:
I have two different images and I want to select one specifi band from each one of them and calculate regression between them. Each one of those bands is a result of calculation I did (in the first image , band' no. 4 is the results of NDVI calculation, and the other band is SAR db values).
The problem is that I haven't found a way to do it.

I used  "read" in order to try to select one band but it didn't work when I tried to do things with this.

My end goal is to actually calculate linear regression between those two bands ,which I also still don't know how to do in rasterio library.

my code:

import pandas as pd
import rasterio
import numpy as np
from rasterio.plot import show
import seaborn as sns
%matplotlib inline

#NDVI image

NDVI = rasterio.open(r"C:\Users\PATH\TO\IMAGE\python\NDVI.tif")

show(NDVI)

#SAR image

SAR = rasterio.open(r"C:\Users\path\to\IMAGE\python\SAR.tif")

NDVIband=NDVI.read(4)
SARband=SAR.read(2)

#from here nothing work and I feel lost.
#I thought to use this library but seems like it works only with pandas (?)
from sklearn.linear_model import LinearRegression

as i'm very new to python and image processing, any clue will be valuable for me ...


texts: Copying only selected with arcpy.CopyFeature_management?

texts:
I am trying to test if each feature of a feature class list intersect each feature of a shapefile (clip layer). I'm working with arcpy in ArcGIS 10.0

The arcgis desktop help says that if the feature comes from a geodatabase, all the features will be copied but if it's a layer, only selected will be copied (http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//001700000035000000).

In my code, I'm working with a feature class list from a geodatabase.In order to copy only selected, I first convert each feature class to a feature layer :

for fc in fcList:
    for row in rows:
        where_clause = "secteur = '{0}'".format(secteur)
        arcpy.MakeFeatureLayer_management(pochoir, "ptmp", where_clause)#pochoir is #the clip feature. For each row I make a tmp corresponding to the geometry that I'll use #after.
        arcpy.MakeFeatureLayer_management(fc, "toselect")
        arcpy.SelectLayerByLocation_management("toselect", "INTERSECT", "ptmp")
        arcpy.CopyFeatures_management("toselect", "selected") ## Here is the problem : #"selected" should be empty when not intersecting with "ptmp" but the whole "toselect" (so #the whole fc) is copied even if not intersecting.
        nbrow = arcpy.GetCount_management("selected")
        outName = "{0}_{1}_{2}".format(fc, secteur, count)

        if nbrow > 0: ## So all the layers goes in the if even if they do not intersect.
            do something
        else:
            do something else


texts: Switch language of labels in Mapbox studio based on users locale

texts: Is there an easy way to change all the map labels to use the users locale instead of a fixed locale. In Mapbox studio the default layers are set to {name_en}. Does that mean that each langugage needs its own Map-Style, or is there an API to change the language?


texts: Getting Esri Brushes to appear as dynamic resources in Visual Studio?

texts:
I used the ArcGIS Pro SDK to create my Add-In project and to create the DockPane but the Esri Brushes don't show up as Dynamic Resources.

According to Esri (https://github.com/Esri/arcgis-pro-sdk/wiki/proguide-style-guide)

"ArcGIS Pro XAML styles should be used as DynamicResources while styling the UI elements in an add-in. Add this code snippet to your xaml so that these Dynamic resources will load in Design mode in Visual Studio."

<UserControl ...
      xmlns:extensions="clr-namespace:ArcGIS.Desktop.Extensions;assembly=ArcGIS.Desktop.Extensions"
      ...>
<UserControl.Resources>
  <ResourceDictionary>
    <ResourceDictionary.MergedDictionaries>
     <extensions:DesignOnlyResourceDictionary Source="pack://application:,,,/ArcGIS.Desktop.Framework;component\Themes\Default.xaml"/>
            </ResourceDictionary.MergedDictionaries>
      </ResourceDictionary>
 </UserControl.Resources>

However, the SDK already puts this reference in the .xaml file but no dynamic resources show up in the Designer options. And manually typing in the Esri Brush name doesn't do anything.


texts: Seeking QGIS equivalent to ArcMap Editor Menu Clip?

texts:
In ArcMap I can start an edit session, drop down the editor menu and choose "Clip" then decide whether to preserve or discard the area(s) in any other displayed layer that intersect with the clipping layer

Is there a QGIS equivalent for this?


texts: Using the Wizard in Overpass Turbo to download all layers within the Bounding-Box

texts:
I am using the Bounding-Box and the wizard query to download data from Overpass Turbo.
I am using the id openStreetMap editor to find the right query word, and it is working great.
However, I was wondering is there a way do download all the data from the bounding box without the necessity of run each layer separately?
I was looking in the Overpass API/Language Guide, but I did not find an answer.


texts: Linking to a specific point in an arcgis online map

texts: I have a map in arcgis online with point data. Each of these points shows houses that have been destroyed by time. We have a wiki that tells the history about these houses, so I was wondering if there is a way to either link to the points location in the map (so I could just drop the link in the wiki), or extract a png file of an extent around the point and show that image in the wiki?


texts: Is there a Filesize Limit for Reading GDBs using Fiona?

texts:
I have python code that currently works fine with a GDB whose size is about ~ 82.7MB.  However when I try to use the same code against a GDB whose size ~ 612.8MB the code seems to hang eventually leading to the script terminating.

Specifically the code where the hang takes place is (first line):

buildingFeatures = [pt for pt in fiona.open(sourceDatabaseFile, 'r', layer=layer1)]

waterFeatures = [pol for pol in fiona.open(sourceDatabaseFile, 'r', layer=layer2)]

    for building in buildingFeatures:
        buildingShape = shape(building['geometry'])

        #if the geometry intersects the spatial index
        for id in idx.intersection(buildingShape.bounds):
            water = waterFeatures[id]

            #Check to see if current hydrographySrf feature intersection is a perennial water
            if (('ZI024_HYP' in water['properties']) and (water['properties']['ZI024_HYP'] == 1)):
                waterShapeGeometry = shape(water['geometry'])

                #Check to see if perennial water actually contains building point
                if waterShapeGeometry.contains(buildingShape):
                    discoveredBuildingsContainedInWaters.append(building['id'])

Is there a better way to open a large GDB using Fiona or is there some sort of size limit for GDBs used with Fiona?


texts: From raster histogram/raster make frequency tables in Python

texts:
Im new in programming... I would like to know, how I can get a frequency table from a rasters histograms? Or how can I generate the table from a raster?

I would like a table of the histogram showed below


texts: Relative Path issue within Model Builder

texts: I keep getting a 000229 Error that my model cannot open the output from the first task and use it as an input for my second task.  I am using relative paths because the model is intended to be published as a geoprocessing service for a webapp.  I have narrowed the issue down to relative paths because the model exicutes perfectly when i use a full path name that points to the same output folder that the relative path poitns to.  I am also stumped because the first output is created and stored propertly using a relative path but the next task cannot open it that output.  Any help is apprecaited.


texts: What can be used to show a 'contour' of equal cost across a raster cost surface, radiating from a single point/line/polygon?

texts:
My challenge is to create a map showing the propagation of waves, similar to sound waves, from their source point, by marking 'contours' of certain intensity at predetermined distances. I have managed this by simply creating buffers at these distances. However, it has additionally been requested that the 'contours' now account for surrounding buildings and infrastructure which will attenuate the waves more quickly. I have created a cost surface based on the level of attenuation provided by the surrounding features. However, I would like to know if there is a way to use the distance/cost path tools to force a process whereby the raster calculates radially out from a point and stops when a certain value has been cumulatively attained by all the cells it has crossed.
Otherwise, is there any software which could be used for a similar purpose which is compatible with GIS and could be used and/or any other solutions within arctoolbox.


texts: Optimizing/reducing geolocation data for speeding up visualization

texts:
I have an application that receives geolocation data from clients (cars) every minute and shows this on a web page (like owntracks). The data is stored simply in an SQL database - 1 point is 1 row, and the display is also similarly simple, all points are shown for the queried period.

My main objective is simple: I need to collect and visualize the location of the clients. The current 1 minute data collection frequency is sufficient.

However, as more and more data gets collected, this is becoming frustratingly slow on the query side. Due to the amount of data collected, it is not even possible to show data spanning multiple days, therefore I am looking at some optimization measures, that would reduce possibly the granularity of the data, but enable visualization of longer periods quicker.

I am not sure however which of the following measures I should follow:

Optimize data as it is received, i.e. do not store data points that are not (significantly) different from previous data points. Some of the clients are quite stationary yet transmit data every minute, therefore this could eliminate the storage of this redundant data, albeit at the cost of losing some valuable information as well (e.g. if I only store the first point for a given location, I cannot retrieve later how much time was spent by the client at that location).
Store all incoming data, optimize regularly, e.g. by running a map-reduce job every hour that would eliminate duplicate data points and possibly optimize others e.g. with the Douglas-Peucker algorithm.
Store and keep all incoming data and optimize upon query, similarly as to point (2). I suspect this would slow down the queries to an unacceptable extent.

Can you suggest the best approach and/or alternative measures? I would be also happy to consider moving to a different storage mechanism instead of SQL, if that could help with this.


texts: Copy geometry from one layer onto attributes of another?

texts: I have two polygon layers, A and B, that have the same features but totally different attributes. I want to use the attributes from layer A, but I want them to have the more accurate geometry from layer B. How can I copy the geometry from a feature in layer B and paste it onto the attributes of the same feature in layer A?


texts: Iterating through columns in table using ModelBuilder?

texts:
I am unable to find any info on how to use ModelBuilder to run a certain spatial analysis function on each column in a table.

I have a table to sample locations and at each location there are a series of measures (e.g. temperature, flow rate, number of certain species).  For each of those attributes (columns in the table), I was to run the IDW function from the surface tools.

This would be easy if I had each attribute in a separate shp/feature class, but not very efficient.  I tried Iterate Feature Selection, but was unable to pass the field name as the Z value in the IDW function.

This is for ArcGIS 10.1 (Basic and Spatial Analyst license)


texts: Get all the Tile layers added to the map using OpenLayers?

texts:
I want to get all the Tile layers which are added in the base_map. When I am using  map.getLayers() I get the list in the array_.  How can I identify all the charts and remove the individual charts.
My process to add the code is:
var layer = new ol.layer.Tile({
     source: new ol.source.TileWMS({
     url:url,
     serverType: 'geoserver'
    })
});
map.addLayer(layer);

Any suggestions will be highly helpful.


texts:
On one of our projects, I have been creating .mxd Files that reference a PostgreSQL Database. For some corporate reasons, the database had to be relocated to a new server, leaving all the mxd connections broken.

This was supposed to be a simple correction with python and the arcpy library. I have been running into some issues, though. Here is the code I've been using to achieve said purpose:

import arcpy
import os
mxdPath = r"..." #Path where the files are stored
originalRef = "Database Connections\Connection to xxx-180.sde"
newRef = "Database Connections\Connection to new-server.sde"

for root, dirs, filenames in os.walk(mxdPath):
    for f in filenames:
        mxd_path = os.path.join(root, f)
        mxd = arcpy.mapping.MapDocument(mxd_path)
        mxd.findAndReplaceWorkspacePaths(originalRef, newRef) #Point where execution halts.
        mxd.save()
        del mxd

The problem is: This script takes forever to run (it's working, but taking over 2 hours on a single file!). I have tried debugging it, and apparently, the findAndReplaceWorkspacePaths takes just too long to update (I would guess that as the layers start to be updated, the program tries to redefine the query layers to the connection that no longer exists, but doesn't recognize it immediately. I am running into the same problem while opening the files per se with ArcMap - a long wait until the document manages to load).

What could be done?

I was hoping there was a way to pause the reloading of the query layers and then reenabling them.


texts: Prevent SAGA in QGIS 3.0 from creating certain outputs?

texts:
I'm wondering if there is a way to prevent SAGA from generating temporary files. For example, if I just want to generate a slope map using SAGA I have to wait for all 11 other grids to get processed. Is there a variable I can set or enter that would prevent them from even being created?

The only option we seem to be given is to not add them, but that doesn't prevent having to wait for them to process. I'm processing huge LIDAR raster files, so it would be great to avoid all this unnecessary wait time.


texts: Deleting inner holes/rings/borders of a spatial polygon in R?

texts:
I am relatively new to GIS. I have a shapefile with spatial polygon.

> class(Park)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"

I use the following code to plot my data

tm_shape(Park)+tm_borders()

I would like to get rid of all the holes in the middle of my polygon, e.g. to keep only the outer border.

How could I do that in R (hopefully in tmap package)?


texts: Export arcpy.ListFeatureClasses to txt for later use

texts:
I have to re-project all Feature classes of a very large .sde. In order to iterate through the whole .sde I'm using arcpy.ListFeatureClasses. Getting all Feature classes into the Feature class list takes about 6 hours. As the Support window is open for a limited time only I would like to generate the Feature class list in advance and save it into a txt in order to be able to directly start the re-projection in the next Support window.
I'm exporting and importing the list with pickle (Export and Import are working) but when I try to iterate through the imported pickle list, trying to apply arcpy.describe an error appears saying, that this Feature class does not exist.

How do I have to Export/Import the ListFeatureClasses Output in order to be able to reuse it later on with arcpy.describe?

I'm working with ArcGIS Desktop 10.6.

My Code in the first file to generate the fcList roughly looks like this:

In a Definition:

fcList = arcpy.ListFeatureClasses("%s.*" % (schema))
return fcList

In the Main section:

with open(fcDictPath, "wb") as FCs:
     pickle.dump(fcList, FCs)

My Code in the file to perform the re-projection Looks like this:

with open(fcDictPath, "rb") as FCs:
 fcList = pickle.load(FCs)

for fc in fcList:
    dsc = arcpy.Describe(fcWithoutSchema)


texts:
I am trying to process MODIS files in R, which are in HDF4 format (link to data). I have experience reading in HDF5 and netCDF files in R, but HDF4 seems to be a common problem in R. I am using a MAC (Mountain Lion OS X 10.8.5), and have found the following webpages (link1, link2) helpful in setting up GDAL and the r package, rgdal. My problem is that is doesn't look like I have any HDF4 drivers, and I'm wondering what step I have missed?

library(rgdal)
gdalDrivers()

...
...
44              GTX                                 NOAA Vertical Datum .GTX   TRUE FALSE
45              GXF                             GeoSoft Grid Exchange Format  FALSE FALSE
46             HDF5                       Hierarchical Data Format Release 5  FALSE FALSE
47        HDF5Image                                             HDF5 Dataset  FALSE FALSE
48              HF2                               HF2/HFZ heightfield raster  FALSE  TRUE
49              HFA                              Erdas Imagine Images (.img)   TRUE  TRUE
...
...

So, nothing listed for HDF4 as shown in point 6.2 of the example. Does anyone know what my problem might be? The issue of reading HDF4 data into R seems to be a common one. I would have no problem doing some sort of batch conversion to HDF5 or netCDF if someone has a suggestion for this.


texts: Purpose of "myPlugin_dialog.py" file created by QGIS Plugin Builder

texts:
I usually use QGIS Plugin Builder to create working templates for QGIS plugins.

While a generic Python plugin needs the files below;

__init__.py
metadata.txt
resources.py (translated from resources.rcc via pyrcc4)
myPlugin.py
myPlugin_UI.py (translated from myPlugin_UI.ui via pyuic4)

Plugin Builder also creates the following file;

myPlugin_dialog.py

I can understand and use it but why do we need this file myPlugin_dialog.py?

If I import "myPlugin_UI.py" at the beginning, I can already access and interact with the form objects directly from the main module "myPlugin.py"

Why should I need an additional file? What is the advantage of this structure?

(Please explain in basic word since I am not fully fluent with PyQt)


texts: Mapbox Geocoder with a Leaflet map. Zoom to result, how?

texts:
I have a Lealfet map together with the Mapbox geocoder outside the map. The Geocoder and the map works. But i need to zoom to the result from the geocoder.

Some code:

          var map = L.map('map').setView([35.936050, 10.541736], 5);

          var tiles = L.tileLayer('http://server.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer/tile/{z}/{y}/{x}.png', {
                maxZoom: 18,
                attribution: 'Map data © OpenStreetMap contributors',
            }).addTo(map);

          var geocoder = new MapboxGeocoder({
           accessToken: mapboxgl.accessToken,
           mapboxgl: mapboxgl,
           marker: false
           });

It just zooms into a random place, when i click a result from the dropdown. In below Example it zooms to a place just outside Somalia

Map result:

Can I somehow log the result from the geocoder, so i can see the data?


texts: ArcGIS JavaScript API Sandbox doesn't properly save html files (they doesn't work)

texts:
I'm learning html/JavaScript code from ArcGIS API Sandbox.
(http://help.arcgis.com/en/webapi/javascript/arcgis../sandbox/sandbox.html?sample=map_webmap)
The code in the sandbox works fine and creates a map in the 'results pane'. However, If I download and save the code from the sandbox to my desktop in html format, and then I open the html it won't generate any map.

I've tried changing parts of the code (like some parameters in the tags, even the '<'!DOCTYPE> line) but I still can make it work. Any advice on how to make that code run?


texts: Polygon layer became point layer when re-opening the project

texts:
I´m making a geological map on QGIS 3.6 and yesterday I finished mapping all the lithology in a categorized polygon layer. Today I opened the same project, to continue working on it, and I found that the lithology layer had become a point layer, without any points so nothing is showing. The attribute table is fine and so are the categories.
Is there anything to be done or do I just start over? And is there a way to avoid this in the future?


texts: Some point features in feature class not displaying in ArcMap?

texts:
I am trying to display a point feature class, but a few features from the same data set are not displayed.

Originally I have heaps of point features in a point feature class, which are well displayed in my map extent (ArcGIS 10.2 / ArcGIS 10.3.1).

Incremental updates are provided by the supplier which are handled and updated (added/deleted) to the main feature class through SQL.

main feature class

When imported to ArcMap 10.2

In ArcMap using SQL table I run a query and 4 points are returned (2 previous and 2 updated) but when exported only 2 records (previous) are displayed.

I am struggling to find out what is wrong with the updated records, so are not displayed.


texts:
I just uploaded a new photo but it has yet to show up.  I applied the query again in cartodb but still does not show.

Shouldn't it be syncing automatically?  And how long should it take using a free account?


texts: Dealing with data of regions situated in two UTM zones using GDAL

texts:
I must work with a projection, in order to execute the intended geometric operations. I'm looking for a projection that could cover an area of up to 200km X 200km. I'm aware of some local projections, but I don't have an EPSG for them.
The ideal solution to me would be a UTM projection  centered in my region of study, but I don't have an EPSG for it.
I'm working in Python (GDAL library). The program is intended to be used everywhere in the globe (Possibly not the polar region)


texts:
I have been unable to run R from within QGIS on my computer. I have R version 3.03 installed and QGIS 2.4 on Windows 7. R is located at C:\Program Files\R\R-3.0.3\bin\and this is the path e specified in Processing Options -> Providers -> R Scripts -> R Folder in QGIS.

The error message I receive is

Missing dependency.This algorithm cannot be run :-(
This algorithm requires R to be run.Unfortunately, it seems that R is not installed
in your system, or it is not correctly configured to be used from QGIS
Click here to know more about how to install and configure R to be used with QGIS

I have tried changing the R folder to other R executables, but nothing works.


texts: Lots of layers in one or multiple services? (and why)

texts:
I have a conundrum that I am getting mixed advice on how to go about.  Therefore id like to put it to GIS-SE for some justified answers.

Scenario:

Client has a web mapping application.  Does not want to split into multiple smaller applications.  Although this goes against what the modern day approach is for maps on the web,  (i.e. lots of focussed web-map apps over one main web-map) I strongly believe that for some users, trying to replicate a GIS application on the web is ok (sometimes).
Client has cached as much of his basemap layers into separate services.
Client still requires an additional 600-700 layers in a dynamic map service...
The service will be published with all of these layers turned off.
It is not anticipated that users will turn on more than 10-40 layers at a time.

I imagine your initial reaction to this is similar to mine (600+?! WTF?!)

However - the requirement is set in stone and why not?  Their prior ArcIMS application had similar functionality, so why can't this newer ArcGIS Server product do the same?  The users potentially need to be able to cross compare and perform analysis on the entire range of layers, even if the layers belong to other departments.

Before you jump to conclusions, the client is a clued up ArcGIS Server admin.
They have administered the 600 layers by all of the best practice rules:
e.g. scale ranges combined with definition queries; annotation over labelling; generalising complex layers at small scales; publish as MSDs; etc

Problem:

What is the better approach here?

Publish all 600 layers into one dynamic map service
Split the layers into logical groupings (hydrology, planning, ecology, utilities, etc)

If you go with #1, and you have a few complex layers turned on.  If you want to turn on a simple points layer, then ArcGIS Server will still have to render the entire layers being displayed all over again.

If you go with #2, each time your making a request, potentially, the web application could have to make several GET requests for ExportMaps from the individual map services (is this bad, or does it create additional load to ArcGIS Server over #1?)

And then this leads to configuration & tuning to ensure everything is as quick as possible.
We can scale the back end of ArcGIS Server to multiple hosts and have some good hardware to sit it on.

If you go with #1, you can throw the maximum # of instances that you AGS can handle.

If you go with #2, I assume you evaluate the performance of the map services (load testing and look at wait times) and address the min/max instances accordingly to ensure there is not one service that is a 'weak link'.

Im currently leaning towards the #2 approach, as my head is still telling me that having 600 layers in one service is madness, but if they are all turned off by default, there really is no issue.

Would love to hear your thoughts.  Let me know if you need more info via the comments, but not looking for answers like 'use a desktop application' or 'educate them to do things differently'

From discussions in comments, I failed to mention another consideration.  The application that the service will be consumed within, has the ability for layer level security (at the application level).  Therefore the group of users (that is fairly large) are assigned to a particular role, and that role will have access to the full 600 layers.  Other roles will be limited.


texts:
This seems silly to me, but I haven't hit this before. I have multiple (30+) point tables, 100 data points each, and want to start designing symbology and doing things like point distance. However, my data source is not finalized. (I have been asked to start now, but will receive updated coordinates in a week). Is there a table format (access table, gdb, etc) where I can calculate pt distance (multiple times per table), design a multi-part symbology based on attribute, throught the table (not making it a shp), and when the data is finalized, have the table auto-update? (and then, in theory, create a shp from that 'final' data that retains symbology, etc)?

This must be possible- but i know i can't edit/replace a dbf that's part of a shp (without starting over), and it seems that edits can only be made in either ArcGIS   OR  in Access, but not both.

thanks for ideas,


texts: Iterating through z-values in ModelBuilder

texts:
I have a series of fields that I am trying to run the IDW tool against. I'm looking for a way iterate through the feature class, using a different field as the z value each time.

This is my data: http://imageshack.com/a/img849/115/zqkc.png

I'm looking to use each field following Village as a z value. Typically I would just do this all manually but I have a ton a data that I need to run the IDW for so I need to automate it.

I've looked through the forums already and I haven't been able to find any responses that include all the necessary details for me to make this happen.

Any insight would be incredibly helpful.


texts: Automating ArcGIS map production using Python?

texts:
I have this layout; it's a map that I divided it to 25 parts using FISHNET from arctoolbox. I want to prepare only one layout and then using a Python script I want to be able to automate the production of the other 24.

The problem is when I use the layout of one part from the fishnet I got 24 empty PDFs and only one with the layout I already started with. But when I use the layout of the whole map it kind of works, and I got 25 PDFs with the 25 parts of the Fishnet but it look like this:

So who can help me to get PDFs with those parts but on the whole page!

My code:

gp=arcgisscripting.create()
arcpy.env.overwriteOutput=True
mxd=arcpy.mapping.MapDocument(r"D:\PFE\data\carte.mxd")
df=arcpy.mapping.ListDataFrames(mxd, "*")[0]
print df

source_geod= arcpy.mapping.Layer(r"D:\PFE\data\PFE\traitement\geod.lyr")
source_RR= arcpy.mapping.Layer(r"D:\PFE\data\PFE\traitement\RR.lyr")
source_oueds= arcpy.mapping.Layer(r"D:\PFE\data\PFE\traitement\oueds.lyr")
source_SBKA= arcpy.mapping.Layer(r"D:\PFE\data\PFE\traitement\SBKA.lyr")
source_Veg= arcpy.mapping.Layer(r"D:\PFE\data\PFE\traitement\Veg.lyr")

i=0
fc=(r"D:\PFE\data\fishnet1.shp")
rows=arcpy.SearchCursor(fc)

for row in rows:
    print i
    print row

    layerFishnet_src=gp.MakeFeatureLayer_management(fc, "lyr")
    layerFishnet=gp.SelectLayerByAttribute_management("lyr", "NEW_SELECTION",'("FID" = ' + str(i) + ')' )

    geod_i= arcpy.Clip_analysis(r"D:\PFE\data\Geod.shp", layerFishnet,
    r"D:\PFE\data\Geod_"+str(row.getValue("FID"))+".shp")

    rr_i= arcpy.Clip_analysis(r"D:\PFE\data\RR.shp", layerFishnet,
    r"D:\PFE\data\RR_"+str(row.getValue("FID"))+".shp")

    oueds_i= arcpy.Clip_analysis(r"D:\PFE\data\oueds.shp", layerFishnet,
    r"D:\PFE\data\oueds_"+str(row.getValue("FID"))+".shp")

    SBKA_i= arcpy.Clip_analysis(r"D:\PFE\data\SBKA.shp", layerFishnet,
    r"D:\PFE\data\SBKA_"+str(row.getValue("FID"))+".shp")

    Veg_i= arcpy.Clip_analysis(r"D:\PFE\data\Veg.shp", layerFishnet,
    r"D:\PFE\data\Veg_"+str(row.getValue("FID"))+".shp")

    geoidlyr=gp.MakeFeatureLayer_management(geod_i, "geoidlyr")
    rrlyr=gp.MakeFeatureLayer_management(rr_i, "rrlyr")
    ouedslyr=gp.MakeFeatureLayer_management(oueds_i, "ouedslyr")
    SBKAlyr=gp.MakeFeatureLayer_management(SBKA_i, "SBKAlyr")
    Veglyr=gp.MakeFeatureLayer_management(Veg_i, "Veglyr")

    geoid=arcpy.mapping.Layer(geoidlyr)
    rr=arcpy.mapping.Layer(rrlyr)
    oueds=arcpy.mapping.Layer(ouedslyr)
    SBKA=arcpy.mapping.Layer(SBKAlyr)
    Veg=arcpy.mapping.Layer(Veglyr)

    lyrr = arcpy.mapping.ListLayers(mxd)[3]
    print lyrr

    lygeoid = arcpy.mapping.ListLayers(mxd)[2]
    print lygeoid

    lyoueds = arcpy.mapping.ListLayers(mxd)[4]
    print lyoueds

    lySBKA = arcpy.mapping.ListLayers(mxd)[6]
    print lySBKA

    lyVeg = arcpy.mapping.ListLayers(mxd)[5]
    print lyVeg

    arcpy.mapping.UpdateLayer( df, lygeoid ,geoid, false )
    arcpy.mapping.UpdateLayer( df, lyrr ,rr,  False)
    arcpy.mapping.UpdateLayer( df, lyoueds ,oueds,  False)
    arcpy.mapping.UpdateLayer( df, lySBKA ,SBKA,  False)
    arcpy.mapping.UpdateLayer( df, lyVeg ,Veg,  False)

    titleTxt= arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT", "title")
    titleTxt=row.getValue("name")
    print titleTxt
    arcpy.mapping.ExportToPDF(mxd,r"D:\PFE" + row.getValue("name") + ".pdf")
    i = i + 1


texts: How to set Vertical Coordinate system when creating a 3D shapefile?

texts:
I'm working with GDAL v1.11 and am seeing the vertical CS stripped when exporting a 3D shapefile (from oracle spatial). This is the command I'm running:

ogr2ogr -f "ESRI Shapefile" B_POINT.shp (oci_url) -sql(sql_select_statement) -nlt POINT25D -a_srs "EPSG:4326+3855" -nln B_POINT

I've tried specifying the a_srs with several different options

EPSG Code: "EPSG:4326+3855"

PROJ.4 Definition: "+proj=longlat +datum=WGS84 +geoidgrids=egm08_25.gtx +vunits=m +no_defs"

Filename: /apps/gdal-1.11.0/compound_wkt.wkt

Where file>cat /apps/gdal-1.11.0/compound_wkt.wkt

COMPD_CS["WGS 84 + EGM2008 geoid height",
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]],
VERT_CS["EGM2008 geoid height",
    VERT_DATUM["EGM2008 geoid",2005,
        AUTHORITY["EPSG","1027"],
        EXTENSION["PROJ4_GRIDS","egm08_25.gtx"]],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Up",UP],
    AUTHORITY["EPSG","3855"]]]

I haven't had any luck in exporting the vertical CS with any of these ogr2ogr a_srs options. I've even installed PROJ.4 and tried setting the s_srs and t_srs with no luck.

The egm08_25.gtx file is downloaded and copied to /apps/proj-4.8.0/nad/ directory, per this post: http://osgeo-org.1560.x6.nabble.com/GTX-format-td3840948.html

Is there any way to export the vertical CS for a 3D shapefile, or is this still a bug? I have the same issue trying to export a FileGDB.

EDIT: I can view the WKT of the compound CS here:

gdal-1.11.0>testepsg epsg:4326+3855
Validate Succeeds.
WKT[epsg:4326+3855] =
COMPD_CS["WGS 84 + EGM2008 geoid height",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    VERT_CS["EGM2008 geoid height",
        VERT_DATUM["EGM2008 geoid",2005,
            AUTHORITY["EPSG","1027"],
            EXTENSION["PROJ4_GRIDS","egm08_25.gtx"]],
        UNIT["metre",1,
            AUTHORITY["EPSG","9001"]],
        AXIS["Up",UP],
        AUTHORITY["EPSG","3855"]]]

Simplified WKT[epsg:4326+3855] =
COMPD_CS["WGS 84 + EGM2008 geoid height",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    VERT_CS["EGM2008 geoid height",
        VERT_DATUM["EGM2008 geoid",2005],
        UNIT["metre",1]]]

Old Style WKT[epsg:4326+3855] = GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]
ESRI'ified WKT[epsg:4326+3855] =
GEOGCS["GCS_WGS_1984",
    DATUM["D_WGS_1984",
        SPHEROID["WGS_1984",6378137,298.257223563]],
    PRIMEM["Greenwich",0],
    UNIT["Degree",0.017453292519943295]]
PROJ.4 rendering of [epsg:4326+3855] = +proj=longlat +datum=WGS84 +geoidgrids=egm08_25.gtx +vunits=m +no_defs

The "ESRI'ified WKT" is what I see in my exported shapefile and filegeodb:

shapefile>ogrinfo -al B_POINT.shp | more
INFO: Open of `B_POINT.shp'
      using driver `ESRI Shapefile' successful.

Layer name: B_POINT
Geometry: 3D Point
Feature Count: 1
Extent: (39.212616, -6.130917) - (39.212616, -6.130917)
Layer SRS WKT:
GEOGCS["GCS_WGS_1984",
    DATUM["WGS_1984",
        SPHEROID["WGS_84",6378137,298.257223563]],
    PRIMEM["Greenwich",0],
    UNIT["Degree",0.017453292519943295]]

I know that shapefiles support vertical coordinate systems (e.g., page 26 here: http://www.geo.hunter.cuny.edu/gis/docs/Whats_New_In_ArcGIS_92.pdf). There must be a configuration in GDAL to get the shapefile ESRI WKT to pickup the compound CS.


texts: i created grid panel using ExtJS with bbar button, the buttton appears to be similar with the background of the panel which makes it to resemble like plain text not really button. so i want to add icon to the button or give it different color. is that possible? how to do that?


texts: Accessing x and y coordinates of a POLYGON in MySQL

texts:
In my database. I have a spatial column called coordinates which is of datatype POLYGON. so for example. The lunch room has coordinates POLYGON(0 2, 1 2, 3 2, 5 2, 0 2). That's just a hypothetical example.
How do I access the individual x and y coordinates for each polygon in the database?


texts: Why my Geoserver on JBOSS/RHEL is slower than on Tomcat/Windows?

texts:
I have my Geoservers in 2 different environments: Tomcat on Windows and JBOSS on RHEL. I was guessing that JBOSS on RHEL will be faster, but to make sure I made some test plans with JMeter.

Both:

are VMware servers
have same number of processor and same amount of memory
are using exactly same SHP file on each disk (no database)
are on the same network segment, so network distance to my JMeter (on my Workstation PC) are the same
are dealing with similar WMS requests and style from my JMeter

My observation shows that Geoserver on JBOSS/RHEL is noticeably slower than on Tomcat/Windows. This brought me to early conclusion that this is because the JBOSS/RHEL is not optimally configured.

But to be sure I made another test with both, but instead of using Geoserver, now only with static content (simply file) with format similar to Geoserver response:

8-bit and 32-bit 256x256 PNG
GML file
geojson file

Those static content are hosted on both JBOSS/RHEL and Tomcat/Windows. Interestingly now JBOSS/RHEL is faster than Tomcat/Windows.

I returned to the Geoserver test plan an played with Image Processing and Image Encoding options (image below) in Geoserver JBOSS/RHEL.

This is because I suspected the java/image processing is the source of problem. But the result is the same: Geoserver on JBOSS/RHEL is slower.

The speed difference of the two Geoserver is not only measured in JMeter, but even visually noticeable in OpenLayers(!).

Here is the summary of the 2 servers (server #20 vs #21):

Now the suspect is clear: the Geoserver on my JBOSS is not working optimal.

Where do i start to try and fix this?

ADDED INFOS:
Java version is added in the above comparison.

I checked this as suggested by @bugmenot. I think the most relevant item is the JVM Optimization in the Container Considerations. The link suggests to set several JVM parameters (1-7 in the table below). Another site suggests few other parameter (8-11).

Therefore i checked  both server:

The JVM parameter in JBOSS does not differ too much with Tomcat.

But still, I tried to increase (in JBOSS server #21) the InitialHeapSize to 128MB and MaxHeapSize to 2048MB, with the following command:

java -Xms128m -Xmx2048m -version

No error. But when i check with the command:

java -Xms128m -Xmx756m -XX:+PrintFlagsFinal -version | grep -iE 'InitialHeapSize|MaxHeapSize'

Both heap size doesn't change to the value I wanted.

So my question now:

Both server doesn't differ too much. So why?
Does Java version affects ?
How to set JVM memory amount (eg. increase heap size) for
JBOSS?

Below are images showing the test plans.

The above shows Summary Report of Geoserver test plan: JBOSS is slower than Tomcat.

The above shows Response Time of Geoserver test plan: JBOSS is slower than Tomcat.

The above shows Summary Report of Static Content (no Geoserver) test plan: JBOSS is faster than Tomcat.

The above shows Response Time of Static Content (no Geoserver) test plan: JBOSS is faster than Tomcat.


texts: Google Earth Engine - tile glitch when exporting GeoTIFF

texts:
I used GEE to export a GeoTIFF of WWF HydroSHEDS DEM (dataset here https://developers.google.com/earth-engine/datasets/catalog/WWF_HydroSHEDS_03VFDEM, GEE script below).

I realized that when I import the DEM GeoTIFF to ArcGIS, a glitch shows at the edge of (what I believe are) the GEE tiles. The glitch then disrupts the creation of Flow rasters (like Flow Accumulation here, only showing accumulations values 0 and 1).Edit: the "tiles" seem to have a dimension of 0.025° x 0.025°.

Has anyone encountered this issue and found a way to avoid it?

GEE code:

// WWF HydroSHEDS Void filled DEM
// Global
// 3 arc-seconds (90m)
var scale = 90

var collection = 'WWF/HydroSHEDS/03VFDEM';
var band = 'b1';
var region = ee.Geometry.Polygon([[[-75, 17.5],[-68, 17.5],[-68, 20.5],[-73.5, 20.5],[-75, 19],[-75, 17.5]]]);

var dataset = ee.Image(collection).select(band);
                //  .filterBounds(region);
print (dataset);

var elevationVis = {
  min: -50.0,
  max: 3000.0,
  gamma: 2.0,
};
Map.setCenter(-71, 18.5, 9);
Map.addLayer(
    dataset, elevationVis, collection);

// Exporting as geoTIFF, to personal Google Drive
var image = dataset
print('crs', dataset.select(0).projection()) // Use to get native projection
Export.image.toDrive({
  image : image,
  region : region,
  description : "HydroSHEDS", // Task name
  folder : 'ee',
  fileNamePrefix : 'WWF_HydroSHEDS_03VFDEM_EPSG3857', // File name
  scale : scale, // Size of cells in exported raster, in m - use same scale as native resolution
  crs : 'EPSG:3857'
  });


texts: Temporal subset and mosaic from a Landsat 8 image collection for each image in Google Earth Engine

texts:
What I need are image mosaics from a cloud-masked, chronologically ordered Landsat 8 image collection for each position in the collection. So, the first mosaic should be built from all images 1 through n, the second from the images 2 through n and so forth. The result should be a new image collection with n mosaics where the last mosaic is not actually a mosaic, but image n itself.
Is there a way I can build this into one function? So far, I have tried somehow iterating through the collection after turning it into a list, making mosaics and appending them to a previously-created image collection which is based on the first mosaic made from the whole image collection.
However, from the code below, I am getting the following error (my list/collection contains 15 images so this is definitely related):
List (Error)
 List.get: List index must be between 0 and 14, or -15 and -1. Found 15.

The code snippet looks like this (full code including standard L8 cloud masking here):
var first = l8Masked.mosaic(); // This is the first mosaic over the entire image collection (1 through n)
print(first);
Map.addLayer(first, visParams, 'Cloud masked');

var newCol = ee.ImageCollection.fromImages([first]); // This is the new image collection all mosaics should be merged to, starting with the first mosaic
var l8List = l8Masked.toList(l8Masked.size()); // Turn the image collection into a list to iterate over it

var makeMosaicCol = l8List.map(function(image) {
  var index = l8List.indexOf(image);
  var nextIndex = index.add(1); // Since the first mosaic was created in line 56, I am starting with the second one by adding 1 to the initial index
  var nextImage = ee.Image(l8List.get(nextIndex));
  var date = ee.Feature(null, {'date': nextImage.date().format('YYYY-MM-dd')}); // Extract the image's date to use it for an image collection subset
  var newStart = ee.Date(date);

  var l8new = ee.ImageCollection(l8Masked) // Filter the initial image collection by the current start date and the initial end date to subset it
      .filterDate(newStart, end);

  var mosaic = l8new.mosaic(); // Create a mosaic based on the current image up to the nth image
  var mosaicCol = ee.ImageCollection(mosaic); // Turn the mosaic image into an image collection to be able to merge it with the newCol collection
  var mosaicCollection = newCol.merge(mosaicCol); // Merge each newly created mosaic to the previously-created newCol image collection

  return mosaicCollection;
});

print(makeMosaicCol);


texts: Error when register database with ArcGIS server, failure to access the DBMS

texts:
The connection property set was missing a required property or the property value was unrecognized. Failure to access the DBMS server:

Machine: DESKTOP-TFF79J4 [db: The connection property set was missing a required property or the property value was unrecognized. Failure to access the DBMS server]


texts: Can an entry in Esri GeoPortal Server contain both ArcGIS Server and *.zip links?

texts:
Does anyone know whether a dataset in Esri GeoPortal Server can support both Preview (in ArcGIS Server), and Download (of a zip file) on the same dataset?

Eg this dataset has a link to Preview in ArcGIS Server. Could this same dataset also contain a Download option, which linked to a zip file? If so, how would that be set up?

(Open takes you to the REST endpoint. Preview, KML, NMF and LYR are auto-generated by ArcGIS Server. Can I also include a Zip file of related data, which I would create manually?)


texts: QGIS 3.0 - qConsolidate or similar plugin

texts: It appears that the qConsolidate has not been updated for Qgis 3.0. Are there any alternative plugins/functions available for exporting all layers and saving the project to a new folder?


texts: Leaflet control window - how to get image & text to appear together in the content of the modal window?

texts:
I have implemented leaflet-control-window plugin into my map and it works. Modal opens on marker click and shows text.
I need image to go below the text in the window but can't get it to work. When I put image url inside of content, it either shows text or image (depending what I put first in the brackets).
How to write it so that both img and text are visible? Here is an example of what I'm trying to achieve, Leaflet postcards, modal with short text description and an img gallery (2-3 photos).
Defining a marker and initializing:
   const marker1 = L.marker([45.899491, 16.844165],{icon:leafletIcon}).addTo(map);

   marker1.on('click',function(){

    L.control.window(map,{
        title: 'Test drive'

    })
        .content('./images/open-mailbox-with-raised-flag.ico'+'Loorem ipsum dolor sit amet.')
        .show()
   });


texts: How to build a model with multivalue input and output in dynamic folders?

texts:
I am trying to build a model that gets a multivalue as an Input which is located in different folders within a main folder, and stores the Output in new folders that have the same structure like in the input, but with different pathway.

I tried to select the data with a "multivalue iterator", then simultaneously create a folder with "parse path" (=x) that has a c:...\main_folder\%x%.

there are two problems with that:

I get just a single output (I'm running the data through alot of processes) from my multivalue input. How can I fix that? I want that every vector\raster will have a new name (%oldname%_%newname%), and each will be exported to a different folder.

So, for example, I want to take a 2 files:

c:...\main_folder\folder1\vector1

c:...\main_folder\folder2\vector2

and I want that each of them will be exported like this:

x:\outputs\main_folder\folder1\vector1_cliped

x:\outputs\main_folder\folder2\vector2_cliped

I want a dynamic pathway, so the tool could create the folders in the main folder by itself, and because there are a lot of sub-folders within the sub-folders, I want to create only the folders that contain the selected data.

Also, the folders structure may change, so I do not want to create a mother-folder and then create more folders with preconditioning.

I should mention I do not have any experience with python and scripting, but I believe that if you are able to give me a script to merge into my modelbuilder, then I'll learn how to edit it a little to match my data.


texts: ArcMap issue - Table text in data view area

texts:
I had the oddest of errors in ArcMap today. After zooming out after doing some union and polygon to raster operations, the system become very slow and one of the layer's table information is drawn in the data view window, column by column, line by line. There are 3700+ items in the table, and it takes some time to draw it, slowing down each scaling action with approx 30 sec.

Image showing how the table is drawn in the data view section of ArcMap (the whole table is not shown here, only some columns. However, the whole table is drawn if zoomed out enough. The text stop stacking if zoomed out.

The list still remains when removing all associated layers (containing the information drawn), it remains after save+quit (sadly I had already saved before recognizing the error). It does not remain when creating a new document and adding the shapefiles. Labels are turned off for all layers with associated information.

Does anyone have any idea what is going on?

Do I have to create a whole new document or can I save my mxd?


texts: IndexError when adding multiple fields to shapefile

texts:
As part of a larger project i'm looking to add multiple fields to a shapefile using an index. This is the snippet in question:

listFields = [['ITEM', 'TEXT'], ['SYSTEM', 'TEXT'], ['WATERCOURCE_BANK', 'TEXT'], ['AIMS_REF', 'TEXT'],
          ['CHAINAGE', 'TEXT'], ['GRID_REF', 'TEXT'], ['CUT', 'TEXT'], ['LENGTH_KM', 'TEXT'], ['COMMENTS', 'TEXT']]

print("Adding fields")
for addField in listFields:
print("   ...{}".format(addField[0]))

arcpy.AddField_management('table_stamp', addField[0], addField[1], addField[2], addField[3], addField[4],
                          addField[5], addField[6], addField[7], addField[8])

I'm getting the following error, pointing at the line of code above:

IndexError: list index out of range

As far as I can tell I have got my Indexing right. Nine fields to add so the indexes should be 0 to 8.

Anybody have any idea what the issue might be?


texts: Conversion of X,Y in Decimal Degree form

texts:
I found a DBF file from some govt site, It has X value as 6265001.698875606
and Y value as 2287922.0458578467 and so on.

I want to convert these values into Decimal Degree format so that I could create points and could plot in my ArcGIS software, can anyone suggest me how to convert these values to decimal degree format?

Few more examples are:-

      X          ,          Y
6189755.067873105,  2289006.0368001014
6203424.773145199   2289971.025355682
6257970.097453281   2297143.056615263
6360789.086368695   2298715.5160687715
6342105.368158773   2299252.8617788553
6186809.299486443   2300760.374511853


texts: Using SQL View as ImageMosaic Store index in GeoServer

texts:
I have an application server creating rasters from the vector data it receives periodically. The data contains either "real" measurements or forecasts about a specific weather phenomenon, and each case is mapped to the relevant database tables, like the following simplified schema

As I understand it, I can create a SQL View that contain all the needed information for a valid mosaic schema, getting the union of MeasurementReportRaster and ForecastReportRaster for a specified phenomenon, for example
the_geom      (the raster extent polygon as in the documentation, BoundingBox in my schema)
location      (the LocationAttribute)
ingestion     (the TimeAttribute, ReportDtm for measurements, ForecastDtm for forecasts)

I was planning on not using the GeoServer API to add/edit/remove granules, but make the neccessary changes in the database from outside GeoServer, and let GeoServer use the derived View for indexing purposes, only reading and never writing in the database.
I did not yet have any success setting up the ImageMosaic store, is this even possible?
If yes I seek some help with the configuration in datastore.properties and indexer.properties, since I still don't clearly understand how everything works, even after reading the documentation multiple times already.

I was adding more info about the configuration, and while I was trying some more things it finally worked. I still don't get it completely but I will probably ask some more specific questions in different Q&As.
I am using the following view:
                   View "public.rain_mosaic_index_v"
  Column   |            Type             | Collation | Nullable | Default
-----------+-----------------------------+-----------+----------+---------
 the_geom  | geometry(Polygon,4326)      |           |          |
 location  | character varying(255)      |           |          |
 ingestion | timestamp without time zone |           |          |

The location is the absolute path to a raster file.
In the folder where I store the raster GeoTiffs, I also have the two configuration files:
indexer.properties
Schema=*the_geom:Polygon,location:String,ingestion:java.util.Date
TimeAttribute=ingestion
Name=rain_mosaic_index_v
AbsolutePath=true
Caching=false
Recursive=false
CanBeEmpty=true
UseExistingSchema=true

datastore.properties
SPI=org.geotools.data.postgis.PostgisNGDataStoreFactory
host=localhost
port=5432
database=dbname
schema=public
user=username
passwd=password

Loose\ bbox=true
Estimated\ extends=false
validate\ connections=true
Connection\ timeout=10
preparedStatements=true


texts:
i have trouble working with my pointlayer, because it has too much information and i want to split it in pieces. There is GPS-pointdata recorded with GPS-Collars and all of the data is in one file. I want to have one file per Collar. Each Collar is declared by a value in the attributetable.
Is there a plugin or another way to select each collar and create a shapefile automatically?

I am using Qgis 3.


texts: How does one replace characters in a layer name using Python ("_" with " ")?

texts:
I can rename layers easily enough, but I really need to replace any "_" with a " ".

This code just renames layers, but I would like to do a find and replace on each layer name

mxd = arcpy.mapping.MapDocument("current")
layers = arcpy.mapping.ListLayers(mxd)

for lyr in layers:
    if lyr.name == "X":
        lyr.name = "Y"

arcpy.RefreshTOC()


texts: Is 'Shortest Route/Distance' analysis is possible in GeoServer using OpenLayers 2.8...?

texts:
I have configured a web based GIS Map using GeoServer and OpenLayers controls. Now I want to add 'Shortest Route' analysis in the same...

Any solution?


texts: Problem opening JP2 image using GDAL (interaction of -stats and -hist flags)

texts:
I came across unexpected GDAL behaviour while working with Sentinel-2 data granules. The Sentinel data should have values between 0 and 10000 when opened correctly.

However, I have noticed that when working with the data, pixel values above 10000 are retrieved. After doing some digging, I found this post which points to similar problems in ESA's Snap software. They refer to a bug in the OpenJpeg library.

This could be the problem in my case as well, as my GDAL binaries use JP2OpenJPEG in the background (see below).

The strange behavior is the following: The values are in the expected range when using the gdalinfo -hist command, but are wrong when using gdalinfo -hist -stats.

I also get values above 10000 when using the GDAL c-api directly.

My questions are

Why is this happening?
Is there a way to open such jp2 images with 15bpp correctly in GDAL?

Here is an example and a link to the corresponding file:

http://sentinel-s2-l1c.s3.amazonaws.com/tiles/48/M/WA/2016/6/29/0/B11.jp2

$ gdalinfo B11.jp2 -hist
Band 1 Block=640x640 Type=UInt16, ColorInterp=Gray
0...10...20...30...40...50...60...70...80...90...100 - done.
  256 buckets from -15.7549 to 8050.75:
  23286152 2037325 ... 11251
  Overviews: 2745x2745, 1372x1372, 686x686, 343x343
  Overviews: arbitrary
  Metadata:
    STATISTICS_MAXIMUM=8035
    STATISTICS_MEAN=274.9824
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=944.69220833573
  Image Structure Metadata:
    COMPRESSION=JPEG2000
    NBITS=15

$ rm B11.jp2.aux.xml
$ gdalinfo B11.jp2 -hist -stats
Band 1 Block=640x640 Type=UInt16, ColorInterp=Gray
  Minimum=0.000, Maximum=12356.000, Mean=272.236, StdDev=901.137
0...10...20...30...40...50...60...70...80...90...100 - done.
  256 buckets from -24.2275 to 12380.2:
  23286928 2372470 ... 0 0 1 0 1
  Overviews: 2745x2745, 1372x1372, 686x686, 343x343
  Overviews: arbitrary
  Metadata:
    STATISTICS_MAXIMUM=12356
    STATISTICS_MEAN=272.236306847
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=901.13698170459
  Image Structure Metadata:
    COMPRESSION=JPEG2000
    NBITS=15

I am using Ubuntu 16.06 with the Ubuntugis gdal version

$ gdalinfo --version
GDAL 2.1.0, released 2016/04/25

$ gdalinfo --formats | grep JP2
JP2OpenJPEG -raster,vector- (rwv): JPEG-2000 driver based on OpenJPEG library


texts: How can I apply a color ramp to a shapefile layer in QGIS 2.8?

texts:
How can I apply a color ramp to a shapefile layer to show something simple like population variance? If this was ArcGIS, I would do it by using the Quantities tab in Symbology but I have been unable to find an equivalent in
QGIS.  Do I have to use Python? If so, can someone point me to a good script for color ramps?

I'm using QGIS 2.8 on a Mac.


texts: Selection Events in ArcObjects using Java

texts:
I am trying to have an ArcMap addin receive ISelectionEvents, but I can't figure out how to register the event handler. I'm using ArcObjects 10.0 for Java. I'm able to catch IDocumentEvents as follows:

public void init(IApplication app) throws java.io.IOException, com.esri.arcgis.interop.AutomationException{
    super.init(app);
    this.app = app;
    MxDoc doc = (MxDocument) this.app.getDocument();
    doc.addIDocumentEventsListener(new EventCatcher()); // EventCatcher implements IDocumentsEvents and ISelectionEvents
}

Any suggestions on where to look would be greatly appreciated.


texts: What are the differences between 81989 and 27700 CRS for the British National Grid?

texts: We are using Oracle Spatial to run point in polygon queries in Great Britain using Ordnance Survey MasterMap and point data. MasterMap is on British National Grid (BNG) and has the Coordinate Reference System (CRS) 27700, while our point data has the Oracle Spatial 'BNG version' 81989.  Will it matter that we are using different systems?