texts: qgis2web, popups with information are displaying for the bottom layer (not the top layer)

texts:
I have a map, the very bottom (behind) layer contains united states (polygons for each state); and the very top layer contains my study locations (points layer).

The layers are displaying correctly, but when I hover over one of my study location points, the information is given as "Montana", which is the information for my "behind" united states layer and not for my "top" study locations layer. (Please see photo 1)

I changed the order of my layers (put the field locations behind the US polygons). I can't see my field locations, because they are hidden behind the polygons, BUT now when I hover over my study site locations (even though I can't see them) the information for my study site appears.

It seems that even though the map is displaying layers in the correct order, the information for popup-on-hover is accessed in reverse order. This is true for qgis2web preview within the QGIS desktop and once I extract the files and look at my Index.html in a web browser.

I added a duplicate of my study sites layer, behind the united states layer (but this is an ugly fix and the "Study Sites copy" layer shows up in the layer switcher, which will confuse the user) (Please see photo 2)

Can someone help me understand the connection between the order of layers which I think is determined in the layers.js:

var layersList = [lyr_USA1,lyr_MajorUSCities2,lyr_MT3,lyr_MTRoads4,lyr_MTCities5,lyr_StudySites6];

--> The USA layer already has 1 attached to it?

and the qgis2web.js code

map.forEachFeatureAtPixel(pixel, function(feature, layer)

How can I order layers in a different way (one layersList to display the layers and reverse layersList to extract hover popup information)?


texts: Versioning error when creating a Topology on a File GDB

texts:
Im having an issue with the 99999 Error. Im creating a Topology dataset for each of the featureclasses within a GDB. For some reason my script works for some feature classes and errors on others. I was having an issue with the "add rules" tool, and switched background processing off, which seemed to have fixed that error. However, I now get this unknown error...

  Cannot register as versioned on this database. Failed to execute
  (AddRuleToTopology).

Has anyone had this issue before, before I go and log this with ESRI?

My script:

ws = arcpy.env.workspace = workspace
fc = arcpy.ListFeatureClasses()
log_write("starting to create topologies for {}".format(workspace))

for f in fc:
    fc_name = arcpy.Describe(f).name[4:]
    log_write("creating topology for {}".format(fc_name))

    sr = arcpy.Describe(f).spatialreference
    log_write("spatial ref : {}".format(sr.name))

    input_dataset = arcpy.CreateFeatureDataset_management(ws, "{}_dataset".format(fc_name), sr)
    log_write("dataset created")

    new_fc = arcpy.FeatureClassToFeatureClass_conversion(f, input_dataset, fc_name)
    log_write("added fc to feature dataset")

    topo = arcpy.CreateTopology_management(input_dataset, "Topology", 0.001)
    log_write("topology created")

    arcpy.AddFeatureClassToTopology_management(topo, new_fc, 1)
    log_write("fc added to toplogy")

    arcpy.AddRuleToTopology_management(topo, "Must Not Have Gaps (Area)", new_fc)
    arcpy.AddRuleToTopology_management(topo, "Must Not Overlap (Area)", new_fc)
    log_write("rules added to topology")

    arcpy.ValidateTopology_management(topo)
    log_write("topology validated")


texts: Is the GeoWebCache S3 Blob Store only for use with preseeded tiles?

texts:
I'm thinking about moving over to the S3 blobstore, but it looks like it only works with preseeded layers which isn't an option for me.

Is it possible to use the S3 blobstore and have the cache built on the fly?


texts: How to disable polygon Label on the vertices (OpenLayers) in Edit Mode

texts:
I added the style on the polygon layer like so.

var style = $.extend(true, {}, OpenLayers.Feature.Vector.style['default']); // get a copy of the default style
style.label = "1233456 ha"
style.fillOpacity = 0.1
style.strokeWidth = 3
var styleMap = new OpenLayers.StyleMap({"default": style});
polygonLayer = new OpenLayers.Layer.Vector("Polygon Layer", {styleMap: styleMap});

The polygon looks fine on normal mode.

However, on ModifyFeature, labels are showing up on every vertex.

So, I test this http://trac.osgeo.org/openlayers/ticket/2176.

var styleMap = new OpenLayers.StyleMap(new OpenLayers.Style({
        label: "${getLabel}"
        // your other symbolizer properties here
    }, {context: {
        getLabel: function(feature) {
            if(mycontrolIsNotInEditMode) {
                return feature.attributes.label;
            }
        }
    }}
));

Then the labels on the vertices are disabled, BUT the label between vertices are still showing up.

how can I disable the labels between the vertices?
Thank you..


texts: Getting Transaction Support enabled in WFS on GeoServer (Boundless OpenGeoSuite)?

texts:
I am trying to publish my feature service using the Boundless GeoServer. It appears that in order to have edit capabilities (I want to be able to add new points to an existing point layer tied to my PostGIS db and edit attributes in existing points if necessary), I have to use the WFS (Web Features Service) in the dashboard. I followed the instructions and set it up with transaction level capabilities connected to my workspace.

However, when I go to preview the layers in GeoExplorer, I get the "transaction support not enabled" error when I use edit to create or edit. I then tried checking in Demos (as directed by this link: http://suite.opengeo.org/opengeo-docs/sysadmin/wfst/index.html) to test and see if my transaction support was enabled. It came back that for my user and password it was not enabled.

Is this somehow tied to a role or password?


texts: Viewshed analysis: how to assess the amount of a building that would be visible?

texts:
I am performing a Viewshed analysis in ArcGIS (10.1), and I use an OFFSETA (observer) of 1.7m and an OFFSETB (target) of 6m. The latter is the height of a typical building for which I wish to estimate the visibility across the landscape of a given study area. I would like to know:
if a location (actually, a pixel) in the resulting binary Viewshed Raster is flagged as  visible, I assume that a 6m-tall building placed at that location would be visible. But I also understand that we would not know if it is the whole  building that would be entirely visible, or just (say) its top 1 meters.

In the context of a Viewshed analysis, is there a way to assess the amount of the building that would be actually visible?


texts: How to compute average slope and area-weighted average slope?

texts:
I am learning GIS for use in hydrologic analysis.  Related to computation of overland flow there are two types of slope that I would like to be able to obtain from GIS.  The first is average slope for a given polygon area.  The second for the same polygon area is an AREA-Weighted average slope.  Presumably one would need to divide the polygon into subareas based on a slope-range (say 0-2%, 2%-4%, 4% -10%, etc) then compute a weighted average.

Is there a plugin that will provide this information, . . . or am I at the point where I need to learn how to write a python script?


texts: Select Data won't work... Parameters are not valid in ArcGIS ModelBuilder?

texts:
I'm creating a model in the ArcGIS 10.3 ModelBuilder and have problems to get the select data tool to work.
At first the complete model: IHS-Klassifizierung
With the Select Data tool (Input Data Element: Raster, Child Data Element: %layerband%_1) I'm selecting Band_1 or Band_3 from a Rasterfile, for this I need to check with Calculate Value if the Bands are called "Layer_x" or "Band_x" (In my raster data are different names and I mustn't change them), with the following code:

def x(Rasterfile):
import arcpy            # I think this is not necessary, but I'm not sure
desc = arcpy.Describe(Rasterfile)       # Get the necessary information
bandlayer = desc.children[0].name       # select the first band
ergebnis = str(bandlayer.split("_",1)[0])   # Remove the _ and the number from the Name
return ergebnis

as Expression I use x("%Rasterfile%") and Date type String.

When I run the model on my sample data (just one raster file) I get the following error:

Executing: IHS-Klassifizierung W:\arctest W:\arctest\output
Start Time: Tue Nov 03 16:39:06 2015
Executing (Iterate Rasters): IterateRasters W:\arctest # # NOT_RECURSIVE
Start Time: Tue Nov 03 16:39:06 2015
Succeeded at Tue Nov 03 16:39:07 2015 (Elapsed Time: 0,52 seconds)
Executing (Parse Path (2)): ParsePath W:\arctest\orthoclip.tif FILE
Start Time: Tue Nov 03 16:39:07 2015
Succeeded at Tue Nov 03 16:39:07 2015 (Elapsed Time: 0,00 seconds)
Executing (Calculate Value): CalculateValue x("orthoclip.tif")
"def x(Rasterfile):\n
import arcpy
desc = arcpy.Describe(Rasterfile)
bandlayer = desc.children[0].name
ergebnis = str(bandlayer.split("_",1)[0])
return ergebnis\n"
String
Start Time: Tue Nov 03 16:39:07 2015
Value = Band
Succeeded at Tue Nov 03 16:39:07 2015 (Elapsed Time: 0,04 seconds)
Executing (Select Data): SelectData W:\arctest\orthoclip.tif Band_1
Start Time: Tue Nov 03 16:39:07 2015
**Failed to execute. Parameters are not valid.**
**ERROR 000800: The value is not a member of Band_1 | Band_2 | Band_3 | Band_4.**
Failed to execute (Select Data).
Failed to execute (IHS-Klassifizierung).

If I try it without the Calculate Value Tool and put in the Select Data- Child Data Element directly Band_1 or Band_3 it works.

Sometimes it happens (I can't figure out why and I don't connect it to changes witch I made) that the following error appears:

Executing: IHS-Klassifizierung W:\arctest W:\arctest\output
Start Time: Wed Nov 04 13:23:47 2015
Executing (Iterate Rasters): IterateRasters W:\arctest # # NOT_RECURSIVE
Start Time: Wed Nov 04 13:23:47 2015
Succeeded at Wed Nov 04 13:23:47 2015 (Elapsed Time: 0,35 seconds)
Executing (Parse Path (2)): ParsePath W:\arctest\orthoclip.tif FILE
Start Time: Wed Nov 04 13:23:47 2015
Succeeded at Wed Nov 04 13:23:47 2015 (Elapsed Time: 0,00 seconds)
Executing (Calculate Value): CalculateValue x("orthoclip.tif") "def x(Rasterfile):\n
import arcpy\n
desc = arcpy.Describe(Rasterfile)       \n
bandlayer = desc.children[0].name       \n
ergebnis = bandlayer.split("_",1)[0]        \n
return ergebnis\n"
String
Start Time: Wed Nov 04 13:23:47 2015
ERROR 000539: Error running expression: x("orthoclip.tif")
Traceback (most recent call last):
  File "<expression>", line 1, in <module>
  File "<string>", line 3, in x
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\__init__.py", line 1246, in Describe
    return gp.describe(value)
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\geoprocessing\_base.py", line 374, in describe
    self._gp.Describe(*gp_fixargs(args, True)))
IOError: "orthoclip.tif" does not exist

Failed to execute (Calculate Value).
Failed at Wed Nov 04 13:23:48 2015 (Elapsed Time: 0,28 seconds)
Failed to execute (IHS-Klassifizierung).
Failed at Wed Nov 04 13:23:48 2015 (Elapsed Time: 0,93 seconds)


texts: How to update datasource of handle bad layer with Python console in QGIS 3.4

texts:
How do I update the datasource of file project in QGIS, for example myexample.qgz from Python console?

What is the command needed, or scripts standalone?

For example:

layer = iface.activeLayer()
newDatasource = "/home/usuario/EJES2020.shp"
newProvider = "ogr"
?.setDataSource('EJES2018',newProvider, newDatasource, extent)


texts: Warp a shapefile from transverse Mercator to Lambert projection

