texts: Accessing geojson feature from pointtomarker in Leaflet?

texts:
I'm trying to my geojson feature 'type_id' to add my json data to different layers.

OnEachPlant = (latlngTemp) ->
  switch ''.feature.properties.type_id
    when '2'
      myIcon = L.divIcon(
        className: 'sIcon'
        html: 'S' )
      shrubMarkers.addLayer L.marker(latlngTemp, {icon: myIcon})
    when '3'
      myIcon = L.divIcon(
        className: 'tIcon'
        html: 'T' )
      treeMarkers.addLayer L.marker(latlngTemp, {icon: myIcon})

$.ajax
    dataType: 'text'
    url: 'map.json'
    success: (data) ->
      L.geoJSON(JSON.parse(data),
        pointToLayer: (feature, latlng) ->
          OnEachPlant(latlng))


texts: Label points at regular intervals using ArcGIS for Desktop?

texts:
I'm working with lines on roads represented by a large number of sequencial points.
Since labeling all of them would make the points impossible to identify, I'd like to know if there is a way to label them in a interval of each 25.

I've tried an SQL query for showing anything ending with 25 but no records were returned.


texts: Drawing the Tropics on D3.js geojson using Linestring

texts:
I am using the geoMercator projection on D3.js and want to draw the Tropic of Cancer and Tropic of Capricorn on my map.

Right now I added the Equator using

plot.append("path")
  .datum({type: "LineString", coordinates: [[-180, 0], [-90, 0], [0, 0], [90, 0], [180, 0]]})
  .attr("class", "equator")
  .attr("d", path)

I tried adding this code for adding the data to draw the Tropic of Cancer at 23.43687 degrees

{type: "LineString", coordinates: [[-180, 23.43687],[180, 23.43687]]}

It works on Geojosnlint.com but not in D3. Also tried adding more points in the coordinates but the line was curved.

{type: "LineString", coordinates: [[-180, 23.43687],[90, 23.43687],[180, 23.43687]]}

How would I add a line that goes straight across?


texts: I have NDVI layer in Earth Engine. I want to assign and calculate area of good and bad ndvi area

texts:
var mergedCollection = s2.map(function(image) {
    return image.select(['B8', 'B4'], ['NIR', 'red'])
    })

   var mergedAllFunction = function(image) {
   var ndvi = image.normalizedDifference(['NIR', 'red']).rename('NDVI');

   var thres1 = ndvi.gte(0).rename('thres1')
   var thres2=ndvi.gt(0.1) && ndvi.lt(0.3).rename('thres2')
   var thres3=ndvi.gt(0.4) && ndvi.lt(0.7).rename('thres3')
   var thres4=ndvi.gt(0.8) && ndvi.lte(1).rename('thres4')

   return image.addBands(ndvi).addBands(thres1,thres2,thres3,thres4);
   print(ndvi)
   }

How do I get area (count of pixels under each threshold)


texts: How to keep the cell size fixed when clipping a NetCDF/Tiff with a shapefile using rioxarray?

texts:
I am trying to clip a Tiff raster with a shapefile. I follow the examples provided in rioxarray (here)
The link to raster file (here) and shp file (here)
import xarray as xr
import rioxarray
import geopandas as gpd
from shapely.geometry import box, mapping

#read the shp file
geodf = gpd.read_file('shp_file.shp')

# Assign the crs
prj = [l.strip() for l in open('shp_file.prj','r')][0]
geodf.crs = prj

# Select one of the polygons
tmp = geodf[geodf['OBJECTID']==1]

da = xr.open_rasterio('raster_file.tif')

file_clipped = da.rio.clip(tmp.geometry.apply(mapping),tmp.crs)
file_clipped.rio.to_raster('file_clipped.tif')

However the resulted raster file (file_clipped.tif) has a different pixel size than the original one. The pixel size for the original raster is (463.3127165m, 463.3127165m) and for the clipped is (830.1019504m, 506.502885m). I am not sure what I do wrong. Is there any way to fix the raster pixel size when clipping?
EDIT1: Updated the codes and also uploaded the shp and raster files.
EDIT2: clipping the TIFF file with the rasterio.mask (in cases where polygons are separated, refer to the answers).
def tif_clip(tif_file, shp_file):
    from rasterio.mask import mask

    def getFeatures(gdf):
        """Function to parse features from GeoDataFrame in such a manner that rasterio wants them"""
        import json
        return [json.loads(gdf.to_json())["features"][0]["geometry"]]

    coords = getFeatures(shp_file)

    out_img, out_transform = mask(tif_file, shapes=coords, crop=True)

    out_meta = tif_file.meta.copy()
    #epsg_code = int(data.crs.data["init"][5:])

    out_meta.update(
        {
            "driver": "GTiff",
            "height": out_img.shape[1],
            "width": out_img.shape[2],
            "transform": out_transform,
            #"crs": pycrs.parse.from_epsg_code(epsg_code).to_proj4(),
            "crs": tif_file.crs,
        }
    )
    out_file = 'clip_raster.tif'
    with rasterio.open(out_file, "w", **out_meta) as dest:
        dest.write(out_img)

tif_file = rasterio.open('raster_file.tif')
geodf = gpd.read_file('shp_file.shp')
tmp = geodf[geodf['OBJECTID']==1]
tif_clip(tif_file,tmp)

This approach keeps the pixel size unchanged and also crop the extend to the extend of the shapefile. However, it still would be nice to be able to do the same thing with rioxarray for NetCDF files.


texts: Line vs. Polygon Intersection Coordinates

texts: I am working with Python, Shapely and Fiona... Is there any way where I can give in an input of a line shapefile and a polygon shapefile, and get a result of the intersection points and their coordinates? The illustration provides an accurate description....


texts: Comparing geometries using Search Cursors

texts:
I'm working on a little arcpy script to compare the geometry of features (linked by a matching ID field

I've got a couple of search cursors populating dictionaries for the two feature classes which I then am trying to compare to find where features (that should be identical shapes) don't match.

I have run into a snag in my concept though in that the output from SHAPE@ is returning False when I compare the features I know are identical (should return True).

import arcpy
fields = ['LineID', 'SHAPE@']
x = 'Line1' # Feature Class #1
y = 'Line2' # Feature Class #2
xDict = dict()
yDict = dict()
with arcpy.da.SearchCursor(x, fields) as cursor:
    for row in cursor:
        xDict[row[0]] = row[1]
with arcpy.da.SearchCursor(y, fields) as cursor:
    for row in cursor:
        yDict[row[0]] = row[1]

When I test I get False returned:

xDict[1] == yDict[1]
> False

This returns False even though I know that record #1 is identical in both feature classes.

Is there another way I can use my search cursors to compare the shape of the two features?  I know there are .WKB, .WKT, .JSON options for Shape - are these valid options for a comparison?  My brief testing suggests they are, but I'm unsure whether there are occasions where they would also return False.

xDict[1].WKB == yDict[1].WKB
> True

xDict[1].equals(yDict[1]) # Should be True
> True
xDict[2].equals(yDict[2]) # Should be False
> True
xDict[4].equals(yDict[4]) # Should be False
> False


texts: Activate Map on the click of a button ArcGIS Pro SDK

texts:
In this case, when an ArcGIS Pro project is opened, non map panes are active. I would like the map pane to be activated on the click of the button as proceeding code requires the MapView.Active to not be null. Inside the OnClick function, I am attempting to activate the map pane but it is not activated in time. It's activated later on which causes errors as the active map view returns null for code that runs after this.

internal class MyDockpane_ShowButton : Button
{
    internal static void activatePane()
    {
        var mapPane = FrameworkApplication.Panes.OfType<IMapPane>().First() as Pane;
        mapPane.Activate();
    }
    protected override void OnClick()
    {
        activatePane()
        MyDockpaneViewModel.Show();
    }
}

I believe that I need to await the activatePane() function so I re-wrote the code as below:

internal class MyDockpane_ShowButton : Button
{
    internal static async Task activatePane()
    {
        await QueuedTask.Run(() =>
        {
            var mapPane = FrameworkApplication.Panes.OfType<IMapPane>().First() as Pane;
            mapPane.Activate();
        });
    }
    protected override async void OnClick()
    {
        await activatePane();
        MyDockpaneViewModel.Show();
    }
}

With the new code, I get System.InvalidOperationException: "The calling thread cannot access this object because a different thread owns it". What am I doing wrong here?


texts: Convert Text to Mtext and rearrange with QGIS

texts:
I have multiple names and parcel data in text form in AutoCAD.  I want to make them in the format as follows:

N/F
First and Last Name (Same as Shown)
T.A #Tax Number (Add T.A. # to each then keep the numbers the same after that)

I feel like there is a fast way I am missing instead of manually changing all of them by hand.

The text has no attributes with them that I need to keep.

I tried to export them to QGIS and have them all in but each text has multiple dots with information associated with it.  I need to only keep the top two in each circumstance and then flip them somehow.


texts:
Is there a specification for the protocol OpenLayers uses to communicate with the server?
I'm trying to reverse engineer what happens on GET, POST, PUT and DELETE, but I don't like to do it this way.

Example:

var vectorLayer = new OpenLayers.Layer.Vector("Vector Layer", {
    protocol: new OpenLayers.Protocol.HTTP({
        url: '...',
        format: new OpenLayers.Format.GeoJSON({})
    }),
});


texts: How to use the replace function to replace "specific" values to another? QGIS

texts:
This question is actually similar to another question I had asked, but this is more a response question, since I am limited to only one question per topic. SO PLEASE READ THIS FULLY.

"Using replace function to replace value to one of value map using QGIS field calculator?."

Also this problem I underfound is not related to the situation I asked that question in. The problem is that the replace function also changes the text that is has in "common" of the field value I want to replace.

I discovered a new problem with the replace function and that is it also changes the text of value's that is "common text".

Let say I have I want to change the text '2015', '2016' and '2017' to 'n.v.t.'. This works great with the replace function, however there is another text called VJ2015. When I execute the 2015 replace command it also changes this to  'vjNULL'. Again not what I want, since I only wanted that "specefic text value" to be changed not the text that is "common" with other value's. So here is my case expresion':

replace("VJ/NJ (LC)", '2013' and '2014' and '2015', NULL)

I know that programs like Excel and others have "Identical letter option", I suspect that there is a way for QGIS to turn this option on.

I also understand that "" stands for column and '' stands for the value, so I suspect I need to type some specific special character to "Indenitical on" but I do not know which. So what do I need to do?


texts: How to create an OD road distance matrix in QGIS?

texts:
I'm looking for the steps required to produce an origin-destination (OD) cost matrix.

What I'm looking for is a matrix of road distances. The distance matrix tool only produces euclidian distances. The QGIS Roadgraph plugin produces the shortest road distance from one origin to N destinations. What I need is a file containing distances from all origins to the N closest destinations. I have a lot of origins so repeating the analyses with Roadgraph is not an option. Thanks again.

I'm looking for a step-by-step description, starting from scratch (ie, with my 3 shapefiles: origins, destinations, road network).


texts: How can I extract flat surfaces from a DEM with GRASS GIS?

texts: How can I use GRASS GIS to get the flat surfaces in this DEM into polygons? The picture is used to show the clearly defined flat surfaces.


texts:
Using pyQGIS, is there a way to hide the loading gml dialog when making WFS http requests via a plugin (for QGIS versions >2.4):
After running this statement the load dialog shows:
QgsVectorLayer(url, "location", "WFS")

