texts:
Using the code below in this jsfiddle example:

tree = new Ext.tree.TreePanel({
        renderTo: "tree",
        root: new GeoExt.tree.LayerContainer({
            text: 'Map Layers',
            layerStore: panel.layers,
            leaf: false,
            expanded: true
            /*,
            loader: {
                filter: function (record) {
                    return record.get("layer").name.indexOf("tasmania") !== -1
                },
                createNode: function (attr) {
                    // add a WMS legend to each node created
                    attr.component = {
                        xtype: "gx_wmslegend",
                        layerRecord: mapPanel.layers.getByLayer(attr.layer),
                        showTitle: false,
                        // custom class for css positioning
                        // see tree-legend.html
                        cls: "legend"
                    }
                    return GeoExt.tree.LayerLoader.prototype.createNode.call(this, attr);
                }*/
        }),
        enableDD: true,
        rootVisible: false,
        width: 170,
        height: 300,
        floating: true,
        x: 380,
        y: 0
    });

I would like to get the legend of the added layers to the tree panel (for instance, "og:archsites"). The jsfiddle code does not work if I uncomment the loader, I think it should work but it doesn't. So, is this possible this way? I use geoext v1.1 in this jsfiddle example.


texts: Solution for Test Automation for web gis/mapping application?

texts:
I've to find a solution for test automation for web gis/mapping application that using OpenLayers 2.x or Leaflet as javascript libraries.

I'm oriented on evaluate Selenium WebDriver (or other similar tools, suggestion? ...): are there somewhere examples / tutorials / experiences or similar?


texts: What's the best way to display a huge aerial image using GeoServer

texts:
I have an aerial image - almost 1.5 GB , GeoTIFF - and I have been trying to add it as a layer in GeoServer...

Adding the layer went successful, but the performance was very bad.

I tried to cache the layer using GeoWebCache, and the performance is still bad.

Any tips?


texts:
I need to add interactivity to several WMS layers using leaflet, the problem i have is how i can identify what layer the user is clicking,  i usually use geojson for this and is easy, but for WMS layers how i can achieve this?

I have seen this question, i guess i can iterate over the layers but then again how i can know which one is the one i need?


texts: Ordnance Survey Town Plan plan vs other scales and which reference system they use

texts:
I am trying to layer up several ordnance survey maps on QGIS which include historical town plans (1st ed 1856-1908 1:500)  and terrain maps (OS Terrain 5 contours).

Although they are for the exact same area. The two sets of maps end up in different area on QGIS (ie quite a few miles apart) They seem to be on two different grid reference systems

Has anyone else had this problem, and how did they deal with it?


texts: Merge tiffs into a single multi banded tiff

texts:
I have two geoTIFFS. The files cover the exact same area, both have 1 band(Block = 256x256) and both have a no data value of 0. They also have the same coordinate system.

I want to merge the files into a single tiff with two bands.

I tried this command:

gdal_merge.bat -n 0 -a_nodata 0 -separate -of GTiff -o merged.tif first.tif second.tif

But I kept getting a Memory error. After researching online, I then used the same command with the -createonly flag:

gdal_merge.bat -createonly -n 0 -a_nodata 0 -separate -of GTiff -o merged.tif first.tif second.tif

Followed By:

gdalwarp --config GDAL_CACHEMAX 500 -wm 500 merged.tiff mergedFINAL.tiff

However opening mergedFINAL.tiff in qgis just yields a blank screen. Can anyone see something Im doing wrong and offer me some advice?


texts: Merging multiple point shapefiles that have same coordinates but different years

texts:
I am using ArcGIS Pro.
I have 35 shapefiles of point data, with each layer representing a different year of data collection.  Each layer has points with the same coordinates.  I want to bring these layers together in order to export a single attribute table with the time series data.
Would a spatial join or a merge be the appropriate tool?
Is there a way to bring all 35 layers together at once?


texts: Detecting ring direction of a polygon using ArcPy

texts:
I am using ArcPy in ArcGIS 10.1 to process some geometry using the arcpy.da data access module.

However I have run into a problem with complex Polygons that contain inner rings/holes.

From the documentation in ArcGIS and this question I can detect rings via the None objects in the Point array.

But now I need to know whether these rings are clockwise or counter-clockwise so that I can tag them as inner or outer.

Is there an ArcPy function I can use to detect the rings' directions?


texts: QGis Grid statistics for polygons error - Module caused access violation!

texts:
I'm trying to calculate statistics for a raster layer using an overlaying polygon layer, I'm using SAGA "Grid Statistics for polygons" function, but it returns an error, sadly.. So in the sextante log it says :

Parameters

Grids: No objects
Files: "Y:/.../clipper_grassexport_0to11"
loading: Y:/.../clipper_grassexport_0to11
Driver: Y:/.../clipper_grassexport_0to11/
Cells: x 44069, y 42309
Bands: 1
Transformation x' = 2635889.298915 + x * 100.001118 + y * 0.000000
Transformation y' = 5421705.687421 + x * 0.000000 + y * -100.000407

Module caused access violation!

GDAL: Import Raster: Module caused access violation!
error: executing module [GDAL: Import Raster]

***

go...
Load grid: C:\Users\...\AppData\Local\Temp\sextante\1368631882.798.sgrd...
failed
error: Grid file could not be opened.
error: input file [C:\Users\...\AppData\Local\Temp\sextante\1368631882.798.sgrd]
error: empty input list [GRIDS]
Usage: Grid Statistics for Polygons -GRIDS <str> -POLYGONS <str> [-RESULT <str>] [-COUNT] [-MIN] [-MAX] [-RANGE] [-SUM] [-MEAN] [-VAR] [-STDDEV] [-QUANTILE <num>]
-GRIDS:<str>    Grids
Grid list (input)
-POLYGONS:<str> Polygons
Shapes (input)
-RESULT:<str>   Statistics
Shapes (optional output)
-COUNT          Number of Cells
Boolean
-MIN            Minimum
Boolean
-MAX            Maximum
Boolean
-RANGE          Range
Boolean
-SUM            Sum
Boolean
-MEAN           Mean
Boolean
-VAR            Variance
Boolean
-STDDEV         Standard Deviation
Boolean
-QUANTILE:<num> Quantiles
Integer
Value Range: 0.000000 - 50.000000

error: executing module [Grid Statistics for Polygons]

I have tried to load my raster and vector layers into SAGA directly, but the raster layer wouldn't load, and same error appears - "Module caused access violation!"

Checked if SAGA is properly installed, and is 2.0.8 version (QGis 1.8.0)

Raster is an GEOTIFF exported from GRASS Tools in QGis in Byte format (Values 0-255), CRS's also match, I have run out of ideas, what could possibly be wrong?

Update - June 17 2013

So to sum it all up - with a large dataset it still gives an error, but if I want to divide the raster in say 4 smaller parts, raster Clipper Tool does it for ages and freezes with an approximate filesize of 100-140 MB, that means I have found no solution yet, that does the job efficiently.

Suggestions?


texts: Changing color for Shapefile in Google Earth?

texts:
I have a shapefile that desrbibes the post/prior routes from/to a certain venue. The white lines are the routes and each segment has a different foot traffic value. The segment in the image has a value of 5923 for example.

Is there any way to edit this shapefile in Google Earth in such a way that each value will be of a different color?

Like in a gradient, or a certain color would fit a range, for example <6k visits is blue, 6k-10k visits is green and so on.


texts: How to show combine results or export them with PyQGIS

texts:
I'm quite new in PyQGIS and I think PyQGIS could make my work easier. I first want to try with a simple code, here is what I want:
merge some features and then export them in a shapefiles.

For the combine parts, it's OK, I have this code:

layer=qgis.utils.iface.activeLayer()
geomfus = QgsGeometry.fromWkt('GEOMETRYCOLLECTION EMPTY')
for elem in layer.getFeatures():
      geomfus = geom.combine(elem.geometry())

it is working, no problem. But now I don't know what to do to see this results in QGIS or to export the combine result in a shapefile.


texts: Inconsistent map exporting using the composer tool in QGIS

texts:
I am using QGIS Desktop 2.8.5 on a windows machine and I am using the composer tool to create images with a scalebar and legend. The following is a screenshot of the my working environment:

The image above appears exactly how I would like.

Now, when I export as a PDF, the resulting image is:

And when I export as an image (tiff, jpg, png), the resulting image is similar (my reputation is <10 so I cannot upload a third image).

In both cases, the exported image is distorted from the original. I've tried to export the image several times and the amount of distortion varies each time. I have also tried to use a different computer (also Windows with 2.8.5) and I get the same results.

Why is this happening and how can I remedy it?


texts:
I am creating a tool which the user will select point or line features from the map and it will display their names in a list then after clicking on any of these features from the list, it will show the intersected polygons in another list.

The problem that I don't know how to do a SpatialFilter to do this intersecting functionality, which it will be done like the one in the Select By Location tool.


texts: Real-Time GPS Tracker with CartoDB and firebase

texts:
is it possible to track the gps tracker in cartodb in real time?

I saw that cartodb has table sync but this sync only every hour but i need something that can sync every 30 sec.


texts: Watersheds in SAGA Fill Sinks Wang & Liu tool

texts:
I'm delineating some coastal subbasins using open source software. I have asked a similar question before (Delineating coastal subbasins using SAGA GIS or GRASS in QGIS?). I found that using SAGA Fill Sinks (Wang & Liu) in QGIS creates watersheds according to what I expect. the problem is I don't find any documentation about how this tool creates the watersheds.

I looked for in the tool documentation (Module Fill Sinks (Wang & Liu)), but it doesn't say how the tool works. Most of the tools used in SAGA, GRASS and other softwares create watersheds based on a previously created channel network, but here the only input is a DEM.

Does anyone know how this tool create watersheds?


texts: Concatenating strings from two fields using Python in ArcMap Field Calculator?

texts:
I'm new to Python scripting and I don't understand what is wrong with this code. I want to concatenate the string values from GRID1MIL and GRID100K. I've created a text column MGRS to begin with.

I also tried writing the following but it did not work:

Using VB Script I was able to do it easily but I want to learn how it's done with Python.


texts: CreateTable_management as part of a script tool fail

texts:
I am trying to create a table as part of a script tool, and it always fails with error 999999.

The general idea of the script is that it is calculating proportional coverage of things in a series of buffers, and the tables are needed to compile the results.

The code I'm using creates the tables and works fine in PyScripter and also works fine in the Python window within ArcMap, but will not run as a tool and fails at the create tables part.

I'm not sure if part of the problem is that I'm trying to access the path names and files via the GetParamtersAsText function or not. I have tried countless ways of naming the outputs in the CreateTable_management function. I'm not sure what I am missing here (maybe I need output parameters).

Here is my code (the relevant parts of creating the tables):

import arcpy
import os

#set up workspace location
Working = arcpy.GetParametersAsText(0) #user provided geodatabase location
arcpy.env.workspace = Working
arcpy.env.overwriteOutput = True

#set up link to buffer files.
Buff50_R = arcpy.GetParameterAsText(3) #user provided: let's pretend file is called 'RightBufferfile'
Buff50_L = arcpy.GetParameterAsText(4) #user provided: let's pretend file is called 'LeftBufferfile'