texts:
I'm quite a newbie in QGIS.
I imported in QGIS a shapefile of a European country (as well as a TIFF representing the contour). These files are referenced in EPGS 3035.

I would like to

transform a copy to a transverse Mercator projection (e.g. EPGS 3043)
transform another copy into a Lambert projection (e.g. EPGS 3034)
then
overlay the 2 results and see how much they differ.

I tried to use (with the raster layer) Raster > projections > reprojection but I couldn't obtain any result: the error message says that

  ERROR 1: The transformation is already "north up" or a transformation
  between pixel/line and georeferenced coordinates cannot be computed

My final target is to obtain a graphical confrontation between Mercator and Lambert projections.


texts: Autocad 2018 not compatible with ArcGIS desktop?

texts:
Does anyone have any experience with Autocad 2018? The latest release doesn't seem to be compatible with ArcMap Desktop 10.5. I've tried running the CAD to GDB tool and adding the files in the table of contents - both methods are coming up blank, when they would work fine with the previous Autocad 2017.
I'm not very familiar with Autocad, but the drawing comes up with no issues in the cad environment.
The only work around I can think of would be to save with a previous version of cad, but that may come with some loss of functionality.


texts: Automatically create new layers from a single dataset using definiton query for unique values in a field

texts:
I would like to automatically display a featureclass as different layers based on a specific field. I've seen a few tools online and on XToolsPro under the category of Split by Attributes, but they create multiple NEW shapefiles or feature classes.

What I would like is for a way to make multiple "layers" that show different views based on definition queries, but I would like all the definition queries to be created automatically based on the change in a specific field.

So for example, for a field "State", there would be multiple copies of the same layer created, but the definition query in each layer says State="TX", State="AK", etc.


texts: Trimble Data Transfer won't convert files

texts:
I was importing files from Trimble Recon for over 3 years and this never happened with me.

When I conect my controller and select the file to import, the file load 100% and when start to convert appear this error.
I tried in 5 computers and 4 controllers recon too.

Asked to local support but they can't find the solution.


texts: How to create data driven pages at a 1":100' scale covering entire network

texts: I am trying to create data driven pages at a 1":100' scale (1:1200 in Arcmap) that would fill several 11x17 sheets of paper. I am trying to use the Grid Index tool to create them, but it seems no matter what I do I always end up with just one giant rectangle despite specifying the scale, width, and height each grid feature should be. Not sure what else I can change in the settings to get this to work


texts:
Would anyone have any input as to where I could find a table/shapefile/feature class that describes the intersection type of roads? I've looked into TIGER roads, StreetMap premium, OpenStreetMap tables. Specifically looking for "street-light", "4 way stop", etc.

If such a database even exists.


texts: Where is the Merge selected features button in QGIS 2.2?

texts:
I can't seem to find the "Merge selected features" button. I need to merge several polygons inside a shapefile and can't find how to do it in this version. (I was running on Wroclaw until 2.2).

Is the Merge selected features button still included on QGIS 2.2? Where can I find it? What's the current icon for it?


texts:
I have a problem with understanding of the ST_Union (PostgreSQL) and the union function in QGIS. Basically, I have as an example two layers A and B that overlap. The result when using union in QGIS is one single layer with three rows: A without B, A & B, B without A. This is in principal the result I would like to have. Because the tables are quite big, I would like to select some columns and attributes; therefore I created a query in PostgreSQL. Following my approach:

SELECT row_number() over (order by a.id desc) as gid, a.nsg as nsg, b.ffh as ffh,
       ST_Union(a.geom,b.geom) as geom, st_area(St_union(a.geom,b.geom)) as area
FROM   p_gaiac.shape1 as a, p_gaiac.shape2 as b

The result is a single layer with only one row (attributes A & B). I am not able to get separate rows with separate attributes like in QGIS. What is wrong?


texts: How to merge two or more rasters which overlap?

texts: How to merge two or more rasters with values 0 and 1, where the value 1 does not overlap with other rasters. Finally, I want one raster file with pixel values 1 from all raster files as 1 and the rest 0.


texts:
I am trying to load a BigTIFF file into R, so that I can extract some data from a set of coordinates. The file I am working with can be found here: ftp://ftp.soilgrids.org/data/recent/CLYPPT_M_sl3_250m_ll.tif

I generally load rasters and extract points using the following code:

require(raster)
require(rgdal)
require(sp)

#some sample lat lon for extracting data
points <- structure(c(-89.426, -89.441, -72.18, 43.046, 43.038, 42.48), .Dim = c(3L, 2L), .Dimnames = list(NULL, c("longitude", "latitude")))

#load a raster layer
layer <- raster('/path/to/raster.tif')

#extract data from the layer for the points of interest.
output <- extract(layer, points)

When I run this code with the raster linked above I receive the following error:

/fs/data3/caverill/SoilGrids/CLYPPT_M_sl4_250m_ll.tif: This is a BigTIFF file.
  This format not supported by this version of libtiff..


texts: Using ArcGIS with dual/multiple monitors?

texts: ArcGIS doesn't seem to have much support for dual monitors. Maybe I'm missing something, but as far as I can tell you can't have the data view and layout view open at the same time on the two monitors, nor can you expand the map view to fill one monitor and have the toolbars and such on the other. How do you take advantage of dual monitors?


texts:
I am having trouble splitting a BigTIFF. I am orginally using @Ivan code from this question. I had to modify it a bit since I am using a .tif that was larger than 4GB. Pillow does not support BigTIFFs. So I ended up using tifffile. The code is running, though the image is not being clipped to the box parameter. I am thinking it has to do that tifffile reads the image in as a numpy array and is not actually being clipped by anything???
I've also noticed that when this code is running my memory is just about maxing out. I tried assigning the data type to a unassigned 8-bit and that drastically helped. Should I compress it as well? I don't want to change the data too much since I will be classifying them and do not want to lose/change data.
import os
from itertools import product
import tifffile

def tile(filename, dir_in, dir_out, d):
    name, ext = os.path.splitext(filename)
    img = tifffile.imread(os.path.join(dir_in, filename))
    print(type(img))
    w = img.shape[0]
    h = img.shape[1]

    grid = list(product(range(0, h - h % d, d), range(0, w - w % d, d)))
    for i, j in grid:
        box = (j, i, j + d, i + d)
        out = os.path.join(dir_out, f'{name}_{i}_{j}{ext}')
        img = img.clip(box)
        img = img.astype('uint8')
        tifffile.imsave(out, img)

tile('Orthomosaic_export_MonFeb01193821460106.tif',
     r'D:\ortho',
     r'D:\model_images',
     1000)


texts: How to use ST_VoronojDiagram for points that are almost aligned?

texts:
I am trying to use the spatialite ST_VoronojDiagram function on a table containing four 4 points that are almost aligned: (1.01 2.02, 2.02 3.03, 10.01 10.02, 11.01 11.02).
I wrote the following script in Python:

import sqlite3

#connexion to databse
conn = sqlite3.connect("D:\\my_path\\test.sqlite")
conn.enable_load_extension(True)
conn.execute("SELECT load_extension('mod_spatialite')")
cur = conn.cursor()
conn.text_factory = str

#creating points table
cur.execute("CREATE TABLE test_geom (id INTEGER NOT NULL PRIMARY KEY AUTOINCREMENT)")
cur.execute("SELECT AddGeometryColumn('test_geom', 'geometry', 4326, 'POINT', 'XY')")
conn.commit()

#inserting points
cur.execute("INSERT INTO test_geom VALUES (NULL, ST_GeomFromText('POINT(1.01 2.02)'))")
cur.execute("INSERT INTO test_geom VALUES (NULL, ST_GeomFromText('POINT(2.02 3.03)'))")
cur.execute("INSERT INTO test_geom VALUES (NULL, ST_GeomFromText ('POINT(10.01 10.02)'))")
cur.execute("INSERT INTO test_geom VALUES (NULL, ST_GeomFromText('POINT(11.01 11.02)'))")
conn.commit()

#creating voronoi polygons table
cur.execute("""DROP TABLE IF EXISTS 'test_vor'""")
cur.execute("""CREATE TABLE test_vor (id INTEGER NOT NULL PRIMARY KEY)""")
cur.execute("""SELECT AddGeometryColumn('test_vor', 'geometry', 4326, 'MULTIPOLYGON', 'XY')""")
conn.commit()

#inserting polygons
cur.execute("""
    INSERT INTO test_vor
    SELECT NULL, CastToMultiPolygon(ST_VoronojDiagram(ST_Collect(geometry),0,0.5))
    FROM test_geom
    """)
conn.commit()

#getting voronoi polygons geometry as text
cur.execute("""SELECT AsText(CastToMultiPolygon(ST_VoronojDiagram(ST_Collect(geometry))))FROM test_geom""")
conn.commit()
diag = cur.fetchall()
print diag

I get [(None,)] as a result.

When I add a random fifth point (18 -4):

cur.execute("INSERT INTO test_geom VALUES (NULL, ST_GeomFromText('POINT(18 -4)'))")

I get this (which is the form of result I want with the four first points):

[('MULTIPOLYGON(((8.334226 -4.294226, 10.725264 1.140878, 16.562488 4.467512, 21
.932716 6.966707, 21.932716 -8.226943, 6.940762 -8.226943, 8.334226 -4.294226)),
 ((10.725264 1.140878, 8.334226 -4.294226, -57.262816 61.302816, -57.262816 78.2
92816, -56.825362 78.292816, -53.3301 74.3601, 10.725264 1.140878)), ((8.334226
-4.294226, 6.940762 -8.226943, -57.262816 -8.226943, -57.262816 61.302816, 8.334
226 -4.294226)), ((16.562488 4.467512, 10.725264 1.140878, -53.3301 74.3601, 16.
562488 4.467512)), ((16.562488 4.467512, -53.3301 74.3601, -56.825362 78.292816,
 21.932716 78.292816, 21.932716 6.966707, 16.562488 4.467512)))',)]

I tried creating the Voronoi polygons in QGIS using the toolbox and I get this:

left: first 4 points; right: 5 points

blue: created with QGIS toolbox Voronoi polygons; purple: created with my python script

QGIS manages to create those polygons even for almost aligned points, so there must be a way to do it with a python script.

Is there any option/configuration/function that I'm missing?


texts: Truncating table values using ModelBuilder?

texts: I currently have this table with the "ACT_DATE" field displaying a long string of numbers and characters.  I only need the first 4 digits, for the year, and everything else I am attempting to discard.  Is there anyway of doing this in Model Builder, or am i stuck having to clean up the table first in Excel?  I have been attempting to work with Calculate Field and Calculate Value, but I just can't seem to get anything to work out my way.


texts: Does France have addressing products similar to British AddressBase?

texts:
In Britain people can license the Ordnance Survey's AddressBase product, which provides a list of almost all the country's addresses, with good building location  accuracy even for remote country buildings. This is obviously useful for address suggestion and geocoding applications!

Like all geodata products it's got some errors and omissions, but it's generally good (if expensive to license)

I'd like to broaden my horizon beyond just the UK. What's the best data available for France, and how good is it?


texts: I am using GRASS and connecting it to another application. The other application needs an FLT file of the raster map which I should export from GRASS. The other application only reads FLT files. Can GRASS export an FLT file


texts: Upload content to ArcGIS Portal from Python

texts:
Is there a way to upload new, custom data to an ArcGIS Portal instance from Python?