I would like to hide this dialog in favor of using the QtGui.QProgressDialog instead, is this possible?


texts: Editing 2007 MS Access *.accdb using OLE DB connection?

texts:
Can you edit a 2007 MS Access .accdb file format using the OLE DB connection?

I have made the connection using the following:

Connecting to a 2007 Microsoft Access database (.accdb) in ArcGIS

I can add and view the tables from this file in ArcMap, but when I try to start an edit session on this table, it errors out with the follwoing 2 errors:

No editable layers
Check to see that you have the appropriate permissions

I created a .accdb file from scratch, so it is not read only.

Under the Advanced tab in the Data Link Properties for the OLE DB connection, I have set the ReadWrite option check mark and unchecked the rest and I am still unable to edit this table in ArcMap.

Click image for larger version Name: OLE_DB_Advanced.jpg Views: 14 Size: 44.7 KB ID: 25411

There isn't much documentation on the OLE DB connection in the help file other than how to connect as posted above. Can anyone point out what I may be doing wrong?


texts: Google Earth Engine: Evapotranspiration Linear Regression MODIS

texts:
I am trying to do a Linear Regression about Evapotranspiration using MODIS and I get this error when I try to inspect the value of the pixels in the result layer: "Error Array index 1 out of bounds.  Expected value between 0 and -1, found, 0."

The code is the next:

// Simple regression of year versus NDVI.

// Define the start date and position to get images covering TDPS,
// Bolivia, from 2000-2010.
var start = '2000-01-01';
var end = '2014-12-31';

var region = table2

// Filter to Collection
// time of year to avoid seasonal affects, and for each image create the bands
// we will regress on:
// 1. A 1, so the resulting array has a column of ones to capture the offset.
// 2. Fractional year past 2000-01-01.
// 3. SAVI.
// 4. Mask Clouds
var col = ee.ImageCollection("MODIS/NTSG/MOD16A2/105")
  .filterDate(start, end)
  .filterBounds(region)
  .map(function(image) {
    var date = ee.Date(image.get('system:time_start'));
    var yearOffset = date.difference(ee.Date(start), 'year');
    return ee.Image(1).addBands(yearOffset).toDouble();
  });

print(col)
// Convert to an array. Give the axes names for more readable code.
var array = col.toArray();
var imageAxis = 0;
var bandAxis = 1;

// Slice off the year and ndvi, and solve for the coefficients.
var x = array.arraySlice(bandAxis, 0, 2);
var y = array.arraySlice(bandAxis, 2);
var fit = x.matrixSolve(y).clip(region);

// Get the coefficient for the year, effectively the slope of the long-term
// NDVI trend.
var slope = fit.arrayGet([1, 0]).clip(region);

Map.addLayer(slope);

// Export the image, specifying scale and region.
Export.image.toDrive({
  image: slope,
  description: 'ET_slope_MODIS',
  scale: 1000,
  maxPixels: 1e9,
  region: geometry
});

Does anyone know which is the main issue?


texts: Using GeoPandas to convert Linestring UTM data to LatLong

texts:
I have a few shapely.linestring objects that are read into geopandas with the call:

piped = gpd.read_file("Pipes.dbf")

In the geometry column the values are in UTM. I have identified the sector as 44Q and I can use utm to extract the lat/long for a single value. I'm pretty new to using Geopandas but was hoping there was an .apply where I could duplicate and convert the linestring objects into a list of tuples and store them in another column.

Example row value: LINESTRING(376753.3276983045 2042051.6105498213, 376764.03744285746 2042051.728052152)

I'd like to have a corresponding column called say lat_long that contains the same information from the linestring column but in a the format [(18.465193133286725, 73.83275226395054),(18.465194819587577, 73.83285367143624)]

def utm_to_latlon(vals, utm_lon, utm_lat):
"""Function to convert UTM TO lat long
   returns a tuple of latlong
"""
    return utm.to_latlon(vals[0], vals[1], zone_number=utm_lon, zone_letter=utm_lat)


texts: Altering graticule intervals logically in varying map scales using ArcPy and/or Data Driven Pages?

texts:
This question appears on many sites and I was wondering if there was a script that is now available that resolves it?

Please also upvote the ArcGIS Idea that Reference grids/graticules should be accessible to arcpy.mapping.

My workarounds have been to create multiple mxd's that just hold the same scale map sets and then export only the relevant pages to that scale but this adds a lot of time to the process when changes are required as they have to be done on each mxd.

The other related issue is with the centering of elements where when a map is 1:10000 and another is 1:100000 one of them will be off centered from other elements such as the "N" Arrow.


texts: ESRI GDB Support in GDAL / ogr2ogr on Mac OS X 10.8

texts: How do you incorporate support for ESRI's fgdb format into GDAL / ogr2ogr on Mac OS X 10.8?


texts: Combining several layers in ArcGIS Desktop using Intersect gives WARNING 000117: Warning empty output generated?

texts:
I'm using ArcGIS 10.3

I have the following 4 layers (all are polygons)

You can find these layers in this link

DATA

Sites: 5 subcatchments (Site1 to site5) and site 5 combine all
subcatchments upstream
Nloss: nitrogen loss for each polygon
Geology: Rock types
Soils: Soil types

What I want?

I want to create a polygon shapefile. In this layer, each polygon will have ALL the attributes from all the 4 layers.

I used intersect to do so. I think union could be used as well to do the same.

I named the resulting shapefile intersect_4_layers.shp. The result of this intersection could be found in this link.

Adding slope

I have slope shapefile that has two categories (0 and 1).

zero: if slope is less than or equal 10
1: if slope is more than 10

I wanted to add slope information to the shapefile resulting from intersecting the 4 layers above (Sites/Soils/Nloss/Rocks) named "intersect_4_layers.shp".

I tried the following:

1- to intersect slope.shp with intersect_4_layers.shp

2- to intersect all 5 layers (Sites, Nloss, Geology, Soils and Slope)

Both trials failed and resulted in this warning

  WARNING 000117: Warning empty output generated

QUESTION

Why intersection worked the 1st time (intersecting 4 layers) and didn't work while trying to intersect them with the slope?

Any suggestions how to get all the attributes of the 5 layers into one shapefile?


texts: Some zones disappear when running Zonal Statistics as Table in ArcGIS

texts:
I'm new to GIS. I want the sum of nighttime lights (DMSP) by district for some country. I'm running this arcpy.gp.ZonalStatisticsAsTable_sa(Limites_distritales, "IDDIST", rc, temp_table, "DATA", "") and then converting to .csv. However, when I check the output some districts (identified by IDDIST) are missing. I've checked all my script and they shouldn't be dropped; my guess is that they are being assigned NoData because they overlap multiple grids (and no single grid is in the polygon) and thus no statistics are being calculated for them.

One of the problematic districts is shown in the picture below.

Is there a workaround for this?


texts:
Tiles are rendered using PostGIS as the store. Is it possible to output the tiles using a raster like style with 16bit values instead of 8bit values?
The goal is to provide more depth in the grayscale. The colouring is done on client side using WebGL. Because of the 8bit png output only 255 bins are available and this is limiting.


texts: Change annotation symbol properties in bulk?

texts:
How can a person change annotation properties in bulk?

For example: take the 4500 records in "Waterbody text" annotation feature class and change all occurrences of "Times New Roman", "Cyan" to "Revival", "Cretan Blue".

I don't want to use layer symbol substitution if it can be avoided (and Esri recommends against it for this anyway: "...is not designed to be a wholesale replacement of annotation symbol editing and symbol management")


texts:
I'm trying to unpack shapefiles in C# using GDAL-OGR library, but I'm confused about the structure of multipolygon with holes in OGR. As far as I know, the models of shapefiles and OGR geometry (the latter following simple feature specifications) are different. Shapefiles have a "polygon" type that might be in OGR term either a polygon or a multipolygon. So I must handle both cases. When reading a shapefile, OGR does some topological analysis to figure out the spatial relationship between rings and determine if they must be outer or inner rings.

An OGR polygon is a collection of rings, the first one being the outer ring and the following ones the inner rings. So, this code will work fine.

Feature feat;
while ((feat = layer.GetNextFeature()) != null)
{
     var geom = feat.GetGeometryRef();
     for (var i = 0; i < geom.GetGeometryCount(); i++)
     {
          //First one is the outer ring and, all the next ones are the inner rings
          var ring = geom.GetGeometryRef(i);
          ...
     }
}

But, if the geometry is multipolygon, e.g 2 polygons with holes. How do you determine which one is the outer ring of the second polygon? Because if I understand it correctly, geom.GetGeometryRef() will only return the first polygon. Shouldn't Feature have something like feat.GetGeometryCount() (to determine the number of polygons - the outer rings only) and feat.GetGeometryRef(i) (to get the i outer ring of i polygon) or something? I checked, and there is none. How do you handle such case? Can anyone help me with it, please?

Also, how many nested levels that OGR support for polygons (i.e, holes inside holes)? How do handle such cases?

Sorry for my English! I'm not a native speaker.


texts: Extracting raster values to point feature class failing in ArcMap?

texts:
I need to extract values to point features from Mosaic dataset made out of  DEMs (Whole world SRTM) . My feature class has 18,000,00 points.

When I ran the below tools, only <20K points receive elevation values ,rest all are null though DEM at those location is not null.

Spatial Analyst

Extract Values to Points
Extract Multi Values to Points

3D Analyst

Add Surface Information

Per suggestions in https://community.esri.com/  , I ensured below. None seems to work .

CRS of DEM and Feature Class is same
Do not open the fc and moasic  in ArcMap
Use Spatial Analyst or 3D Analyst tools in Catalog
Put FC and Mosaic Dataset in a common Geodatabase
Run using ArcPy

Often I see claims of of this tool working after a relevant patch to 10.1 . I am using latest version 10.6 and seems to have no solution to this problem.

My production system is built on ArcMap so using ArcGIS Pro is not an option for me.

How the tool is executed ?

Tool Boxes > System Tool Boxes > Spatial Analyst Tools > Extraction > Extract Values to Points

Input point features --> D:\Data\Terrain.gdb\points
Input raster  --> D:\Data\Terrain.gdb\SRTM3
Output point features --> D:\Data\Terrain.gdb\points_elev

Environment of the ArcMap

OS : Windows 10
Virtual Technology : None
Webserver: None (using only desktop)
Service Pack: None
DBMS : File geodatabase (for storing Mosaic  dataset and feature class)
Virtual Technology : None
DBMS OS: Windows 10
DBMS Version: File Geodatabase

ArcPy Script

import arcpy
from arcpy import env
from arcpy.sa import *

# Set environment settings
arcpy.env.workspace = r"in_memory"
arcpy.env.overwriteOutput=True
arcpy.CheckOutExtension('Spatial')

# Set local variables
inPointFeatures = "D:\Data\Terrain.gdb\points"
inRaster = "D:\Data\Terrain.gdb\SRTM3"
outPointFeatures = "D:\Data\Terrain.gdb\points_elev"

# Execute ExtractValuesToPoints
ExtractValuesToPoints(inPointFeatures, inRaster, outPointFeatures,
                      "INTERPOLATE", "VALUE_ONLY")

What happens with subset of points ?