#Loop through and make table based on the name of each of the buffers for each metric of interest
#When I use PyScripter or the Python window in ArcMap the tables are created in the working
#geodatabase and are named 'InC1RightBufferfile', InC1LeftBufferfile', etc...
fcList =[Buff50_R, Buff50_L]
for eachFC in fcList:
    NatVeg_Table = "InC1" + eachFC ##Blank file for any natural veg within each buffer
    arcpy.CreateTable_management(Working, NatVeg_Table)
    arcpy.AddField_management(NatVeg_Table, "BID", "TEXT", 12)
    arcpy.AddField_management(NatVeg_Table, "SUM_Per_Cov","DOUBLE", 20, 4)

    NatForest_Table = "InC2" + eachFC ##Blank file for any forest within each buffer scale
    arcpy.CreateTable_management(Working, NatForest_Table)
    arcpy.AddField_management(NatForest_Table, "BID", "TEXT", 12)
    arcpy.AddField_management(NatForest_Table, "SUM_Per_Cov","DOUBLE", 20, 4)

    HumanDev_Table = "InC3" + eachFC ##Blank file for any Human Disturbance within each buffer scale
    arcpy.CreateTable_management(Working, HumanDev_Table)
    arcpy.AddField_management(HumanDev_Table, "BID", "TEXT", 12)
    arcpy.AddField_management(HumanDev_Table, "Per_Cov","DOUBLE", 20, 4)

I've tried using Delete_management to delete the table ahead of creating the table with no luck, and have tried commenting out the first table, but get then get the same error starting with the NatForest_Table.

Here is the print out from the Results window for extra info:

Here are the parameter configurations; the buffer files correspond to parameters 3 and [4], and all the input feature classes are configured the same way.


texts: Avoid a for loop in raster::extract(rst,shp)

texts:
I am working on R to extract the mean and maximum value of a raster within a 3 meter buffer of some buildings.
For this, I have created a for loop that iterates iterates through each building to extract this two values. My current code looks as follows:
for (b in c(1:nrow(buildings_shp))){

    building <- buildings_shp[b,]

    buffered <- st_buffer(building, 3)

    raster_cropped <- crop(raster, extent(buffered))

    mean <- extract(depths_cropped, buffered, fun = mean, na.rm = TRUE)
    max <- extract(depths_cropped, buffered, fun = max, na.rm = TRUE)

    buildings_shp[b,"mean"] <- mean
    buildings_shp[b,"max"] <- max

  }

This loop, however, takes a considerable amount of time (~17 minutes for 1500 buildings), and the step that seems to take most time is the two extract lines. I would like to know if there are ways to speed up this process by:
a) avoiding the use of a loop - the reason for this loop is that I fear that if I use st_buffer on the entire dataset, then when buildings are closer than 3 meters I would generate overlapping geometries, which may cause an error. UPDATE - having tested this, I confirmed that the results produced are different when using overlapping geometries than when iterating through each one of them
b) parallelizing the for loop (i have tried the raster clustering feature, but it did not speed up the process, probably because it did not parallelize the loop itself but the extract function)
c) using other function than raster::extract. I have seen some posts recommending the velox package, but it seems like this package has been removed from CRAN.
Some dummy data (copied from the referenced question above)
library(raster)
library(sf)

raster <- raster(ncol=1000, nrow=1000, xmn=2001476, xmx=11519096, ymn=9087279, ymx=17080719)
raster []=rtruncnorm(n=ncell(raster ),a=0, b=10, mean=5, sd=2)
crs(raster ) <- "+proj=utm +zone=51 ellps=WGS84"

x1 <- runif(100,2001476,11519096)
y1 <- runif(100, 9087279,17080719)

buildings_shp <- st_buffer(st_sfc(st_point(c(x1[1],y1[1]), dim="XY"),crs=32651),200000)


texts: Approaches to delete Polygon overlap in PostGIS?

texts:
Is there function to delete the areas of overlap between polygon features to achieve visualized results?

Id like to check for each geometry (row) if a geometry is intersected, or in other words if a geometry is overlapped by another. If this is the case Id like "erase" the overlap from the other geometry without creating a hole. When there is no direct way to delete the overlap, is it possible to create new features from the overlaps, and then filter them out by a certain size treshold?

Maybe there is a simple way I did not see but please excuse PostGIS is new to me.


texts:
I have installed QGIS via osgeo4w 32 bit yesterday, to take pyqgis api in standalone script, but I have some problems.

I have python 2.7.13 32bit and pyqt4 installed from Gohlke. After, I created a file (qgis.pth) and I put it in python27/Lib/site-packages to get the path for QGIS api.

Paths are :

C:\OSGeo4W\apps\Python27
C:\OSGeo4W\apps\Python27\libs
C:\OSGeo4W\apps\Python27\Lib
C:\OSGeo4W\apps\Python27\DLLs
C:\OSGeo4W\apps\qgis\bin
C:\OSGeo4W\apps\qgis\grass\bin
C:\OSGeo4W\apps\qgis\grass\modules
C:\OSGeo4W\apps\qgis\grass\scripts
C:\OSGeo4W\apps\qgis\i18n
C:\OSGeo4W\apps\qgis\python
C:\OSGeo4W\apps\qgis\python\console
C:\OSGeo4W\apps\qgis\resources
C:\OSGeo4W\bin
C:\OSGeo4W\bin\gdalplugins
C:\OSGeo4W\bin\osgPlugins-3.1.7
C:\OSGeo4W\etc\ini
C:\OSGeo4W\etc\postinstall
C:\OSGeo4W\etc\preremove
C:\OSGeo4W\include\geos
C:\OSGeo4W\lib

(and if I add the paths from QGIS pyqt4 and sip I had some errors)

The problems are:

First test :

import qgis.core

Error :

import qgis.core
  File "C:\OSGeo4W\apps\qgis\python\qgis\__init__.py", line 26, in <module>
    from qgis.PyQt import QtCore
  File "C:\OSGeo4W\apps\qgis\python\qgis\PyQt\QtCore.py", line 38, in <module>
    from PyQt4.QtCore import QPyNullVariant
ImportError: cannot import name QPyNullVariant

Second test :

from PyQt4.QtCore import *
from PyQt4.QtGui import *
#from PyQt4.QtSvg import *
#from PyQt4.QtCore import pyqtSignal as Signal
#from PyQt4.QtCore import pyqtSlot as Slot
import qgis.core

error :

sip.setapi(api, 1)
ValueError: API 'QDate' has already been set to version 2

and finally, I tried this code :

import sip
API_NAMES = ["QDate", "QDateTime", "QString", "QTextStream", "QTime", "QUrl", "QVariant"]
API_VERSION = 2
for name in API_NAMES:
    sip.setapi(name, API_VERSION)
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from PyQt4.QtSvg import *
from PyQt4.QtCore import pyqtSignal as Signal
from PyQt4.QtCore import pyqtSlot as Slot
import qgis.core

error if version is 1:

from PyQt4.QtCore import QPyNullVariant
ImportError: cannot import name QPyNullVariant

error if version is 2:

from PyQt4.QtCore import QPyNullVariant
ImportError: cannot import name QPyNullVariant

Any idea? If I remove the import for QGIS then I don't take some error with imports for sip ot pyqt. It's only if I try to import PyQGIS.


texts: Why does the OpenLayers plugin override the predefined projection?

texts:
What I am trying to say is:

You import a shapefile into an empty QGIS and the QGIS environment gets the projection system of the imported shapefile.
You import another shapefile into QGIS. The projection system assigned earlier by the previous shapefile remains in place for the QGIS environment.
But, when you import a basemap using the OpenLayers plugin, the environment gets the projection system of the imported basemap.

Why doesn't QGIS just follow one standard by keeping the projection system that was assigned by the first layer imported into its environment?

Note: I have tried this with QGIS 2.10.1 only.


texts: symbols for different polygon layers in ggplot2

texts:
I'm trying to plot the boundary of different spatialPolygonsDataFrames in ggplot2 or in rasterVis gplot(), but I just can't find the way, they tend to overwrite the previous one; this is something quite common in mapping, having different symbols for each layer, it's something like this what I'm trying to achieve (done with base graphics)

library(raster)
mex <- raster::getData(country="MEX", level=1)
munis <- raster::getData(country="MEX", level=2)
plot(munis[munis$NAME_1=="Veracruz"|munis$NAME_1=="Oaxaca",], border="green")
plot(mex[mex$NAME_1=="Veracruz"|mex$NAME_1=="Oaxaca",], add=T, lty=1)
legend("topright", c("state boundary", "munic. boundary"), pch=14, col=c("black", "green"), bty="n")


texts:
I have >500 very large KMLs that contains only point data and a reference to a .png pop-up

Each KML looks like this

<kml xmlns:xsd="http://schemas.opengis.net/kml/2.2.0/ogckml22.xsd" xmlns:xmlns="http://www.opengis.net/kml/2.2/" version="1.0">
  <Document>
    <name>test.kml</name>
    <open>1</open>
    <Folder>
      <name>TS</name>

      <Placemark>
        <name>1</name>

        <description><![CDATA[]]></description>
        <Point>
          <extrude>1</extrude>
          <altitudeMode>clampToGround</altitudeMode>
          <coordinates>-117.766574236,56.204770715,0</coordinates>
        </Point>
      </Placemark>

    </Folder>
  </Document>
</kml>

The file above contains only one point, but most of my files contain more than 1 million equally spaced points.

I can do the regionation in Google Earth Pro, but for >500 KMLs, it will take very long time and effort.

Is there a way to do similar regionation in R?

How to define the dimensions of a regions and how to cascade them in R?

I used R to write the KMLs in the first place using the following code:

  write_header <- function (outfile) {

  cat('<kml xmlns:xsd="http://schemas.opengis.net/kml/2.2.0/ogckml22.xsd" xmlns:xmlns="http://www.opengis.net/kml/2.2/" version="1.0">\n')
  cat("  <Document>\n")
  cat(paste0(
      "    <name>",outfile,"</name>\n"))
  cat("    <open>1</open>\n")
  cat("    <Folder>\n")
  cat("      <name>TS</name>\n")
}