In essence, I want to programmatically execute the steps involved in adding new GeoJSON content via the web-interface, but do everything from within a python script.

The context is that I'm doing analysis in Python on one machine, generating results, and want to upload that data as a feature layer to the GIS system.

I cannot use arcpy.


texts: Measuring distances from each pixel center (in raster) to each point (in vector) in PyQGIS?

texts:
I have one raster layer and point vector layer. I want to measure distance between center of each pixel in raster to each point in vector layer. I extracted x,y center of each cell and added them to a list, but I can't make a relation between this list and vector points to calculate distances. How can I measure distances?

from qgis.core import *
from osgeo import gdal, ogr, osr

pntLayer = QgsVectorLayer("/Data/Villages.shp","pointLayer")

feat = [ feat for feat in pntLayer.getFeatures() ]

# Open tif file
ds = QgsRasterLayer("/Data/Study.tif","Study")

pixelWidth = ds.rasterUnitsPerPixelX()
pixelHeight = ds.rasterUnitsPerPixelY()

# extent of the layer
ext = ds.extent()

originX ,originY = (ext.xMinimum(),ext.yMinimum())
src_cols = ds.width()
src_rows = ds.height()

outDs = drive.Create("/Data/pop_Rst.tif", src_cols, src_rows, 1, gdal.GDT_Float32)
outBand = outDs.GetRasterBand(1)

outData = numpy.zeros((src_cols, src_rows), numpy.float32)

def pixel2coord(x, y):
    xp = (pixelWidth * x) + originX + (pixelWidth/2)
    yp = (pixelHeight * y) + originY + (pixelHeight /2)
    return(xp, yp)

pntRstList = []

for i in range(0, src_cols):
    for j in range(0, src_rows):
        rspnt = pixel2coord(i,j)
        pntRstList.append(rspnt)

print pntRstList


texts: I'm primarily looking for a GIS data source for postal codes in Japan as well as census data for different scales of geography.  Has anyone found a good source of these data sets that is in English?


texts: How to draw rectangle from the ReactBounds using OpenLayers

texts:
I am looking for an solution where I can draw the rectangle with the reactbonds I have found the solution one but its a legacy version of the OL. I have react bonds as
 ["-1535320.98225989", "6134670.12008471", "379695.245336351", "9059640.13281387"] and what I am looking for something similar to this:
http://dev.openlayers.org/examples/boxes.html
Here is the code for the following :
 var box_extents = [
                [-10, 50, 5, 60],
                [-75, 41, -71, 44],
                [-122.6, 37.6, -122.3, 37.9],
                [10, 10, 20, 20]
            ];
            var map;
            function init(){
                map = new OpenLayers.Map('map');

                var ol_wms = new OpenLayers.Layer.WMS( "OpenLayers WMS",
                    "http://vmap0.tiles.osgeo.org/wms/vmap0?", {layers: 'basic'} );

                var boxes  = new OpenLayers.Layer.Boxes( "Boxes" );

                for (var i = 0; i < box_extents.length; i++) {
                    ext = box_extents[i];
                    bounds = OpenLayers.Bounds.fromArray(ext);
                    box = new OpenLayers.Marker.Box(bounds);
                    box.events.register("click", box, function (e) {
                        this.setBorder("yellow");
                    });
                    boxes.addMarker(box);
                }

                map.addLayers([ol_wms, boxes]);
                map.addControl(new OpenLayers.Control.LayerSwitcher());
                map.zoomToMaxExtent();
            }

EDIT:
After some search I have found I can draw the rectangle from ol.interaction.Draw.createBox(). But still no luck with the reactBonds Here is the complete code:
rectangle = new ol.interaction.Draw({
                    source: new ol.source.Vector({wrapX: false}),
                    type: "box",
                    geometryFunction: ol.interaction.Draw.createBox(),
                });
                map.addInteraction(rectangle);

Any suggestions?


texts: Allowing users to add data to marker that they create on Leaflet map?

texts:
I have created a basic web map and implemented leaflet.draw functionality to it. What I would like to do is allow anyone that visits the web map to use the draw tools to add a marker (point) and when they add that marker a popup would appear allowing them to fill in per-determined fields that pertain to that point, this data would then be saved to an existing geojson file.

Is there a leaflet plugin / some other js wizardy out there that can make this happen?


texts: Way to export and import TIFF file symbols

texts:
I've read several reference to import saved 'style reference' for shapefile but couldn't find same works for TIFF files (or any image files).
I'm user of ArcGIS 10.3.1 and now I have TIFF files that have almost 20 category symbols (which L2_CODE in the picture below). And I want to save these symbols and apply to other TIFF files. I've already tried saving as .ri f(style reference), but this file doesn't show up when apply to new TIFF file

Is there anyone who knows how to use a style reference for TIFF file,
or is it impossible in ArcGIS?


texts:
I would like to know if my import of some osm2po results into a PostgreSQL database have been successful or not.

Is there something like: Import successful or a simple Done?

All I got was:

public.eu_2po_4pgr.geom_way SRID:4326 TYPE:LINESTRING DIMS:2
1(row)


texts:
I have a delimited text later and I need to only display every tenth label from a series of field values using an expression.
What’s the easiest way to do this?

Field values are a range of ascending numbers.


texts: Configuration of GeoServer and OpenLayers

texts:
I did the following:

I installed GeoServer by the installer for Windows as a service and standard folders .
Open http://localhost:8080/geoserver/www/ol-demo.html to test my GeoServer WMS layers and not appear the Google, Yahoo, Earth and WMS layers from my GeoServer service.
I copied into the www folder with examples of OpenLayers beginners, this folder contains the img folder, theme, metadata file, openlayers.js OpenLayers and my code examples from the book.

So I opened and ran http://localhost:8080/geoserver/www/code/chapter_5_ex3_google_layer_mashups.html Google Layer, the overlay of WMS example, not my local layer) .

Questions:

I used the Windows installer, it installed the GeoServer software is it alright? Thus no Tomcat or other container?
What can be wrong in the settings of my GeoServer server to not be serving the layers?

Installed according to the instructions from the official website, to run manually as a service, with personalized folder in the default folder and no attempt worked.

In this example http://localhost:8080/geoserver/www/ol-demo.html it not only calls the WMS layer (which is on my localhost), but also calls Google, Yahoo and Earth these should appear on the map, but don't.

My initial goal is to make the overlay of my WMS layers provided by PostGIS on Google Maps, where I did not succeed. So I decided to use the demo to see if my GeoServer server is working.

Any help is welcome, as I am with this problem for days and I do not know that attempts to do.


texts: Unable to perform join raster-table ArcGIS 10.3