When we experimented with subset say - 600 points or 20,000 points - the same issue is seen. Some gets elevation values(correct values), the others have null values.
Neither the toolbox nor the Arcpy script throw any error.
We noticed ,using one raster(Geotiffs on disk) at a time is giving us correct values. This is however taking 12 hours or more - not suitable for task at hand.


texts: Are there limitations for QGIS Composer images larger than 32768 pixels?

texts:
We are outputting a map file with composer : size 45,312 pixels width and 4,352 pixels height in TIFF and PNG format.

The map is calculated but in the file we have only 32,768 pixels width show, the rest of the image is blank (actually transparent due to the alpha chanel).
The output file size is 45,312x4,352 as demanded.

There's a limitation for QGIS on image calculation width (and/or height) ?

We are running on QGIS 2.8.2 / Ubuntu 14.04 LTS 64bit / 32GB RAM

Another limitation : look like the file size cannot exeed 1.5 GB...

Does those limitations still exist in QGIS 2.10 ?


texts: Splitting polygons and getting areas to split also

texts: I have been provided with a shapefile with certain land areas to split. How do I split the polygons so that the new areas are calculated and are listed in the attribute table? Currently I've been able to split the polygons but the area for each part is still listed at the same area of the original polygon in the attribute table before it was split.


texts: QGIS 2.6: GPX with 1546 records only reads 230

texts: I have a gpx with 1546 waypoints. But: only 230 waypoints are read in QGIS 2.6 (using "add vector layer"), the rest is skipped. Can't understand why. File structure seems ok since the first 230 waypoints are read.


texts:
I installed mapnik and created XML files for creating map. Using a python I created the image corresponding to the XML File point to a postgis enabled postgresql database. Then I downloaded the tilecache and configured the .cfg file as follows,

[Outer]
type=Mapnik
mapfile=/var/www/html/mf/outer.xml
bbox=291627,904686,958569,1426831
srs=EPSG:32643
debug=no

[cache]
type=Disk
base=/var/www/html/cachedata

I configured in apache 000-default.conf file as,

<Directory /var/www/html/tilecache>
   AddHandler cgi-script .cgi
   Options +ExecCGI

</Directory>

Then I run the url.

http://localhost/tilecache/tilecache.cgi/1.0.0/Outer/0/0/0.png

got 404 Page not found error

http://localhost/tilecache/tilecache.cgi shows the text code itself.

No tiles creating and  given the cache folder full permission.  Please help.


texts: Converting GRIB to CSV or to binary data

texts:
The challenge is to extract multidimensional data (weather and ocean data in this case) from .grib files, without having to

download and install any new software on the computer by running installation programs,
download and compile any new source code,
rely on pre-installed software on the computer (like java, python or installation packages)

so that

Windows 10 is supported (without Linux emulators),
float64 is supported,
it does it "bulk", ie. not only one record ("band") at time (which can take a full day for a large dataset, at speed of for example 1 band/second). Instead, it shall do the full extract with one single command (or only a handful of commands), and
it outputs the data from the grib to a named file, as binary or as readable ASCII.

Ie. needing a downloadable command-line tool(set) that (after virus testing) extracts the data from .grib files, for further processing!
The GRIB format

GRIB (GRIdded Binary or General Regularly-distributed Information in
Binary form1) is a concise data format commonly used in meteorology
to store historical and forecast weather data.

https://en.wikipedia.org/wiki/GRIB


texts: Expose the Coordinate System Parameters in FME

texts: I would like to take fields in the Coordinate System Properties. (Data Inspector) to use as user attributes.


texts: Shouldn't QGIS read the proj definitions from the spatial_ref_sys table?

texts:
I have several layers in PostGIS that use a custom projection. I added the projection SRID and it's definition to the spatial_ref_sys table but, when I load the layer in QGIS, I get the "Coordinate reference System Selector" dialog. Shouldn't QGIS automagically read the proj definitions from the spatial_ref_sys table?

The spatial_ref_sys table is in the public schema and the layers are on a different schema. The custom SRID that I added with code 90001 is:

  'PROJCS["Brazil / Albers Equal Area Conic (WGS84)",GEOGCS["WGS
  84",DATUM["WGS_1984",SPHEROID["WGS
  84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["longitude_of_center",-48.5],PARAMETER["standard_parallel_1",-20.5],PARAMETER["standard_parallel_2",-24.5],PARAMETER["latitude_of_center",-25.0],UNIT["Meter",1.0]]'


texts: Accessing older Bing imagery in ArcGIS Desktop?

texts:
I've been using the Bing basemaps in ArcGIS 10. I'm studying bird ecology and have been overlaying bird GPS locations on to Bing imagery to see where the birds go and what habitats they chose.

The Bing imagery has been super clear, however, Bing has just recently updated the imagery for my study location in New Zealand and now the imagery is very dark and almost unusable.

I was wondering if there is a way to access previous Bing imagery? I previously downloaded the Bing aerial layer but this has automatically updated to the newer images too.


texts: Calculate field : find the minimum value from other fields, but exclude the zero

texts:
Using solution from post : https://stackoverflow.com/questions/29284929/rejecting-zero-values-when-creating-a-list-of-minimum-values-python-field-calc

My goal is to populate my "least" column with the minimum value of 4 "time" columns and to exclude zeroes when calculating.The "least" column format is double and my "time" column format is long integer.

This code seems to work for others but not for my table, I've spent a whole day on this and can't find the culprit.

Here is my pre-logic script code:

def my_min(time, time_1):
  lst = [time, time_1]
  return min([x for x in lst if x !=0])


texts:
I am running the nearest neighbour algorithm in QGIS 3.4.3.

I have tried various parameters but the problem persists. I have also tried other algorithms (simple Kriging) and had the same problem. However, I did succed running the IDW interpolation algorithm from the QGIS plugin itself.

this is the message log :


texts: JPG file all black when added as raster in QGIS

texts:
I add a new raster layer in the latest QGIS and select a JPG file, but it is all black. I tried adding the TIF version of the same image and that worked.

What is it with the JPG, and how should I import it correctly? Thank!


texts: Processing large amounts of Rasters in R

texts: I have 6500 processed and clipped rasters that indicate snow levels on a daily basis in the northern hemisphere. I need to create a raster that represents a sum of the snow for the season. The winter season would then go from September 1999 to march 2000. I was just wondering how I would go about recalling the daily rasters into R using the loop function/ how i would go about calculating the seasonal sums. Keep in mind that the rasters span two different years and that I don't want to include all 1999 snow data because that would include the previous winter(98/99) as well.


texts: Data type gets converted from "qlonglong" to "int" when creating Virtual Layer

texts:
I have a source layer with ID column as qlonglong. When I execute a query in DB Manager the types seems correctly preserved.
The problem is that if I "Load as New Layer" the resulting Virtual Layer has qlonglong converted to int.
As a consequence I cannot use this ID field to join to other tables.
It seems a very basic problem, but I've spent more than two hours and cannot find any better workaround than Refactoring the resulting Virtual Layer. I would love some help fixing the root problem.
Here is a link to example file
Here is the Virtual Layer dialog, trying the cast suggestion. But I think type conversion happens not on source read, but on layer write.


texts: Spatial Join In PostGIS where centroid of a polygon is inside another polygon

texts:
I have two overlapping polygon layers as shown in figure below. One with small grid polygons and another with a large irregular polygon. I want to join them such that attributes of the large polygon is joined to the each grid polygon whose centroid falls inside the large polygon. So in most cases it does but at the boundary it may not be so.

This can be accomplished for example in ArcGIS with spatial join where match criteria is Input feature "HAVE THEIR CENTER IN" join feature. However, since the data is huge and not manageable with ArcGIS I prefer to use PostGIS.

Is there any function in PostGIS which fulfills the match critera 'HAVE THEIR CENTER IN'. I could not find one so I tried the below query but it does not return any join result.

select g.*, p.ags, p.rs
from grid_polygon as g
left join large_polygon as p
    on st_intersects(st_centroid(g.shape),p.shape)

What can I do to correct this?


texts: Change shapefile coordinate system using Python

texts:
I am struggling with Python and Kartograph.py. I'd like to change coordinate system of the whole shapefile from EPSG:5514 to EPSG:4326. Found some code in here working for points, but I don't know how to cope with shapefile. Code for point goes like this:
from kartograph import Kartograph
from pyproj import Proj, transform

K = Kartograph()
# S-JTSK EPSG:5514
# Kartograph uses EPSG:4326

inProj = Proj(init='epsg:5514', preserve_units=True)
outProj = Proj(init='epsg:4326')

x1,y1 = -599411.949672, -599411.949672
x2,y2 = transform(inProj,outProj,x1,y1)

Do I need to do anything like this for shp?
from kartograph import Kartograph
from pyproj import Proj, transform

K = Kartograph()
# S-JTSK EPSG:5514
# Kartograph uses EPSG:4326

inProj = Proj(init='epsg:5514', preserve_units=True)
outProj = Proj(init='epsg:4326')

for point in shapefile:
    SX = #coordinateX in shp attribute table
    SY = #coordinateY
    x1,y1 = -599411.949672, -599411.949672 #do something to transform every point
    x2,y2 = transform(inProj,outProj,x1,y1)


texts: Change attributes' values based on geometry relation

texts:
MY GOAL IN A NUTSHELL

I am trying to calculate the angles between lines. For doing this, I calculated the lines' slop respect to x axis and now I am trying to subtract them to find the angle inside. I implemented the equation in Python (and for now I don't want to deal with modules like shapely or fiona but I am looking
for simple Python scripts to get confidence).

PROBLEM

In the attribute table, data only for one feature is filled.

SCRIPT TRIAL

import itertools
layer.startEditing()
layer.dataProvider().addAttributes([QgsField("ang_vertex",QVariant.Double)])
layer.updateFields()
layer.commitChanges()

layer.startEditing()
field=layer.dataProvider().fieldNameIndex('ang_vertex')
print field
for geom1,geom2 in itertools.permutations(layer.getFeatures(),r=2):
    if geom1.geometry().touches(geom2.geometry()):
        equation=geom1[3]-geom2[3]
        print feat[2], equation
        layer.changeAttributeValue(feat.id(), field, equation)
layer.updateFields()
layer.commitChanges()


texts: The raster or shapefile exported from Ecognition does not align with the original raster when I open it in QGIS


texts: I'm using QGIS and I created a grid and joined a column with observation points to it.  I'm trying to use symbology to shade the cells of the grid based on the values that I joined.  But it looks like it just punched a hole in the grid at each of cells instead of coloring them.  Does symbology not work normally with grids?


texts:
I have the same question like the one aked here:
How to get the “absolute value” of a raster cell in QGIS?

I came up with the same solution for the problem some time ago. However I showed this solution to a software developer some days ago, and he was very amused. Is there a build in function in QGIS to acquire the modulus of a real number?


texts: Create Table (Data Management), Unable to Save to Feature Dataset

texts:
I'm trying to create a table in a feature dataset and get the error: Parameters are not valid. The code works when I write only to the gdb. Am I not allowed to write to a feature dataset or am I doing something wrong?

import arcpy
import datetime

now = datetime.date.today()
exoticPlant = "Kudzu"
outPath = r"V:\chat_gis\gis_projects\Vegetation_Managment\Data\Derived_Data\Derived_Data.gdb\Report_01_14_2013"
exoticReportTemplate = "\MASTER_Data\VegetationManagment_Chatt.gdb\Exotic_Report_Template"
exoticReportName = exoticPlant + "_UnitCoverage_5years_" + now.strftime("%m%d%Y")