write_styles <- function(df){

  for (j in 1: nrow(df)){
    cat(paste0(
        '      \n")
  }
}

write_placeholders <- function(df){

  for (k in 1: nrow(df)){
    row=df[k,]
    x=row$x
    y=row$y
    cat("      <Placemark>\n")
    cat(paste0(
        "        <name>", k, "</name>\n"))
    cat(paste0(
        "        \n"))
    cat(paste0(
        '        <description><![CDATA[]]></description>\n'))
    cat("        <Point>\n")
    cat(paste0(
        "          <extrude>",k,"</extrude>\n"))
    cat("          <altitudeMode>clampToGround</altitudeMode>\n")
    cat(paste0(
        "          <coordinates>",x,",",y,",",0,"</coordinates>\n"))
    cat("        </Point>\n")
    cat("      </Placemark>\n")
  }
}

write_end <- function(){

  cat("    </Folder>\n")
  cat("  </Document>\n")
  cat("</kml>\n")
}

writ_TS2kml <- function(outDir,df, outfile) {

  setwd(outDir)
  tryCatch({
    sink(outfile)
    write_header(outfile)
    write_styles(df)
    write_placeholders(df)
    write_end()
    sink()
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})}
}


texts: Google Maps, HERE maps data in one database

texts: I want to create a database from HERE and Google Maps to combine with my data. Then I create an application that can use this database. Basically I want to host the server with Google and here map data combine with my data. How should I proceed? Is it possible to combine my data in to these maps data.


texts: Obtaining longitude and latitude from GeoTIFF with geotools

texts:
Continuing what is written in this post, I have another question:

for each point (X, Y) of my GEOTIFF file I extract the coordinates as follows:

Envelope2D pixelEnvelop = geometry.gridToWorld(new GridEnvelope2D(x, y, 1, 1));

double[] coords = new double[]{pixelEnvelop.getCenterX(), pixelEnvelop.getCenterY()};

For example I obtain

-577460.1192321777, 640293.6331481934
-577460.1192321777, 639293.9553527832
... (and so on)

However, I would like to have latitude and longitude available. Is it possible to convert them?

For information, I print also the content of pixelEnvelop.getCoordinateReferenceSystem() method

PROJCS["unnamed",
  GEOGCS["WGS 84",
    DATUM["World Geodetic System 1984",
      SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]],
      AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
    UNIT["degree", 0.017453292519943295],
    AXIS["Geodetic longitude", EAST],
    AXIS["Geodetic latitude", NORTH],
    AUTHORITY["EPSG","4326"]],
  PROJECTION["Transverse_Mercator"],
  PARAMETER["central_meridian", 12.5],
  PARAMETER["latitude_of_origin", 42.0],
  PARAMETER["scale_factor", �],
  PARAMETER["false_easting", 0.0],
  PARAMETER["false_northing", 0.0],
  UNIT["m", 1.0],
  AXIS["Easting", EAST],
  AXIS["Northing", NORTH]]

As you can see, scale_factor parameter does not contain a readable value.


texts: Associating Inset Map with Layer using ArcPy?

texts:
I'd like to only display an inset map, inset map scale bar, frame and leader if an associated layer is displayed. Similar to Associate Point Layer, Annotation Group and Line Layer in Map Legend but for an inset map rather than annotation.

I can't put up an image through the firewall but the inset is the same as the main map at a different scale in ArcGIS 10.2.1 for Desktop under Win7. It has a few line and point layers, same labels and the scalebar.

Following on from https://gis.stackexchange.com/a/311066/121893 and using a framework from https://community.esri.com/thread/178874. This arcpy code has no affect on inset or scalebar. How can I get it to work (including where to put code in project)?

def display_layer_inset(glayername):
    mxd = arcpy.mapping.MapDocument('CURRENT')
    df = mxd.activeDataFrame
    df2 = arcpy.mapping.ListDataFrames(mxd, "inset")[0]
    layers = arcpy.mapping.ListLayers(df)
    for lyr in layers:
        if lyr.isGroupLayer and lyr.name == glayername:
            print "Layer: " + glayername
            if lyr.visible == True:
                print "visible"
                df2.visible = True
            else:
                print "not visible"
                df2.visible = False

# use as
display_layer_inset("MyGroupLayer")

Output:

Layer: MyGroupLayer
not visible


texts: How to resample a raster to higher resolution in QGIS?

texts:
I'm fitting a SDM using R. A requirement is that all environmental predictor layers have the same extent, origin and resolution.

I have high resolution bathymetry data (25m) that I will resample to 100m (precision of occurrence records). However, some environmental predictors pulled from a hydrodynamic model are at a much lower resolution (~7km) e.g. temperature / salinity.

I am aware of the introduction of error in resampling from low -> high resolutions. However, resampling the from 25m -> 7km is not suitable for my requirements.

Is it permissible, using QGIS, to define the resolution of the low resolution layers by saving the resolution to that of the higher res layers? Please note I am not resampling - just changing the number of columns & rows in the saved tiff so that R sees the layer as having an equal resolution to the other layers.


texts: Spatial Clustering With Roughly Equal Sizes

texts:
I'm not sure which plugins or tools I'd want to further research to achieve my goal.

I have a layer which breaks my geography into regional blocks and an attributes table with a column for number of people per regional block.  I want to automatically assign those regions into subsets of approximately equal sized number of people, and have those subset groups be geographically near each other.

Example:
Springfield has 100 regional blocks, each ranging with people counts from 200 to 600 or so per regional block.  I'd want to say "give me sets of several regional blocks each with around 1000 people total per set, where the regional blocks in the set that are near each other." So some sets may only have 2 regions, other sets may have 6 regions, etc.  I can then assign those sets to certain employees to drive out to that side of town, work in those few regions, etc.  So if "Set A" were to have one region in the far west of town and another region on the far east of town, that isn't very helpful because of all the travel time.  So I'd want "Set A" to be regions kind of near each other, ideally touching.

Again, I realize this is probably a multi-part project, but I don't even know what this type of work is called or what QGIS tool I'd want to learn how to use to achieve my goal.

Where should I start my research?


texts: Looking for shapefiles for Canadian shale plays with labels and/or meta data about the plays

texts: I work for a company in the oil and gas industry.  We are making a geospatial map with mapbox and leaflet.  Our task is to show shale plays in relation to other industry items like well sites.  We have completed the work in the U.S. and gathered most of our data from RigData and EIA.gov; these are both US based and do not contain other countries in the dataset. I have searched the BC Ministry of Natural Gas; BC Ministry of Energy and Mines; Alberta Geological Survey; http://www.nrcan.gc.ca/ and googled several combinations of "Canada/ian" and "Shale play/s".  So far, I cannot find a shapefile that says [shale play name] shapefile. The shape files should be able to look like this map.  Can someone help point me to a dataset that can reproduce the shale plays and names?


texts:
I'm making a graduated style map and I would like to add a class that would represent null values.  I can add a class for zero (0) value... meaning some of the areas reported zero report, but some areas don't have reports (i.e. null values).

How can I add a class for null values?


texts: Loading intersections consolidated graphs

texts:
To get the right connected components, I projected a graph then use osmnx.simplification.consolidate_intersections to consolidate intersections.
import osmnx as ox
import pandas as pd
import networkx as nx
%matplotlib inline
cities = {'Copenhagen':'Copenhagen Municipality, Denmark'}

G1=ox.graph.graph_from_place('Copenhagen Municipality, Denmark', simplify=False, retain_all=True, truncate_by_edge=True, which_result=None, buffer_dist=None,clean_periphery=True, custom_filter='["highway"~"cycleway"]')

G11=ox.simplification.consolidate_intersections(ox.project_graph(G1, to_crs='EPSG:3812'), tolerance=15, rebuild_graph=True, dead_ends=True, reconnect_edges=True)
G2==ox.project_graph(G11, to_crs='EPSG:4326')

ox.io.save_graphml(G2, filepath='data/Copenhagen_bike.graphml', gephi=True, encoding='utf-8')
print('Copenhagen bike_net saved')

I saved G1 as GraphML file, but when I tried to load Copenhagen_bike.graphml, it would come error:
ValueError: invalid literal for int() with base 10: '14-1' ".

I wonder how to load the file to get a graph with no error.


texts: QGIS 3.4. - page numbers in print layout

texts:
I'm using QGIS print layout to create an atlas with around 180 features. I'm trying to find a way to add page number in print layout.

I tried adding text label with [% @atlas_featurenumber %] to show page number, but I have multiple pages for each feature (and they are added dynamically depending on the amount of data in tables), so it does not work well (it shows feature number instead of page number because features have multiple pages).

Is there a direct way to insert page number to print layout?


texts: Dissolve by group and get the lowest or highest in other field

texts:
I've got the next table of attributes:

I want to perform a dissolve by the first column (date) and get:
- Earliest value by group in the second column.
- Latest value by group in the third column.

There is any easy way to do so ? The aim is to do so into the modeler so I would try to avoid complex code, or at least an easy way to implement it.


texts: Trouble changing 'unknown' CRS of vector layer in QGIS

texts:
I'm new to GIS. I've been given a shapefile that I've added as a new vector layer in QGIS and it displays OK, but it says it has an Unknown CRS. The CRS from the prj file is:

PROJCS["Transverse_Mercator",GEOGCS["GCS_Airy 1830",DATUM["D_unknown",SPHEROID["airy",6377563.396,299.3249753150316]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",49],PARAMETER["central_meridian",-2],PARAMETER["scale_factor",0.9996012717],PARAMETER["false_easting",400000],PARAMETER["false_northing",-100000],UNIT["Meter",1]]

and the one displayed in the details window of the CRS selector is:

Unknown CRS
WKT
PROJCRS["Transverse_Mercator",
    BASEGEOGCRS["GCS_Airy 1830",
        DATUM["D_unknown",
            ELLIPSOID["airy",6377563.396,299.324975315032,
                LENGTHUNIT["metre",1,
                    ID["EPSG",9001]]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",49,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-2,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996012717,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",400000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",-100000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
Proj4
+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +a=6377563.396 +rf=299.324975315032 +units=m +no_defs +type=crs
Extent
Extent not known

I'd like to convert it to a well known CRS within QGIS but when I do select a different CRS from the options, all the features disappear, with the exception of EPSG:3994 which displays the right map but inverted N/S. Does anyone have any advice on what I'm doing wrong or what CRS I should be changing it to?


texts: ArcPy script extremely slow with Search and Update Cursor

texts:
I have developed a Python script to automatically retrieve the maximum height in a certain distance around my input polygon feature class. The height should be calculated for each 5 meters.
For a small input with around 50 polygon features, the script takes around 5 minutes. Whenever I input a really large dataset (10.000 features), the script gets extremely slow and after hours interrupts.
My (pseudo) code is as follows:
with arcpy.da.SearchCursor(fc, ['FID']) as rows:
        for row in rows:
            my_fid = row[0]
            expression = '"FID" = ' + str(my_fid)
            fc_tmp = r"fc.shp"
            arcpy.Select_analysis(fc, fc_tmp, expression)

            output_tmp = r"output_tmp.shp"
            arcpy.management.CreateFeatureclass(r"", output_tmp)

            -- pseudo code
            # 1) Buffer  feature class
            # 2) Clip DGM to raster
            # 3) convert feature class to polyline
            # 4) Make points along polyline
            # 5) make tmp result feature class
            # 6) Buffer points
            # 7) get max height in buffer around points using zonal stats
            # 8) join max height with tmp result feature class

            with arcpy.da.UpdateCursor(output_tmp, ["FID"]) as recs:
                for rec in recs:
                    rec[0] = my_fid
                    recs.updateRow(rec)
                del rec, recs

            arcpy.management.Append(output_tmp, r"myResult.shp", "NO_TEST")
            arcpy.management.Delete(output_tmp)
            arcpy.management.Delete(fc_tmp)
        del row, rows

Can it have to do something with the cursors?
I thought maybe the selection takes long because the feature has to be selected from 10.000 features, but the selection is done very fast.
I do not understand that when I input 10.000 features, also the first 50 features are not done within 5 minutes, although when I input onyl 50 features, it does not take longer than 5 minutes. I just tested it, the first 3 features are done within 20 minutes. Thats absolutely disproportionate.
Any ideas?


texts: query syntax to select every multiple of nn record (shapefile)

texts:
I have data like multiples of 40.
40
80
120
160
200
and so on.

I would like to label every record that is a multiple of 100.  I have tried a couple of methods and can't get them to work either as a query in the label or as a select statement.

"STATION" =  %00
RIGHT ("STATION", 2) =  00

I get invalid SQL statement.
Does anyone see what I am doing wrong?

Respectfully I did search and did not see any of the example questions you presented.  However it still does not answer my question, re-stated, what am I doing wrong?


texts: Displaying measured grids in ArcGIS Pro 1.4?

texts:
Is it really not possible to display reference grids in an ArcGIS Pro 1.4 layout view ie Eastings & Northings?

At the moment there only seems to be an option for graticules. The ESRI support webpage (http://support.esri.com/technical-article%5C000012529) mentions that it's not been implemented yet, it seems strange to not include this as an option.


texts:
I have a very large image which I want to display using OpenLayers so visitors can zoom in and pan. I have tiled it (1 image for the first layer, 4 for the second, 9 for the third, etc.) and I am trying to display it using a TMS layer.

Following this information I have it working more or less. However, on some zoom levels, I can not pan all the way to the right and the image seems cut off (even though the tiles exist).

I have set my image height and width in px as the layer's maxExtent.

Is there a better way to do this?

Thanks.


texts: ArcScene can't symbolize LAS points by classification code or return number

texts:
I'm using ArcScene 10.3 and loaded a LAS file into a LAS Dataset. When I look at the LAS files statistics I can see that there are several classification codes and 7 different returns.

Everything looks OK when I symbolize the LAS dataset points by Elevation. One is able to identify trees, buildings, silos a windmill and other things:

It is also possible to set the Filters to "Ground" and fade out the buildings, trees and the windmill.

The strange thing is, that I only get yellow points when I symbolize the LAS dataset points by Classification (Filters set back to "All"):

When I symbolize the LAS dataset points by Return I even end up with no points at all.

Why is ArcScene able to filter ground points but can't visualize the different returns and classifications?


texts: Summarizing GPS coordinates inside polygon in Google My Maps?

texts: I need to set up the most effective routes to visit different locations. I have several teams to visit the locations depending on where they are. I am given Latitude and Longitude of all the locations. What I am currently doing is upload all the locations from an Excel file to Google My Maps and give them a number sequence. The number sequence is in the same order as the list in the excel file and the GPS points are shown in the map as circles with the number inside, for example if there are 100 locations they are 100 circles with numbers from 1 to 100 each one representing an item from the list in Excel. After I have all the points represented in the map I decide how to organize them. If I have 100 locations and 5 teams, each team has to visit 20 locations. I then draw a polygon on top of the map and the GPS points like a fence to divide the points according to proximity, so when I finish I have 5 poligons with 20 locations each. Here comes the tricky part.. I haven´t found a way Google Maps export what GPS points are inside each polygon, so what I do is to write each number in a column in Excel and in the colum next to it I write the polygon it is inside, then with Vlookup I find in the original Excel file which polygon each number belongs to. This is very difficult when there are 500 or more GPS points.


texts: Geolocation of a drawn point on printed map without scanning

texts:
I have a printed map of a district with several points on it that can be used as geolocation reference points.

I have the shapefile of this map in QGIS.
If someone draws a point on the printed map, how can I find the position of the drawn point without scanning the printed map ?

I was thinking of measuring the distance from that point to three other points on the paper, so I can have the zoom used and find the position. Is there a plugin to help with that? Are there other ways to accomplish this?


texts: Calculating distance features within a shapefile - ArcGIS

texts:
I have a file that has a list of postal codes of customers and the main store they frequently visit (e.g., Customer 1, Customer 2, and Customer 3, visit Circle Centre). In ArcGIS I want to be able to calculate the distance between customers and their frequently visited stores. Has anyone ever done this and know what tool to use?

I know there is a tool that assigns points to the nearest feature, however, that's not what I am trying to achieve.


texts:
I am trying to reverse engineer a consultants work. They are using a postgis spatial view as a source layer in MapServer:

LAYER
            NAME 'Parcels2009'
            CONNECTIONTYPE postgis
            DATA "
                    parcel_geometry FROM (
                            SELECT id, year, parcel_working_title, appropriation_title, legal_citation, link_to_accomplishment_plan, project_type, acres, county, township, range, section, description, easement_ohf, fee_ohf, x, y, parcel_geometry
                            FROM map_view_final
                            WHERE year = '2009'
                    ) AS foo USING UNIQUE id USING srid=4326
            "

I can view the table in postgres:

And following this question I can add the layer into QGIS:

But the layer has no data/empty attribute table. If I add the layer with the "Spatial Type" set to 'NoGeometry' it will load the table just fine.

I suspect the issue is with my geometry column.

When I try to the following pgsql2shp command I get the following error:

pgsql2shp -f parcels -h localhost -u ####-P #### #### "SELECT * FROM map_view_final where parcel_geometry is not null"

Initializing...
Done (postgis major version: 2).
Output shape: Null Shape
Dumping: XUnknown WKB type (7) for record 0

I should note that I have similar issue trying to use ArcGIS to accomplish my task:

In the end I just need a shapefile to QC the data and for use in a separate project.


texts: Converting PBF to Shapefile with name:xx tags as a column in the attributes

texts:
I have a PBF file containing only rivers information (extracted using osmosis when waterway=river).

A conversion of this PBF to SHP works fine with ogr2ogr when I am selecting fields like "osm_id" and "name". However, I would like to have the name of the rivers in different languages in the attribute table. The available translations appears in the "other_tags" column like this example: "name:en"=>"The Thames","name:fr"=>"La Tamise".

Is it possible to create a SHP (or any other format) where the different translation of the rivers name would appear in different columns (eg. name:en, name:fr)?

This gdal command:

ogr2ogr -f 'ESRI Shapefile' -sql 'select osm_id,name,"name:en" from lines where waterway="river"' rivers.shp rivers.osm.pbf --config ogr_interleaved_reading yes -lco ENCODING=UTF-8

returns the following error: ERROR 1: Unrecognized field name name:en


texts: Finding script location dynamically when using ArcPy?

texts:
I have created a script that uses Esri's arcpy site package and it functions properly in a stand-alone environment. I am now trying to integrate this script into a toolbox to be run within an open ArcMap session and I am having trouble. In integrating this script I am using the tool structure that ESRI suggests... http://goo.gl/eesrV . The script I have created has information in the 'ToolData' folder that I need to access to run the script.

How do I locate the script's directory dynamically when running it in an open ArcMap session?

In my stand-alone script I use 'os.getcwd()'; this works because it returns the proper location of the script and I then know the structure around it. Whereas when incorporated as a script in a tool box on run within an ArcMap session 'os.getcwd()' returns the directory of the map document while the toolbox, script, and tool data are located where ever the user installed them and my script can't find this structure.


texts:
Can someone help me please, I'm searching for a way to highlight a Feature Class in the TOC.
I tried so many snippets and I didn't get the result I wanted

So what I'd like to have is a function written in vb.net in which I pass a string that is the name of the Feature Class that shall be highlighted with this blue rectangle, basically just like when a user would make a selection in the TOC


texts:
I try to measure a polygon with Export/Add geometry columns tool. The units for area measurement are square meters.

The problem is that I get different area results for the polygon:

When the ellipsoid is WGS 84- the area result is 1809 m2 (square meters)
When I choose calculate using Layer CRS in the Export/Add geometry columns tool, then the area result is 1799 m2

I don't understand:

What do I choose in the Ellipsoid field (in General tab)?

What do I choose in the "Add geometry columns" tool?

Which of the results is the right for me?


texts:
Im loading a shapefile using the method QgsVectorLayer() and adding a style to it using the method loadSldStyle() and the path to the sld file. So far so good.

What I want now is to export that style as .qml file. I suppose that maybe I could use the method exportNamedStyle() but this method use an argument of type QDomDocument. I have no idea how to create this this object.

Can anyone tell me how to export the .qml file of this loaded layer using python?

Using the method saveNamedStyle() it only works if I use the python console inside qgis.
In the python console in the qgis desktop I did this:

vlayer = iface.addVectorLayer("/home/inesf/sag/prodgen/SAGautodataset/Cystoseira_compressa.shp", 'mylayer', 'ogr')
vlayer.loadSldStyle("/home/inesf/sag/prodgen/SAGautodataset/Cystoseira_compressa.sld")
vlayer.saveNamedStyle("/home/inesf/sag/prodgen/SAGautodataset/Cystoseira_compressa.qml")

And it worked fine, it exported a .qml file with the same style defined in the SLD file. Perfect.

When I do this in the python console outside the QGIS desktop, it exports a .qml file but empty.

When I run a python program with this code with the same code, does not give errors, but nothing happens.


texts: OSM: Getting the correct Overpass Turbo query using a keyword (and values)

texts:
I'm completely new to Overpass queries for getting OSM data, so sorry for this maybe trivial question.

I'm looking for the correct query using Overpass turbo. I have found some demo queries, but I am not sure if they are correct for my needs, so I would like to be sure.

My aim is to get all Points of Interests, which are using the KEY: "tourism" in a specified region (a bounding box with geographic coordinates, or a specified district, like a country).

I know I can get all data in a specified region with this query:

(node(51.249,7.148,51.251,7.152);<;);out meta;

but how do i include the keyword, that only POIs with they key: tourism are returned?

2nd question: how could i expand this query to receive only specific values of the tourism key, like all apartments, viewpoints and attractions (values: "apartment", "viewpoint", "attraction") in one file?

More information for tourism key:
http://wiki.openstreetmap.org/wiki/Key:tourism


texts:
I'm using QGIS 2.0.1. I want to compute areas (in km^2) of all municipalities in Mexico. I understand that I need to convert my map from WGS84 to UTM to get units in km, but I'm not sure how to account for the fact that Mexico spans six UTM zones (from 11N to 16N), or which CRS I should pick from the hundreds available.

Which CRS should I use for this task?


texts: Promotion radius in miles around venues using PostGIS

texts:
Bear with me as I'm new to PostGIS and GIS in general. I also don't know the domain well enough to create a good title for this post! Here is what I'm trying to accomplish:

I have Venues in my system. Each Venue has a promotion radius in miles. Given an arbitrary point (lat/lng), I want to find all the Venues which encompass the arbitrary point in their promotion radius.

Here's what venues tables look like:

venues
id | name | latitude | longitude | promo_radius

Is this possible with PostGIS? Are there any resources that would help me further understand the query?

Thanks!


texts:
I am trying to categorize a polygon layer with several buildings, by the number of floors.

I want to aggregate by classes like this way:

1 floor = 1 to 2 floors;
2 floors = 1 to 2 floors;
3 floors = 3 to 4 floors;
4 floors = 3 to 4 floors;
5 floors = 5 to 6 floors;
6 floors = 5 to 6 floors

But by just simply entering the classification/legend I want (by hand) I get this result, which as you can see it duplicates the results:

Do I have to use "Ruled-based" and if so... how?


texts: Where I can find real world data for triangulation?

texts:
I want points for practicing triangulation algorithms.

Just points with coordinates [X,Y] and elevation [Z] of real world terrain.

I can use Google Maps Elevation API, but there is no method how to get only significantly terrain points.


texts: How to turn polygon into base layer in Surfer 14?

texts:
How do I turn a polygon I have drawn into a base layer in Surfer 14?

In Surfer 12, I used to be able to select my Base Map --> Digitise --> Select Polygon --> Digitise the Polygon Extents --> Export as Shapefile --> Import as Base Layer. However, in Surfer 14 it seems that when you select Digitise, you cannot select the polygon as the digitise function then switches off.

Does anyone have any tips on how to do this?


texts: Arcpy Batch Spatial Join Error code 000732

texts:
I am relatively new to python so apologies if this is a basic mistake.  I would like to make a python script to batch run the Spatial Join Tool.  The Target shapefiles are in one folder and the join features are in another.  Both folders contain 3 shapefiles with identical names.  The output location is to a third folder.  Script below:

import arcpy
import os

#Set environments

# 1. Target features. Split_Sample Points - SSPFworkspace and List variable
SSPF_1_workspace = r"D:\GIS\Python\Trackplot_py\Test_rev0_1\Split_Sample_Points"
arcpy.env.workspace = SSPF_1_workspace
SSPF_1_List = arcpy.ListFeatureClasses()
print (SSPF_1_List)

# 2. Join Features. Split point folder is SPF_2_workspace
SPF_2_workspace = r"D:\GIS\Python\Trackplot_py\Test_rev0_1\Split_Points"
arcpy.env.workspace = SPF_2_workspace
SPF_2_List = arcpy.ListFeatureClasses()
print (SPF_2_List)

# 3. Spatial Join Output folder.  SJ_Points is SJF_3_workspace
SJF_3_workspace = r"D:\GIS\Python\Trackplot_py\Test_rev0_1\SJ_Points"
arcpy.env.workspace = SJF_3_workspace
SJF_3_List = arcpy.ListFeatureClasses()

for TarF in SSPF_1_List:
    for JoinF in SPF_2_List:
        if TarF == JoinF:
            outputFSJ = os.path.join(SJF_3_workspace, "SJ_" + TarF)

            # Spatial join tool code

            arcpy.SpatialJoin_analysis(target_features= TarF,
                                       join_features= JoinF,
                                       out_feature_class= outputFSJ,
                                       join_operation="JOIN_ONE_TO_ONE",
                                       join_type="KEEP_ALL",
                                       field_mapping='',
                                       match_option="CLOSEST",
                                       search_radius="1",
                                       distance_field_name="SJ_Dist_m")
print (SJF_3_List)
print ("Complete")

The result is:

C:\Python27\ArcGIS10.8\python.exe D:/GIS/python/Trackplots/Spatial_Join.py
[u'SSS_I13152_PRC.shp', u'SSS_M13031_01_pt1_PRC.shp', u'SSS_X13426_PRC.shp']
[u'SSS_I13152_PRC.shp', u'SSS_M13031_01_pt1_PRC.shp', u'SSS_X13426_PRC.shp']
Traceback (most recent call last):
  File "D:/GIS/python/Trackplots/Spatial_Join.py", line 41, in
    distance_field_name="SJ_Dist_m")
  File "C:\Program Files (x86)\ArcGIS\Desktop10.8\ArcPy\arcpy\analysis.py", line 533, in SpatialJoin
    raise e
arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Target Features: Dataset SSS_I13152_PRC.shp does not exist or is not supported
ERROR 000732: Join Features: Dataset SSS_I13152_PRC.shp does not exist or is not supported
Failed to execute (SpatialJoin).

Process finished with exit code 1

As you can see the print lists appear to work specifying the shapefile list for both target and join features.  Yet the error (000732) states that the shapefile does not exist or is supported. Any suggestions?


texts: Leaflet map with jQuery autocomplete: fitBounds (click input) on selection?

texts:
I'm trying to fire a click event on my search button and/or the input field as soon as a user selects something from the jQuery autocomplete <li>:

https://jsfiddle.net/grh2dLyw/1/

I'm not sure where to put this (#zoominmap is my button, #columnsearch the input field):

document.getElementById('zoominmap').click();

This is my code:

        function initAutocomplete() {
        function log( message ) {
          $( "<header>" ).text( message ).prependTo( "#log" );
          $( "#log" ).scrollTop( 0 );
        }

        $( "#columnsearch" ).autocomplete({
          source: function( request, response ) {
                SQL_CLIENT.request({
                params: {
                    q: "SELECT name FROM " + dataSet + " where name Ilike '%" + request.term + "%'"},}).then(function(data) {
               response($.map(data.data.rows, function (item) {
                  return {
                    label: item.name,
                    value: item.name
                  }
                })
              )
            })
          },
          minLength: 3,
          select: function( event, ui ) {
            log( "Selected: " + ui.item.value );
                }
            });
        };

        // start jQuery autocomplete
        initAutocomplete();

        // search with AJAX
        $('#zoominmap').click(function () {
            if(previousSelection){
                map.removeLayer( previousSelection );   // if you're using Leaflet and map variable is named 'map'
            }
            input = $("#columnsearch").val();
            $.ajax({
            type: "GET",
            url: `https://schwaebische.carto.com/api/v2/sql?format=geojson&q=SELECT the_geom FROM ${dataSet} where name Ilike '${input}'`,
            dataType: "json",
            success: function (response) {
                geojsonLayer = L.geoJson(response)
                map.fitBounds(geojsonLayer.getBounds(), {maxZoom: maximumZoom, pan: {animate: true, duration: 0.5, easeLinearity: 0.25}});
            L.geoJson(response, {
                style: function (feature) {
                    return {color: "#FFF", opacity: 1, weight: 3};
                },
                onEachFeature: function(feature, layer){
                    // remove this layer after another polygon is selected
                    previousSelection = layer;
                    }
                }).addTo(map);
        }});
    })


texts: "Clip Raster by Mask Layer" produces zero values (and not Nodata) outside mask boundaries

texts:
This is my first question here. I have a big raster of the world's nighttime lights and want to clip it down to the Adm0 of a country:

In order to do so, in QGIS 2.18 what I do is using the Clipper function and checking the option "Crop the extent of the target dataset o the extent of the cutline":

I expected to get a clipped raster with Nodata values for the pixels outside the boundaries of my mask vector. What I get instead is a zero value for those pixels:

Apparently, the algorithm sets a zero value for all the pixels outside my Adm0 boundaries.
Specifying "No data value" to zero is not an option, since some original pixels within the Adm0 boundaries actually have a value equal to zero, so I would miss some significant information.

I haven't had this problem with previous clipping operations on other raster layers.

I get the same problem in QGIS 3.2.


texts:
When I use quickOSM in QGIS to download OSM data, I check the "relations" option, but then?
My need is the following :

Download OSM data for a specific area
In the table for the "point" geometry, I add a field : "StreetName"
For all nodes having a non empty addr:housenumber attribute, I want to get the street name from the relation "associatedStreet > Name" and fill the field created in 2.

The thing is, I can't find the relation anywhere.
Example : https://www.openstreetmap.org/node/685883404 I would have "Rue des vinaigriers" in my field.
This is my query
<osm-script output="xml" timeout="25">
    <union>
        <query type="node">
            <bbox-query {{bbox}}/>
        </query>
        <query type="way">
            <bbox-query {{bbox}}/>
        </query>
        <query type="relation">
           <has-kv k="type" v="associatedStreet"/>
            <bbox-query {{bbox}}/>
        </query>
    </union>
    <union>
        <item/>
        <recurse type="down"/>
    </union>
    <print mode="body"/>
</osm-script>

I tried with recurse type = relation-node, relation-way, with both of them...


texts: Format and value of Esri Technical Certification?

texts:
Starting mid-January you will be able to register for a series of Technical Certifications from Esri and I was interested in doing this, but have not been able to find out too much information about it online.

Does anyone know the format of the test?  From the Skills Measured section I get the impression that you will be working on a computer but I could not find anything definitive.

My only other question is if anyone has even heard much about this and how much weight it is going to carry?  I got a notification about this once a couple months back but otherwise I have not heard anything.  Just curious if anyone in their jobs or GIS circles are being told/recommended that they need to get this or if people are just different for now.


texts: QGIS, faulty Zonal Statistics calculations in case of empty polygons

texts:
I run QGIS 2.10 on a 64-bit Debian system and am currently trying to do some basic raster vector statistics.
However, I found that when applying the Zonal Statistics plugin, I end up with completely wrong values (instead of NULL) in case a zone does not have any raster cells within its boundaries. Apart from this, the calculations appear to be accurate.
Surfing SE for this error doesn't provide me with much insight except for knowing that the Zonal Statistics plugin seems to have a bit of an error history.

Any ideas how to resolve this?

In the meantime, I'll try my best using GRASS and SAGA...

After giving the original data a thorough inspection (there was no initial reason to be sceptical about its quality since it's supposedly official data), invalid geometries might be the problem. I'll double-check by drawing all polygons anew and get back to this later.

Using clean vector data didn't change anything. Thus, here's the bug report: https://hub.qgis.org/issues/13081


texts: OpenLayers not displaying GeoJSON layer with features added after the fact

texts:
I have a vector layer with a GeoJSON source layer. OpenLayers isn't showing the layer when I try to add features using the addFeatures method:

myGeoJSONObject = { /* some GeoJSON */ };

editLayer = new ol.layer.Vector({
    source : new ol.source.GeoJSON({
        style: styleFunction,
        defaultProjection: 'EPSG:4326',
        projection: "EPSG:3857"
    })
});

editLayer.getSource().addFeatures(myGeoJSONObject);

map.addLayer(editLayer);


texts: Implement system:time_start for NDVI trend analysis in GEE

texts:
I am developing a yearly NDVI trend analysis using Landsat data in Earth Engine over the last decades. My process is to add a NDVI band to every scene for each summer (green) month, then to reduce the collection by median to one image per month and finally to reduce all monthly images by max to generate one image for one year. My plan is to create a time series to see the positive/negative trends of NDVI over the years, but somewhere along the way the system:time_start property gets lost while reducing the images. I want to generate a chart but I always end up with: Error generating chart: No features contain non-null values of "system:time_start". I know you can copyProperties(image, ['system:time_start']) but I cannot find out where and how.
To reduce the length of the code, I cut it to two months:
//cloud mask & UNIX time
var timeField = 'system:time_start';

function maskL8sr(image) {
  var cloudShadowBitMask = (1 << 3);
  var cloudBitMask = (1 << 5);
  var qa = image.select('pixel_qa');
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
              .and(qa.bitwiseAnd(cloudBitMask).eq(0));
  return image.updateMask(mask)
    .select('B[0-9]*')
    .copyProperties(image, ['system:time_start']);
}

var geometryexpl = ee.Geometry.Polygon([[155.23045236462508,69.80628971253647],
                                      [159.03172189587508,70.77679819617978],
                                      [155.23045236462508,70.77679819617978],
                                      [155.23045236462508,69.80628971253647]]);

var landsat = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR");

//compute NDVI of biweekly plots
var addVariables = function(image) {
  //compute time in fractional years since the epoch
  var date = ee.Date(image.get(timeField));
  var years = date.difference(ee.Date('1970-01-01'), 'year');
  //return image with added bands
  return image
    //add NDVI band
    .addBands(image.normalizedDifference(['B5', 'B4']).rename('NDVI'))
    //add time band
    .addBands(ee.Image(years).rename('t'))
    .float()
    //add a constant band
    .addBands(ee.Image.constant(1));
};

//filter landsat8 data by county boundary, date, apply cloud mask, add variables
var NDVImap_may_2014 = (landsat
  .filterBounds(geometryexpl)
  .filterDate('2014-05-01', '2014-05-31')
  .map(maskL8sr)
  .map(addVariables));

var NDVImap_june_2014 = (landsat
  .filterBounds(geometryexpl)
  .filterDate('2014-06-01', '2014-06-30')
  .map(maskL8sr)
  .map(addVariables));

//reduce collection to one image using median of the rasters
var NDVI_may_2014 = NDVImap_may_2014.reduce(ee.Reducer.median());
var may_2014_NDVIbound = NDVI_may_2014.reduceRegions({
  collection: geometryexpl,
  reducer: ee.Reducer.median(),
  scale: 30
});

var NDVI_june_2014 = NDVImap_june_2014.reduce(ee.Reducer.median());
var june_2014_NDVIbound = NDVI_june_2014.reduceRegions({
  collection: geometryexpl,
  reducer: ee.Reducer.median(),
  scale: 30
});

var collection_2014_median = ee.ImageCollection.fromImages(
  [ee.Image(NDVI_may_2014), ee.Image(NDVI_june_2014)]);

var max_2014 = collection_2014_median.reduce(ee.Reducer.max());
var max_2014_bound = max_2014.reduceRegions({
  collection: geometryexpl,
  reducer: ee.Reducer.max(),
  scale: 30
});

var palettes = require('users/gena/packages:palettes');
print(palettes, 'palettes');
var palette_ndvi = palettes.colorbrewer.BuGn[5];

var ndviParams = {min: -1, max: 1, palette: palette_ndvi};

Map.addLayer(max_2014.select('NDVI_median_max').clip(geometryexpl), ndviParams, 'NDVI');

var chart = ui.Chart.image.series(
  collection_2014_median.select('NDVI_median'),
  geometryexpl,
  ee.Reducer.mean(), // default
  400, // nominal scale Landsat imagery
  'system:time_start')
    .setOptions({title: 'NDVI 1‐Year Time Series',
                 vAxis: {title: 'NDVI'}, });

print(chart);


texts: Storing Circles and run efficient query using PostGIS

texts:
In my database I have two tables where:

Table Pin contains a column of type point representing lat and long.
Table Fence contains (1) a column of type double representing a radius (2) and a foreign key reference to a Pin ID.

Given a lat, long and radius I would like to find any Pin+Fence that intersects (ST_Intersects) with the given circle.

Since I can't represent/save a circle in PostGIS as is and some answers on Stackoverflow refer to using ST_Buffer. Can I store the output of ST_Buffer or ST_Envelope in a new column in the Fence table and maybe use a trigger to update it (when Pin row get's updated)?


texts: Setting ArcGIS Server 10.2 Data Store Folder as Local vs. Network (UNC) path

texts:
Is there any difference?

How setting ArcGIS Server Data Store Folder as local vs. network path will affect the publishing process?

I saw warnings in ArcMap 10.1  when I Analyze the MXD before I publish, complaining that network path is used. (Assuming this might slow down the REST service if network access to MXD data is slow)

Example:
My file server has designated shared folder for ArcGIS data called "GisDataStore". It is mapped as local drive on my ArcGIS server, so I can access it both ways:

\MyFileServer\GisDataStore  and  Z:\GisDataStore

I wonder which I should use to Register the "GisDataStore" folder as ArcGIS 10.2 Server datastore.

Here's the warning I am getting:

SEVERITY    STATUS  CODE    DESCRIPTION

Medium      Unresolved  10027   Layer's data source is referenced via a UNC path


texts: Open Data Cube how to continue after the creation of the database?

texts:
I have some questions about using Open Data Cube. First of all I have successfully installed the virtual environment with miniconda and the database with PostgreSQLPortable-10.1.1.
Now I have already used the ingest command to initialize the database with .yaml files (always following the online guide) such as bom_rainfall_grids.yaml etc ..
Now I don't really know how to continue a UI exists? How do I get images processed?


texts:
I'm using: QGIS version 1.8.0-Lisboa on Windows XP, installed using OSGeo4W

I am behind a proxy but I have configured it in Options > Network.

When I go to Plugins > Fetch Python Plugins..

All of the the Repositories are "Status: Connected" however on the Plugins tab I can only see the 6 installed plugins (Repository: only locally available) - See images below:

Any ideas why I cannot see any plugins from the "connected" repositories?

Many thanks.

Update - As suggested by vinayan I can install plugins by manually downloading them and putting them in the relevant folder and using the QGIS Plugin Manager.  However still no joy using the Fetch Python Plugins.


texts: Show transparency after adding transparent MapBox layer to OpenLayers

texts:
I'm new to tile set maps. I am trying to add MapBox layers that I made in TileMill (and have sitting in my MapBox Hosting account) into OpenLayers. I have successfully added a non-transparent layer. I have also successfully added a transparent layer, however, that layer does not keep its transparency in OpenLayers. It seems to put a white background behind it, thereby hiding anything underneath it.

I've seen this example which I can get to work fine, but does not address transparency.

I've tried adding every transparency option I can think of in OpenLayers, but no luck. I'm wondering if maybe that XYZ format doesn't support transparency, because the OpenLayers Getting Started tutorial Ex. 5 shows a working example using the WMS format (which has a transparent: "true", parameter.

Anyway, here is the code I am using:

    <title>OpenLayers MapBox Example</title>


texts: QGIS Web Client printing error with some text added to template

texts:
I have a problem using printing. If I put a text label object (I want to add some text to the template) then pdf would never generate (it is error generating pdf shown). It only generates without any text added. My template works well with added elements like scale, image, legend. I don't know what to do to add also text.

Do you have any idea or any examples of working print layouts with some text added?


texts: Are there any examples or tutorials showing Branch Version editing with ArcGIS Server JavaScript API?

texts:
I am trying to implement branch versioning in an application that I'm building using the ArcGIS Server JavaScript API.

I found a blog post showing how to edit branch-versioned data in ArcGIS Pro, and the Version Management Service documentation shows the required syntax for each individual REST operation, but I was hoping to find a working demo or tutorial/blog post illustrating how to use branch versioning from start to finish.

For example, how to start the editing session and specify a new or existing version, then check in the edits and reconcile/post them, and deal with any conflicts.

Is there any example code or tutorial which illustrates how to use branch versioning while editing in a JavaScript application?


texts:
I want to extract a feature info from a WMS automatically (behind the WMS there is a GML Database). Is there a way to that easily in QGIS? As a coordinate that requests the feature info I want to use the centroid of a shape. With the cursor I can select and copy the information...
I read that queries like that are not possible since the WMS is exactly made for keeping users from getting all the information?


texts: Splitting single shapefile into many to use as geoprocessing extents in ModelBuilder?

texts:
I am trying to do the following things using ArcGIS model builder.

select by expression to output feature class.
use output feature class as spatial extent
raster to point

The main goal is to reduce loading of computer, because the original raster is super large. If I convert all raster to point, the point shape file clashes.

The problem is: how can I iterate the select function? In my model, it just exports one file. I have 15 shapes in a file, and I want to export each shape automatically.
How can I revise my model?


texts: Why project is not imported by QWC Admin GUI? (qwc-services)

texts:
I am testing QGIS's qwc-services docker container stack deployment.
I have cloned the 'qwc-docker' GitHub repo and customized the docker-compose.yml to create only the services needed to implement a basic WMS server.
'Docker-compose up' goes successful and the deployment works OK with the available qwc-demo QGIS project.
Next is trying to import a simple .qgs project I created with QGIS desktop 3.16 with 1 single layer showing data in a .csv file.
The project works OK in QGIS desktop and I have configured both the map and the layer for the QGIS-Server publishing according to the docs (I think). But my project is not shown in the qwc Config GUI.
I can see that the project is NOT imported into the QGIS server /data path and, looking into few log messages, I discovered that if the layer does not have set a custom property convert_categorized_layer to 'true' (in the .qgs xml file) it will be silently skipped and, being the only layer in the project, the whole project is silently skipped.
How can I configure a layer to have such custom property convert_categorized_layer set to true using QGIS-3.16 desktop? And, in general, is there a proven procedure available helping to publish QGIS 3.16 desktop projects into qwc-services stack deployment?


texts: Setting value of field in all selected rows of table using ArcGIS Desktop?

texts:
I have an attribute table in ArcGIS 10 (it's large - around 300,000 elements). I've selected some rows of the table using 'Select by attributes' and I want to set the value of a field for all of those rows.

How should I go about doing this? I'm hoping there will be a nice, easy, GUI-based way - as I'm teaching someone with very little GIS experience.


texts: FME 2019: Generate a new PK based on a concatenation of a string and a three digit increment

texts:
I'm transforming a number of geographic tables from one MCD to another, so I'm using FME (I'm new to FME).
I can map 80% of my fields with out any problem (MCD_A.fld_1 = MCD_B.fld_1) but in the in the second MCD I need to generate a new PK based on a concatenation of a string and a three digit increment something like 'BBBSYN_'||Increm.
I've tried generating the value directly in the FME writer (Feature Type -> User Attributes -> Value -> text Editor (see below) using @Count() - but this only generates simple increments '1','2','3' - I need '001', '002' etc. Also I already use @Count in another field and FME seems to confuse them both, i.e. in fld_1 and fld_2 of the same feature, count results in fld_1 = 1 and fld_2 = 2, not ideal.

I've looked through the safe website and the youtube channel - which is full of useful stuff but it's like finding a needle in a haystack - I also found this on SE : How to create autoincrementing primary key with dynamic writer in FME
Which looked like it would solve at least part of the problem - but it's a couple of years old and I don't have same options available in my feature type dialog box (I can't set type to identity or index to primary key for example)  - I'm using FME 2019 and this post is over two years old so i'm guessing they're not using FME 2019...
This is probably incredibly easy (it really should be...), but I can't find the right function ?


texts: Using pyshp to read a file-like object from a zipped archive

texts:
Presume I have an archive, foo.zip. It contains the parts of a shapefile: foo.shp, foo.dbf, foo.shx, foo.prj, and foo.cpg. I want to work with this shapefile in python, using pyshp, AKA the shapefile library.

Right now, I can extract the archive to disk, then read in the shapefile:

r = shapefile.Reader('foo')

I have quite a few such archives, though, and it would be nice to be able to read in these files without having to first extract each archive. According to pyshp's documentation, you can read in a file-like object with syntax like this:

myshp = open('foo.shp', 'rb')
mydbf = open('foo.dbf', 'rb')
r = shapefile.Reader(shp=myshp, dbf=mydbf)

And I quote: "This feature is very powerful and allows you to load shapefiles from a url, from a zip file, serialized object, or in some cases a database."

Thus I have tried creating a file-like object from my archive. First, creating the archive:

archive = zipfile.ZipFile('foo.zip')

Works fine. But if I then try

r = shapefile.Reader(archive.open('foo.shp')

I get the error, "shapefile.ShapefileException: Shapefile Reader requires a shapefile or file-like object."

That's legit; I should be passing it a file-like object. So I tried modifying the documentation's canonical example:

myshp = archive.open('foo.shp', 'r')
mydbf = archive.open('foo.dbf', 'r')
r = shapefile.Reader(shp=myshp, dbf=mydbf)

But this produces the error, "io.UnsupportedOperation: seek".

For what it is worth, I cannot use 'rb' as an option with archive.open(). If I do, I am told, "RuntimeError: open() requires mode 'r', 'U', or 'rU'". Also, using "archive.read()" instead of "archive.open()" just produces the "...requires a shapefile or file-like object" error.

I am stumped. As far as I can tell, I have created archive correctly, I have correctly used archive.open() to grab files as file-like objects, and I have used the correct syntax to have pyshp read from a file-like object. And yet, these errors. Does anyone have any idea what I'm doing wrong here? Has anyone had success passing file-like objects from zip archives to pyshp?


texts:
I have a request to provide GPX files of certain trails (around 100) that are present in OpenStreetMap. Given a local PostGIS install with OSM data, I think the process would be something like:

Define query to select the correct features, with appropriate metadata
Extract those features as shapefiles
Convert the shapefiles to GPX

Is this on the right track? Would a different format (KML, GeoJSON...) be better?

If this is the right way:

Step 1: Straightforward enough, I think.

Step 2: Using ogr2ogr with the postgis option? (eg http://www.bostongis.com/PrinterFriendly.aspx?content_name=ogr_cheatsheet)

Step 3: The only command line tool I could find, shp2gpx, doesn't seem to be available for download. Is there something else that could do the job?

Just wanted to check whether this is a sensible approach, and that there isn't an easier way, before I dive in.


texts:
Is it safe to use ray casting algorithm to check if a coordinate is inside a polygon specified as coordinates? Purpose of my question is because the longitudes are not parallel in a Cartesian sense, so is that a problem? If so, how big of a problem is it (for how big polygons it can matter for example), and how can I correct it? For example, can I map lat/lon to another domain and use ray casting safely?

I don't know much about GIS, but willing to learn anything related to my question.


texts: Tiling high resolution image using GDAL tools

texts:
I have a .tif image that I georeferenced using QGIS utility. It was originally in UTM with NAD83 as datum.
I need to make tiles out of it and get a kml file to overlay onto google earth.
Here's how I am doing it:

gdal2tiles.py -p geodetic -k filename

My image is slightly off where it should have been on google earth.What is going wrong here? Also do we need to warp the image using gdalwarp as my image though georeferenced is still a projection, it's not for a round globe. If yes, what command/parameters I need?


texts: Converting Mars Global DTM to Spherical Mesh

texts:
I am working on a video game that is attempting to provide a semi-realistic model of Mars. My starting point has been to find and use reliable geographic information to create the "model" of Mars.

I have found this global 465m Digital Terrain Model of mars. It is provided in a CUB file which I believe just represents an elevation as a height value in a large 2D image.

I was hoping someone here knows of either a formula or methodology I could follow in order to convert from a raster "heightmap" into a spherical mesh. I know how to mesh a point cloud so even if I could get help in going from a raster heightmap to a point cloud I can take it from there.


texts: Putting a label beside QgsVertexMarker object

texts:
Recently I learned how to create a point in the map independently from any layer, using QgsVertexMarker.
Now, I would like to put a label beside it, using Python. However, I do not see an option 'label' or something like that belonging to the class QgsVertexMarker.

How to do this?


texts: How to know when a widget from CARTO Airship is loaded

texts:
I'm getting some information from a server using fetch and using jQuery append to dinamically add the widgets to my page, but I need to know when the widget is done loading, there's any trigger that shows me that the widget is already loaded? here's my piece of code

  function drawLayers(sqlSource, style, ctrl, id_layer, widgetsObj, visibility){
    fetch(sqlSource)
    .then(response => response.json())
    .then(function(data){
    source[ctrl] = new carto.source.GeoJSON(data);
    vizLayer[ctrl] = new carto.Layer('layer', source[ctrl], style);
    vizLayer[ctrl].id = id_layer;

    layerBridge[ctrl] = new AsBridge.VLBridge({
      carto: carto,
      map: map,
      layer: vizLayer[ctrl],
      source: source[ctrl]
    });
    drawWidgets(widgetsObj, vizLayer[ctrl], layerBridge[ctrl]);
    vizLayer[ctrl].addTo(map, 'watername_ocean');
    if(!visibility){
      vizLayer[ctrl].hide();
    }

    layerBridge[ctrl].build();

 });
}

function drawWidgets(widgetsObj, vizLay, bridgeCtrl){
  for(var i = 0; i<widgetsObj.length;i++){
    if(vizLay.id == widgetsObj[i].source.id){
      if(widgetsObj[i].type == 'category'){
        $('.as-sidebar--right').append('
<as-category-widget visible-categories="5" show-clear heading='+widgetsObj[i].title+'  id='+widgetsObj[i].options.column+'></as-category-widget>
');
    bridgeCtrl.category('#'+widgetsObj[i].options.column, widgetsObj[i].options.column, {
      readOnly: false
    });
    vizmod = new carto.Viz('@widgets: viewportFeatures($'+widgetsObj[i].options.column+')');
    vizLay.blendToViz(vizmod);

    }

I need to use jquery to mod some properties from the widget after is loaded, I've tried to use the code below, just for an example, but it doesn't work

$('#'+widgetsObj[i].options.column).on('loaded', example_function);


texts: Add Values to New Column in Attributes Table using Python Console in QGIS

texts:
I've added a new column to my table called "STRNAME" using:

path = "/Users/shapefiles/roads/"
metro_roads = QgsMapLayerRegistry.instance().mapLayersByName('Metro_roads')[0]
processing.runalg('qgis:addfieldtoattributestable', metro_roads, 'STRNAME', 2, 20, 0, path+"Metro_roads_new.shp")

Now I want to populate the new column with values.  But I've looked for a processing algorithm (similar to "qgis.addfieldtoattributestable') and can't find one.

How can I add values to a table column from the Python Console?


texts: Exporting shapefiles from QGIS to PostGIS?

texts:
The tutorials I've been seeing such as this one https://saarapakarinen.wordpress.com/2015/06/13/tutorial-creating-a-postgis-database-for-qgis-part-2/

They mention a SPIT plugin that I can't find anywhere even after installing the WIEN version. DB Manager doesn't quite seem to do the same thing as described in the link above

Is there some other plugin to do the job?


texts: I created a dynamic map service in Arcmap and published to my server. The spatial reference used: GCS_WGS_1984. I used the javascript without basemap to display in a html. I cannot do a centerandzoom(). I suspect it's the spatial reference because if I add a basemap from esri the centerandzoom() works. How Can I solve this?


texts: Is it possible to print a map in PyQgis3 of a database table without rendering the layer in the project gui?

texts: There is a way to print a map using pyqgis3 without rendering the vector layer in the project interface? (by giving just the connection to the table, the map extent and the style )


texts: Fill holes with an area bigger of a specific value in a raster coverage

texts:
I want to fill some areas of a raster coverage with a value of my choice. Those areas now contain NODATA values. So, they look like holes in my coverage.

I know I could use r.clump function of GRASS, but I don't know how to specify the maximum area. Because if a hole is bigger than that area, I don't want to consider it as hole. Is that even possible with GRASS?


texts: Parameter values to create a vector grid with 1km x 1km cells?

texts:
I created a shapefile from WKT in QGIS using the values from this link.

I want to create a vector grid with 1kmx1km cells but I can't figure out what values need to be set for the X and Y parameters. How do I calculate those?


texts: Which monitor for working on printable maps?

texts:
I have a 24" monitor with a resolution of 1920x1080. Mostly I'm working on printable maps in QGIS. The size of some labels is sometimes smaller then 5mm on the final printed map. At my monitor this labels are not readable (as you can (not) see on the screenshot below).

Screenshot from canvas in scale 1:75000:

Is there a (software-)way to enlarge the view on my monitor or do I have to buy a monitor with a better resolution? Which one would you recommend?


texts: Statistical comparison between two maps with aggregated average observations and average estimates from a model on a 1 x 1 km² grid

texts:
I've been searching for methods for the past two days, but I couldn't find something usefull, so I hope someone could help me out:
I've got two datasets in python containing PM2.5 values:

Dataset which consists of hourly estimated values done by a model for a geographical location on a 1 x 1 km² spatial resolution
Dataset that contains low cost sensor observations for the same geographical location as point values

I want to analyze the difference between te observed values and the estimated values from the model on a daily scale (daily averages per 1 x 1 km² grid).
Thus far I created a vector grid of 1 x 1 km², aggregated all daily observations within the grids and assigned the daily averages to the grids (full code found here: https://github.com/ibatoukhtine/PM2.5/blob/main/PM25%20analysis.ipynb). I've done this for both the estimated values and the sensor observations.
Now I want to do an analysis between the two, but I'm stuck in what is possible for this comparison. Is it just as easy as creating a map which shows the mean difference per grid cell? For example, Grid_estimates - Grid_observations = Grid_difference? This will tell me if the model is overestimating or underestimating at certain locations, however I also still need to say something about the uncertainty of the grid cells depending on the amount of unique sensor contributed values to the daily average for the grid cells (maybe a CI?)


texts: How to Select Attributes from a Selection of Attributes

texts:
I have a land-cover map for an area in which I want to identify certain attributes that fall within a category. I want to select the shapes that represent wetlands that are in a grassland area. To do this I have tried any different methods such as separating the attributes into their own layers and subtract from the other, difference or minus, and tried to put one on top of another using the mosaic tool but neither worked. I am now thinking it may be best to write a specific query that would allow me to select an attribute within a selection of attributes. If I selected grassland, what can I do next to select wetlands are are within that selection of grassland? There are not options on the drop down menu for selection type that would allow me to do so. What can I do?

I am using ArcGIS Pro 2.3.2 and some screenshots are included with my question.


texts: How to register to event in OpenLayers 3?

texts:
I have a web application in OpenLayers and I am starting to rewrite it to OpenLayers 3.
In the application I have some features that I need to delete/move/add/edit  then to save the changes to the server.

In the previous version I
added a new control, and an event to the control.
These are all the controls in the previous version:

    controls = {
    point: new OpenLayers.Control.DrawFeature(editable,
                        OpenLayers.Handler.Point, { **featureAdded: onCompletAdd** }),
    line: new OpenLayers.Control.DrawFeature(editable,
                        OpenLayers.Handler.Path, { **featureAdded: onCompletAdd** }),
    polygon: new OpenLayers.Control.DrawFeature(editable,
                        OpenLayers.Handler.Polygon, { **featureAdded: onCompletAdd** }),
    regular: new OpenLayers.Control.DrawFeature(editable,
                        OpenLayers.Handler.RegularPolygon,
                        { handlerOptions: { sides: 5} }, { featureAdded: onCompletAdd }),
    select: new OpenLayers.Control.SelectFeature(
                editable,
                {
                    clickout: true, toggle: false,
                    multiple: false, hover: false,
                    box: false,
                    toggleKey: "ctrlKey", // ctrl key removes from selection
                    multipleKey: "shiftKey" // shift key adds to selection

                }
            ),
    modify: new OpenLayers.Control.ModifyFeature(editable, { **onModificationEnd: onCompleteMove** }),
    deletef: new OpenLayers.Control.SelectFeature(editable, {
        clickout: false,
        box: false,
        toggle: false,
        title: "Delete"
    })
};

It's calling the events featureAdded,onModificationEnd etc...
and for deleting I register to the event 'featurehighlighted' like this:

    controls.deletef.events.register("featurehighlighted", this, function (e) {
    editable.removeFeatures([e.feature]);
});

So how can i do it in OpenLayers 3?
I didn't find any helpful examples in the API.


texts:
I've download the Leaflet package and there is "leaflet-src.map"
What is it exactly? Do I need it to use Leaflet correctly?


texts: IFeature.Store() method throws COM exception

texts: I am developing a custom ArcFM attribute autoUpdater to update a couple of field values within a feature class. The autoUpdater will fire while creating as well as updating any feature with the feature class. Both ArcFM Objects and ArcObjects are used for the development. During feature creation, while trying to execute IFeature.Store() method using ArcObjects, the following exception occurs : "Error HRESULT E_FAIL has been returned from a call to a COM component". The feature class I am updating is participating with a Geometric network and also have relationship class associated with that.


texts: Finding coordinates of given raster value?

texts:
I would like to get the coordinates of a raster cell with a given value. I suppose many cells can take the same value, in which case several coordinates should be returned. In my first problem I'm interested in the location of the maximum value(s), but also seeking other values in planned future analysis.

I will work with very small rasters, which I cut based on polygons before processing, so speed is not so important in this case. Although fast solutions are welcome, also for wider applicability for others.

I do batch processing using python and QGIS, and can of course use other tools such as SAGA, GRASS, numpy, scipy to do this. The problem seems simple, but I have been completely stuck at this for a while now. I'm quite new to using python and QGIS.


texts: Avoiding errors when adding rows to blank table using ArcPy?

texts:
I am having problems when using ArcPy with ArcGIS Desktop to add in a set number of rows into a .dbf table.

What should I do to prevent these errors?

Here is a sample code that gives me errors:

##########setup###########
import arcpy, os
from arcpy import env

##########parameters##########

ws =  r'C:\Users\ddsouth\Documents\Python_tests'
env.workspace = ws

newFile = "adding_table_rows_test_1"
newFileFull = os.path.join(ws, newFile) + ".dbf"

###########code###############

#create the table
arcpy.CreateTable_management(ws, newFile)
arcpy.AddField_management(newFileFull, "RouteName", "TEXT", "", "", 15)
arcpy.AddField_management(newFileFull, "DepotName", "TEXT", "", "", 15)

#edit the table
cursor = arcpy.da.InsertCursor(newFile, ["ROUTENAME"])

for q in range (0,10):
    cursor.insertRow((str(1),))
del cursor

print "Completed."

The creation of the table goes fine the first time, then when it comes to the edit portion I get the following error message:

Runtime error
Traceback (most recent call last):
  File "<string>", line 26, in <module>
SystemError: error return without exception set

Then if I run the program again, I get a different error, this time during the creation of the table, even if I have the "Overwrite the outputs of geoprocessing operations" option selected:

Runtime error  Traceback (most recent call last):   File "<string>", line 18, in <module>   File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\management.py", line 15715, in CreateTable     raise e ExecuteError: ERROR 001143: Background server threw an exception.

I also noticed that while the tables I created appear in ArcCatalog, they do not appear in Windows Explorer, there is only a LOCK file while ArcMap is open.

What should be done?


texts:
I want to divide linestring into separate segments via Postgis, like it "Explode lines" tool in QGIS does, e.g. from linestring
I want to get separate segments

I have some thoughts about take it by using (st_dumppoints(Linestring)).geom then split linestring by their nodes with ST_split and dump again received geometry collection, but I got an unexpected result, because it returns me segments like

(not pair of coordinates as I expected!)

Maybe anybody have some ideas?


texts:
I made a map in QGIS, and I now want to render it into a web interactive product. I've tried to use Leaflet but I'm not able to reproduce the same pathways.

So far, I managed to produce simple lines as per code below:

        map.createPane('pane_Onthogenetic_Juveniles_Migration_13');
        map.getPane('pane_Onthogenetic_Juveniles_Migration_13').style.zIndex = 413;
        map.getPane('pane_Onthogenetic_Juveniles_Migration_13').style['mix-blend-mode'] = 'normal';
        var layer_Onthogenetic_Juveniles_Migration_13 = new L.geoJson(json_Onthogenetic_Juveniles_Migration_13, {
            attribution: '',
            interactive: false,
            dataVar: 'json_Onthogenetic_Juveniles_Migration_13',
            layerName: 'layer_Onthogenetic_Juveniles_Migration_13',
            pane: 'pane_Onthogenetic_Juveniles_Migration_13',
            onEachFeature: pop_Onthogenetic_Juveniles_Migration_13,
            // style: style_Onthogenetic_Juveniles_Migration_13_0,
        });
        bounds_group.addLayer(layer_Onthogenetic_Juveniles_Migration_13);

I tried:

        var polyline = L.polyline(layer_Onthogenetic_Juveniles_Migration_13, {}).addTo(map);
        var arrowHead = L.polylineDecorator(polyline, {
            patterns: [
                {offset: '100%', repeat: 0, symbol: L.Symbol.arrowHead({pixelSize: 60, polygon: false, pathOptions: {stroke: true}})}
            ]
        }).addTo(map);

And this doesn't work.

Example of what I would like to get in figure below (bot the red arrows and the arrow-shaped pathways).

Is there a way to do it?


texts: SelectByAttribute is not working (ArcGIS, Arcpy. Python script) -

texts:
My script is below.
I use to file from filegeodatabase to select all of rows (in a loop, one by one) to create a simple map and export it to jpeg.

For some reason, once I updated my script to run outside of arcmap (so I am not using reference to layer in arcmap now), the exported jpegs are focused on just a random spot rather than the feature location.

Any ideas what I need to update for it to run again?

import arcpy
arcpy.env.workspace = r"H:\NetworkAnalyst\TopologyVer4\mb_outputs.gdb"
arcpy.env.overwriteOutput = True
mxd = arcpy.mapping.MapDocument(r"H:\NetworkAnalyst\RoadSegments.mxd")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
arcpy.MakeFeatureLayer_management(r'H:\NetworkAnalyst\TopologyVer4\mb_outputs.gdb\T_Err_BridTunn',"lyrBridTun")

#Function returns unique values for specified field (all works fine)
def unique_values(table, field):
       with arcpy.da.SearchCursor(table, [field]) as cursor:
            return sorted({row[0] for row in cursor})

lstBT = list()
lstBT = unique_values(r'lyrBridTun', 'OriginObjectID')

for each_f in lstBT:
    whereClause = 'OriginObjectID = %s' %each_f
    arcpy.SelectLayerByAttribute_management(r'lyrBridTun', "NEW_SELECTION", whereClause)
    df.zoomToSelectedFeatures()
    df.scale = 1000
    arcpy.RefreshActiveView()
    str_name = str(each_f)
    arcpy.mapping.ExportToJPEG(mxd, r"H:\NetworkAnalyst\BridgeTunn_" + str_name, data_frame = "PAGE_LAYOUT", resolution = 300)
    print each_f, "Bridge Tunnel done"


texts:
I'm trying to figure out how to provide WMS - Services with QGIS and in my Research I encountered QGIS Server.

But I'm still uncertain about how it actually works. I install Apache Webserver and make the configurations like in the tutorial https://docs.qgis.org/testing/en/docs/training_manual/qgis_server/install.html

Than I open my QGIS Project-> Project Settings -> QGIS Server and in the following window I make all the necessary Settings.

What I should get is a WMS URL like "https://localhost/cgi-bin?SERVICE=WMS&REQUEST=GetCapabilities"   right?

But my client can't Import this WMS Service yet or can he? I think the moment I turn off my computer he won't have access to the URL anymore right? My Question is, after installing QGIS Server and creating a get capabilities request, what else do I need to do so that the URL is permanently available for users with permission.


texts: How to create rectangular buffers around points in QGIS with Python?

texts:
Using the following code from a similar question (How to create square buffers around points in QGIS with Python?) I am able to create square buffers (5x5 degree cells) around my points.

However, I also need to make 5x10 degree cells and 10x20, is there a way I can modify the following code to get this result?

Code:

layer = iface.activeLayer()

feats = [ feat for feat in layer.getFeatures() ]

epsg = layer.crs().postgisSrid()

uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer&field=x:real&field=y:real&field=point_id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'square_buffer',
                           'memory')

prov = mem_layer.dataProvider()

for i, feat in enumerate(feats):
    point = feat.geometry().asPoint()
    new_feat = QgsFeature()
    new_feat.setAttributes([i, point[0], point[1], feat.id()])
    tmp_feat = feat.geometry().buffer(1000, -1).boundingBox().asWktPolygon()
    new_feat.setGeometry(QgsGeometry.fromWkt(tmp_feat))
    prov.addFeatures([new_feat])

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)


texts:
I'm trying to build a GUI that allows replacement / update of the PostGIS database that holds my OSM extract (which I use for rendering, typical mapnik / tile server arrangement).

I'd like to show the user when the data that is already loaded was "current" to. That is, if my PBF extract (from http://download.geofabrik.de/) "contains all OSM data up to 2016-07-06T19:29:02Z", then I'd like to show the "2016-07-06T19:29:02Z" part in the GUI.

However I don't know how to get that from the database.

From http://wiki.openstreetmap.org/wiki/PBF_Format I see there is a osmosis_replication_timestamp field in the header, but I don't see how to retrieve that with just osm2pgsql.

At this stage, I'm thinking about either to "shell out" to something like osmosis, or trying to parse the HeaderBlock myself. I'd then have to store that (perhaps in the database, perhaps in my local state, perhaps as a file like the mod_tile approach). Both options suffer from lack of integration with the import process (e.g. if someone does a command line import using osm2pgsql, my separate state won't reflect that import), so I end up having to track import date / timestamp pairs. Both osmosis and a parser are pretty heavy for one timestamp...

Is there a way to get the replication timestamp directly from the results of an osm2pgsql import?


texts: Exporting region of interest (ROI) from scatter plot in ENVI

texts: I have a scatter plot of LST (y-axis) amd NDVI (x-axis) on ENVI 5.3. I selected some points from the distribution as a region of interest. When I export the ROI, I found only the value of Y. However, I need to get both the value of X and Y for deriving a linear equation. So, can any one assist me the steps to do.