texts:
I am trying to replicate some results from a publication. There, one of the datasets they use is WISE30sec map (available for download here: https://www.isric.org/explore/wise-databases). My problem is that even tough they state the following:

"The zip file includes a raster GIS file (WISE30sec) as well as a range of soil attribute data files (see Appendix 3 of the documentation). Using ArcGIS® or similar, users may join the raster data to the derived soil properties files, as managed the MSAccess® database (WISE30sec.mdb). Linkage is through the map unit code or grid cell identifier (NEWSUID) of the raster set and the NEWSUID of the various soil attribute data files."

I am unable to perform the above-mentioned join as I fail to find the grid cell identifier (NEWSUID) when using the Add Join function. ArcMap only shows the fields "ObjectID", "Value" and "Count" and none of them coincides with values shown at the tables I am trying to access. I have tried to copy raster and build attribute table but it didn't work. Also, when I put it in R, I cannot see these attributes that the raster is supposed to have.

EDIT: I am adding a some info, just in case they might be useful:

Raster properties:

  43201,
  1675310,0083333333,
  0,00833333331,35
  GBAIGGenericsigned
  integer16
  Bit-32768absentlevel:
  7, resampling: Nearest
  NeighborNoneBasicPermanent83,6249994416-180180,008331893-55,9833333333GCS_WGS_1984Degree (0,0174532925199433)D_WGS_1984skipped
  columns:1, rows:1, ignored value(s):
  0321560-10

Add Join function:


texts: Selecting one feature among intersections using ModelBuilder?

texts:
I have a model which check the location of objects and select those which have the same position on the map.

I need to have a selection of only one object in intersections but not to doi it manually. Do any tool in ArcToolbox or process in ModelBuilder can complete such a task?

For instance, here I have six point objects which have intersections. How to choose only three of them, which have to be moved?


texts: Extracting OpenStreetMap Bus Routes using ArcGIS Desktop?

texts:
I am wanting to extract bus routes from OpenStreetMap to use as shapefiles in ArcGIS.

Is there a simple way to go about this?

As i am trying to extract them directly from OSM however i keep getting the error message:

  This XML file does not appear to have any style information associated
  with it. The document tree is shown below.


texts: Choosing coordinate system for cost distance route selection analysis?

texts:
I'm performing a cost distance route selection analysis for a type of linear infrastructure in Australia.  The process I will use will be similar to this ESRI tutorial. The route will traverse several states and therefore several UTM zones.  Normally, when calculating long distances like this in Australia, I would use the Geoscience Australia Lambert projection.  Indeed, I began converting my data to this projection in order to begin calculating things like slope (which I also resampled to a lower resolution to match Land Use).

I'm now beginning to wonder if I am using the correct coordinate system for the route analysis.  When dealing with Cost Distances, I'm assuming this type of coordinate system is appropriate.  However, I'll be rasterising many polygons and I wonder if I should be using an Equal Area projection such as Australia Albers.

In short, what type of coordinate system is appropriate for performing a cost distance route analysis?


texts: Merging road segments into one road line in QGIS

texts:
I have a Shapefile with a road network (lines) which includes several road segments (lines).

Is there any way to merge all these lines into one line using QGIS?


texts: Restrict layers in attribute table widget in JavaScript?

texts:
I'm building a webapp, where you can change from several map services. In every map service there are several layers with same naming, but not with the same number of layers. I cannot set the attribute table widget up in Web AppBuilder, because the JSON too hardcoded. I'm really new in JavaScript.

I would make a list in js, what contains the allowed layer's  name:

layersToShow: [Layer1, Layer2, Layer3]

in the widget's _ResourceManager.js

I think here is the section from the manager, where I could repoint the reader. But I don't know how. I would add all the possible allowed layers to the list, and the js may fail to get the layerinfo. How would you accoplish this?

      if (updateConfig) {
        if (this.config.layerInfos.length === 0) {
          // if no config only display visible layers  ,egy tombot csinalok es amikor itt filterel, akkor egy consistal kerem le a tombbol
          var configLayerInfos = attrUtils.getConfigInfosFromLayerInfos(layerInfos);
          this.config.layerInfos = array.filter(configLayerInfos, function(layer) {
            return layer.show;
          });


texts: Generating equal sized polygons along line with PyQGIS

texts:
I would like to create polygons along a line to use them for AtlasCreator in a next step.
ArcMap has a tool called Strip Map Index Features.
With this tool I can choose the height and width of my polygons (say 8km x 4km) and produce/rotate them along the line automatically.
One of the generated attributes of each polygon is the rotation angle that I need to rotate my north arrows in Atlas Generator afterwards.

Does anyone have an idea how to solve this task in QGIS / with PyQGIS?
GRASS- or SAGA-algorithms or a prossessing-toolbox-model which could be used inside a custom plugin would be fine, too.
I need not only the print extents but also the polygons itself as I want to print a map with all polygons/extents as some sort of overview map. I am looking for a PyQGIS-solution that can be used in a QGIS-Plugin without the need for installing a software aside from QGIS (no RDBMS like PostGIS / Oracle).


texts: Converting PyQt5.QtCore.QVariant to string

texts:
I am using QGIS with the in-built Python console. I want to replace some values in a shapefile with some other values. But when I try to change the attribute values, I get an error about the data types (see the code below):

Error: AttributeError: 'QVariant' object has no attribute 'split'

However, I find no way to convert the QVariant to a String. If I run the code without the editing turned on, everything works fine, but the changes are not saved. How can I convert the QVariant to a string and are there ways around this issue?
shp_layer = QgsVectorLayer(shp_path, "", "ogr")

with edit(shp_layer):
    for feature in shp_layer.getFeatures():
        height_str_list = feature[height_list_field].split(":", 1)[1][:-1].split(",")  # some fancy string editing in here, but i swear it works
        if '...' in height_str_list:
            height_str_list.remove('...')
        height_flt_list = [float(i) for i in height_str_list]
        feature[height_field] = sum(height_flt_list)/len(height_flt_list)  # because QGIS pretends not to know mean()
        shp_layer.updateFeature(feature)
    shp_layer.setSubsetString("")

I also discovered if I say height_str_list = feature[height_list_field].value().split(....)
I get another error:

"AttributeError: 'str' object has no attribute 'value'"


texts: QgsRasterCalculator fails in standalone Python application

texts:
I am trying to write a simple Python program which shall be a wrapper for QgsRasterCalculator() function, to use it in some batch processes. In both QGIS 2.14 (on Linux) and QGIS 3.10 (on Windows 10) I have problems to find an expression that works.

Copying a raster, when the expression holds just the band reference works. But including an operator or a second band will fail.

E.g.:

layer = QgsRasterLayer(path_to_dted_tile, 'n45_dt1')
entries = []
band = QgsRasterCalculatorEntry()
band.ref = layer.name() + '@1'
band.raster = layer
band.bandNumber = 1
entries.append(band)

expression = '"n45_dt1@1"'

context = QgsCoordinateTransformContext()  # QGIS 3.10 only
calc = QgsRasterCalculator(expression,
                           path_to_output_file,
                           'GTiff',
                           layer.extent(),
                           layer.width(),
                           layer.height(),
                           entries,
                           context)   # QGIS 3.10 only
result = calc.processCalculation()

works as expected. When I use: expression = '"n45_dt1@1" + 1.0'

QgsRasterCalculator crashes with the following output on the console:

  Process finished with exit code -1073740940 (0xC0000374)

If the expression is wrong, how should it look like and, if not, what may cause this failure?


texts:
I have a shapefile British Columbia with postal codes. The attribute table gives the codes with no spaces, eg V6P2A8. I am trying to join a csv file by postal code (the key I am using). My csv file has postal codes with a consistent space, eg V6P 2A8. Naturally the join doesn't work because of that space.

How can I either remove the space?


texts: Can anyone list out great sources of GDAL training? Books, videos, companies and/or people offering tutoring?


texts: How to pull data from an existing PostGIS database into Carto?

texts: I have an existing PostGIS database. I want to be able to access my tables from that DB on my Carto enterprise account. Any suggestions?


texts:
I am having a problem with leaflet markers. The problem is that when I drag a marker with a lower z-index under another marker with a higher z-index and let it go, it makes another marker. I have a click event that adds a marker to the map when I click it.

This happens in Chrome and Internet Explorer 11.
This does not happend in Firefox.

Another problem I noticed was that when I would drag a marker outside the browser window and relase in Firefox, it would mess up things up.

JSFiddle
http://jsfiddle.net/eaflores/x9JzP/


texts: How to add a rasterised line feature to raster via raster calculator in ArcGIS

texts:
I want to sum different layers (Population Density, Land Use, Rivers, Roads) to derive an Exposure map for landslides using the raster calculator in ArcGIS.

First I converted the two line shapes (Rivers & Roads) into raster. Then I aggregated both to be able to reclassify.
Afterwards I reclassified(standardised) the raster files in order to be able to weight them.

Summing raster LU and PopDens works flawless, but as soon as I add one of the rastered line features, the result looks like it uses them as mask --> only cells, where river or road occurs are considered, instead of the whole study area.

Does anyone have an idea?


texts: Plugins crashing with "No module named qgsprojectionselectionwidget" in Windows?

texts:
My Python QGIS plugin works on Ubuntu (at least for me, and all fatal bugs reported so far are solved), but when I run it on QGIS installed from the official 2.10.1 64-bit installer for Windows (on Windows 7 64-bit), I can't install it due to following error:

This plugin is broken
No module named qgsprojectionselectionwidget

Another user reported it this problem as well. but with another error message (I don't know which operation system and architecture does he have, but I'd guess Windows and 64-bit EDIT: I guessed wrong, he has Mac OS X):

This plugin is broken
No module named qgsmaplayercombobox

Both refer to widgets in the GUI form, which I built through Qt Designer in Ubuntu. Qt Designer is packed with QGIS installer for Windows as well and claims to contain the same QGIS custom widgets, so the problem shouldn't be in those widgets missing, and QGIS should recognize them - but it doesn't.

So how to fix it? My plugin is still experimental, so I can just add "sorry, doesn't work under Windows" to the metadata, but it's not a real solution.

EDIT: the same behavior confirmed on 32-bit version too.

Now I found a good hint in comment to this answer: use os.path.join to make the paths cross-platform. This is easy for the few hardcoded paths I have, but the errors seems to suggest the problem is deeper, perhaps even in Qt. How can I handle this problem in my plugin without tweaking anything else?


texts: Getting z-Value from xy-Values within IPolygon using ArcObjects?

texts:
How can I get the z-Value of given xy-Values within an IPolygon?

I have an IPolygon, which is zAware. For given xy-Values I want to get the z-Value of the Polygon-Shape. What I want to do is something like intersecting the Polygon with an vertical line and getting the Intersecting Point.

Has anyone an idea how to do something like this?


texts: Displaying attributes for features within current map extent in ArcGIS Desktop?

texts:
I have found Landsat 8 PanSharpened WMS service here. It is a mosaic made by Esri.

I can open attribute table of this mosaic dataset.

Can I sort rows of attribute table to see only rows related to rasters in my current view extent using ArcGIS for Desktop?

The option I need is similar to Attribute Table widget's Filter by Map Extent:

  Only displays attributes for features within the current map extent


texts: joining tables with different number of entries

texts:
I work in ArcGIS 10.3 and I have this problem:

I have a dbf table of about 300 restaurants, each restaurant has a different values like address, number of seats etc.
I have another dbf table of beers that are tapped in the above mentioned restaurants. For each record of a restaurant in the first table, there are one or many records in the second table, for example:

Table 1

Restaurant "A", Elm Street 15
Restaurant "B", Ash Road 4

Table 2

Restaurant "A", Elm Street 15, Beer_type1, Beer_type1_price
Restaurant "A", Elm Street 15, Beer_type2, Beer_type2_price
Restaurant "A", Elm Street 15, Beer_type3, Beer_type3_price
Restaurant "B", Ash Road 4,    Beer_type2, Beer_type2_price
Restaurant "B", Ash Road 4,    Beer_type3, Beer_type3_price

I'm trying to join (or relate?) the 2nd table to the 1st one, so I could see  (hopefully in a single row) which beers are available in which restaurant. It would be something as "join one to many". I researched some but didn't find a suitable solution. I found very useful tool called "Add Relate" which is available only in ArcGIS Pro, so that's not an option.
When done, I will add points (X,Y coords) of the restaurant to the map and each point should have attributes of the 1st table AND of the 2nd. Any tips how to handle this?


texts: Select Iterate feature, use 2 Iterate rasters

texts:
For multiple polygon data, I want to save the values of multiple raster data that exist for each record for each record.
I have tried various things, and I have not been solved, so I ask a question.

I converted the polygon data into points and saved the raster values using'Extract value to points' and'Extract multi values to points'. The data is output horizontally, so I tried using'transpose fields', but the field I wanted was not created.
We tried to convert all raster data into pointer data, select one record of polygon data, and then save the points included in the area. In this case, the point data was made into one shp through'merge'. It is possible for a small area, but data was not generated when shp exceeded 2Gb.

The process I want to do is as follows.
I have raster data and and vector data.
1000 raster data are stored in the'img' folder, and the format is *.img.
The vector data is polygon and there are 1000 pieces.
The model builder was planned by writing the following scenario. Image reference

I select the first record (ex. eco1) from the grid data. And it creates the'aaa' field.
Read the raster data stored in the'img' folder one by one using'iterate rasetr'.
Read the raster data of 2. included in the record area selected in 1. and save it in the'aaa' field of 1.
After reading all raster files, eco1 saves them as shp files.
Next, repeat steps 2-4 after selecting eco2.

Questions
Various attempts are made, but failure occurs.
I want to get help on how to approach it.


texts: Finding Highest Percentage Overlap between two polygons?

texts:
I'm currently working with electoral and census information.
Quick summary of the two geographical units I'm working with:

Voting Precinct --> These polygons provide me with electoral results from a specific area.
US Census Tract --> These polygons provide me with information from the US Census Bureau.

And you guessed it! These two type of polygons do not match.

There is also no rule to determine if one type is bigger than the other, as there might be some Voting Precinct polygons larger than some census tracts polygons but also shorter.

Using R, Python or QGIS, is there a way I can get a table that tells me for each Precinct which is the US tract that overlaps the most?

I'm expecting to know the following information
Precinct 1 --> Matches the most with Tract 43
Precinct 2 --> Matches the most with Tract 45
Precinct 3 --> Matches the most with Tract 43 (yes, the same tract can match the most with two or more precincts.

I know it's a very broad question and any comment, solution, library, code will definitely help me!


texts: Feeding gdalwarp with binary data instead of ascii, how to specify header information

texts:
I have a working pipeline where I download a tar.gz file which contains multiple binary files in the german RADOLAN format (Translate PNG obtained through RADOLAN from DWD's projection to EPSG:3857).

I parse these binary files and save them as ASCII files. For each binary file I generate two ASCII files, one with the precipitation data to generate the colored clouds and one with the "valid" radar data, which I then use with an edge detect algorithm to generate the radar boundaries (grey line).

I generate the image via gdalwarp, because I need to reproject this RADOLAN projection into Web Mercator (EPSG:3857)

For this I use the following command:

gdalwarp
  -overwrite
  -s_srs "+proj=stere +lat_0=90.0 +lon_0=10.0 +lat_ts=60.0 +a=6370040 +b=6370040 +units=m"
  -t_srs "EPSG:3857"
  source-pixels.asc
  target-pixels.tif

gdalwarp
  -overwrite
  -s_srs "+proj=stere +lat_0=90.0 +lon_0=10.0 +lat_ts=60.0 +a=6370040 +b=6370040 +units=m"
  -t_srs "EPSG:3857"
  source-edges.asc
  target-edges.tif

I use that command twice for each binary file, since I have two ASCII files.

Finally I read those tif files wit PIL (Python Image Library) where I colorize the clouds, detect the edges, and merge the images.

My problem is with the ASCII files.

Writing these ASCII files is very slow. I use numpy's savetxt method to write the arrays into the file. This takes about 0.33 seconds per file, and I'm saving 48 ASCII files (from 24 bin files). I'm waiting 16 seconds just to get the files written.

In contrast, when I use numpy.save, each file takes a bit less than 0.03 seconds.

So I'm thinking about not using the ASCII files, but to feed binary files into gradwarp.

The issue is that these ASCII files contain additional information in a header.

This is the content of one of the edge files:

ncols         900
nrows         900
xllcorner     -523462
yllcorner     -4658645
cellsize      1000
NODATA_value  -1
668 668 668 668 668 668 668 668 668 668 668 668 668 668 ...

668 indicates that this pixel is out of range, that no radar was able to collect information about that position. Those pixels (positions) which do contain valid radar information, even if this means that nothing is there, are indicated by a 0.

The pixel data ASCII files have the same header, but the values vary according to the precipitation intensity.

The important thing here is the header information:

ncols         900
nrows         900
xllcorner     -523462
yllcorner     -4658645
cellsize      1000
NODATA_value  -1

These ASCII files apparently follow some gis-related specification which I don't know about. I found those values in a PDF file.

Now, if I want to feed binary data into gdalwarp, this header won't be in the file. How can I pass that header information to gdalwarp? Also, how do the pixels in the binary file need to be specified? How do I specify gdalwarp how it should interpret that binary data?

I'm a hobbyist at this, so don't expect any background knowledge from my part.

The most I came to think of was to create GeoTIFF files instead of ASCII files, as it apparently is possible to add some geoinformation to the TIF file, so that gdalwarp would be doing a TIF to TIF reprojection, which I guess sounds sensible.

How would I structure such a GeoTIFF file? Would that be what comes the closest to a binary format as I'm intending to use it? Are there any libraries for generating those files with NumPy? Or should I take another approach?


texts: How to make QGIS label only the shapes that are showing up?

texts:
I have a layer where shapes are defined by rules (country) and secondly by categories (a type). I then select (check box) the country to show the shapes only for that country.  This works very well, as I change them on each composer and display data relevant to that country.
The problem is with the labels: they show up for all countries all the time and I cannot hide them.

Does anyone know how to make QGIS label only the shapes that are showing up?


texts:
I have two layers with roads of different types
I wish to have pretty intersections, as here

Symbol levels work only for features in single layer.


texts: Preventing diagonal flow direction in ArcGIS Desktop

texts:
I would like to prevent diagonal flow direction/accumulation, as it leads to frequent "river shortcuts". As shown on the picture, instead of following the river in blue, because of the coarse raster resolution, a shortcut is created between the two diagonally adjacent cells (indicated by the red arrows). Here, the grey raster is the filled DEM.
Is there any option in either the flow direction or the flow accumulation function that satisfies this constraint?


texts: Remove unwanted black area around clipped raster in QGIS 2.18

texts:
Following on from Remove white background in Raster in QGIS

When I choose the unwanted value - in my case I want to remove black area around the clipped raster - I write 0 as the unwanted value, but I also lose a lot of important detail in the clipped raster. So there is a lot of differences between the original raster to the clipped raster

How can I choose to remove only the black area around the clipped raster?


texts: Unable to use (and even to list) pure attributes tables (with no geometry in them). I need to do a flow using one containing {from, to} cities codes

texts:
I have to do a flow of home-work trips.
The base of my work is a geographical layer of cities
an OpenStreetMap map of French cities that has their shapes and cities codes (named : INSEE) :

The second source, an attribute table with the list of the trips
it's a normal SGBD table named deplacement_domicile_travail_2017 that has no geometry.
This table tells from what and to what cities are the trips, and will be used to ask QGIS to draw the arrows between them.

But tables with no geometry are not listed anywhere under QGIS
However, I see no way under QGIS 3.14 to perform my flow layer.

In explorer, PostGIS folder doesn't list me the deplacement_domicile_travail_2017, and in fact, no other non-geographical table (even if they have a primary key).
If I edit my layer and go the Joins menu, only geometrical tables or view are offered for targets, and not pure attributes tables.
If I was to create a view on deplacement_domicile_travail_2017 in order to add geometry on it (and join it would have to link it to my OSM table), this view would have two geometries on it : one for the starting city, one for the destination city. It would cause other troubles and I want to avoid that.

Question :
How are used pure attributes tables, those containing only data columns and no geometry, on QGIS ?


texts: How can I get the average elevation for an area in QGIS?

texts: I have got .tif file for a country and a dbf file which divides the country area into 42 subunit. I want to get the average elevation data for each subunit in meters.


texts:
Following this example and according to Leaflet WMS documentation I extended L.WMS.Source class and overwrote hooks (e.g. showFeatureInfo).

My example can be found here.

When I use the code as given in the example above

name = info.split(/'/)[3];
info = "Name: " + name;

the popup displays

  Name: border:none

which seems to be a part of HTML inside my popup window.

What do I have to do to only get the name of the feature displayed?


texts: CORS issue loading an external local GeoJSON into OpenLlayers 5

texts: I have a local folder on my network with my HTML file and two subfolders, one with the downloaded folder from OpenLayers (renamed ol) and the other with my GeoJSON file (data). When I try to load the HTML file in Chrome, I get a CORS error. Why do I get this even though the document is loading everything from the same file? Is there any way to load it with a JS file like you can with Leaflet? Here is my code:


texts: Fill Sinks of multiple rasters in a Python expression

texts:
I have 31 rasters that represent different watersheds and I would like to use the Fill (Spatial Analysis) tool to fill the sinks of all rasters at once in Python. All the rasters have the same name followed by a number (0-30).

Here is my code so far:

Import arcpy
From arcpy import env
From arcpy.sa import *
Arcpy.env.workspace = (r“Workspace\Path”)

For I in range (31):

outFill = Fill(“raster_name”+str(i)+“.tif”)
outFill.save (“Workspace\Path”)

So far this code only fills sinks of the raster_name0 and does not put it in the folder location I specify. Any ideas on how to modify this code to work on all 31 rasters and save them in a specficed location?

I have ArcGIS 10.1 with an advanced arcinfo license.


texts: Interpolation of holes resulting from DTM groundfiltering of SfM pointcloud

texts: I am wondering if there is an efficient  interpolation method to close the holes in a dtm resulting from groundfiltering. At the moment I am directly interpolating the point cloud to a raster using idw or kriging but finally there are dents in the dtm where for example trees where located before. In general the result is sufficient but it would be nice to know if someone found another way in R or SAGA not creating a dent or small mount. I am thinking of something like closing the gaps from outside to the inside, in extending the edges.


texts:
I'm working with ArcGIS Desktop 10 and I have different polygon features in different layers with their attributes like shape_area. I would like to create a unique graphic where it is possible to see the respective percentage of each one.

Can I do this if I have different layers?


texts: Identify incorrect direction of travel of Polyline

texts:
I was trying to work on a simple way to identify incorrect direction of travel of polyline as compared to base map. However was unsuccessful. I tried few plugins but didnt get anything to help me with it.
As per the image below the polyline on right has incorrect direction of travel

Did someone face the same issue before ? Any recommendations on how to address this ?


texts: reading raster without making QGIS apply the scale factor to the raster

texts:
Since QGIS multiplies the image with the scale factor when reading, it messes up the image if the scale factor is not properly defined to the data as QGIS expected. For example, I am reading MODIS NDVI into QGIS and array is int16. Int16 range is between (-32768, 32767). Range of input data values are between -2000 and 10000 and the scale factor is 10000. QGIS multiplies the array with 10000 and values go beyond the acceptable range of int16 and QGIS rescales the pixel values and make them between (-32768, 32767). In the end, you can not fix the problem by dividing it by the scale factor because values are rescaled once they go beyond int16 range.
Is there any possibility make QGIS not apply the scale factor to the data?

download the data and extract NDVI sublayer
http://e4ftl01.cr.usgs.gov/MOLT/MOD13A2.005/2009.07.28/MOD13A2.A2009209.h10v06.005.2009228232821.hdf

My recommendation to QGIS guys. You shouldn't worry about scale factor because you don't know how it is defined (either multiplication or division).


texts: AddIndex_Management throws maximum open cursors exceeded

texts:
The following script is run to create indexes across all the feature classes within a particular dataset:

import arcpy
import os.path

def listFcsInGDB():
    for fds in arcpy.ListDatasets("FRED.GS*", "Feature") + ['']:
        for fc in arcpy.ListFeatureClasses('', '', fds):
            yield os.path.join(arcpy.env.workspace, fds, fc)

def listIndexNamesInFeature(f):
    return [i.name for i in arcpy.ListIndexes(f)]

def listFieldNamesInFeature(f):
    return [i.name for i in arcpy.ListFields(f)]

arcpy.env.workspace = "C:\Users\rmorgan\AppData\Roaming\ESRI\Desktop10.1\ArcCatalog\ds08.sde"
idx_stub = "IXGS_"
attributes = [
                "QASTATUS", "POINTUSAGE",
                #... several more here
            ]

fc = listFcsInGDB()
for f in fc:
    indexes = listIndexNamesInFeature(f)
    fieldList = listFieldNamesInFeature(f)
    for a in attributes:
        if a not in fieldList:
            print a + " is not in " + f
            continue
        if any(a[:4] in s for s in indexes):
            print a + " has an index already in " + f
            continue
        arcpy.AddIndex_management(f, a, idx_stub +  a)
        print a + " is now indexed in " + f

After some time of successfully creating around a hundred indexes, the following error is thrown:

Traceback (most recent call last):
  File "E:\PyScripts\AddIndexes.py", line 57, in <module>
    arcpy.AddIndex_management(f, a, idx_stub + str(counter) + a)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\management.py", line 5271, in AddIndex
    raise e
ExecuteError: ERROR 999999: Error executing function.
Underlying DBMS error [ORA-01000: maximum open cursors exceeded] [FRED.GS_FC1]
....
Underlying DBMS error [Error executing stored procedure sde.version_util.close_state::ORA-01000: maximum open cursors exceeded] [SDE.DEFAULT][STATE_ID = 2561192]
Wrong column type [ORA-01000: maximum open cursors exceeded]
Underlying DBMS error [ORA-01000: maximum open cursors exceeded]
Underlying DBMS error
Underlying DBMS error
Underlying DBMS error
Underlying DBMS error
Underlying DBMS error
Underlying DBMS error
Failed to execute (AddIndex).

This is ArcMap 10.1 running on top on Oracle 11g SDE.

Running the Oracle command show parameter open_curs; results in 2000 cursors.  Certainly, this script doesn't need 2000 cursors.  There are about 35 feature classes and 15 fields needing indexed.


texts: Using query in a while loop causing android application to freeze

texts:
I am having some troubles with the while loop + query. What I'm trying to do is track the user location (using onlocationChanged) and keep querying a FeatureLayer until the query returns a result and then stop the while loop is using a boolean.
The problem once the "while" loop is launched is that the application freezes, here's my code:

boolean noResults=true
while(noResults)

{
    pointGeometry=MyPosition;

//myPosition is changing it's value in "onLocationChanged" using the listener, see code below

Graphic pointGraphic = new Graphic(pointGeometry, simpleMarker);
SpatialReference spatialRef = MapFragment.mMapView.getSpatialReference();
Unit unit = spatialRef.getUnit();
double adjustedAccuracy = b_dist;
Polygon p = GeometryEngine.buffer(pointGeometry, spatialRef,adjustedAccuracy, unit);
Query q = new Query();
q.setGeometry(p);
q.setSpatialRelationship(SpatialRelationship.INTERSECTS);
fLayer.selectFeatures(q, SELECTION_METHOD.NEW, callback4);
}

public void onLocationChanged(Location loc) {

                double locy = loc.getLatitude();
            double locx = loc.getLongitude();
            Point wgspoint = new Point(locx, locy);
            mLocation = (Point) GeometryEngine.project(wgspoint,
                    SpatialReference.create(4326),
                    mMapView.getSpatialReference());

//Here's where myPosition changes its value.

            myPosition = mLocation;
}

CallbackListener<FeatureSet> callback4 = new CallbackListener<FeatureSet>() {

        public void onCallback(FeatureSet fSet) {
            Graphic[] grh=fSet.getGraphics();

            if(grh.length>0){
                        noResults=false;
                        showTaost("there's a result");
                    }
});

Is there an alternative to the while loop to fix this? If not, how can I use something like a delay (query the Feature Layer every 5 secondes or something like that)?


texts: Selecting Polygon Feature using ArcGIS API for JavaScript?

texts:
I used ArcGIS JavaScript API 3.25 for developing.  I want to select polygons within 1km radius. so I used a sample code that esri provides. But it works only for point features. I want to use that code for select polygons.

Can you help me to correct my coding?

var featureLayer = new FeatureLayer("http://203.189.68.219:6080/arcgis/rest/services/Sites/inginimitiya_Tank/MapServer/11"
/*
  ,{
    mode: FeatureLayer.MODE_ONDEMAND,
    infoTemplate: template,
    outFields: ["Name", "OBJECTID"]
  }
*/
);

var symbol = new SimpleLineSymbol(
  SimpleLineSymbol.STYLE_SOLID,
  12,
  new SimpleLineSymbol(
    SimpleLineSymbol.STYLE_SOLID,
    new Color([247, 34, 101, 0.9]),
    1
  ),
  new Color([207, 34, 171, 0.5])
);
featureLayer.setSelectionSymbol(symbol);

// Make unselected features invisible
//var nullSymbol = new SimpleMarkerSymbol().setSize(0);
//featureLayer.setRenderer(new SimpleRenderer(nullSymbol));

map.addLayer(featureLayer);
var circleSymb = new SimpleFillSymbol(
  SimpleFillSymbol.STYLE_NULL,
  new SimpleLineSymbol(
    SimpleLineSymbol.STYLE_SHORTDASHDOTDOT,
    new Color([105, 105, 105]),
    2
  ), new Color([255, 255, 0, 0.25])
);
var circle;

// When the map is clicked create a buffer around the click point of the specified distance
map.on("click", function (evt) {
  circle = new Circle({
    center: evt.mapPoint,
    geodesic: true,
    radius: 1,
    radiusUnit: "esriMiles"
  });
  map.graphics.clear();
  var graphic = new Graphic(circle, circleSymb);
  map.graphics.add(graphic);

  var query = new Query();
  query.geometry = circle.getExtent();
  // Use a fast bounding box query. It will only go to the server if bounding box is outside of the visible map.
  featureLayer.queryFeatures(query, selectInBuffer);
});

function selectInBuffer(response) {
  var feature;
  var features = response.features;
  var inBuffer = [];
  // Filter out features that are not actually in buffer, since we got all points in the buffer's bounding box
  for (var i = 0; i < features.length; i++) {
    feature = features[i];
    if (circle.contains(feature.geometry)) {
      inBuffer.push(feature.attributes[featureLayer.objectIdField]);
    }
  }
  var query = new Query();
  query.objectIds = inBuffer;
  // Use an objectIds selection query (should not need to go to the server)
  featureLayer.selectFeatures(query, FeatureLayer.SELECTION_NEW, function (results) {
    var totalPopulation = sumPopulation(results);
    var r = "";
    r = "The total Census Block population within the buffer is " + totalPopulation + ".";
    dom.byId("messages").innerHTML = r;
  });
}


texts: On a map of the globe that is 1920x1080px, how big would each pixel be?

texts: I heard from one source 21km square, and from another 245km square...


texts: How to create a color ramp to change the color of a svg image in arcgis javascript api

texts:
I have an SVG image to create a PictureMarkerSymbol for a feature. How to create a color ramp to change the SVG color based on a criteria, like a attribute on my feature?

The idea is like that:

var mySymbol = new esri.symbol.PictureMarkerSymbol("../resources/img/myimage.svg", 32, 32).setOffset(0, 8);
var renderer = new esri.renderers.SimpleRenderer(mySymbol);

renderer.changeColor ( field, colors, criteria )


texts: Interpreting Historical Aerial Imagery - Agricultural Midwest

texts:
I was recently looking at some aerial imagery of a Midwest farm and wondering how to interpret a couple things. First, in the 1938 imagery, there are several fields that appear to have weird geometrical lines in them... Could these be related to the installation of crude tile drainage systems? Just tillage? They seem too frequent to be haybails, as the field is on the order of hundreds of feet wide.

Also, in a ravine draining into the river, there appears to be a relatively abrupt elevation change as if there is a structure or human-made mound of some kind within it, but I am not sure how to interpret this? Any idea what these features could be?

Pictures of the 1938 photograph as well as modern Google Earth snapshots also attached.


texts:
I have an arcpy script that uses a XTools Pro geoprocessing tool, which works if I run it from ArcMap (with the XTools Pro extension enabled), but if I run it outside of ArcMap I get an error.

Tool I'm running:

arcpy.ImportToolbox("C:/Program Files (x86)/DataEast/XTools Pro/Toolbox/XTools Pro.tbx")
arcpy.XToolsGP_DispersePoints_xtp(InputLayer, OutputLayer, "0.25 Meters", "1 Meters", "BY_CLUSTERS", "RING", "", "")

Error message:

arcgisscripting.ExecuteError was unhandled by user code
Message: Failed to execute. Parameters are not valid.
ERROR 000824: The tool is not licensed.
Failed to execute (XToolsGP_DispersePoints).

I have attempted to enabled this extension from arcpy with no success, and have found the XTools Pro help pages don't mention arcpy at all (so now I wonder if it's even supported).

How can I enable the XTools Pro extension from within arcpy?  I've attempted to use the following, but am unsure what arcpy would call the XTools Pro extension, so may be using the incorrect Extension name.  (I have tried XTools, XToolsPro, XTools-Pro, XTools Pro, xtp, XToolsGP)

if arcpy.CheckExtension("XToolsPro") == "Available":
    arcpy.CheckOutExtension("XToolsPro")
else:
    print "XTools Pro Not Available"


texts: What is source of "TabletPC inking error code. RtpEnabled called multiple times" error from within ArcMap?

texts:
I am attempting to remotely troubleshoot an issue with an ArcMap AddIn (.NET 4) that I wrote where the users (only a couple of them) are receiving this error:

  TabletPC inking error code. RtpEnabled called multiple times
  (Exception from HRESULT:0x80040239).

But I believe this exception is be thrown within a Stop Operation and that my AddIn is picking up the error.

Unfortunately, it is at a place in the code where I don't log the exception w/stack (which I will fix in the next build). This is occurring on only 2 machines that are both running 10.2.1. One is Windows 7 Pro and the other is Windows 8.1 Pro.

The issue is intermittent and recently resolved (temporarily) by clearing the users local profile. But the issue will return.

Has anyone seen this before within ArcMap and are there any known causes (ArcGIS Desktop itself, popular AddIns, extensions, specific system configs, user profile configs)?

As I am readying a new build with more logging, is there any diagnostic information that I can be collecting that would help isolate the cause if this error?


texts:
I've been playing with Mapbox's blog post about converting heat maps to contour lines, and I'm stuck at the extraction part of this process.

I used the CSV plugin to create a vector layer of points, then a heat map raster layer based on that, but I can't seem to be able to see the contour lines based on that raster layer after that. When I looked at the properties and selected categorized, it didn't seem to have any symbols listed.

I'm guessing this is probably some kind of type error since I had to create the raster layer based off of map units and not meters, but I don't know how to correct it. What am I missing?

Here's a screenshot. It's a map of racially charged geographic locations based on their name. The main difference being that I've stripped any instances of Arroyo Negro, but left in more explicit cases.


texts: Quick question regarding geocoding limits on google's API

texts:
We have a google-api business license, so we can do 100k requests in a day, but at the moment, we are hitting our limits as we have a bunch of initial data to geocode. Instead of constantly hitting google, I want to set a option to turn the API off as soon as we get the QUERY_OVER_LIMIT status back. I also want to turn it back on as soon as it's available. I was going to do this by setting a timestamp in the DB the first time we get the status.

What I'm wondering is, at what time does google reset the keys so they can do more call? If it's at midnight, what timezone is it (CST/PDT/EST/GMT/etc)?


texts: Prevent Leaflet from requesting non-existing custom tiles

texts:
I have created custom tiles using TileMill. extracted tiles as png images and adding them on base layer. using following code

var venu_map = new L.tileLayer('static/b18/{z}/{x}/{y}.png', {
               maxZoom:22,
              attribution: '&copy; OpenStreetMap contributors'
              }).addTo(map);

It works fine. But for every zoomlevel change leaflet makes a number of(around 20 -30) requests to non existing tiles resulting in "NetworkError: 404 NOT FOUND. Complete error msg is "NetworkError: 404 NOT FOUND - http://127.0.0.1:8000/static/b18/16/47372/30375.png". The particular image does not exist. Even after setting maxbounds this error show up. Any solutions?


texts: Filling gaps in raster to zero value cells?

texts:
I am trying to add some forest elevations to a DEM using an ArcView (Basic) level license of ArcGIS Desktop.  I have converted a polygon feature class of a forest into a raster with the same resolution as my DEM.  Then I reclassed the forest raster to 15 (the height of the forest).  This has resulted in a patchy raster with gaps of nothing between.

To get a raster that covers the entire study area I created a constant raster with the extent of the forest polygon and containing the value of 0.  when I add the two rasters together I end up with a replica of my forest raster rather than a raster with forest cells of value 15 and 0 values everywhere else.

Any Ideas of how to fill in the gaps between my forest cells?


texts: Help excel table not transferring to Arc correctly!

texts: New to GIS and Ive been trying to figure this out for a while. I have an excel sheet with CT names such as 0128.03, 0128.04, starting from 0001.00 and up. Not all CTs have .01, .02 at the end just .00. anyways when I transfer it to a dbf and open it up in arc the field shows with no decimals. I did a custom format for 0000.00 but for some reason 0128.03 shows up twice and I was missing 0128.04. I realized the duplicate of 0128.03 actually was 0128.03999999 and should be 0128.04 theres others doing the same thing as well. I have no clue how to get 128.04 or stop the continous decimal. Do I format the cells a certain way in Excel? I need these values to be correct so i can join the CT field with another table.


texts: How many and which versions of QGIS should I use to test a plugin?

texts:
Greetings I'm working on a plugin. I'm developing it in QGis 3.4 and now testing it in QGis3.0 (with several bugs).
With these two versions is enough if I want it to be approved in the repository? I don't want to uninstall and install new versions, it's a risk.
What are the requirements? Minimum versions?


texts: How do edit a .shp file geometry manually?

texts:
In the program that I've designed, I am loading a .shp file with custom .dbf attributes. Once a week the national institute releases an updated database in the form of a .shp file. Because my .shp includes custom .dbf attribute fields, I would have to manually find all the newly added objects and import them to my .shp file and then fill out the attribute fields (open attributes table tool).

Is there a method to edit a .shp file geometry in a "text" editor, similar to how you can edit the attribute table in QGIS?

I've managed to convert the .shp and .dbf file to SQL geometry table, edit the fields, but I am unable to find a method to convert SQL geometry back to a .shp/.dbf file.


texts: ArcGIS Pro Network Analyst service area not calculating distance correctly

texts:
Using ArcGIS Pro version 2.7* updated from 2.6.0 (which did not solve this issue)
I built my own network based on street center-lines data. I'm not concerned with looking at vehicular speeds or direction along the network at this point so I'm using length as my 'cost'.
I input the facilities and set the distance to feet, established the breaks at 150, 300, 450, and 600 ft. Looking at the generated polygons critically I'm seeing some issues. The red line in this picture is measured out to 150 ft but the polygon is far larger than what I would expect.

These screenshots show my settings in ArcGIS Pro which I imagine illustrate the proper setup of this service area.

What could be the issue I am overlooking here?


texts: Obtain address from PostGIS OSM data like Nominatim does?

texts:
We currently use the Nominatim service to do address lookups. We want to store the OSM data in our own servers as our usage gets higher.

How do I calculate the address of a given coordinate, like how Nominatim does?


texts: How do I create a layer with popups that has hundreds of unique markers/symbols?

texts:
I am creating a layer in ArcGIS JavaScript API.

I am looking for the best way to have a feature layer with popups and unique markers.

This layer has hundreds of markers at a time.

The markers are unique.

I have attempted to use a uniqueValueRenderer, but discovered that the unique values are assumed to be predefined.

My markers are customized based on the values in the data, but there are too many variations to create predefined unique values.


texts: QGIS create point data attribute as the name of the polygon it falls in

texts:
I'm working in QGIS with a point data set overlaying a farm with 140 fields. While I can analyse the points from each field individually, I would rather use an awk script to calculate statistics I need.

The problem is that I am unable to import the containing polygon_id as an attribute into the point table.

Any suggestions would be greatly appreciated.


texts: Make permanent layer in ArcMap, Image Analysis Functions to display hillshade

texts:
I'm using ArcGIS 10.2.2.

I want to make relief map, so that it looks good at all scales. I used the hillshade function but after exporting the data, the relief looks good only at one scale.

I'm using this method:
Using Image Analysis Functions to Display Layer Tints on Hillshades

How can I save the result without loss of quality?

P.S. I already tried to use mosaic.


texts: Re-centering a map layer using manually input coordinates in GeoExt

texts:
I am trying to create some code in GeoExt that will allow me to centre that map layer I am viewing to coordinates that have been submitted through a formPanel.  I would have thought that this would have been fairly straightforward but this is proving to be quite troublesome and I am wondering if any of you guys could help?

So far I have created the the panel with 2 boxes to accept the eastings and northings but can't work out how to get the values that have been submitted to actually move the image.

The code I have cobbled together this far is:

formPanelCoord = new GeoExt.form.FormPanel({
    border: true,
    layout: 'anchor',
    bodyStyle:'padding: 20px',
        items: [{
        xtype: "textfield",
        width: 100,
        fieldLabel: "Easting",
        name: "search_lc",
        value: "",
        emptyText:'Easting',
        allowblank: false
    },

{
        xtype: "textfield",
        width: 100,
        fieldLabel: "Northing",
        name: "search_lc",
        value: "",
        emptyText:'Northing',
        allowblank: false
    }],

    listeners: {
       actioncomplete: function(form, action) {

        features = action.response.features;
        console.log(features.length)
     }

  }

});

formPanelCoord.addButton({
        text:"Go to Coord",
        handler:function() {
        console.log(formPanelCoord);
            //alert(formPanelCoord.items.items[0].getValue());
            //alert(formPanelCoord.items.items[1].getValue());
            formPanelCoord.items.items[0].getValue();
            formPanelCoord.items.items[1].getValue();
                       //function recentreMap (form) {
                       var easting = from.Easting.value;
                       var northing = form.Northing.value;
                       //var easting = formPanelCoord.items.items[0].value;
                       //var northing = formPanelCoord.items.items[1].value;
                       map.setCenter(new OpenLayers.LonLat(easting,northing));

formPanelCoord.addButton({
        text:"Go to Coord",
        handler:function() {
            this.search();
     },
     scope:formPanelCoord
});

            //this.search();
     },
     scope:formPanelCoord
});

Any help gratefully received!


texts: Calculating Topographic Exposure with ArcGIS Raster Calculator?

texts:
I need you help trying to calculate Topographic Exposure (TOPEX) on a DEM in ArcGIS.  I am trying to figure out wind exposure on a DEM with 25m cells.
There is a great blog post where someone explains how to do it but he uses GRASS. Jamie Popkin's GIS Blog I thought I would be able to put the same equation into ArcGIS Raster Calculator but it says the equation is invalid. I've tried playing with it as much as I can but I can't seem to make it work.

His first equation is:

r.mapcalc incl4north = "atan(((dem - dem[4,0])) / ((25 * 4)))"

And then he goes on to string a bunch of them together

Does anyone know how the equation should appear in Raster Calculator in ArcGIS Desktop?


texts: Splitting OSM roads and adding additional data

texts:
I am using OSM road data as a basis for my project. For some of the roads I have additional info. Pink roads are OSM roads and green is the road data I got elsewhere.

What I need to do is to find which roads are identical and add the extra info I got to the OSM road layer.

The problem is that the line segments are not the same and the info I want is connected to each of the line segments. Are there any tools that would split the osm line segments so that if the line segment of the other data ends then the OSM segment would also end.

Also are there any tools that compare if the segments are the same? I've been using qgis, gdal, python and postgis for this project.


texts: Transfer raster values into point shapefile in Google Earth Engine

texts:
I have a shapefile of points and I want to transfer the pixel values of a collection in GEE (e.g. Sentinel-2 max NDVI) to the corresponding points.
With polygons it seems that this can be done with a reduction, does it work also for points?


texts:
I am just wondering what I am doing wrong when trying to use the GetGeomType. I need to focus on separating layers from a pix file. Two categories are the bitmaps(which are vectors) and then the raster bands. I have imported gdal and ogr already.
Below is what I have tried so far and the errors I get when trying to use them.

    raster_image = gdal.Open(in_raster)
    band3 = raster_image.GetGeomType(3)
    print  band3

Above is how I am calling it but it gives me the error:

  AttributeError: 'Dataset' object has no attribute 'GetGeomType'

Then trying it by first getting the raster band and then checking the geom type I get:

    band3 = raster_image.GetRasterBand(3).GetGeomType()
    print band3

  AttributeError: 'Band' object has no attribute 'GetGeomType'


texts:
I've been reviewing the theory behind the UTM projection and I think I've gotten a handle on the approach, including the difference between transverse secant cylinders and tangent cylinders, standard lines and central meridians.

For any given single UTM projection, the zone only seems to apply to one side of the secant cylinder.  Does each UTM zone have its own secant cylinder, or does one cylinder include opposite sides of the reference sphere?

I've been through many of the available sources on the net, and looked through a couple GIS text books, without a clear answer.  I found the Penn State GIS textbook helpful https://www.e-education.psu.edu/natureofgeoinfo/c2_p23.html as well as the USGS Map Projections Working Handbook https://pubs.usgs.gov/pp/1395/report.pdf  but I've not yet found a source that claims a solid answer.

I'd be happy with a "Let me Google that for you" answer... I've had no luck.


texts: How to use a raster dataset to store and retrieve "vectors" using Python?

texts:
For a rectangular area of about 10 by 10 km, I have to generate "displacement vectors" with a 1 meter resolution (i.e. a delta-x and delta-y value for every cell of 1 by 1 meter). I will end up with 100 million tuples (x, y, delta_x, delta_y). In a plain text file, this would look like:

x       y       delta_x  delta_y
------  ------  -------  -------
228682  184554  0.429    1.221
228682  184555  0.428    1.222
228682  184556  0.428    1.223
228682  184557  0.427    1.224
228682  184558  0.427    1.225
...
229526  184822  -0.084   0.593
229526  184823  0.185    0.002
229526  184824  0.185    0.001
229526  184825  0.186    -0.001
229526  184826  0.186    -0.002
...

I was wondering if there is an efficient way to store these values as a dataset of type raster/coverage...? The x, y would be the pixel or cell position, and the delta_x and delta_y would be the values. Is HDF5 the way to go? I have no experience using HDF5.

I need to store the values in Python and retrieve them in Python in an efficient way (more efficient then first reading 100 million tuples from a text file, then storing them in a R-Tree index and finally using the index doing intersects).

I am looking for some guidance, tips, best practices, small examples or use cases where this has been done before.


texts: OpenLayers: disable selection from label of feature

texts:
I have a feature with a circle and a text attribute. The text attribute is displayed via a text-style above the circle center.

The current behavior:
I can select feature from text style label.

The wanted behavior:
I want to disable selection of feature from text style label

Is it possible to create such behavior?


texts: Extracting border between contiguous polygons using ArcGIS for Desktop?

texts:
I have a shapefile with polygons representing unique populations. I want to determine contiguous polygons in my shapefile, and then extract the border separating contiguous polygons as a separate polyline file. Importantly, I require the created polyline file to indicate the two polygons it separates.

I am aware of Polygon Neighbor to produce a table of all contiguous pairs, but this doesn't help me extract borders as polylines.

Consider group A in the following example map:

Group A has 8 contiguous pairs: AB, AC, AD, AE, AF, AG, AH, AI. I would like the output file to separate the aggregate border for group A into 8 separate lines, each line referencing only the shared border between each contiguous pair; e.g., the border of AH in red would be a single line in the output data. As noted before it is important that the output file also identifies the two groups each line separates.

To clarify, I require a general solution for all possible contiguous pairs in my dataset. I do not have a single group such as the example group A that I am interested in.

This question is similar to Extracting polyline of border between differing polygons using R or GRASS?, except I'm interested in a solution using ArcGIS 10.2 for Desktop. I have some knowledge of Python if an ArcPy solution is easiest.


texts:
I've written a processing script using pyqgis 3.4 which works.  The next step I would like to achieve is to load openstreetmaps as the base map from within the processing script.  I can do so easily from the Python Console, but that is not how I want to go about this task.  I have tried three different ways to load the raster layers, and can't find or think of any others.  The code doesn't give errors, but also doesn't load openstreetmap.

from qgis.PyQt.QtCore import QCoreApplication, QVariant
from qgis.core import (QgsProcessing, QgsProcessingAlgorithm,
QgsProcessingParameterRasterLayer,QgsProcessingParameterNumber,
QgsProcessingParameterRasterDestination, QgsRasterLayer, QgsProject)

import processing

class RasterAlg(QgsProcessingAlgorithm):
    OUTPUT_RASTER_A = 'OUTPUT_RASTER_A'

    def __init__(self):
        super().__init__()

    def name(self):
        return "RasterAlg"

    def tr(self, text):
        return QCoreApplication.translate("RasterAlg", text)

    def displayName(self):
        return self.tr("RasterAlg script")

    def group(self):
        return self.tr("RasterAlgs")

    def groupId(self):
        return "RasterAlgs"

    def shortHelpString(self):
        return self.tr("RasterAlg script without logic")

    def helpUrl(self):
        return "https://qgis.org"

    def createInstance(self):
        return type(self)()

    def initAlgorithm(self, config=None):

        self.addParameter(QgsProcessingParameterRasterDestination(
            self.OUTPUT_RASTER_A,
            self.tr("Output Raster A")
            ))

    def processAlgorithm(self, parameters, context, feedback):

        urlWithParams = 'type=xyz&url=https://a.tile.openstreetmap.org/%7Bz%7D/%7Bx%7D/%7By%7D.png&zmax=19&zmin=0&crs=EPSG4326'
        rlayer = QgsRasterLayer(urlWithParams, 'OpenStreetMap', 'wms')
        QgsProject.instance().addMapLayer(rlayer)

        view = QgsRasterLayer(urlWithParams, 'OpenStreetMap2','gdal')
        QgsProject.instance().addMapLayer(view)

        results = {}
        results[self.OUTPUT_RASTER_A] = rlayer

        return results```


texts: Select Layer by Attribute: selection using loops

texts:
I'm quite new with python and i'm currently stuck trying to make a loop... I have a shape file of a grid with several fields of different species (more than 150). The values of each record are 0 or 1 depending on the presence of the species on that grid. What i need to do is to make a selection of grids for each specie. So far this is my attempt:

import arcpy
import os, sys, datetime, time
arcpy.CheckOutExtension("Spatial")
from arcpy.sa import *
from arcpy    import env
arcpy.env.overwriteOutput = 1

arcpy.env.workspace = "G:\\Sara_hab_sp\\test"
arcpy.env.scratchWorkspace = "G:\\temp\\default.gdb"

arcpy.MakeFeatureLayer_management("G:\\Sara_hab_sp\\Confocc_pergrid.shp", "G:\\Sara_hab_sp\\test\\grid50")
gridlist = arcpy.ListFields ("G:\\Sara_hab_sp\\test\\grid50")

for grid in gridlist:
    print(grid.name)
    query = 'grid = 1'
    selectgrid50 = arcpy.SelectLayerByAttribute_management ("G:\\Sara_hab_sp\\test\\grid50", "NEW_SELECTION", query )

I always get an error message about the syntax of the query but i tried different things and it seems to me the error it's more related to the fact that maybe it's not possible to enter a list of fields in a query, although i have seen it's possible to do it with a list of different values of the same field...

Any suggestions on how to deal with this?


texts:
I drew a bus map of San Juan, Puerto Rico in leaflet.js.

Using the layers control I let users toggle the bus routes they want to see.  Unfortunately, the list of 45 or so bus routes is too long to display on the page.  Route 21 is hidden behind the Leaflet link, and there are more routes beyond that.

How can I let users see all their options in the layer control popup?


texts: Multiplying rasters in standalone PyQGIS

texts:
I'm trying to multiply some rasters in a standalone python environment. I have done it successfully in the QGIS python environment, and manually in QGIS, but when I run the following code, nothing occurs, nor do I receive any indication of an error.
alg_params_cdsm_filt1 = {
'INPUT_A':home_folder+'\\CDSM.tif',
'BAND_A':1,
'INPUT_B':home_folder+'\\buildings_buffered_raster.tif',
'BAND_B':None,
'INPUT_C':None,'BAND_C':None,'INPUT_D':None,'BAND_D':None,'INPUT_E':None,'BAND_E':None,'INPUT_F':None,'BAND_F':None,
'FORMULA':'A*B',
'NO_DATA':None,
'RTYPE':5,
'OPTIONS':'',
'EXTRA':'',
'OUTPUT':home_folder+'\\CDSM_filt1.tif'}

 cdsm_filt1 = processing.run("gdal:rastercalculator", alg_params_cdsm_filt1)

I took the code from the QGIS processing history when it had run successfully.
Why does nothing happen?


texts: Programmatically determine ArcGIS Layer(*.lyr) file version

texts:
I've been through IStorage, Microsoft Compound File Spec, ArcObjects deserialization etc. trying to detect the .lyr file version.

Using ArcObjects, I can get at the Layer properties, data source, extent, etc, but the documentVersion property of the LayerFileClass always reports back 'esriArcGISVersion10' on 9.3 Layer files.

ILayerFile l = new LayerFileClass();
l.Open(fileName);
LayerFileClass t = (LayerFileClass)l;
string version = t.DocumentVersion.ToString();

Any help is appreciated.

Note this crusty forum post from 2007.


texts:
I am confused regarding the OpenLayers installation. In the official page it says it is highly recommended to download using npm install ol while in this page you can download just the zip file.

So which one is the best approach ?

I used npm and it creates a node-modules folder and I am getting a message that says npm WARN enoent ENOENT: no such file or directory, open 'my-installation-path/package.json'. While going through the node modules folder I can't find the ol.js file


texts:
I have been using and experimenting with Intra-Maps Roam to collect data in the field for just over a year now.

However, I would like to build from source to take advantage of some of the recent changes that have been added. I tried with instructions from http://roam-docs.readthedocs.org/en/latest/dev-guide/Building/ , but have been unsuccessful. Currently using 2.4.2.

Any recommendations for a non-programmer?


texts:
I try to clip raster data with GDAL (in fact I am using QGIS, which produces me a command line for GDAL)
The raster has been transformed to EPSG:3857, which is also the projection of the shapefile which contains the boundaries.

gdalwarp -dstnodata 0 -q -cutline D:\GIS\data\cutline_shape.shp
 -crop_to_cutline -of GTiff D:\GIS\data\raster_3857.tif D:/GIS/data/clipped_raster.tif

The clipping processes without error, but the GDAL clips a wrong portion of the raster, it has the correct shape but is about 20 km to the north of the shapefile.


texts: Which Esri software is needed to open MXD files?

texts:
I received a group of MXD, AHD, and XML files on a CD from a client.

What is the least expensive Esri software that would be needed to view the information contained in the MXD files?

I only need to view and print.


texts: How to detect that QGIS is shutting down in a plugin?

texts: I built a plugin in QGIS and needs to stop it when user wants to exit to remove signals and close properly. How can I get an event when QGIS is exiting ?


texts: Canopy Height Model of LiDAR data in QGIS?

texts: I have a DTM and DSM files. I want to develop a Canopy Height Model (CHM) in QGIS.  Is it possible?


texts:
The word 'mosaic' is often used for quite different things within the worlds of ArcGIS and GDAL.

With GDAL and it's downstream applications like Qgis it is common to see mosaic used in the context of "assemble all these tiled images so they appear as one big image", with gdalbuildvrt and gdalmerge put forth as a convenient means of getting there.

ArcGIS has the same starting point, "assemble all these tiled images so they appear as one big image", and then goes on to add "while doing so [blend, min, max, mean, sum, ...] the overlapping pixels" and a few other things.

To put in uncharitably, in this regard ArcGIS is smart and GDAL is stupid.

However it's pretty evident GDAL is not stupid, so my question is: in the GDAL and/or Qgis tool chain what is the equivalent workflow to ArcGIS mosaic?

(note, I'm not talking about the Mosaic Dataset, which is a different beast altogether)


texts: What type of coordinates are these and how to transform them to standard WGS84 (EPSG:4326)?

texts:
I have the coordinates :
North = 861597.4129 , East = 1063313.808
All I know is that they should belong somewhere in the city of Cali, Colombia. I have tried different Colombian EPSG codes but I still do not know which format they are in.


texts: I am having problems with the zoom in QGIS 1.8.0 where a slight move of the scroll wheel on the mouse zooms way way out. It will re-render at each zoom level and with large datasets it takes over the computer for sometimes up to a minute. Is there a way to disable the mouse wheel zoom in QGIS?


texts: Styles embedded in Spatialite database are lost while copying database

texts:
I have created a Spatialite database with several layers. I have customised the style for each layer and have saved the style as the default style in the database. All works fine when I connect to this Spatialite database and load the layers.

But after I have copied the Spatialite file to a new location, loading the layers results in reversion to the default system styles, and the customised styles are lost.


texts:
How can I put x,y values in an Excel spreadsheet into Arcmap?

They do not show on the map when I try.

I think the problem is the projection. I tried coordinate system: GDA_1994_South_Australia_Lambert and the projection: Lambert_Conformal_Conic. It still doesnt work :(.

y               x
-34.9167    138.5956
-34.9211    138.6216
-34.9473    138.5827
-34.9464    138.6519
-34.9726    138.579


texts: Search for the syntax to display an HTML link in an OpenLayers popup

texts:
I use a code found on the site : http://www.untourenvelo.ch/Carte_Plein_Ecran.html . It uses the Api Openlayers 2.13.1 .
In the development of this website it is possible to display GPX data that are in XML format, for points of interest called Waypoints, when clicking on one of them appears in a popup window the name of the point.
In my use I would like to display an additional data which is included between <desc> tags and which includes a hypertext link allowing to reach either photos or an HTML page.
The display of this popup is processed in the file : "Carte_Plein_Ecran.js" : http://www.untourenvelo.ch/Carte_Plein_Ecran.js
function onFeatureSelect(event) {
        var feature = event.feature;
        // Since KML is user-generated, do naive protection against Javascript
        //var content = "
"+feature.attributes.name + "
" + feature.attributes.description;
        var content = feature.attributes.name;
        if (content.search("


texts: Converting from oblique and non-cartesian custom coordinates system with PostGIS?

texts:
I have a custom coordinates system that is:

Not cartesian: The X-axis unit is 25m & Y-axis 6.25m
Oblique: None of both axis' azimuth is North-South (0°). X-axis goes 54.2° & Y-axis 324°.
Orthogonal: There 90° between both axis.

For instance if I have those 3 points (The first column is long/lat, the second is in the custom coordinates system):
SRID=4326;POINT(-85.647 21.318) <=> POINT(0 0)
SRID=4326;POINT(-85.599 21.353) <=> POINT(0 1000)
SRID=4326;POINT(-85.794 21.497) <=> POINT(1000 0)

Can I write a custom PROJ.4 text to use with ST_Transform?
Inspired from the PostGIS documentation's example, would something like that work?:
 SELECT ST_AsText(
   ST_Transform(
     ST_GeomFromText('POINT(1000 1000)'), custom_proj4, 4326))
 FROM (
   SELECT
     '?????????????????????????????????'::text AS custom_proj4
 ) AS data;

                                      st_astext
 --------------------------------------------------------------------------------
  POINT(-85.7467905609375 21.5321592956322)

I saw this good answer for the oblique part but it is with a cartesian metric grid.
Edit:
I simplified for the example, but the real spatial reference is:
SRID=3795;POINT( 7217.670 193658.490) <=> POINT( 901   501)
SRID=3795;POINT(-7052.270 213450.630) <=> POINT(1877   501)
SRID=3795;POINT(66458.500 266451.180) <=> POINT(1877 15001)
SRID=3795;POINT(80728.440 246659.040) <=> POINT( 901 15001)

The most eastern corner is POINT(901 501). It is located in the caribbean, between Mexico and Cuba.

I already know how to convert it in my application but I'd like to do it in PostGIS by storing the whole spatial reference in spatial_ref_sys.proj4text and using ST_Transform. It could be a transformation from the custom grid to SRID 4326 or from the custom grid to geographical coordinates.