arcpy.CreateTable_management(outPath, exoticReportName, exoticReportTemplate)


texts: GeoExplorer: "Could not create layer record for layer group"

texts:
I installed GeoServer via OpenGeo Suite, and still working on it.

I created a "Layer Group" and add layers from local PostGIS. I tried to preview it from "Layer Preview" and show perfectly through the "OpenLayers" option, but unable to show in "GeoExplorer".

I got the javascript warning message says: "Could not create layer record for layer 'testgroup'. Check if the layer is found in the WMS GetCapabilities response." on GeoExplorer.js:1441

To view the source code shows like:

1441: {name:"fixed",type:"boolean"},{name:"selected",type:"boolean"},{name:"restUrl",type:"string"},{name:"infoFormat",type:"string"},{name:"getFeatureInfo"}];c.fields.each(function(a){i.push(a)});b=new (GeoExt.data.LayerRecord.create(i))(h,b.id);b.json=a}else window.console&&0<this.store.getCount()&&void 0!==a.name&&console.warn("Could not create layer record for layer '"+a.name+"'. Check if the layer is found in the WMS GetCapabilities response.");return b},getProjection:function(a){var b=this.getMapProjection(),
    "Could not create layer record for layer 'testgroup'. Check if the layer is found in the WMS GetCapabilities response."
1442: c=b,a=a.get("srs");if(!a[b.getCode()]){var c=null,d,e;for(e in a)if((d=new OpenLayers.Projection(e)).equals(b)){c=d;break}}return c},initDescribeLayerStore:function(){var a=this.store.reader.raw;this.lazy&&(a={capability:{request:{describelayer:{href:this.url}}},version:this.version||"1.1.1"});var b=a.capability.request.describelayer;if(b)a=a.version,1.1<parseFloat(a)&&(a="1.1.1"),a={SERVICE:"WMS",VERSION:a,REQUEST:"DescribeLayer"},this.describeLayerStore=new GeoExt.data.WMSDescribeLayerStore({url:this.trimUrl(b.href,

I am sure that I can preview them separately in GeoExplore from "Layer Preview" page.

It seems everything goes well when I checked it from http://localhost:8080/geoserver/wms?SERVICE=WMS&REQUEST=GetCapabilities, without any  error massages or exceptions shown, and the layer group description show at the bottom of the page.

Does anyone run into the same problem?

//UPDATE

I found that I missed a javascript error massage when I view the REQUEST=GetCapabilities page.

It says: "Uncaught TypeError: Cannot set property 'boxSizing' of null" on ContentVeil.js:11

To view the source code shows like:

11: veil.style.boxSizing = "border-box";
    "Uncaught TypeError: Cannot set property 'boxSizing' of null"
12: veil.style.borderStyle = "solid";

I am wondering where the error comes from.

Furthermore,

I cannot load the "field list" when editting limit condition in style's rules in GeoExplorer.


texts:
I can't figure out how to wrap the text on the leaflet tooltip. I dont want to display the marker so its hidden, but i want to the text to wrap the text in the tooltip if possible

.myTooltip {
 background: green;
 border: 2px solid cyan
}

.leaflet-tooltip-left.myTooltip::before {
 border-left-color: cyan;
 }
.leaflet-tooltip-right.myTooltip::before {
  border-right-color: red;
 }

then adding to the map, i've change the properties of the no wrap but doens't seem to effect anything. Also, for some reason the background color doens't seem to be changing, but thats a different questions

  var marker = L.marker(new L.LatLng(item.Latitude, item.Longitude), { opacity: 0 }); //can be set to zero
                marker.bindTooltip(item.Text, { permanent: true, noWrap: true, className: "myTooltip", offset: [0, 0] });
                marker.addTo(mymap);


texts: How do you provide a magnified view within layout view?

texts:
I want to label each land parcel and realize some groups of parcels are too small and close together to label effectively. Therefore, I want to create a magnified view to include on the sheet for this region.

Is there a way to do this?


texts: Easy way to map & share some spreadsheet points on a web map

texts:
What is an easy solution for placing some points on a web map and sharing the map and info with others from a spreadsheet I have?  I want to be able to change the spreadsheet info and have the map reflect the changes.

The area of interest is Vietnam.  I have been able to export the points to a Google Earth file and maybe that is the easiest solution.


texts: Problem creating polygon features using ArcObjects / ArcGIS 10

texts:
I need to create a grid of polygons. Each polygon is rectangle, thus has 4 corners (vertices). When I create polygon features using ArcObjects they seem to have wrong geometry, which results in wrong label placement (outside polygons).

I use the following approach:

Dim pPntColl As IPointCollection
pPntColl = New Polygon

Dim pPoint As IPoint
pPoint = New Point
pPoint.PutCoords(dX, dY)
pPntColl.AddPoint(pPoint)

...

When I use 4 vertices I get a polygon which has one segment missing (visible in Edit session).

When I add the first point at the end of polygon (as a 5th point), geometry looks good, but the label is still off.

I know I can use ITopologicalOperator2 and Simplify() method but:

This affects the performance (I need to create many thousands of polygons).
I don't understand, what is the reason of the problem. I have 4 points, they should create a topologically correct polygons...

So, my question is: what is the best way of creating polygon features within a polygon feature class using ArcObjects?

I already checked some ESRI recommended approaches, e.g.:
http://help.arcgis.com/en/sdk/10.0/arcobjects_net/conceptualhelp/index.html#//0001000002wm000000
but the problem is still the same: I need to use Simplify() to get polygons right.

Any help and hints how to optimize this will be appreciated.


texts:
I am trying to show map based on the gridset defined in geoserver, in such a manner that map zooms to that Scale only.

Scale is as follow : 1:1000, 1:500, 1:400, 1:300

        view = new ol.View({
            center: [564878.4195196485, 4928488.1371189915],
            zoom: 16,
            projection: projection
        })

        map = new ol.Map({
            controls: ol.control.defaults({
                attribution: false
            }).extend([mousePositionControl, scale]),
            layers: [L1, L2, L3, L4],
            target: 'map',
            view: view
        });


texts: Flip orientation of Font Markers on line when upside down using QGIS

texts:
I would like to have font markers flip when they are upside down but otherwise follow the line angle.
Here is an image of how they look now:

I would like them to look like this:

I achieved this by switching the line direction but I want to do it through the style and I would like to do it at the marker level not the line level.
I have tried a geometry generator to switch the line.
if(y_max( $geometry)>y_min($geometry),$geometry,reverse($geometry))
But this only switches the issue to the bottom line.

I also tried a data defined override on the font marker rotation with these expressions.
if( @symbol_angle < 180,0,180)
if(y_max( $geometry)<  y_min($geometry),0,180)


texts: Expand or Collapse Layer Groups in PyQGIS3

texts:
PyQGIS2 has a function to expand/collapse a layers with a subset of styles in the Layers panel as follows:

self.iface.legendInterface().setLayerExpanded(self.segments_layer, False)

I know that in PyQGIS3 self.iface.legendInterface() has been replaced by gsProject.instance().layerTreeRoot() but making that switch in this case produces an error:

  AttributeError: 'QgsLayerTree' object has no attribute 'setLayerExpanded'

What is the appropriate replacement for setLayerExpanded in PyQGIS3?


texts: How to create an SVG marker symbol layer in QGIS?

texts:
I am using Qgis 2.14.3-Essen on a MacBook Pro running OSX 10.8.5.

I have a number of SVG files stored in a directory named QgisSymbols.

I have set the path to the Qgis directory in Settings/Options/System/SVG paths.

I have also added the svg file directories to the svg directory within the QGIS application (Contents/Resouces/svg).

On the properties/Style tab when I select the symbol layer type as 'SVG marker' I don't see any SVG groups? I just get the message 'This layer doesn't have any editable properties' - even when the layer I am trying to change the symbol for is editable?

I have tried this with a line layer and a point layer and get the same result for both. I want to use SVG symbols for lines and points but can't.


texts: QGIS raster calculator ^2 operation gives negativs results

texts:
I want to compare 2 same tiles at different dates.

I use the raster calculator like this :

(IR_band_frb - IR_band-aug)^2

But I get this result :

The data range goes from -6049 to 4486 ... Why I get negativ results ?

(Same when I write "abs" instead of "^2" symbol)

I tried

sqrt((im1-im2)**2)

But I get missing values (transparent "colors" in the output raster). So negatives values are passed to the square root operation.

Maybe do you have a better way for comparing 2 rasters in QGIS, do not hesitate !

Ty !


texts: Geoserver WFS GetFeature returns no features when reprojecting PostGIS layer

texts:
I may be misunderstanding the GeoServer documents, but when accessing a GeoServer 2.2 WFS 1.1 layer, should I be able to reproject to an arbitrary SRS via the SRSNAME parameter?  I'm generating a bounding box from the location of a click on a WMS map where I know there is a feature.  The map is projected into UTM zone 15 (EPSG:32615).  So I generate the following URL:

http://server.domain.com/geoserver/MyWorkspace/wms?SERVICE=WFS&VERSION=1.1.0&REQUEST=GetFeature&TYPENAME=MyWorkspace%3AMyLayer&MAXFEATURES=9&SRSNAME=EPSG%3a32615&BBOX=382138.16774714%2c4587006.11489167%2c382138.16774714%2c4587006.11489167&

This query returns 0 features, and I don't know why.  It's essentially the same query I was using for another layer, also served from GeoServer, that was sourced from a shapefile in EPSG:32615, where the query worked.  So here's what's different about this layer:

I'm using PostGIS for a backend on the new layer; the old one was a shapefile
The shapefile was in EPSG:32615; the PostGIS table is in EPSG:4326
The PostGIS table is using the "geography" data type rather than the "geometry" data type, if that's relevant.

Any ideas?


texts: Labeling expression to label one value in a attribute table column

texts:
I am using QGIS 2.0 on Windows. I have a shapefile for a point layer that has the the following attribute table:

I am trying to label using the Status column and only want the value that equals Office Location to show on the map. I tried the expression Status = 'Office Location' and I get either a 1 or a 0 show on the map. I cannot figure out the appropriate expression to do what I want and I haven't been able to find anything specific to this on here or through a Google search.

I could just make another column an only provide data for the specific values that I need, but that seem counterintuitive.

How do I label a single value and not show the other values?


texts: qml style differs from the same sld style saved in QGIS

texts:
I am trying to save style that I made in QGIS via layer properties in *.SLD format, but when opening my layer and import this saved *.SLD style - it differs from the initial one!

As you can see:

my style in QGIS:

my layer with saved *.SLD style:

All labels are disappeared and all contour lines become very thick! (width 0.1 becomes 1.0)

Why is that styles are different?


texts: Which geodatabase type is most appropriate for both open source and ESRI use?

texts:
Main Question:  Is there a geodatabase type that seamlessly works with both ESRI and open source products?

Details:  My company currently passes around geodata in various vector and raster file formats and isn't particularly organized.  Most processing is done using open source tools, but some customers demand that we host our data in ESRI Map service format.  We have an ArcGIS server instance already that hosts these services, but we still do our processing upstream using open source tools before dumping the results in an ESRI arcsde geodatabase.  From what I can tell, ogr/gdal have somewhat limited abilities to read/write from an arcsde geodatabase, but not enough that I could reasonably store all of our data in an arcsde geodatabase and expect our processing tools to be able to use the arcsde geodatabase as the primary data store for the company.  Am I wrong?  If I switched to an open source geodatabase such as postgis, our processing tools would work fine but I don't think that I could use a postgis database as the data store for hosting ESRI services.  Am I wrong?  Is there a geodatabase type out there that can service both open source and ESRI seamlessly?  If not is it common practice to maintain parallel geodatabases, one for open source and one for ESRI?

Additional wrinkle:
The answer I got was to use PostgreSQL/PostGIS .  The only potential wrinkle with that is that I'd like to host this database on Amazon Web Services (AWS) and if possible I'd like to use Amazon Relational Database Service (RDS) to minimize our database management overhead (backups, load balancing, etc).  However, RDS only currently supports MySQL, Oracle and SQL Server as database engines.  I'm thinking that if RDS starts supporting PostgreSQL that this would be the ideal scenario, but just to double check, PostgreSQL is my only option for ArcGIS + gdal/ogr-based open source interoperability correct?  If so I'll likely create a PostgreSQL database on AWS, I'll just have to spend more resources managing it.


texts: Viewing MapInfo .MRR rasters in ArcGIS Desktop?

texts:
I'd like to view MapInfo Multi Resolution Raster (.MRR) in ArcMap or ArcGIS Pro.

Esri Supported Rasters in ArcGIS Pro

I've found zero methods to view mrr. I'm not keen to export to a different format. I'm going to submit an enhancement request to Esri Inc.


texts: How to remove broken linked layers from ArcMap in ArcObjects?

texts:
In Python there is a function to list all broken data sources in a MapDocument:

brknList = arcpy.mapping.ListBrokenDataSources(mxd)

You could then use this function to remove such broken links:

arcpy.mapping.RemoveLayer(df, lyr)

How can I do a similar operation in ArcObjects 10.0?


texts: ArcGIS Network analyst - categories of facilities

texts:
Using the network analyst in ArcMap I am finding 5 closest facilities of each incident.

Is there any chance to work with categories of facilities and find 5 closest facilities that are not in the same category?

For example:  I am trying to find 5 closest pharmacies from a village that is close to a big city "Ahoj". All of the 5 closest pharmacies are in that city "Ahoj" but I want to find only 1 pharmacy from city "Ahoj", then 1 pharmacy from the second closest city, then 1 pharmacy from the third closest city etc.

I can add a printscreen if necessary.


texts: Displaying Google Maps and OpenStreetMap map in Java?

texts:
Is there any Java API that is like OpenLayers ?

I would like to display google Maps and Openlayers Maps in a java application


texts: Why is this GRASS script (from GRASS WiKi) causing an error?

texts:
I'm trying to run a script from the GRASS WiKi to create a simple map using ps.map. I'm confused because it's generating an error that the output was not set, when it appears to me that it is set at the beginning of the script. The script and error follow. The first code box contains the contents of the file spearfish.map, which is referenced as the input value in the second code box.

# Simple ps.map example using the Spearfish dataset
g.region rast=elevation.10m

ps.map output=spearfish.ps << EOF
paper us-letter
  end
scale 1:120000
raster elevation.10m
colortable
  where 2 6.75
  end
text 50% -7% Spearfish County, ND (Mount Rushmore)
  fontsize 16
  end
mapinfo
  where 4.5 7.25
  end
vlines roads
  where label ~ 'highway' OR label = 'interstate'
  color grey
  end
vpoints archsites
  symbol basic/triangle
  end
end
EOF

The command for running the script and the error it generates are below:

> ps.map input=spearfish.map
ERROR: Required parameter <output> not set:
    (PostScript output file)

UPDATE: Based on ShaunLangley's comment, I tried ps.map output=out.map. This seems to have started ps.map in interactive mode and I was able to enter the rest of the script one line at a time. That did create a map as a postscript file, so I'm headed in the right direction. I'll close this thread as solved and open a new one more targeted on the question of how to run the commands as a script.


texts:
Is it possible to find all Sentinel 2 tile ids by country? For example, for a country such as China, is it possible to get a list of all Sentinel 2 tiles id?

Things I've tried so far include looking at Sentinel 2 KML in Google Earth. There is also a matlab code snippet that takes input as lat,long and returns tile ids.

Matlab code snippet: https://www.mathworks.com/matlabcentral/fileexchange/61640-iskandervito-s2_tilefinder


texts: Is SpatiaLite the only single file/db spatially enabled exchange format?

texts: I am trying to figure out if there are other viable exchange formats for spatially enabled data. So far, it seems SpatiaLite is the only one out there, but it has yet to be adopted by the industry.


texts: gdal_polygonize error with Arc/Info binary grid raster file

texts:
I am trying to convert an Arc/Info binary grid raster file to a shapefile using the gdal_polygonize.py utility. But the operation fails with the following error

  Warning 1: Unsupported Arc/Info Binary Grid tile of type 0x80
  encountered. This and subsequent unsupported tile types set to no data
  value.

Here is a screenshot of the error

I am able to open and view the binary grid file just fine in ArcMAP etc. Here are the results when i run gdalinfo on the binary grid file. It seems to be fine

Driver: AIG/Arc/Info Binary Grid
Files: /home/vish/Sample_Data/nh_priority
       /home/vish/Sample_Data/nh_priority.ovr
       /home/vish/Sample_Data/nh_priority/vat.adf
       /home/vish/Sample_Data/nh_priority/w001001x.adf
       /home/vish/Sample_Data/nh_priority/w001001.adf
       /home/vish/Sample_Data/nh_priority/hdr.adf
       /home/vish/Sample_Data/nh_priority/dblbnd.adf
       /home/vish/Sample_Data/nh_priority/prj.adf
       /home/vish/Sample_Data/nh_priority/sta.adf
Size is 4987, 9645
Coordinate System is:
PROJCS["NAD83 / New Hampshire",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.01745329251994328,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4269"]],
    UNIT["METERS",1],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",42.5],
    PARAMETER["central_meridian",-71.66666666666667],
    PARAMETER["scale_factor",0.999966667],
    PARAMETER["false_easting",300000],
    PARAMETER["false_northing",0],
    AUTHORITY["EPSG","32110"],
    AXIS["X",EAST],
    AXIS["Y",NORTH]]
Origin = (228040.838957945175935,311318.187724410032388)
Pixel Size = (30.000000000000000,-30.000000000000000)
Corner Coordinates:
Upper Left  (  228040.839,  311318.188) ( 72d35'2.81"W, 45d17'53.85"N)
Lower Left  (  228040.839,   21968.188) ( 72d32'41.44"W, 42d41'39.84"N)
Upper Right (  377650.839,  311318.188) ( 70d40'35.98"W, 45d17'51.67"N)
Lower Right (  377650.839,   21968.188) ( 70d43'8.53"W, 42d41'37.85"N)
Center      (  302845.839,  166643.188) ( 71d37'52.26"W, 44d 0'0.06"N)
Band 1 Block=256x16 Type=Byte, ColorInterp=Palette
  Min=0.000 Max=2.000
  NoData Value=255
  Overviews: 2494x4823, 1247x2412, 624x1206, 312x603, 156x302
  Color Table (RGB with 3 entries)
    0: 255,255,255,255
    1: 211,255,190,255
    2: 0,84,3,255
<GDALRasterAttributeTable>
  <FieldDefn index="0">
    <Name>VALUE</Name>
    <Type>0</Type>
    <Usage>5</Usage>
  </FieldDefn>
  <FieldDefn index="1">
    <Name>COUNT</Name>
    <Type>0</Type>
    <Usage>1</Usage>
  </FieldDefn>
  <FieldDefn index="2">
    <Name>STEW_STATUS</Name>
    <Type>2</Type>
    <Usage>0</Usage>
  </FieldDefn>
  <Row index="0">
    <F>0</F>
    <F>7209043</F>
    <F>non-eligible</F>
  </Row>
  <Row index="1">
    <F>1</F>
    <F>7672024</F>
    <F>eligible, non-priority</F>
  </Row>
  <Row index="2">
    <F>2</F>
    <F>11779591</F>
    <F>priority</F>
  </Row>
</GDALRasterAttributeTable>


texts: Is there a way to distinguish between digitized and pasted features?

texts:
Is there a way to tell if an added feature in QGIS is in fact a newly digitized one or a pasted feature?

I think QGIS is treating both cases in the same way, emitting the same signal:

QgsVectorLayer::featureAdded(QgsFeatureId fid)

for a digitized feature as for a pasted one!

Maybe there's a way to tell the difference?


texts: Extract symbology labels to attribute table

texts:
In Arcmap 10.5, I have a feature layer with soil types where there has been made labels under symbology. These labels do not appear in the attribute table of the feature layer although they are based on one of the columns, the column "class", in the attribute table. E.g. a "class" value of "1" has been made into the symbology label "Sand".
How do I transfer these symbology labels to a column of its own in the attribute table?
I could of course make a statement in field calculator, but considering the number of levels of "class", that would be very time consuming.
I have tried to search for if the question was already posted, but the only thing I have been able to find was Export symbology to attribute table and there the question is not answered.


texts: Calculate float value in ModelBuilder and use as input for "Minus" tool

texts:
I'm using ModelBuilder and have three variables (var1, var2 and var3) which I want to use to calculate a new value. This value should be used as an input for the "Minus" tool. The calculation to be performed looks like this:

NewValue = var1 + (var2 / 2 / var3)

When calculating this manually for var1 = 2, var2 = 2 and var3 = 3 I get 2.333... (2 + (2 / 2 / 3))

When using the "Calculate Value" tool I get: 2

I was forced to use the data type "Formulated Raster" as output for "Calculate Value" because it seems to be impossible to use the output as an input for the "Minus" tool otherwise. Maybe this is the reason for the rounding.

So I tried using the "Raster Calculator" instead but this just leads me to an error message I don't understand:

How can I get my Model to calculate and use the correct value of 2.333...?

EDIT:

Changing the Rastercalculator Expression like suggested in the comments doesn't fix the problem:


texts: Is there a way within QGIS to hide/display (freeze/thaw) an individual node feature within a layer or does this control only exist at the layer level.


texts: ERROR - Loading Data using Geocode Addresses to ArcGIS

texts:
I am very new to ArcGIS and I am facing an error while using ArcMap 10.5 to load the below attached excel sheet through Geocode Address option.

Note that Multiple field mapping option and Mapping ADDRESS, CITY, STATE and FULLZIP column


texts:
I'm writing a script which takes, for its first two parameters, feature classes (or tables). For the third parameter, I want to provide the user with a field mapping option similar to what ArcMap uses for Append or Spatial Join, but I haven't been able to figure out how to add this functionality.

I found the data type "Field Mapping" in the tool properties, but I don't know how to get it to work. Help?

Edit: The original purpose of my tool is to copy the most recent inspection data from a related table into the parent table (the child table has multiple records for each record in the parent). I need the field mapping to help me know which fields from the child table to pull into equivalent fields in the parent table.

This is what I have:

# Based off code from: http://support.esri.com/technical-article/000011948
import arcpy

def updateParent(fc1, fc1_fields, fc2, fc2_fields):
    '''
    @param fc1 : feature class : The parent feature class, which we will update
    @param fc1_fields : list of field names : First name is date field, second is field to join, could modify code and add more fields to copy over
    @param fc2 : feature class or table : The child feature class, from which we will pull data
    @param fc2_fields : list of field names : See comments on second param, fields must be in the same order
    '''

    # We want to update values in our first feature class
    cursor = arcpy.da.UpdateCursor(fc1, fc1_fields)

    for row in cursor:
        # Create a new date list for each value in fc1
        date_list = []

        # Search all the values in fc2
        cursor2 = arcpy.da.SearchCursor(fc2, fc2_fields)
        for row2 in cursor2:
            # If the value matches the currently selected on in fc1, add date to list
            if row[1] == row2[1]:
                date_list.append(row2[0])
        # Get the latest (maximum) date
        print max(date_list)
        most_recent_date = max(date_list)

        # Update row in fc1 based on that value
        row[0] = most_recent_date
        cursor.updateRow(row)

    del row, cursor, row2, cursor2

# Define the workspace that will contain the two feature classes.
arcpy.env.workspace = r"C:\GIS_Temp\Temp.gdb"

#fc1 = arcpy.GetParameterAsText(0) # Parent feature class
#fc2 = arcpy.GetParameterAsText(1) # Child feature class / table

# This is where the problem is!!!!
# I don't know how to get a field map as an input and I don't know how to extract the information once I have it.
#fldmapping = arcpy.GetParameterAsText(2)

# If I hardcode my script, I just give it lists with the fields in the correct order so that I know DATE = DATE and Name = PostName
updateParent("TestRelateFC", ["DATE","Name"], "RelateTestTable2", ["DATE","PostName"])

Right now I'm only copying over DATE, but I would also like to be able to pull over the data that goes with that inspection date (like condition, etc.) The fields in the parent and child may not have the same name, which is why I need field mapping.

I assume I need to use updateParameters() in tool validation, but I don't know how it works for field mapping.


texts:
After making a check in from an arcpad AXF file, i notice that there's duplicate ObjectID records.

This ArcSDE PostgreSQL table is versioned and edited by Axf checkins and a .net engine application.

The weird behaviour is that the ObjectIDs from checkin start from the last ObjectID checked in from AXF in the last time, but between these times i added several records using the .net app so the ObjectIDs produced here are duplicated by the last check in, the .net app adds record directly through PostgreSQL, so i set the ObjectID  to (Max(ObjectID) + 1).

i'm using ArcSDE 10.1 and PostgreSQL 9.0.8.

i'm not sure but i suppose that the problem comes from generating ObjectID in .net app level.

How could I generate an ObjectID and inform ArcSDE that this ObjectID has been used ?


texts: OSM Python Utils: Modifying generate_tiles.py to create 512 x 512 tiles?

texts:
[Edited to include example image outputs]

SE folks, I'm wanting to modify generate_tiles.py to create 512 x 512 tiles.

At om_henners' recommendation, I tried replacing 12 of the 13 instances of 256 in the Python file (the 13th specifies the .PNG image format). I've since tried a test run, and I can already tell something is wacky...

Replacing all the relevant 256's with 512's merely produces a larger tile, when the intended result is increased mapped-area per tile (I think it amounts to 4x coverage).

Posting the code isn't worthwhile becuase you can see it here, but it is worth noting this is the only line for which I did not change 256 to 512 between test runs:

im.save(tile_uri, 'png256')    <-- the only 256 not changed to 512

So I'm getting these two outputs, and I have two curious observations:

1) The file address is the same. I would expect---if the mapped-area was changed between runs---the same location should fall into a different file, no?

2) Even though the mapped-area remained the same, the 512 tiles gained labels. Very interesting.

256x256 (%/8/61/97.png)

512x512 (%/8/61/97.png)

(And, I apologize for double-posting, but this question isn't getting any interest at help.OSM, and it seems at least one person at GIS.SE is okay with some of that traffic appearing here.)

Thanks in advance for the help!


texts: Data points collected by GPS on Arc 1960 UTM 37S do not correlate when projected on Google Earth

texts:
The following points were collected using a handheld GPS Garmin 76 set to Arc 1960 UTM UPS zone 37.

Beacon, easting, northing
FD8,274179,9135718
FD9,277586,9136635
FD10,282122,9137245

On entering the points into Google Earth using the easting and northing data, the points come out to approximately 300-350m north and 55-70m west of where they should be.

When I convert the points to WGS 84 before loading (and changing the setting on Google Earth), I have the same result.

Is this a projection issue?

Is it possible that the GPS were set on a different projection to what was recorded?

Here is a GE image of the differential: the green pins are what the gps data shows and the yellow pins are where they should be.


texts: Automatically Panning to Point Created with Absolute X,Y?

texts:
I'm currently working on a project where I have a many sets of XY coordinates along with corresponding marked aerial imagery (all on paper) from which I'm creating a feature class for a database.

After I enter each set of coordinates, I would like for ArcMap to automatically pan to the newly-created point, which would make it easier for me to reference it against the aerial imagery. (I'm currently right-clicking the feature on the Attribute Table and selecting "Pan.")

I feel like this should be easy/obvious, but I'm having no luck.

How would I go about doing this?


texts: How to do geographic operations (intersect) on OSM XML?

texts:
I would like to create bicycle router which routes more through shade. For that I would like to use OSM. Most routers use OSM xml so my question is:

How to add tags to a street that goes through a forest and write that back to OSM xml?
Getting part of a street through forest is easy with POSTGIS. Some intersect operation probably. I can insert OSM to Postgresql with osm2psql but I have found no way to get it back in usable xml again.

My current idea is parsing OSM XML and saving forest polygons somewhere. Then in another pass intersect each road with a forest polygon if bounding boxes intersects. If road is fully in forest tag it as in_forest. If not split road and tag forest part as in forest and other part as is.

I would use python libraries (Shapely) or C++ if Python isn't fast enough. But this would probably take ages. My country is only 350MB unzipped OSM XML.

Is there any better way? I looked into osmosis tag-transform but this is only for regex changes.


texts: Splitting two large TIFs covering the same area (or very similar areas) into smaller chips

texts:
I downloaded a bunch of images from Maxar's Open Data program covering pre- and post-images of recent disasters. See https://www.maxar.com/open-data
The pre-event and post-event images cover the same (or very similar) areas of interest. In other words, the intersection of a pair of pre- and post-event images is almost the same as the union of the pair, but not necessarily the same.
What I would like to do is to use QGIS to take a pair of pre- and post-images covering a given location and then crop the images into, say 512 by 512 pixels, chips such that in the intersection of the pair of images, a chip from the pre-image would have a matching chip covering exactly the same coordinates from the post-image.
Any ideas on how I can do that?


texts: Export Selected records from multiple feature classes to one feature class in ArcGIS Pro

texts:
I have used a select by polygon tool to select features in an area of interest from multiple feature classes.

How do I export the selection into one feature class in ArcGIS Pro?

At the minute the only way I can see is to export the feature from each layer seperately and then combine. I hope to do this in one command.


texts:
While working with the ModelBuilder, I added 2 columns to the data. After you have done "Summary Statistics" according to the GridID, you know the output product in tabular form. I wanted to print the table with the join based on the common column GridID's. But there is a problem here. During the Add Join process, the model wants the input product to be in Table Of Contents.
So it want to data as a "parameter". When I run the model from the inside, I say "Add to Display" to the latest product. Thus, the model sees it and accepts it as input. But when I run the model from the outside, it gives error because there is no input / output product in "Table of Contents". How can i solve it?


texts:
I have about 1.5 billion points in a geography column which represent topography.
I'll be retrieving all the points that are within a certain distance of a given long/lat. I was able to create a GIST index and query the table using ST_DWithin. These queries run in 40-50 ms.
The data in question is static. I believe a BRIN index would be quite a bit smaller and possibly faster but I've been having trouble getting it to run in anywhere near as as fast as the GIST index. I've let it run for about 5 minutes before canceling the query.
Here's how the table and index are setup. I got this method from the PostGIS docs:
CREATE TABLE topo
(
    location geography(Point,4326) NOT NULL,
    elevation smallint NOT NULL
)

#Create a temp index
CREATE INDEX idx_temp_topo ON topo2 USING btree(ST_GeoHash(location::geometry,4326));

#Order rows based on temp index
CLUSTER topo2 USING idx_temp_topo;

#Create BRIN index
CREATE INDEX topo_location_ix ON topo USING BRIN( location );

#Drop temp index
DROP INDEX idx_temp_topo;

#Vacuum
VACUUM ANALYZE topo2;

Query Plan
SELECT * FROM topo WHERE ST_DWithin(location, ST_MakePoint(-124.99875, 30.001250199900003), 100);

    #   Node
1.  Gather
2.  Bitmap Heap Scan on topo as topo
Filter: st_dwithin(location, '0101000020E610000052B81E85EB3F5FC019AFDFEE51003E40'::geography, '100'::double precision, true)
3.  Bitmap Index Scan using topo_location_ix
Index Cond: (location && _st_expand('0101000020E610000052B81E85EB3F5FC019AFDFEE51003E40'::geography, '100'::double precision))


texts: use Spatial filter tool running time take more

texts:
I want to check any 3 way junction or 4 way junction features snapping or not. I have used spatial filer it is taking more time for running could you please tell me any way to improve tool running speed.

    IWorkspaceFactory pfwf = new AccessWorkspaceFactoryClass();
    IWorkspace pwork = pfwf.OpenFromFile(PGDBpath, 0);
    IFeatureWorkspace pfwspace = pwork as IFeatureWorkspace;
    IFeatureClass cenline = pfwspace.OpenFeatureClass("Centerlines");
    IFeatureCursor cen_cur = cenline.Search(null, false);
    IFeature cen_F = cen_cur.NextFeature();
    int cnt = 0;
    int count = cenline.FeatureCount(null);
    while (cen_F != null)
    {
        cnt++;
        label3.Text = cen_F.Class.AliasName + " " + cnt.ToString() + "/" + count.ToString();
        Application.DoEvents();
        IPolyline pline = cen_F.Shape as IPolyline;
        IPointCollection ppcol = pline as IPointCollection;
        for (int i = 0; i < ppcol.PointCount; i++)
        {
            IPoint ppoint = ppcol.get_Point(i);
            ITopologicalOperator ptopo = ppoint as ITopologicalOperator;
            IGeometry pgeo = ptopo.Buffer(0.1);
            ISpatialFilter pslt = new SpatialFilterClass();
            pslt.SpatialRel = esriSpatialRelEnum.esriSpatialRelIntersects;
            pslt.Geometry = pgeo;
            IFeatureCursor pfcur = cenline.Search(pslt, false);
            IFeature cen_F_1 = pfcur.NextFeature();
            double dist = 0.0;
            bool chk = false;
            IProximityOperator prox = cen_F.Shape as IProximityOperator;
            while (cen_F_1 != null)
            {
                if (cen_F.OID != cen_F_1.OID)
                {
                    dist = prox.ReturnDistance(cen_F_1.Shape);
                    if (dist != 0.0)
                    {
                        chk = true;
                    }
                }
                cen_F_1 = pfcur.NextFeature();
            }
            if (chk == true)
            {
                wr.WriteLine(cen_F.OID + "," + cen_F.Class.AliasName + "," + dist + "," + "CenterLine Not snapped");
            }
            System.Runtime.InteropServices.Marshal.ReleaseComObject(pfcur);
        }
        cen_F = cen_cur.NextFeature();
    }
    System.Runtime.InteropServices.Marshal.ReleaseComObject(cen_cur);
    wr.Close();
    MessageBox.Show("Completed");
}
catch (Exception ex)
{
    MessageBox.Show(ex.Message);
}


texts: Web mapping software for a huge raster time series?

texts:
I'm an ArcGIS Desktop user who is stepping into the world of web mapping for the first time.  Following advice I found here on GIS.SE, I started with the free tutorials from OpenGeo.

However, I began to realize that most of the demos and examples were geared towards vector data.  My main project involves the display of a 300 frame time series of 5000 x 5000 pixel rasters.  On my desktop computer these are stored in a single 5000x5000x300 16-bit integer BIP binary, which totals around 8GB.  My goal is to be able to click on a cell (of a single raster in the time series) and have a graph pop up showing the values of that pixel in the 300-element time series.  The raster data used for each time series graph must be stored losslessly, though the overlay maps that are clicked on can be lossy caches.

Is there anything that may be better suited for this project (for a novice web GIS developer) than OpenGeo?  Or should I just keep going with this setup?

For reference, I have programming experience in Python, Java, and PHP.  I don't have much experience with SQL.  This is an open-ended project so I have plenty of time to learn new languages if need be.  I already have a web server, though I'll probably have to switch hosts because GoDaddy doesn't seem to support PostgreSQL without getting a VPS.

Thanks for your help!

EDIT: (Jan 13)  I'm still seeking information on exactly how best to store a 3 dimensional 16-bit integer BIP raster and be able to efficiently query a single z-axis "column" of data.  I don't want to convert it to a 32-bit format (because that would double its file size from its current 16-bit form).


texts:
Can python read the value of a cell in a spreadsheet?

The simplest example would be a script that ran a buffer tool on a given shapefile.

For the distance parameter, instead of just using a number like '1000', the script would point to a value in a cell of a spreadsheet (libre or open office preferred).

If there was then a way to trigger the script from the spreadsheet by way of a button, that would be the next step (then the next step would be to have a map control inside the spreadsheet to see the updated results!)

Just to give some insight into where I'm going with this: I'd like to use a spreadsheet as an analysis 'dashboard' where users could run analysis with different parameters - what would proximity around parks (grocery stores, etc.) be at 1/2 mi vs 1/4 mi...then another sheet in the spreadsheet would have a breakdown of the demographics within that proximity.

Thank you!


texts: Automatically Merge intersecting river Segments based on their strahler Classes

texts:
I have delineated a river network for a large area using a DEM. I have tried different approaches through SAGA GIS and/or ARC GIS.
The resulting file in all cases has a simple problem. The river network needs some merging in order to get acceptable results.
For example in the image the order 2 segment is divided into 3 different parts. I want it to be merged into 1 record in the attribute table.
In the database I have a column  "ORDER" stating the order of each segment. I need to find a way to merge the intersecting segments that have the same ORDER number. i.e. in the picture provided I need to merge the "ORDER 2" segments that intersect one to another but keep separate the order 1 segments. This has to be done for all the area (about 11.000 segments so it is very time consuming to do this manually).
Also I have 2 other columns: NODE A and NODE B, representing for each line segment where it starts and where it ends.

Left: What I have, Right : What I need to do.


texts: Add time stamp to an imageCollection of Google Earth Engine

texts:
On the Mapping over an ImageCollection documentation page, in order for adding time stamp the function is written as

    // This function adds a band representing the image timestamp.
var addTime = function(image) {
  return image.addBands(image.metadata('system:time_start'));
};

However, on the very next page of the documentation, Reducing an ImageCollection, the function is written as:

// This function adds a band representing the image timestamp.
var addTime = function(image) {
  return image.addBands(image.metadata('system:time_start')
    .divide(1000 * 60 * 60 * 24 * 365));
};

Question:

So, why the division does not exist in the first function? What is the difference?


texts: Not able to find out if polygons in one shapefile intersect with polygons in other shapefile using GeoPandas

texts:
Santa Barbara contains some places from CA_Places. But when I use intersects, I get noting but all False. I have spent hours on this. Both shapefiles are here for reference.
CA_Counties = geopandas.read_file("CA_Counties\\CA_Counties_TIGER2016.shp")
CA_places = geopandas.read_file("CA_Places\\CA_Places_TIGER2016.shp")
Santa_Barbara = CA_Counties[CA_Counties['NAME'] == 'Santa Barbara']
CA_places.intersects(Santa_Barbara)

All I'm getting is this warning:
C:\Users\LENOVO\.conda\envs\gdal_home\lib\site-packages\geopandas\base.py:39: UserWarning: The indices of the two GeoSeries are different.
warn("The indices of the two GeoSeries are different.")


texts: Converting Surface Model to Elevation Model?

texts:
Is there a way to convert a digital surface model (DSM) to a digital elevation model (DEM) in QGIS?

Surface models contain lots of noise (e.g. trees, buildings) when I only want to show the ground. I currently have to manually delete individual objects using a masking layer, replacing null values with GDAL Fill NoData and then smooth out with a gaussian blur. This is very laborious.


texts: GRASS r.sun returns "ERROR: insol_time and incidout are incompatible options"

texts:
I'm trying to measure irradiance on a DEM using the r.sun tool in QGIS 3.6. I have grids for the elevation (EPSG:32735 - WGS 84 / UTM zone 35S), the slope and aspect (EPSG:4326 - WGS 84). I've tried to ensure they are in the proper coordinate systems asked for here.

However, every time I run it, I get the following list of red errors in the log:

*ERROR: insol_time and incidout are incompatible options

ERROR: Raster map incidoutce7a2d761eb742ac87d7519caff9005a not found

ERROR: Raster map or group incidoutce7a2d761eb742ac87d7519caff9005a not found

ERROR: Raster map beam_radce7a2d761eb742ac87d7519caff9005a not found

ERROR: Raster map or group beam_radce7a2d761eb742ac87d7519caff9005a not found

ERROR: Raster map insol_timece7a2d761eb742ac87d7519caff9005a not found

ERROR: Raster map or group insol_timece7a2d761eb742ac87d7519caff9005a not found

ERROR: Raster map diff_radce7a2d761eb742ac87d7519caff9005a not found

ERROR: Raster map or group diff_radce7a2d761eb742ac87d7519caff9005a not found

ERROR: Raster map refl_radce7a2d761eb742ac87d7519caff9005a not found

ERROR: Raster map or group refl_radce7a2d761eb742ac87d7519caff9005a not found

ERROR: Raster map glob_radce7a2d761eb742ac87d7519caff9005a not found

ERROR: Raster map or group glob_radce7a2d761eb742ac87d7519caff9005a not found*

Then no output map is generated. I've tried running it in an older version of QGIS (2.18), but the whole program just crashes.

It looks like someone has reported a bug about the first line of the error here. But wondering if there's anything I'm doing wrong as I'm new to using GRASS.


texts: Satellite image of agriculture land cover for some country

texts:
I am looking for a free satellite image of agriculture land cover for North Korea. Can someone guide me to any free resource.

Is it available on earth Explorer, if yes, which image should I download.


texts: How to interpolate missing data from neighbours in PostGIS?

texts:
I have a dataset with point data within a given country. Let's say my dataset looks somewhat like this:

tree_id | species | age | geom
------------------------------
   0    |   Ash   | null| ...
   1    |  Beech  |  70 | ...
   2    |   Ash   |  10 | ...
   3    |  Beech  |  70 | ...
   4    |  Beech  | null| ...
   5    |  Beech  |  60 | ...
  ...   |   ...   | ... | ...

As you can see the dataset has some missing data. For instance, tree_id 0 has no age. Therefore I would like to interpolate those missing values from a 100 meter radius.

I am looking for the mean of the species. The result should also include the number of sample trees used. A result table could then look like this:

tree_id | age_avg | samples
---------------------------
   0    |   11.8  |  113
   3    |   12.2  |   97
   5    |   50.7  |  272
  ...   |   ...   |  ...

Could you get me started with some PostgreSQL query code, please?


texts: Is Enterprise GIS license required to use Production Mapping license?

texts: I'm looking to build up a hybrid GIS Enterprise environment whereby multiple users can access data at the same time via Arcmap. If possible, we'd like to not have to spend a fortune doing this. My understanding is that I can edit data in GeoServer from ArcMap (with a Production Mapping license) via WFS-T. But through further research, it appears that I cannot get a production mapping license without first purchasing Enterprise. Is this correct?


texts:
I try to represent point layer with different size and i want to have a color ramp also:

The problem is that when i  choose the color method, the size method disappear:

I can change the color manually, but i seek for automatic way.

I also have red this How to deactivate the Color Ramp in QGIS graduated vector styling? and didn't got my wanted answer.


texts: Which units are used for the Print Composer grid?

texts: I changed the units on the general tab to inches. When I go to set up the grid under snapping, whats the units. From what I have read it seems as though they are map units, it that correct? I am a attempting to get set up a sheet to fit within en existing title block.


texts: How to show QgsDebugMsg in eclipse PyDev

texts: In the source code of QGIS, I see a lot of  QgsDebugMsg line. What should I do in order to show it when debug/run my standalone application in Eclipse PyDev?


texts: Atmospheric correction in ENVI, export to ERDAS

texts:
I am working on a project using Hyperspectral imagery and have ERDAS Imagine on my home PC and have access to ENVI on another PC at University. Is it possible, for example, to perform the atmospheric correction in ENVI, then import it into ERDAS, or vice versa, and perform Principal Component Analysis, Minimum Noise Transformation (MNF), supervised classification etc?

The file is a .l1r file from EO-1 Hyperion downloaded from GloVis.

I read here (Imagine files between ERDAS and ENVI) that you can export the file as a .bil and then import into ERDAS, so I presume it is possible?

I don’t have ENVI on this PC to double check at the moment.


texts: Using rasterio to show RGB image with nodata

texts:
This question is same as https://stackoverflow.com/questions/59132459/use-rasterio-show-rgb-image-with-nodata

I have a multiband image with nodata value (which is Landsat and been set nodata to 0). And I want to use rasterio show the image,But I found that the show function cant deal the nodata value with multiband image. And I had test that singleband image is well.My code like this:

with rio.open(path) as src:
    arr = src.read([1, 2, 3], masked=True)
    show(arr)

the image is like this:

This is the landsat image,You can find that the image is warped by black color,In this place the data is nodata and would be display as transparent.

How do I show nodata multiband img by rasterio?


texts:
Using Fiona for reading my VRT file I'm receiving an error that OGR_VRT driver is not supported:

fiona.errors.DriverError: unsupported driver: 'OGR_VRT'

My test.vrt file looks like:

<?xml version="1.0" encoding="UTF-8"?>
<OGRVRTDataSource>
    <OGRVRTLayer name="test">
        <SrcDataSource>absolute_path_to_test_shp.shp</SrcDataSource>
    </OGRVRTLayer>
</OGRVRTDataSource>

I'm trying to read it using:

>>> import fiona
>>> c = fiona.open('abs_path/test.vrt')

How can Ogr VRT be read by Fiona?


texts: pgr_dijkstra is considering one-way streets as bidirectional

texts:
I imported data from a map extracted from OpenStreetMap. With osm2pgrouting I imported it into PostgreSQL, to be worked with pgRouting and PostGIS. Command used:
osm2pgrouting -f DIRECTORY/map.osm -d DATABASE -p 5432 -U postgres -W PASSWD-c mapconfig.xml --addnodes --attributes -- clean

He created the following table structure:

I am using the pgr_dijkstra function to calculate a minimum cost route between two locations:
SELECT * FROM pgr_dijkstra('
SELECT gid as id, source, target, ST_Length(the_geom::GEOGRAPHY) AS cost,
ST_LENGTH(the_geom::GEOGRAPHY) AS reverse_cost
FROM ways',
24212, 5521, directed := true
) as di
  JOIN ways
ON di.edge = ways.gid
order by seq;

It is working perfectly when I use directed: = false. But when I put directed: = true pgRouting is returning the same result. This is wrong because even with directed: = true he is considering my unidirectional streets as bidirectional. The image below represents some return lines from my query:

oneway = YES and one_way = 1 is a one-way street that the function should not be returning. oneway = UNKNOWN and one_way = 0 is my two-way street. The image below shows this one-way street that is being considered as a route path, which is wrong, is in the opposite direction of traffic. In the drawing it can be a bit confusing, but the street "Florestano de Macedo" is unidirectional, as shown by the return in the database. The figure is a bit confused because it is a street that passes next to a highway, and on the other side of the highway there is another street. But this problem occurs for any route with directed: = true.

What am I doing wrong? Do I have to prepare/change my data in any way for these cases? Is my query wrong? PgRouting workshops only say to use directed as true or false, they do not mention data preparation for this.


texts: Geoserver SLD labels dynamic displacement to resolve overlap conflicts?

texts:
Is it possible to style labels in such a way that if two labels were to overlap each other, they could move out of each other's way and display next to each other instead of their original position?

For example, if by default a label is displayed in the center of a polygon, and two polygons have the same geometric center, only one of the two labels for those polygons would show (as they would have the same priority, one of the two labels is chosen at random). I'd like to style the labels in such a way that if this were to happen, both of the labels will display but they'd be displaced or offset slightly from the center so that they can both be displayed at the same time, without overlapping.

Is this possible? Relevant documentation:
http://docs.geoserver.org/latest/en/user/styling/sld/reference/labeling#pointplacement


texts: Python to generate a different number of random points for features using iteration with tool parameter values in feature attribute table

texts:
I want to generate a certain number of random points in N number of zip codes. The number of points to be generated for each individual zip code is stored the attribute of the feature class of all state zip codes, field type = Long. There is also a field with a unique ID for each zip code. This unique ID is a 5 character string.

This is the python script I wrote to accomplish the task. The code throws multiple errors.

import arcpy # Import ArcPy site-package and os modules

fc = r'C:\ArcData\2018Cases\ZCTA_Cases.shp' # Path to input feature class

fields = ['id2','Num'] # Fields to be searched. id2 is zip code. Num is the number of points to be generated in the zip code

p = r'C:\ArcData\2018Cases\\ZCTAScript_Results' # Output Path

scur = arcpy.da.SearchCursor(fc,fields) # Create a searchcursor object, which consists of rows of data. Each row will be a tuple of (id2,Numcase)

for row in scur: # Iterate scur
    query = ' "id2" = row[0]' # Query used to select current zip
    selection = arcpy.SelectLayerByAttribute_management(fc, 'NEW_SELECTION', query) # Selection of current zip
    arcpy.CreateRandomPoints_management(p, row[0], selection,'', row[1]) #Create random points for each
    # row in scur using id2 as output name and Num as number of points

del scur # Delete the cursor to release schema locks

This code is based off of code posted by a previous user in response to a previous question that I asked on this issue that was placed on HOLD. The code the user kindly provided did not throw errors but it did not produce the desired result; it generated random points in the full feature class (all zip codes in the state) extent not in the extent of a given individual zip code. Here is that code:

import arcpy

fc = r'path/to/feature/class' #path to input feature class
fields = ["ID","PointNum"] #the fields you want to search
p = r'output/path' #output folder or gdb

scur = arcpy.da.SearchCursor(fc,fields) #create a searchcursor object,
#which consists of rows of data.  each row will be a tuple of (ID,PointNum).  there will be N rows.

for row in scur: #iterate scur
    arcpy.CreateRandomPoints_management(p, row[0], fc, '#', row[1]) #create random points for each
    #row in scur using ID as output name and PointNum as number of points

del scur #delete the cursor to release schema locks

My code tries to resolve the issue by selecting the individual zip code that is being read by the search cursor and setting the constraining extent of the random point tool to that zip code. I cannot figure out how to write the code to accomplish the desired task.


texts: Accessing arcgis.com webmap layer using ArcMap

texts:
Background:
Our local newspaper contained this link to some upcoming forest management activity in the area:
https://usfs.maps.arcgis.com/apps/webappviewer/index.html?id=021ff3371f584a54beb15cd481a8fbb5
Here's a screenshot of the resulting argis.com webmap:

Problem:
I would like access the prescribed fire polygon layer displayed in this webmap, either as a service or via download.  However, I'm unable to accomplish this seemingly simple task.  Here's what I've unsuccessfully tried:

Look for some sort of link (service url or layer download) in the webmap itself
Search for the polygon layer via my personal arcgis.com account.  Here my steps were to create a new map and search using Add > Search for Layers, Add > Browse Living Atlas Layers, and Add > Add Layer From Web.
From ArcMap (10.8.1), Add > Add Data from ArcGIS Online and search

Question:
How do I access the polygon layer in this webmap?


texts: Embed Static Carto Map in Google Sheet Cell

texts:
Carto is mapping software that provides me with a url to a dynamic map.

I want to show a static version of the map (which updates daily) in a google sheets cell.

Ideally there would be a way to manipulate the carto url to show as an image and just use the IMAGE function in google sheets.

It doesn't seem possible, so I may have to use google app scripts, I've never used before and have no experience with html, so would be good to get some specific advice.


texts: Line lengths calculation in combination with DEM

texts:
I am working in QGIS 2.18.0.

I have a problem regarding calculation of line lengths in combination with raster files, to be precise with DEM raster.

There is not a problem to calculate the line lengths (vector layer), but this calculation assumed the surface to be flat (or not?). By doing this that way I will calculate the length of the "cathetus", and I need the length of the "hypotenuse".

Is there any add-on or is there any formula to calculate the line length that incorporates the change in the altitude?


texts: How can I convert Horizontal dilution of position to a radius of 68% confidence?

texts:
I am working at mobile app that uses external GPS module. What I need is to calculate accuracy for the fix from GPGGA command to be able to compare it with my mobile platform's GPS accuracy. Mobile platform uses the radius of 68% confidence measured in meters as an accuracy.

My question is: what is  Horizontal dilution of position in GPGGA command? How can I convert it to the radius of 68% confidence? I used this description for my parser implementation.


texts: How to convert a coordinates to indice array (vice versa) in a NetworkX graph?

texts:
I want to use a raster for a A* and bidirectional Dijkstra path analysis in NetworkX. I am using Python for this project.

Raster example (it's a png file converted when uploaded, but the real problem is TIFF):
Here

First I read in the raster with GDAL

input_raster = "raster.tif"
raster = gdal.Open(input_raster)

Next I read the raster as an array

bandraster = raster.GetRasterBand(1)
arr = bandraster.ReadAsArray()

So, I'll transform coords using a function:

def coord2pixelOffset(rasterfn, x, y):
    raster = gdal.Open(rasterfn)
    geotransform = raster.GetGeoTransform()
    originX = geotransform[0]
    originY = geotransform[3]
    pixelWidth = geotransform[1]
    pixelHeight = geotransform[5]
    xOffset = int((x - originX)/pixelWidth)
    yOffset = int((y - originY)/pixelHeight)
    return xOffset, yOffset

CostSurfacefn = 'raster.tif'
source_coord = (-41.1823753163, -13.83393276)
target_coord = (-40.3726182077, -14.2361991946)

# coordinates to array index
source = coord2pixelOffset(CostSurfacefn, source_coord[0], source_coord[1])
target = coord2pixelOffset(CostSurfacefn, target_coord[0], target_coord[1])

Next, the graph is loaded and function is called:

G = nx.DiGraph(arr)
length, path = nx.bidirectional_dijkstra(G, source, target)

But, the function return a error:

raise NetworkXError("The node %s is not in the digraph."%(n,))

  networkx.exception.NetworkXError: The node (2612, -352) is not in the digraph.

How to do this in correct way?


texts: Seeking shapefiles for Indian Constituencies?

texts:
The GIS data officially provided by the government is only available as PDF maps. Some people have shared shapefiles for the previous delimitation but these are outdated. Parliamentary constituency data can be also be viewed here but no source has been provided.

Would anyone know a way to get GIS data for the assembly and parliamentary constituencies?


texts: GEOTIFF explanation of .aux xml and tfw world file

texts:
I want to know how raster images work,

I have a set of geotiffs that has no aux, tfw or xml

when I check the spatial information in arccatalog, they do have an projection. and they shows correctly.

if the projection is embedded in the geotiff, why do I need .tfw ?

I need to know what it does, is it used for quick transformation if the geotiff is in real coordinates wgs84 and I want another projection for fast show? or does my geotiffs lack some information?

and when I define another projection, arccatalog creates aux, do i need to keep this aux?

so. shortening the question, what is aux, tfw good for if the information is already stored inside the geotiff?


texts: Using custom contextMenu in Arcobjects(C# and ArcMap)

texts:
I added a new custom tree to TOC which looks like this:

IObjectClass
    IFeature
    IFeature ...
IObjectClass
    IRow ...

made context menu using this sample http://help.arcgis.com/en/sdk/10.0/arcobjects_net/conceptualhelp/index.html#//00010000026q000000 and it shows up alright, but unfortunatelly its context menu commands are not doing anything. Here is onNodeMouseClick:

private void tvwSelected_NodeMouseClick(object sender, TreeNodeMouseClickEventArgs e)
{
    if (e.Button == MouseButtons.Right)
    {
        m_contextItem = e.Node.Tag;

        UID menuID = new UIDClass();
        if (e.Node.Tag is IFeature)
            menuID.Value = "{d32ae970-abb8-4f02-a592-9181f1307ed9}";
        else
            return;
        ICommandBar cmdBar = (ICommandBar)_application.Document.CommandBars.Find(menuID, false, false);
        cmdBar.Popup(0, 0);
    }
}

and here is context menu:

public ContextMenuTest()
{
    AddItem("esriCore.FeatureFindContextMenuItems", 1); //flash
    AddItem("esriCore.FeatureFindContextMenuItems", 2); //zoom
    AddItem("esriCore.FeatureFindContextMenuItems", 3); //identify
    AddItem("esriCore.FeatureFindContextMenuItems", 6); //unselect
}

I believe the probblem could be that IFeature is not what these commands require, or the whole thing works in totaly different way. Any help would be appreciated.


texts: What is the most efficient way to check if coordinate didn't change for 10 minutes (c# or mssql) from a list of coords

texts:
I have a list of GPS coordinates (lat,lon) with a timestamp that I got from GPS logger.
I want to find the first points that didn't change in a buffer of 100 meter for 10 minutes.
What will be the most efficient way to do this (I can do it in c# or in mssql) ?

for example , I drove my bike 20 k"m , stopped for half an hour (logger keep recording data every 10 seconds) and than continue to drive. I want to identify the first time I stayed in place in area of 100 meter (because logger isn't accurate) for atleast 10 minutes.