texts: How to extract geo data out of a polygon shapefile into CSV KWT

texts:
I have a polygon shapefile
that contains all the neighbourhoods in the Netherlands.

My aim is to export the neighbourhoods and coordinates to CSV file (KWT).

Without success I tried (with QGIS) these 'solutions':

How to export Polygons to CSV with coordinates?
How to get Shapefile point coordinates into .xls?

It got me to the point where I have all the attributes (from my vector layer) in a CSV file. The problem lies in that I can't seem to find any geodata (coordinates) in them.

Where can I find the polygon coordinates to finally export them with the neighbourhoods out my shapefile to CSV?


texts: Extract and Add raster class attributes to vector file attributes in ArcGIS

texts:
I have a shapefile having (say 20 polygons). There are 20 rows in the attribute table of this shapefile (1 corresponding to each polygon). I have a binary raster with two classes 1 and 2. I want to calculate and extract the area corresponding to class 2 of raster in each of these polygons.
How can this area be added to attributes of shapefile and show the area in 20 rows of 20 polygons? How to do this with minimal manual effort in ArcGIS?


texts: I want to transform WGS to MRT (GCS for Malaysia) using QGIS. Is "save as" to a new CRS the correct method? With ArcGIS, I can create transformation method to do transformation but I'm not sure in QGIS. If it is possible, how? And what is the string for xyz rotation, xyz translation, scale factor?


texts: Why does Arcpy CalculateField alter string being written to string field?

texts:
I have a "Model_Version" field in my feature class that is being used for record keeping.  I do not intend to ever do any numeric operations with it, so I am storing it as a string field.  Right now my model version is 1.0, so it is being stored as "1.0".

When I used the Field Calculator manually in ArcMap, I can get my desired value by entering "1.0" or '1.0'.

But when I use ArcPy as below, the resulting string in the feature class is "1", without the decimal and trailing zero:

arcpy.CalculateField_management(myLayer, "Model_Version",
                                "1.0", "PYTHON3")

I have tried the following and get the same result of "1" every time:

arcpy.CalculateField_management(myLayer, "Model_Version",
                                '1.0', "PYTHON3")

arcpy.CalculateField_management(myLayer, "Model_Version",
                                str("1.0"), "PYTHON3")

arcpy.CalculateField_management(myLayer, "Model_Version",
                                1.0, "PYTHON3")

Why is it not writing the string as I have coded it?


texts: In DIVA GIS how can I change the polygon co-ordinates of a .shp file?

texts: I am using polygon layers from various sources to create species distribution maps however some of the layers have different polygon coordinates to others, therefore I cannot use them together as layers.  Is there a simple way to change the coordinates of the polygon?  I am relatively new to DIVA gis so any help is much welcomed.


texts: Round Raster values up to nearest 10 QGIS

texts:
I have a raster with cell data varying between 0.1 and 156. I want to round the data to the nearest 10, i.e any value between 0.1-10 = 10, 10.1-20 = 20 and so on. Theres lots online about rounding to the nearest decimal, but not to the nearest 10...

I have tried using the raster calculator, GDAL raster calculator and endless other things, does anyone know how to do this?


texts: REST Queries for ESRI Map Service Raster Layers

texts:
I'm having trouble getting REST queries to work for any ESRI Map Service raster layers. Right now I've just been working with the query generator as part of the ESRI Service Directory. Here's are two Map Server examples:
https://apps.fs.usda.gov/arcx/rest/services/RDW_Wildfire/MTBS_CONUS/MapServer/5/query
https://basemap.nationalmap.gov/arcgis/rest/services/USGSImageryTopo/MapServer/0/query
No matter what my query looks like I get "Invalid or missing input parameters." If this were a vector Map Service I'd put 1=1 in the WHERE clause and the service would return all the features contained in the layer. For raster layers in a Map Service I've tried every combination of inputs I can think of and always get the same "Invalid or missing input parameters.". If I specify JSON as the output format it shows up with error code 400 but I can't find anything on this either.
Edit
To be more specific I'm trying to figure out what is wrong with the REST query I'm trying to use, I know there are numerous other ways to access the information I'm looking for but to integrate this information with the rest of my workflow the data needs to come via an API call. I've tried removing fields and adding text or specifying a geometry, requesting only a count or IDS and I can't get anything to work. Here's an example of what one of my current queries looks like:
https://apps.fs.usda.gov/arcx/rest/services/RDW_Wildfire/MTBS_CONUS/MapServer/5/query?where=1%3D1&text=&objectIds=&time=&geometry=&geometryType=esriGeometryEnvelope&inSR=&spatialRel=esriSpatialRelIntersects&relationParam=&outFields=&returnGeometry=true&returnTrueCurves=false&maxAllowableOffset=&geometryPrecision=&outSR=&having=&returnIdsOnly=false&returnCountOnly=false&orderByFields=&groupByFieldsForStatistics=&outStatistics=&returnZ=false&returnM=false&gdbVersion=&historicMoment=&returnDistinctValues=false&resultOffset=&resultRecordCount=&queryByDistance=&returnExtentOnly=false&datumTransformation=&parameterValues=&rangeValues=&quantizationParameters=&featureEncoding=esriDefault&f=json


texts:
I am trying to put a map into the center of my ArcMap's layout view. I have checked various options (Full Extent; Zoom in and out; File --> Page Print Setup) but neither of them solved the problem.

Currently, my map looks like this:

And a while ago I managed to produce this map. I remember that it involved selecting a pre-defined "A4" sized page in one of the menus. Unfortunately, I am unable to reproduce how I got there.

Now I would like to center the first map in the same way as the second map. I am aware of at least a subset of similar questions that were asked on this site. All solutions that I found there, however, involved manually shifting the map around. Ideally, I would like to reproduce exactly what I did such that the map extent of both maps is the same.
Edit: See comment to answer below


texts: Why do my attribute table headings keep changing?

texts:
I am working on a project which previously had no issues. I recently went in to create another layer and when I created the vector layer the headings that I had given to the entries all changed to things like "Action" "_1" "_2" etc. No matter how many times I tried to add in new columns they would continue to change similar to this pattern. In addition, my previously created layers now seem to be messed up as well. All the data entered are changed to symbols rather than numbers and text.

Any suggestions on what could be the cause or the solution?


texts: Spatial join with field mappings yields zero value when inputs are layers made from feature classes

texts:
this might be more of a bug report than a question, but I encountered this strange behavior in ArcPy and was wondering if anyone else has experienced the same issue.
Summary: I was trying to perform a spatial join between two layers that are created with arcpy.MakeFeatureLayer_management (in order to select a subset of the features for the join). Then I defined some field mappings with fieldmap.addInputField(layer_name, field_name). Finally, I did a spatial join based on the two layers, but the resulting fields defined with the field mappings are all zero values. My workaround for now is to save out the selected features from the layers, and use these feature classes for the spatial join.
Script:
arcpy.env.workspace = 'C:/example_folder'
arcpy.MakeFeatureLayer_management('input_layers.gdb/base_layer', 'base')
arcpy.MakeFeatureLayer_management('input_layers.gdb/Points_of_Interest', 'POI')

arcpy.SelectLayerByAttribute_management('POI', 'NEW_SELECTION', where_clause = '"visits">10')

fm = arcpy.FieldMappings()

nf1 = arcpy.FieldMap()
nf1.mergeRule = 'First'
nf1.addInputField('base','UID')
fm.addFieldMap(nf1)

nf2 = arcpy.FieldMap()
nf2.mergeRule = 'Sum'
nf2.addInputField('POI','visits')
fm.addFieldMap(nf2)

arcpy.SpatialJoin_analysis('base','POI',
                'Output_layers.gdb/base_POI_sjoin',
                'JOIN_ONE_TO_ONE','KEEP_ALL',
                field_mapping = fm,
                match_option = 'WITHIN_A_DISTANCE',
                search_radius = '2 miles')

Result: Both the UID and visits fields are all zero values in the output layer.
Any thoughts as to why this might be happening?


texts: Points on the opposing side of a breakline

texts: If you have a breakline or even a single line segment within a DEM, how can you tell if two points lie on the opposing side of that breakline or line segment? From a computational geometry standpoint is the cross product of two lines the only alternative? Additionally, how would you scale this issue when you have thousands/millions of points and breaklines/segments? I am trying to preserve breaklines by restricting smoothing across them.


texts: LiDAR DSM: difference between Return 1 and First of Many?

texts: I'm making a DSM for an urban LiDAR dataset; I understand these typically are also referred to as "first return" rasters- I don't quite understand the difference though between Return 1 and First of Many, both of which are present in the dataset. Which should be used for a DSM?


texts: "Permission denied" error when trying to create shapefile layer

texts:
I'm not super familiar with QGIS as I am used to ArcGIS. But when I try to create a new shapefile layer, I keep getting an error that says:

Layer creation failed: Failed to create file vacancies.shp: Permission
denied.

Right now I am working off an OSM Standard layer and a shapefile I downloaded from the ArcGIS website. What might be wrong?


texts: Appropriate Interpolation accuracy and method for a particular project

texts: I'm trying to find the best interpolation method for  a project I'm working on.  I've tried the basic settings with spline, IDW & kriging and they all appear ok at first glance, but when I look more closely they don't follow the points exactly like they should.  For example, I will have a point or cluster of points with a value of 59, and a neighboring point or cluster of points with a value of 60.  When I adjust the symbology to 59.5 I would think the line would be drawn evenly between the two points or clusters of points every time.  In some cases the line is drawn between the two, but in other cases it is not.  The resulting map is "generally" correct, but with a lot of little errors here and there.  I'm not real familiar with how interpolation works exactly, and I want to know if one of the methods would be more accurate or have less errors for the particular project I'm working on. I will have 20,000+ points over an area the size of Florida (if that even matters).


texts: Edit large geojson data sets in Openlayers

texts:
I'd like to edit a rather large collection (several MB of polygons and lines) with OpenLayers.

What's the best way working with large data sets, considering client side (rendering) performance and server side scalability?

We're currently using Geoserver on the backend side, though in order to enable scalable editing capabilities we're also looking into other directions.

EDIT

I should add a bit more background. I expect lots of geometries, even when zoomed in. Consider very large venues where there are hundreds of stands / booths, besides also other indoor (floors and others) and outdoor (building) geometries.


texts:
Developing the issue solved in this post:

Trying to pop up image in GeoJSON file?

I would like to have each GeoJSON image thumbnail possible to open in larger window.

I was trying to connect with my local file URL (but it can be obviously some link in the web):

{
  "type": "Feature",
  "properties": {
"Title": "JHG",
    "Head": "Shortstanton sidings",
"Description": "Conor Murphy",
"URL": "JHG",
"Pict": "image.png"
"PicURL": "image.png
},
  "geometry": {
    "type": "Point",
    "coordinates": [
      0.05458831787109375,
      52.29163006501503
    ]
  }

}
                  ]
            };

Then I built up the major code like this:

var sitis = L.geoJSON(sitec, {
pointToLayer: function (feature, latlng) {
feature.properties.myKey = feature.properties.Title + ', ' +
feature.properties.Head
    return L.circleMarker(latlng, geojsonMarkerOptions);
},
    onEachFeature: function (feature, layer) {
        layer.bindPopup("
"+feature.properties.Title+"

Address:
"+feature.properties.Head+"

"+feature.properties.Description+"

Website:"+feature.properties.URL+"

");
}

    }).addTo(map);

And I got a blank screen unfortunately.

Besides the following code works in this case:

var sitis = L.geoJSON(sitec, {
pointToLayer: function (feature, latlng) {
feature.properties.myKey = feature.properties.Title + ', ' +
feature.properties.Head
    return L.circleMarker(latlng, geojsonMarkerOptions);
},
    onEachFeature: function (feature, layer) {
        layer.bindPopup("
"+feature.properties.Title+"

Address:
"+feature.properties.Head+"

"+feature.properties.Description+"

Website:"+feature.properties.URL+"

");

}

    }).addTo(map);

But I got the same link for all GeoJSON placemarks.

Does anyone know how to make the link unique for a single one popup?


texts: How to extract long and lat of start and end points, to seperate columns, from the geometry column in GeoPandas?

texts:
I have a line shapefile in GeoPandas. How can I extract the start and end coordinates as (x,y) pairs in two separate columns?

import geopandas as gpd
gdf = gpd.read_file(pth_line_shapefile)
gdf['start_coordinates_XY'] = gdf.apply(lambda x: [y[0] for y in x['geometry'].coords], axis=1)

The above code extracts X coordinates of all points in the geometry not the first point (in x,y format)

An example line shapefile linked here


texts: Generating composite values ( MAPFID numbers) for previously duplicated vales based on 2 different fields

texts:
Good day everyone!

I have lurked around here for info preiously, but this is my first post, so please be gentle if I commit any faux pas.  Okay, let me try to break this down so we can all understand. (myself included)

I am attempting to edit a file GDB that contains over 1.5 million forest stand polygon records and change those that have a specific duplicated value.  I will most likely be obtaining the information on how to do that from Generating sequential numbering for duplicate values in a field sorted by ObjectID, and since the data is stored on a file GDB, I will also be sourcing some code from Python script for identifying duplicate records (follow up), however that is the part that I have mostly figured out already.

My problem right now is that I am not familiar enough with python terminology to even properly explain the next step so I am going show you my thought process on this as best as I can. Warning: this explanation is not going to be in proper syntax simply because I am not very good at writing it quickly yet.

This is the part I think I understand what to do

IF Field A(called'NP')>0 AND Field B(called 'PREFID')contains a duplicate value THEN generate unique 'PREFID' value counting from +9999

and this is where it gets hinky for me:

Where a new value is created based on this criteria, I need another existing field (called 'MAPFID') to have a new value generated that will be a composite of the values in that row for the fields 'MAP' and the newly created 'PREFID' number. This field is called 'MAPFID'.  Eg: if MAP=1343590 and PREFID=10000 then the value in MAPFID=134359010000

Because the 'MAP' value is based on UTM coordinates, the values for 'MAP' will not necessarily be sequential.

I realize that this may be a green thing to ask about, but this has been stumping me for too long and I can't seem to find anything online about this.  I'm probably just not using the right vocabulary for my searches, but that just goes to show ya how green I am :P
Any help would be greatly appreciated! Thanks guys and gals.


texts:
Thanks in advance for helping me with my issue.  This GIS stuff is using parts of my brain that have remained idle since high school, and some parts are working faster than others.

I have land data stored as Albers projection and I want it converted to Lat, Lon so I can display it on Bing maps.  The prj.adf file in the folder containing my data is this:

Projection    ALBERS
Datum         NAD83
Spheroid      GRS80
Units         METERS
Zunits        NO
Xshift        0.0
Yshift        0.0

Parameters
  29 30  0.0 /* 1st standard parallel
  45 30  0.0 /* 2nd standard parallel
 -96  0  0.0 /* central meridian
  23  0  0.0 /* latitude of projection's origin
  0.0 /* false easting (meters)
  0.0 /* false northing (meters)

I have been searching the web for said transformations, and have come close with the article by James Murtha (http://www.delphi3000.com/articles/article_4479.asp), but this apparently maps Lat Lon to his 1000 scale map, which is not the same scale as my data.  The units described in my prj file show "Meters" - so do I need to calculate the number of meters between the 1st and 2nd standard parallel and assume that this is my unit size in the y direction?

Does anyone have any formulas (or links to information) that I can apply to map between (Lat, Lon) and the Albers (X,Y) projections?  Sorry, but I do not have ARCGIS or any other tools to help me so I have to do this in native code without the conversion libraries in these tools.

Thanks - Steve Tufty


texts:
I have these points on Energy performance and I want to create a heatmap in QGIS following the same logic in the colour scheme.

I would like to create a heatmap ideally a vector or raster based on these values, like the example below.

I know there is a plugin in QGIS but which steps do I need to take so that my raster has the same values as my points?


texts: Finding furthest distance between two UK buildings that share same postcode?

texts:
What is the furthest distance between two UK buildings that share the same postcode?

I'm not sure how to find this out myself.


texts:
Is there any Google API (or any other) which can provide boundary information (i.e. Lat, Lng) for a particular location. Location can be a city, district or a region with in a city. For Example Dammam city (Saudi Arabia), Al-Azizia (District of Dammam) etc

Currently my requirement is, data for, Ta'if. I want Lat, Lng coordinates because I want to draw polygon around region.

Here is another similar question, which dosen't provide complete answer. Seeking neighborhood boundaries for Saudi Arabia?


texts: Find distance between points and borders

texts:
I'd like to ask for your help for the following problem.  I would like to find the distance for each points from the nearest border, represented by the polyline feature here.

I was thinking about using the Near Tool but when I use the latter with the unique polyline of borders as input feature and the points as nearest feature, I get an error 999999.


texts: Using REST Services acquiring limited data from DataBase like LayerDefinition

texts:
I am working with ArcGIS Server 10.2 and MS SQL DB.
I have a requirement in which i need to retrieve limited data using a query from a bulk data using the same Map Service REST URL which can give me filtered Spatial Data while interacting from server to DB. This filtered spatial data will be added to the Map.

Scenario: I have 17 Town data which when added to the map makes the map loading and rendering slow in hosted application. If i need only 4 town data then i should be able to get only data of the same from database using Map Service REST URL without change in URL.

My Understanding:I have some understanding where the REST API can help me in creating a temporary Table in Server side which will collect newly filtered data. This table can be accessed with the same URL but the destination of the table will be changed.

If anyone had came across the same requirement or implementation kindly share with me so that here it will help lot of people who are facing such issues.


texts: How to edit the metadata for individual bands of a multiband raster, preferably with GDAL?

texts:
I created a multiband raster using the gdal_merge.py utility where each band represents a measurement from a different year:

gdal_merge.py -o stack.tif -ot float32 -separate -a_nodata -9999 DEM_2012.tif DEM_2013.tif DEM_2015.tif DEM_2016.tif DEM_2017.tif DEM_2018.tif

However, the information about which year corresponds to which band is lost in the resulting multiband raster. Is there away to retain this information during the creation of the multiband raster? If not, how can I simply access and edit the individual band metadata after the raster is created?


texts: Build a leaflet polygon(s) from GeoJSON feature(collection)

texts:
I find it quite practical to work with leaflet Polygon and it would be great to be able to build such objects directly from GeoJSON content. I did not find any Polygon constructor taking neither a GeoJSON object nor a GeoJSON string as input. Is there some method to build a leaflet polygon directly from a feature or an array of polygons from a featurecollection?

Something like:

geoJ = GeoJSON object returned/parsed from ajax call
pol = L.polygon( geoJ.features[0] )
bounds = pol.getBounds()


texts: landsat 8 images in ERDAS imagine 2014 is not mosaicking properly

texts: I am trying to merge all 4 landsat 8 images in order to create land use land cover map, but its not merging properly.I layerstacked this images before merging. Should i download landsat images which is captured in same day to mosaic it properly. All four images are captured on different dates.


texts: Access qgis.explodelines in Python (outside QGIS)

texts:
I want to access the explodelines function in Python outside of QGIS from a standalone script.

import sip
for api in ["QDate", "QDateTime", "QString", "QTextStream", "QTime", "QUrl", "QVariant"]:
    sip.setapi(api, 2)

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

# Prepare processing framework
sys.path.append('C:/OSGeo4W64/apps/qgis/python/plugins')

from processing.core.Processing import Processing
Processing.initialize()
from processing.tools import *

# Run the algorithm
layerInput = QgsVectorLayer("files/centre/colour/streets.shp", 'test', 'ogr')
general.runalg('qgis:explodelines', layerInput, 'temp.shp')

# Exit applications
QgsApplication.exitQgis()
QApplication.exit()

EDIT:

The result of the algorithm is: 'NoneType' object has no attribute 'error'
I used various shapefiles but none of them worked.


texts: Layers do not align after projection using QGIS?

texts:
I work on QGIS 3.4.4.

My data

I have 4 .tif rasters that will eventually have to be stitched together. The rasters' metadata suggest EPSG:3577 as the initial CRS.
I also have a constriction point shapefile having the following coordinates: 1316205E; -2581235N

My aim

Is to make sure the point shapefile sits exactly where the 4 rasters meet and having the CRS same as the initial ones (EPSG:3577).

My steps

Set the options for default CRS (see pic below):

Created the XY point and saved as a .shp with the EPSG:3577 CRS.

Have used "warp" from GDAL package to project the 4 rasters, with same target CRS as the shp, set the "Output file resolution in target georeferenced units" as 30 (my understanding is that should set the pixels resolution as 30x30m)

At the end, the resulting image is as follows:

This is with google earth image in the background:

What am I doing wrong?


texts: Difference between "hard edge" and "soft edge" in Arc's Add TIN Line?

texts: Is there a technical difference between "hard edge" and "soft edge" in Arc's Add TIN Line on the TIN Editing toolbar? In both you can specify the height of this constraint and they seem to be doing the same thing.


texts: Cartographic Depth-To-Water method using Cost/Path Distance

texts:
I'm having trouble understanding the method for creating a cartographic depth to water index and seek some guidance. I'm using ArcGIS Pro 2.3, but equally happy with ArcMap desktop.

The method is outlined in Murphy, Ogilvie and Arp (2009) "Topographic Modelling of Soil Moisture Conditions: A comparison and verification of two Models" European Journal of Soil Science 60, p. 94-109.

The model takes two inputs: streams and a slope grid.  Excerpt from paper below:

"An iterative function finds the cumulative slope value associated with the least slope path from the cell in the landscape to a surface water cell. That is, the function explores each possible path from the cell in the landscape to a surface water feature across the slope grid and sums the cell slope values along those paths (cumulative slope values). The function takes account of whether the path crosses a cell parallel to the cell boundaries or diagonally by using a multiplier, a. Therefore, the path of least cumulative slope will be identified on the basis of both slope and distance, the cumulative slope value tending to increase with both for any possible path.

The cumulative slope value associated with the least slope path is then multiplied by the grid cell size (resolution) and assigned to the cell in the landscape. This is carried out for each cell. The result is a DTW grid, formally derived from:

where dz/dx is the slope of a cell, i represents a cell along the path, a is 1 when the path crosses the cell parallel to the cell boundaries and 1.414214 when it crosses diagonally and xc is the grid cell size [m].

Because the cumulative values over the whole path of axc and dx will be approximately equal, the index value approximates the elevation difference between the cell in the landscape and the hydrologic source cell (the cumulative value of dz over the whole path).

In this process, all surface water features (lakes, streams, rivers) are assumed to have a DTW value of 0. Lower values indicate wetter soils. Values tend to increase away from surface water features into the landscape, indicating drier soils. DTW increases more rapidly in steeper terrain (greater slope values) and more slowly in flatter terrain (smaller slope values).

So, from this, I understood the streams and slope to be inputs into a cost function and that the output is multiplied by the cell resolution. Firstly, I couldn't decide on whether cost or path distance would be better, and if using Path, in which input would the slope be best (i.e. cost raster?, surface raster? vertical factor?) Both seem to take account of the horizontal or diagonal travel mentioned in the paper.

However, I feel like there's more to it than that. Essentially the method is trying to calculate depth to the water table, and my results are looking nothing like that. The last two paragraphs are not making a lot of sense to me. What am I missing?


texts: Clipping ROI in TerraClimate - Google Earth Engine

texts:
I am trying to analyze a polygon that I imported from Google Earth and get images for that polygon only. How can I do that? I tried clipping the area but it doesn't seem to work. Any ideas?
// Load the TerraClimate image collection and filter it by a date range
var dataset = ee.ImageCollection('IDAHO_EPSCOR/TERRACLIMATE')
                  .filter(ee.Filter.date('1979-01', '2016-01'));

// Load the region of interest
var subwatershed = ee.FeatureCollection('users/cherratc/subwatershed-polygon');
subwatershed = subwatershed.geometry();

**//Clip the region of interest
var clipToCol = function(image){
  return image.clip(subwatershed);
};**


texts: Join selected features from different shapefiles (QGIS)

texts:
I have two different shapefiles, with selected points.
I want to create a new shapefile with only the selected features of those two shapefiles.


texts: Missing 'Expression' option in 'Random Points in Polygon' tool

texts:
I am running QGIS Desktop version 3.8.2
I need to use the 'Random Points in Polygon' tool as part of a process.
The number of points per polygon is defined by one of the fields so I simply add that field name in the 'expression' option and it works without issues.

My colleague who is runnning QGIS Desktop version 3.10.2 needs to run the same process as me however her version of the tool does not have the 'Expression' parameter.

I have tried sharing the tool python script but that errors as the parameter 'EXPRESSION' is not recognised.
Can I locate and share my version of the tool?


texts: Automating GIS slope and aspect algorithm?

texts:
I am currently obtaining slope and aspect of roofs using ArcGIS as follows:

Mosaic several Lidar .asc files to new raster
Extract by mask using Mastermap building outlines to get hold of just the building pixels.
Spatial Analyst Tools – Surface – Slope
Spatial Analyst Tools – Zonal – Zonal Statistics as Table - mean

I have a large amount of data to process.  I could automate this using ModelBuilder in ArcGIS but I’m sure this is not the most efficient method of achieving the results I need.  I’m considering the following:

Script everything in Python
Combination of Python/QGIS/GRASS
R software
Python/Saga combination

Does anyone have any experience of a similar problem and how did they solve it?

Does anyone have an opinion on which software or language I should select?


texts:
I have a sample .csv file (sample download) that I am trying to convert to a gridded raster. So far, I have used two methods, one of the methods is explained in the answer here. But, I am running into problems, how can I fix this problem in R? The code for the first method is as follows:
# CSV to raster
library(raster)
# sample data
sample <- read.csv("Sample_Data.csv", header=TRUE)
x <- raster(xmn=-110.906, xmx=-110.031, ymn=43.03938, ymx=44.3438, res=0.0625, crs="+proj=longlat +datum=WGS84")
rstr <- rasterize(sample[, c('X', 'Y')], x, sample[, 'VIC_SC_04'])
plot(rstr)

This code results the plot nelow and as you can see it's not gridded.

So I tried another method, and the code for this method is:
# CSV to raster
library(raster)
library(sp)
library(rgdal)
library(rgeos)

# sample data
sample <- read.csv("Sample_Data.csv", header=TRUE)
data<-data.frame(sample$X,sample$Y,sample$VIC_SC_04)

# points from scratch
coords = cbind(sample$X,sample$Y)
sp = SpatialPoints(coords)

# make spatial data frame
spdf = SpatialPointsDataFrame(sp, data)
spdf_pi = SpatialPixelsDataFrame(sp,tolerance = 0.016129, data) # Storing incorrect values in Z basically storing X coordinates in z as well**

    raster = raster(spdf_pi[data$VIC_SC_04]) # gives error Error in which(sapply(from@data, is.factor)) :
  argument to 'which' is not logical

# back to data
as.data.frame(data)
plot(data)

# Another try
dfr <- rasterFromXYZ(data)  #Convert first two columns as lon-lat and third as value (throws an error)
plot(dfr)

# coerce to SpatialPixelsDataFrame
gridded(spdf_pi) <- TRUE

# coerce to raster
rasterDF <- raster(spdf_pi)
plot(rasterDF)

When I plot 'data', I get:
plot(data)

But in the following line, I get the following respective error:
dfr <- rasterFromXYZ(data)
Error in rasterFromXYZ(data) : x cell sizes are not regular

The resultant raster should look like this:


texts: Displaying dynamic text from layer other than index on Data Driven Pages?

texts:
For a map project I am working on with Data Driven Pages, I'd like to create a dynamic text using information from a different layer than my index layer.

Is this possible?
Any add ins or extensions available to make this happen?


texts: FME Dynamic schema error - Features With No Schema defined

texts:
I'm trying to update a column in an ArcSDE table using the sde30 reader/writers.

Keep getting an error:

  Features With No Schema defined

and no features are updated.

Scrolling up there is a message in the log:

  readschema resulted in 0 features

The reader and writer are pointing at the same db and table.

I have the writer in update mode, with 'standard feature_type' dynamic schema based on the reader (as I don't want to update all columns, just one per run which i alter dynamically)

I had it working well with the geodb writer, but switched to sde30 for performance reasons.

If I set the 'feature types to read' on the reader, I no longer get the schema missing error, instead I get:

  Error retrieving dbtune keyword list:SDE_ERROR_CODE(-66) Invalid
  parameter value passed to function

I think I am doing something fundamentally wrong!


texts:
Is there an event fired when a kml layer has finished loading in google maps?

Google maps api version 3

To add a layer I have code like:

var layer = new google.maps.KmlLayer(layerUrl, { map: map, preserveViewport : true });

To remove a layer, I have code like:

layer.setMap(null);
layer = null;

The kml are rendered by google maps as images on the clients machine (i.e. the polygons are not created on the client).

I have a number of kml files displayed on a map.  When the user clicks a button, I remove those KML layers, and replace them with other similar KML layers (depending on user input).

The transition is a bit slow.  The user sees the map without any KML layer for about 30 seconds.

I would like to display a loading sign, from the moment the user clicks the button, to the time that the KML layer finally appears.  I need some way to decide when to remove the loading sign.

Thanks

edit:

Looking at stackoverflow (I should have checked), they suggest using:

google.maps.event.addListener(layer, 'metadata_changed', function () {

e.g. here:

https://stackoverflow.com/questions/5668816/google-maps-api-3-loading-icon

On my first draft, that event seems to fire some time before the layer is actually displayed, but I shall investigate further.


texts: Exporting feature collection after being flattened using Google Earth Engine

texts:
I'm trying to export a feature collection of feature collections after being flattened. I cannot get the collection exported with the columns appended (in return value) and just flag an error with 'reduceRegion'.
Can anyone help me export this collection?
The link is here for what I'm trying to do: https://code.earthengine.google.com/c1b2986ccd3019b8bd94220d887a0c32


texts: Geoserver GetLegendGraphic for a layer group shows duplicate styles for all of its layers

texts: I'm trying to get the legend for a Geoserver layer group using a GetLegendGraphic request.  This layer group contains 50+ layers, all of which use the same style.  As expected, the returned image contains the legend for that style, repeated 50 times.  Does anyone know of any way to force geoserver to group identical styles in the legends?


texts: Getting the maximum value of one image in Google Earth Engine

texts: The Google Earth Engine API doc has function to get max value of some bands. But how to get the max value of one image?


texts: GET encoding for static WMTSCapabilities.xml

texts:
Supposedly, one can turn a collection of static tiles, hosted according to the ZXY schema also used by OpenStreetMap tiles and Google Maps tiles, into a standard-conformant OGC WMTS by accompanying it by an (also static) WMTSCapabilities.xml file for satisfying the GetCapabilities WMTS operation.

Amongst others, this document describes how itself may be obtained:

  All WMTS servers operating in a procedure oriented architecture style
  and using HTTP SHALL specify the message encodings that MAY be sent
  using HTTP GET transfer of operation requests. Specifically, an
  ows:Constraint parameter SHOULD be included, with "GetEncoding" as the
  value of the "name" attribute and specifying the values allowed:

  a) The value "KVP" indicates that KVP encoding is allowed, when using
  HTTP GET transfer as specified in clause 8.

If I understood correctly, hosting the static WMTSCapabilities.xml at an URL like https://example.com/mymap/1.0.0/WMTSCapabilities.xml means that the request to it would not be "KVP encoded" but "RESTful encoded". However, the standard only names "KVP" as a possible value.

So, if the WMTSCapabilities.xml is only made available at https://example.com/mymap/1.0.0/WMTSCapabilities.xml, not also at e.g. https://example.com/mymap.cgi?service=WMTS&version=1.0.0&request=GetCapabilities (which would be an URL for a KVP request), should the operations metadata read

<!-- ... -->
<ows:OperationsMetadata>
   <ows:Operation name="GetCapabilities">
      <ows:DCP>
         <ows:HTTP>
            <ows:Get xlink:href="https://example.com/map/1.0.0/WMTSCapabilities.xml">
               <ows:Constraint name="GetEncoding">
                  <ows:AllowedValues>
                     <ows:Value>RESTful</ows:Value>
                  </ows:AllowedValues>
               </ows:Constraint>
            </ows:Get>
         </ows:HTTP>
      </ows:DCP>
   </ows:Operation>
   <!-- ... -->
</ows:OperationsMetadata>
<!-- ... -->

or just

<!-- ... -->
<ows:OperationsMetadata>
   <ows:Operation name="GetCapabilities">
      <ows:DCP>
         <ows:HTTP>
            <ows:Get xlink:href="https://example.com/map/1.0.0/WMTSCapabilities.xml">
               <ows:Constraint name="GetEncoding">
                  <ows:NoValues />
               </ows:Constraint>
            </ows:Get>
         </ows:HTTP>
      </ows:DCP>
   </ows:Operation>
   <!-- ... -->
</ows:OperationsMetadata>
<!-- ... -->

?


texts: WFS - Which values should I enter as `namespace` and `featuretype` - Where can these values be read?

texts:
I like to use WFS - GetFeature request, that should return a selection of features from a GeoServer data source.

The request

http://localhost:8080/geoserver/wfs?service=wfs&version=2.0.0&request=GetFeature&typeNames=namespace:featuretype

shows me the following in the browser Firefox:

<ows:ExceptionReport version="2.0.0" xsi:schemaLocation="http://www.opengis.net/ows/1.1 http://localhost:8080/geoserver/schemas/ows/1.1.0/owsAll.xsd">
<ows:Exception exceptionCode="InvalidParameterValue" locator="namespace">
<ows:ExceptionText>Unknown namespace [namespace]</ows:ExceptionText>
</ows:Exception>
</ows:ExceptionReport>

I do not know exactly where I can see on the server which values I should enter as namespace and featuretype.

I tried a lot.

For example I tried

http://localhost:8080/geoserver/wfs?service=wfs&version=2.0.0&request=GetFeature&typeNames=namespace:featuretype

That showed me the following

http://localhost:8080/geoserver/wfs?service=wfs&version=2.0.0&request=GetFeature&typeNames=%20KRE_ALKIS:sg_flurstueck_f

<ows:ExceptionReport version="2.0.0" xsi:schemaLocation="http://www.opengis.net/ows/1.1 http://localhost:8080/geoserver/schemas/ows/1.1.0/owsAll.xsd">
<ows:Exception exceptionCode="InvalidParameterValue" locator="typeName">
<ows:ExceptionText>Feature type KRE_ALKIS:sg_flurstueck_f unknown</ows:ExceptionText>
</ows:Exception>
</ows:ExceptionReport>

Which values should I enter as namespace and featuretype - Where can these values be read?


texts: How to deserialize geoJon input to locationtech.jts.geom.Geometry

texts:
I have following input to a Spring ReST endpoint:
{
    "featureType": "feature",
    "geometry": {
        "type": "Point",
        "coordinates": [
          65.0390625,
          62.91523303947614
        ]
      }
}

I am trying to map it to a custom class I have defined as:
public class Feature   {
  @JsonProperty("featureType")
  private String featureType = null;

  @JsonProperty("geometry")
  private Geometry geometry = null;
}

Geometry is org.locationtech.jts.geom.Geometry.
I have following pom dependency:
<dependency>
            <groupId>com.graphhopper.external</groupId>
            <artifactId>jackson-datatype-jts</artifactId>
            <version>1.0-2.7</version>
</dependency>

And I have configured JtsModule as bean:
   @Bean
   public JtsModule jtsModule()
   {
      return new JtsModule();
   }

But I get error:
"Type definition error: [simple type, class org.locationtech.jts.geom.Geometry]; nested exception is com.fasterxml.jackson.databind.exc.InvalidDefinitionException: Cannot construct instance of `org.locationtech.jts.geom.Geometry` (no Creators, like default construct, exist)

I am not able to get what part I am missing here.


texts: Opening vector files in QGIS 2.18 causes error and crash

texts:
When trying to open vector files (.shp or .tab) in QGIS 2.18 and 2.18.3 on my desktop using either drag and drop or through Add Layer, I get the "qgis-bin.exe has stopped working" message and crashes. Otherwise QGIS appears to open without problems.
If I open a raster file (.ecw or .tab) I can then open the vector and display it without a problem. Strange.

I have tried:

a) Reinstalling both 2.18 and 2.18.3 on the desktop

b) Using older version 2.14 on desktop - everything appears to work fine with same files

c) Checking 2.18.3 on my laptop - everything appears to work fine

2.18 was working well on the desktop until I upgraded to 2.18.3 however this version seems to be ok on the laptop so doesn't appear to be the programme. Have never had this issue before over the years of installing QGIS as usually installs seamlessly so has got me a bit baffled and annoyed. May be something obvious?

Further to this (64bit version, with or without GRASS):

1) Browser works fine

2) If I remove the raster and leave the vector, it crashes. If I have the raster and add 2 vectors, it's fine. If I remove 1 vector, it's fine, if I remove the other vector leaving the raster, it's fine. If I have 2 vectors and remove raster it's fine. If I then remove 1 vector leaving 1 vector it's fine. If I remove everything it's fine. If I then load vector again, it crashes.

Further comment. In the absence of any other suggestions (eg possible issues with plugins or anything else), that an uninstall/reinstall doesn't seems to fix it and the fact that 2.18.3 works fine on my 32bit laptop and on another available 64bit desktop, seems this will remain in the realm of the unexplained and may need to do a reformat of system to sort it. Not really my idea of an optimum solution but if all else fails....


texts: Customize rules in Random Forest for tree species classification

texts:
Is it possible to customize the rules in Random Forest algorithm for tree species classification?

all the integration or customization of rules in random forest algorithm is done in medical field, crime, drugs.

I have never found it done for tree species.

If it is possible can anyone suggest the suitable software


texts: Network analysis: shortest path to 'visit' all polygons

texts:
I'm totally new to QGIS but I really like it and would like to solve something more elaborated. I found a lot of threads and tutorials about how to solve the traveling salesman problem with QGIS but my 'problem' is a bit different.

I don't want to set a start and endpoint (at least no endpoint) but have multiple polygons and the goals is to 'visit' all of them. In fact the polygons are the districts of a city and I want to know which path to take to visit all of them while taking the shortest route.
I also have a shapefile with the streets of course.

Maybe there are already solutions for this out there and I just didn't use the right wording in my search queries.


texts: Watersheds tool does not create watersheds in ArcMap?

texts:
I`m trying to create watersheds and every time I run the Watershed tool, the result will be very interesting.

I use snap raster and I set the extent of the output to the extent of the flow direction raster, but it does not make any difference.

Does anyone know how to fix this issue?


texts: Python module to delete SHP features (without Desktop GIS installed)

texts:
I have a piece of software (not Arc) that runs nightly on a non-production machine which updates shapefile features from an external database.  Periodically I would like to delete all of the features in the shapefile (not the file itself, which must remain) and let the software "rebuild" the shapefile from scratch.  I would like to automate this process.

I don't have any GIS software currently installed on that machine.  I was hoping that I could script a routine in Python that would automatically delete the features, much like Arc's Delete Features geoprocessing tool.

Are there any Python modules that would allow me to do this?  Preferrably open-source?  I looked at Shapely and PyShp but didn't see anything that would allow me to delete the features in mass or that matched a WHERE clause.  They can write features and analyze them, but haven't seen DELETE FEATURES functions anywhere.

I must certainly be missing something simple...

EDIT : I have 35 folders (different geographic extents, all in their own projection), with 35-65 shapefiles with makes around 1000 shapefiles to deal with.


texts: Adding an inset map to Arcgis Online web app or web map?

texts:
Is it possible to add an inset map frame to an ArcGIS Web app? I would ideally have the mid-east as a main focus, then a small inset map of the united states in the corner.

I know how to do this in a static map using ArcMap desktop, but didn't know if it was possible to do in ArcGIS online.


texts:
Does anyone know how to add labels (using a feature column such as "state") to a CARTO.js layer? This is how I have added the layer to my map:

const boundariesDataset = new carto.source.Dataset(`
      districtmunicipality2016
    `);
    const boundariesStyle = new carto.style.CartoCSS(`
      #layer {
        ::outline {
          line-width: 1;
          line-color: #000000;
          line-opacity: 1;
        }
      }
    `);
    const boundaries = new carto.layer.Layer(boundariesDataset, boundariesStyle, {
        featureOverColumns: ['district_n']
    });

Ps, this line:

featureOverColumns: ['district_n']

Is used to select a feature to display in a pop up, using:

const popup = L.popup({ closeButton: false });
        boundaries.on(carto.layer.events.FEATURE_OVER, featureEvent => {
          popup.setLatLng(featureEvent.latLng);
          if (!popup.isOpen()) {
            popup.setContent(featureEvent.data.district_n);
            popup.openOn(map);
          }
        });


texts: Where is the option to correct for rolling shutter in ESRI's Drone2Map software?

texts: I am looking for the proper way to enable a software fix for rolling shutter in Drone2Map.


texts: Solve OD-matrix for distance, add travel time to outputs

texts:
I have an ArcGIS bike network with two kinds of streets. The first are just regular streets, where a biker has a speed of 6km/h. The second type are dedicated bike lanes where a biker has a speed of 10km/h.

Now i would like to solve a closest facility problem with the network analyst. I would like to solve for closest distance, but I also would like to know the travel time that belongs to this closest distance.

Is there an easy way to do this?


texts:
What are sources where I can learn how to use GeoServer?

I know of these two sites:

http://geoserver.org/ - which in fact is not working (at least at the moment)
http://workshops.opengeo.org/geoserver-intro/ - not only (tutorials) for GeoServer

But I am interested in more complex information, and not to read only documentation.


texts: Why do FusionTable layers sometimes not work / slow down?

texts:
Fusion tables layers are great.  I have uploaded some data and polygons of areas in Europe and most of the time they respond very fast.

First thing in the morning (after prolonged inactivity I assume), the layer will sometimes take some time (about 5 minutes) before it appears.

Occasionally (about once a day), even whilst I am looking at and using the map layer, it drastically slows down, or stops working all together.  This lasts for roughly half an hour.

Are these common and unavoidable problem of using a service I neither control nor paid for, or is there something I might have done wrong?  Is anyone else experiencing issues?

Thanks


texts: importing MapInfo georeferenced ratser on qgis 2.2

texts:
I do not understand why loading mapInfo shapefiles is possible in Qgis but for the georeferenced image file, is there any possible reason for this?

My major concern is that, so far georeferencing on Qgis is problematic to me (I only get to see the final result only after georefencing, which is not always well done as seen by the distortion). Given that MapInfo permits you correct your precision, I will like to georeference on MapInfo and import the file to Qgis, can anyone help me:

Bring in this file, or
Enlighten me more on how to better my georeferencing.


texts:
I am quite new to EE. Currently, I am trying to 'burn' polygons from a feature collection to an image. I tried a few references (link1, link2), but I haven't been able to use them effectively for the scope I am about.
I tried:
dataset = ee.FeatureCollection("WCMC/WDPA/current/polygons") \
    .filterBounds(ee.Geometry.Rectangle(*bbox)) # [xmin,ymin,xmax,ymax]

# Option 1
image = ee.Image().float().paint(dataset, "REP_AREA")

# Option 2
image = dataset.reduceToImage(["REP_AREA"], "mean")

But the output is a GeoTif cropped on the bounding box, with all pixels having the same value. I exported it with geemap:
geemap.ee_export_image(image, 'protected_areas.tif', region=boundary) # ee.Geometry.LineString

If I use the inspector directly into EE's code editor, I clearly see that the polygons are there with the property I'd expect. What am I doing wrong?


texts:
I am currently working with 3 data frames in ArcMap. I exported them to pdf then Adobe Illustrator. I noticed that 2 of the data frames that include rasters come up with striped layers on them (while the third data frame comes up fine and I can edit the features). I deleted the rasters in ArcMap and then imported it to Adobe and I realized that yes, the reason for the problem is because of the rasters.

I was wondering how I can import these rasters to Adobe? I tried converting the raster in ArcMap as a PNG, but when I import it to Adobe, it comes up as a black rectangle.


texts: Overpass XML query with multiple polygons or bounding boxes

texts:
I am having trouble creating an OpenStreeMap Overpass XML combined polygon/bounding box query to return what I want. Can someone provide a simple example or correct the following XML?

I want to construct a query to retrieve data from multiple polygons or bounding boxes. Specifically, I want to retrieve all data within a small polygon and only certain data types, major roads, in a larger containing bounding box and return the combined result.

This XML shows conceptually what I want to do but is returning all data for the larger area rather than all data for the small area and major roads for the larger area:

<osm-script timeout="20" element-limit="1073741824">
    <union into="_">
        <!-- Smaller area, all data -->
        <polygon-query into="_" bounds="54.12 -2.779  54.143 -2.779  54.143 -2.761 54.12 -2.761" />
        <recurse from="_" into="_" type="up"/>
        <recurse from="_" into="_" type="down"/>
    </union>
    <union into="_">
        <!-- Larger area, major roads only -->
        <bbox-query s="54.12" n="54.143" w="-2.779" e="-2.744" />
        <has-kv k="highway" v="motorway"/>
        <recurse from="_" into="_" type="up"/>
        <recurse from="_" into="_" type="down"/>
    </union>
    <print from="_" limit="" order="id"/>
</osm-script>


texts: Make "qgis:delaunaytriangulation" and "qgis:dissolve" save their outputs on a specific predefined layer

texts:
I'm new to python and QGIS. I'm working on modifying a plugin which runs algorithms "qgis:delaunaytriangulation" and "qgis:dissolve" like so

delone_triangles = processing.runalg("qgis:delaunaytriangulation", assistant_layer, None)['OUTPUT']
delaunay_layer = processing.getObject(delone_triangles)
#some code later
dissolved = processing.runalg("qgis:dissolve", delaunay_layer,
                                          True, '', None)['OUTPUT']
dissolved_layer = processing.getObject(dissolved)

and from what I can gather what it does is it creates temporary files, saves the path to those files in both the delone_triangles and dissolved variables respectively and then uses that file to create the delaunay_layer and dissolved_layer objects. So with each iteration it essentially creates more  unnecessary files. According to their documentation both of those algorithms have an output parameter, but I haven't been able to figure out how to properly use it and the documentation is not helpful
https://docs.qgis.org/2.6/en/docs/user_manual/processing_algs/qgis/vector_geometry_tools/delaunaytriangulation.html

https://docs.qgis.org/2.6/en/docs/user_manual/processing_algs/qgis/vector_geometry_tools/delaunaytriangulation.html

What I wanna do is I wanna set up both the delaunay_layer layers before and then have those algorithms save their outputs onto those layers. Is that possible? How do I do that?


texts: Copy field value from a point feature class to another field in a line feature class

texts: I have a point feature  and a line feature in a geometric network. I am trying to find a way ( Tool or script) to copy values from a field in the point feature and paste them in the line feature. Any suggestions?


texts:
I have a problem printing map out of Composer. I need to print Corine Land Cover but when I try to print the map, many of corine layers aren't visible on printed map. They are loaded in project and visible in it, but when I try to print them to PDF, many of the dissapear.

Is this because of too many data?
This is project preview. As you can see, all data is loaded in visible.

And this is the result - visible only few layers, but not all.
I tried to close and open the program, refresh project and composer, print to JPG, but result is the same.

Any idea?


texts:
I use GeoExt2 in my app through online path:

Ext.Loader.setConfig({
    enabled: true,
    paths: {
        GeoExt: "http://geoext.github.com/geoext2/src/GeoExt/"
        ,'Ext.ux': 'http://extjs.cachefly.net/ext-4.0.7-commercial/examples/ux'

    }
});

I need to have maps and popups available offline (without internet), so I think I have to download GeoExt2 source and use it locally. Is that possible/recommended? due to the fact that it is still alpha version.


texts: R - Programmatically convert raster stack in data frame based on polygon extraction

texts:
Sorry for cross-posting. I also posted this question at the R-SIG-GEO discussion list, but since I wanted to get as much feedback as possible I decided to post it here too.

I am trying to extract temperature values from a raster stack for about 400 municipalities in Brazil. My final goal is to create a data frame that is going to be used as a database for an interactive map server - probably using shiny and leaflet.

The final data frame would look like this:

> head(df)
Location        Var   Cut Year Month Freq
Campinas  temperature  10 2010  1    11
Campinas  temperature  10 2010  2    19
Campinas  temperature  10 2010  3    30
Campinas  temperature  10 2010  4    29
Campinas  temperature  10 2010  5    31
Campinas  temperature  10 2010  6    30

I have global raster stacks with daily data and I am counting, for each month in the raster, the number of days above certain temperature threshold. Please see below:

library(raster)
library(zoo)
library(maptools)

# Create a rasterStack similar to my data - same dimensions and layer names
r <- raster(ncol=360, nrow=180)
s <- stack(lapply(1:730, function(x) setValues(r, runif(ncell(r),min=0,max=30))))
idx <- seq(as.Date("2010/1/1"), by = "day", length.out = 730)
s <- setZ(s, idx)
s

# Define functions for 10, 15, 20 and 25 degrees
fun1 <- function(x, na.rm) {
sum(x > 10, na.rm)
}

fun2 <- function(x, na.rm) {
sum(x > 15, na.rm)
}

fun3 <- function(x, na.rm) {
sum(x > 20, na.rm)
}

fun4 <- function(x, na.rm) {
sum(x > 25, na.rm)
}

# Count number of days above the threshold temperature
days.above.10 <- zApply(s, by=as.yearmon, fun = fun1)
days.above.15 <- zApply(s, by=as.yearmon, fun = fun2)
days.above.20 <- zApply(s, by=as.yearmon, fun = fun3)
days.above.25 <- zApply(s, by=as.yearmon, fun = fun4)

Now, what I would like to do is to programmatically extract values for each location on my study area. The locations are defined as a shapefile with municipal contours of the Sao Paulo state in Brazil.

In this example, however, just for reproducibility's sake, I will be using a world polygon. But keep in mind that in my actual data the polygons will be much smaller.

# Import *sample* polygon data and subset only five "locations"
data(wrld_simpl)
locs <- subset(wrld_simpl, wrld_simpl@data$NAME %in% c("Argentina","Bolivia","Brazil","Paraguay","Uruguay"))

# Plot
plot(days.above.10,1)
plot(locs,add=T)

I feel like half of the work is done, but I am just grasping with the conversion to data frames.

Based on this self-contained example I provided, what would be the best strategy to come out with a data frame per location, like this?

> head(Argentina.df)
Location        Var Cut Year Month Freq
Argentina temperature  10 2010    1  11
Argentina temperature  10 2010    2  19
Argentina temperature  10 2010    3  30
Argentina temperature  10 2010    4  12
Argentina temperature  10 2010    5  17
Argentina temperature  10 2010    6  14

> head(Bolivia.df)
Location        Var  Cut Year Month Freq
Bolivia  temperature  10 2010    1  29
Bolivia  temperature  10 2010    2  31
Bolivia  temperature  10 2010    3  30
Bolivia  temperature  10 2010    4  17
Bolivia  temperature  10 2010    5  19
Bolivia  temperature  10 2010    6  12

and so on.

Note that "cut" refers to the temperature thresholds defined in the functions above. Each cut should come from the equivalent raster stack: days.above.10, days.above.15 and so on.


texts: Changing field values selected attributes of shapefile in ArcGIS desktop?

texts:
I have nine fields in a shapefile and want to use Python to change the values of those fields. The value would be given by the user.

The attributes selected would receive the input value each time the script is run. If the first field is non-zero, it would be compared to the value input and if they are not the same, the values would be compared to the second field. If it doesnot find the value in any of the field, it would append the value to the end.

The table looks like this:

If the input is 350, then the table should look like this.

As you can see, I want only the selected rows to be populated with the value. And, the given input may not be received by the same field.

In other words, it should append the given input if it is not present in any of the specified fields.


texts: How to map a dataset with multiple datums

texts:
I have a dataset (in excel) of species with the coordinates in which they were sighted, collected by many different people over decades. Aside from coordinates, the dataset also specifies a datum (Nad_27, Nad_83, WGS_84). I added the data to ArcGIS.

So my fields look like:

The HorizontalDatum column is confusing me. This one excel file says the dataset is composed of 3 different datums depending on the record. Although the original datum is said to be one of those three for any given record, can I set the datum for the entire dataset to Nad_83 in arcGIS and ignore this column? Will that still be accurate? I don't understand datums or transformations enough to be sure.

If not, will I have to split the dataset into 3 shapefiles for records of each datum?

These data points are all within the same province. I just want to keep all the records in the same shapefile.

To summarize and clarify: I will need to set a geographic coordinate system for this shapefile manually, as it's coming from excel data. Can I set the entire shapefile to Nad_83, or should I split it up into 3 shapefiles, and project those with the 3 "Horizontal Datums" stated in the excel field? This would be less than ideal... I'd like to be able to query the whole dataset in one go.


texts: Using mouse-over on raster surface to stream values to text file using ArcGIS Desktop?

texts:
This is a weird one. I have a user who wants to add a statistical surface to ArcMap and then mouse-over the raster surface and stream those values (z) to a text file. MouseTips can show underlying values, but I can't see a way to use it.

I'm not inclined to create features and geoprocess them for the results as the user wanted a real time experience. Several searches in stack and google have not yielded any results.

I am using ArcGIS Desktop 10.3.1.

Is it even possible?


texts: Using ogr2ogr to clip shapefile based on coordinate query

texts:
I have some shapefiles that I would like to convert to GeoJSON and clip to display only Arctic regions (for example, everything above 45°N).

I tried the following, but nothing is clipped.

Any suggestions?

ogr2ogr -f GeoJSON land.json -spat -180 0 180 45 land.shp

and

ogr2ogr -f GeoJSON land.json -clipsrcwhere "y<45" land.shp

Clearly I'm grasping at straws. I haven't been able to find examples of clipsrcwhere in any of my searches. Not sure if it is the right option, but my region is not exactly rectangular.


texts: How does OpenStreetMap specify underspecified addresses?

texts:
I live in Austin, Texas.  If I'm in Google Maps and I type in "2504 E Cesar Chavez", it takes me to an address in Austin and the panel on the left says "2504 E Cesar Chavez St, Austin, TX 78702".  But if I type the same thing into openstreetmap.org, there are 3 different results, none of which are in Austin (one is in New Mexico, one in South Texas, and one in Ecuador).

What's going on here?  I have two big questions:

How does OpenStreetMap figure out what the "real" address is, as compared to Google Maps?
Is there documentation about how to correctly specify street addresses for OpenStreetMap?  Or, does anyone here have tips on how to do it?

(A bit of background: I'm working on a project that will involve putting a few thousand photos on a ViewShare map in Omeka, and am hoping to find a way to standardize location metadata en masse that will play nice with ViewShare, which uses OpenStreetMap.)


texts: Grouping Multiple Raster Layers using ArcGIS Desktop

texts:
How do I create multiple raster layers group, for example when using several aerial photos?

I wish to group the layers into one so that I can edit/select as a batch rather than having to click on individual layers.

I use ArcGIS Desktop 10.2.


texts: Loop Issue: Extracting mean temps for each layer of RasterStack (NetCDF File)

texts:
I am running into a little bit of a problem with a RasterStack loop. I am trying to extract the mean temperature for each county in North Carolina for each hour over the last 10 years. I am able to run some simple code and get this for each layer with relative ease but as I am trying to loop this process to get the mean for every timestamp in the Raster Stack/Brick, I am not having much luck.
library(tidyverse)
library(rgdal)
library(sp)
library(raster)
library(doParallel)
library(foreach)
library(ncdf4)
librarcy(exactextractr)

# Brick the NC File
r<-stack("ERA_Temps.nc")

output <- for (i in 1:length(r)) {

              wgsraster <- r[[i]]

              year <- str_sub(names(wgsraster), 2, 5)
              month <- str_sub(names(wgsraster), 7, 8)
              day <- str_sub(names(wgsraster), 10, 11)
              hour <- str_sub(names(wgsraster), 13, 14)

              NC_Counties <- cbind(shp, exact_extract(wgsraster, shp, c('mean')))

              output_df <- data.frame(County.Name = NC_Counties$NAME, GEOID = NC_Counties$GEOID, temp_mean = NC_Counties$mean, year, month, day, hour)
                    return(output_df)
                  }

Output <- output_df

Examples and More Info:
 r
    class      : RasterStack
    dimensions : 81, 81, 6561, 72  (nrow, ncol, ncell, nlayers)
    resolution : 0.25, 0.25  (x, y)
    extent     : -90.125, -69.875, 24.875, 45.125  (xmin, xmax, ymin, ymax)
    crs        : +proj=longlat +datum=WGS84 +no_defs
    names      : X2018.09.14.00.00.00, X2018.09.14.01.00.00, X2018.09.14.02.00.00, X2018.09.14.03.00.00, X2018.09.14.04.00.00, X2018.09.14.05.00.00, X2018.09.14.06.00.00, X2018.09.14.07.00.00, X2018.09.14.08.00.00, X2018.09.14.09.00.00, X2018.09.14.10.00.00, X2018.09.14.11.00.00, X2018.09.14.12.00.00, X2018.09.14.13.00.00, X2018.09.14.14.00.00, ...

I can fairly easily get what I need for an individual layer:
wgsraster <- r[[5]]

year <- str_sub(names(wgsraster), 2, 5)
month <- str_sub(names(wgsraster), 7, 8)
day <- str_sub(names(wgsraster), 10, 11)
hour <- str_sub(names(wgsraster), 13, 14)

NC_Counties <- cbind(shp, exact_extract(wgsraster, shp, c('mean')))

output_df <- data.frame(County.Name = NC_Counties$NAME, GEOID = NC_Counties$GEOID, temp_mean = NC_Counties$exact_extract.wgsraster..shp..c..mean...,
                                        year, month, day, hour)

output_df
     County.Name GEOID temp_mean year month day hour
1         Saluda 45081  298.3141 2018    09  14   11
2   Williamsburg 45089  300.2301 2018    09  14   11
3     Dorchester 45035  299.3716 2018    09  14   11
4       Berkeley 45015  299.7389 2018    09  14   11
5      Greenwood 45047  298.5175 2018    09  14   11

But when I try and run the loop mentioned above, it either fails and gives no output or gives the following error:

Error in .local(x, i, j, ...) : you must provide an index


texts: Getting WHERE statement to work in v.generalize?

texts:
I have a shape (its topology is terrible but I don't think it's the problem).  I'm trying to simplify it to reduced it's size using v.generalise.  I succeeded at reducing the whole shape but I actually want to do a two steps simplification:

simplifying all the shape with threshold 0.0001 first,
simplifying all the polygons with urban=0 further more with threshold 0.0003

The problem is that my second v.generalize call simplifies everything, like if the WHERE call wasn't working.

For this question, we will just focus on the point 2, so here what I do:

A subset of my data can be found here.

And here is some visualization of what I did:

My shape initially:

My v.generalize call:

The result (below) show the blue line the new polygons border which where simplified even if they were not in urban=0.

Just like if the WHERE statement was ignore.  I don't see any error in the log specific to that SQL statement.  The actual call to GRASS built by QGIS is:

C:\OSGEO4~1\bin>v.generalize input=vector_5be4482a0b8e72 type="line,boundary,area" where="urban=0" method="douglas" threshold=0.0003 look_ahead=7 reduction=50 slide=0.5 angle_thresh=3 degree_thresh=0 closeness_thresh=0 betweeness_thresh=0 alpha=1 beta=1 iterations=1 -l output=output96039149be394bdaa48a3ba7d8b849f7 error=error96039149be394bdaa48a3ba7d8b849f7 --overwrite

I've try multiple WHERE statement (urban="0", urban<1, etc) to try to match the v.generalise tutorial but still no luck.  I also tried to do it directly in GRASS bypassing QGIS and the problem persist meaning is not a QGIS problem more than a GRASS problem.

In this example, the urban field is a real number.

Any idea what the problem is?


texts: Generating points from chainage and offset

texts: I have a polyline layer and an external csv containing road sign data. The csv has unique ID of the polyline, chainage along each unique segment and an offset. I have access to ArcGIS, MapInfo and QGIS, is there any way I can generate point data from this data?


texts: ArcGIS 10.3: Insert rotated image without background

texts: How can I insert this image into an mxd so the white background doesn't show up?  I've tried using photoshop to delete the background and save it as a format that supports transparent backgrounds (gif and png) to no apparent effect.  Creating clipping masks in photoshop and inserting resulting images also doesn't seem to do anything.  I need to fit a bunch of similar rotated signs next to each other, with a bunch of base features in behind them (such as the dark grey lines in the attached image) than need to be visible, so I can't have the backgrounds of each image blocking each other and the base features out.


texts:
I need a contiguous, nationwide USA Road Network with speed limits, distances, etc.  I have been using TIGER but the challenge is at state boundaries - where contiguous road segments are split.

Is there any free or commercially available data set that provides this data, besides ESRI/ArcGIS?


texts: Ensure BaseDynamicLayer is always the most visible layer

texts:
If I add a shapefile to my map and then a dynamic custom layer inherited from a BaseDynamicLayer, then the shapefile is lower on the viewing stack then my custom layer, and it displays correctly.

mapControl.AddShapeFile(directoryPath, shapeFile);
mapControl.AddLayer(dynamicCustomLayer); //this is ok

But if I add my dynamicCustomLayer first, and then add a Shapefile, the shapefile is makes the dynamic layer non viewable because it is higher on the viewing stack...

mapControl.AddLayer(dynamicCustomLayer);
mapControl.AddShapeFile(directoryPath, shapeFile);//now my dynamic layer is invisible

How can I ensure that the dynamic custom layer is always on top of the viewing stack even if my user adds a shapefile after the custom layer was enabled...

I'm using C# ArcEngine 10


texts:
I am trying to copy elements from a shp to a GDB.

I took inspiration from this link, but after execution, my elements do not appear in the database. However, when I look in the table, it says 9 elements are present, but there are no lines in the table.

Also, the OGR log window displays the following error :

  Data source is invalid (Cannot open geodatabase at the moment since it is in 'FID hack mode')

I am new to pyQgis so I'm open to any suggestion about the pseudo code below.

##Conversion to GDB=name
##src1=vector
##src2=vector
##src3=vector

from qgis.core import *
from qgis.gui import *
from qgis.utils import iface
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from time import *
import os
import processing
from datetime import datetime
from PyQt4.QtGui import QDockWidget
import subprocess
import glob
import shutil, errno

#*********************************************************************************************************************************************
#Initialisations
src_layers=[src1, src2, src3]
i=0

print "START____ " + strftime("%Y-%m-%d %H:%M:%S", localtime())

#Groups for output layers
root = QgsProject.instance().layerTreeRoot()
displayGroup = root.insertGroup(0, "GDB_"+str(tm))

#*********************************************************************************************************************************************
#1.Opening

layerlink= #string to open GDB in Qgis(removed from pseudocode), similar path/to/gdb/folder.gdb|layer=dest
GDB_layer = QgsVectorLayer(layerlink, GDB_str, "ogr")

GDB_layer.startEditing()
for srcL in src_layers[i]:
    #Open the source layers
    src_layer= QgsVectorLayer(srcL, srcL, 'ogr')

#*********************************************************************************************************************************************
#2. add source layer fields to the GDB
    src_layer_attrib_names = src_layer.dataProvider().fields()
    src_attributeList = src_layer.dataProvider().fields().toList()
    newattributeList=[]
    for attrib in src_attributeList:
        if GDB_layer.fieldNameIndex(attrib.name())==-1:
            newattributeList.append(QgsField(attrib.name(),attrib.type()))
    GDB_layer.dataProvider().addAttributes(newattributeList)
    GDB_layer.updateFields()

#*********************************************************************************************************************************************
#3. Copy source elements to the GDB

#Copy features over to the the new memory layer
    GDB_layer_attribute_list=GDB_layer.dataProvider().fields().toList()
    cfeatures=[]
    xfeatures = src_layer.getFeatures()
    for xfeature in xfeatures:
        #Create destination feature
        cfeature = QgsFeature()
        #Populate feature with geometry
        xgeometry = xfeature.geometry()
        cfeature.setGeometry(xgeometry)
        #Populate feature with attributes
        cfeature_Attributes=[]
        for GDB_QGSfield in GDB_layer_attribute_list:
            attribute_field=GDB_QGSfield.name()
            #Get old attribute value and append
            idx = src_layer.fieldNameIndex(attribute_field)
            if idx>=0:
                cfeature_Attributes.append(xfeature.attributes()[idx])
        cfeature.setAttributes(cfeature_Attributes)
        cfeatures.append(cfeature)
    GDB_layer.dataProvider().addFeatures(cfeatures)
    GDB_layer.commitChanges()
    GDB_layer.updateExtents()
    print "nb : " + str(GDB_layer.featureCount())
    #here the count should increase but the number of features stays 9

#6. Add layers to display in Qgis

QgsMapLayerRegistry.instance().addMapLayer(GDB_layer, False)
displayGroup.addLayer(GDB_layer)

print "OK_______ " + strftime("%Y-%m-%d %H:%M:%S", localtime())


texts: How to select only areas bigger than x acres?

texts:
I would like to create a new shapefile from another area shapefile. The new shapefile should only contain all areas that are bigger than 5 acres (10.000 m²).

How do I do this?


texts: R: creating an error matrix from two rasters

texts:
I have two land use rasters (6 classes each) and I'd like to create an error matrix (errors of omission and commission) from them both.

The premise is that I have a very fine 'base' raster and have aggregated it to a  new, coarser resolution. I would like to know not just the difference in total area of each land use class but from what and to what land class have changed;

# make a raster to simulate fine res base raster;
require(raster)
ras.fine <- raster(nrows=10, ncols=10,xmn=0, xmx=10, ymn=0, ymx=10)
ras.fine[] <- sample(seq(from = 1, to = 6, by = 1), size = 100, replace = TRUE)

# aggregate it to represent the coarser raster
ras.coarse <- aggregate(ras.fine,fact=2,expand=FALSE,fun=modal,na.rm=T)

Once we have done the above, it's easy enough to establish the total areas of land use classes and establish the difference ;

# use freq to count instances of land class and multiply by square of resolution
data.frame(class = freq(ras.fine)[,1],count = freq(ras.fine)[,2],area = freq(ras.fine)[,2]*(res(ras.fine)[1])^2)

data.frame(class = freq(ras.coarse)[,1],count = freq(ras.coarse)[,2],area = freq(ras.coarse)[,2]*(res(ras.coarse)[1])^2)

but the differences in total area of land class aren't enough; i'd like to know from what land classes the error has occurred from etc, making an error matrix.

As a start, we could directly compare the base raster and the aggregated raster by disaggregating the coarser raster back to the base resolution (is this fair??);

# disaggregate back to original resolution
rc.d <- disaggregate(ras.coarse,fact=2)

# then create a grid of 'disagreement'; that is where the two rasters do not agree
disagree <- ras.fine != rc.d

# and establish what land classes make up those cells of disagreement, from both the fine raster and the coarser raster
fine.cov <- ras.fine * disagree
coarse.cov <- rc.d * disagree

so now I have statistics for how much of each land class is classified incorrectly at a coarser scale, and theoretically what it is going to/coming from.

I'm a bit stuck from here; how do i fashion it into an error matrix? from what land class has another land class commissioned area from and vice versa?

I am essentially analysing the aggregation method but this allows quantification in uncertainty mapped against cost of effort.

Further: Qu on Cross Validation with regards to the validity of error matrices and raster aggregation


texts: GML, KML, GeoJSON - Speed rendering 3109 polygons?

texts:
I am working with Geoserver, serving U.S. Lower 48 counties to openlayers (3109 polygons - lots more vertices).  The counties are loaded into a postgis database.  I am curious about developer experience when trying to push that quantity of vertices to the client.

What WFS format have you attained the best results with?  Has additional tuning to Geoserver been used?

I realize that tiled WMS would be faster, but I want to allow for dynamic changes in a choropleth map using openLayers, ie. the user submits a form, a Python script is called, and new data bins are returned for openlayers to reload the map div.  I also want to try this in full resolution form before reducing polygon complexity in openlayers.


texts:
I have two datasets:

parcel layers
an overlaid boundary layer

When the boundary was initially drawn it was never snapped to the parcel boundaries. We now need to conduct reviews around the entire perimeter to ensure it aligns with features/boundaries.

Is there a quick way to run a boundary investigation?

Even symbolizing the overlap or identifying where it should be snapped to the boundary.


texts: Adding more than 200 codes to domain in geodatabase?

texts:
I want to add a domain in geodatabase as country code and country name (description) (more than 200 codes).

What is the  best way to do it? How can copy/write manually/python the list of coded domain into my geodatabase?


texts: Styling WMS layers in GeoServer with fields using "%" character

texts:
I have an Oracle Spatial 11g geodatabase connected to GeoServer 2.8 through Oracle plugin. Some feature classes have field names with % (percent) characters, e.g. NAME%EN or AREA%KM2.

The layer shows up correctly and all its attributes are passed with GetFeatureInfo request. However, I cannot use SLD styling that is based on those fields (e.g. creating thematic map based on AREA%KM2 field).

Does anyone have experience with tables using field names with % characters? They are accepted by Oracle naming standards so I hoped GeoServer will handle them too.


texts: When I calculate an area with area measure in QGis, it gives as a result values of area 44% more area than in reality, any idea why?


texts: Have one line and I want 1 from that one

texts: I have one line layer which is made up of multiple segments. I want each segment to be a new line feature so that I can select them individually using select features by area, as in the picture.


texts: Do jagged contours indicate more spatial variability in a data set?

texts:
A general question when working with contours, do jagged lines indicate more spatial variability in a data set?

A short background of my problem and objective: I would like to delineate neighborhoods based on historical home prices to create homogenous markets for actively listed homes. I have used the GA extension and chosen the covariance variable and circular model. My data set contains nearly 600 points and returned a regression function of 0.81 with a mean prediction of 1137 or so for home price. I am exporting my filled contour prediction map to GA Contours.


texts: Is there a downloadable data source for mapping out rivers with names for the whole world?

texts:
I'd like to be able to extract rivers with their names and the flow size (not direction) for the whole world. The scale I'll be drawing it at is about 2 centimeters per 10 real world km - nothing too precise! At this point I've looked at the hydrosheds dataset available for download and it seems great except for missing the river names...do you guys know of any other that is available for download?

I've also looked at the HYDRO1K (unable to download) and the Global Runoff World Rivers...I can't get a copy (asked already). Or - do you know how I could automate the naming of rivers after I have drawn them from hydrosheds? In other words - some online DB or API that I can use lat-long to ask for river name at this or that lat-long coordinate? I use Java and am expecting to have to code this out by using lat-long pairs...

Thanks!


texts: Dissolve everything except for one particular polygon

texts:
I want to dissolve these polygons (column FIRST_N03_003) except for FID 34 and 77.
Essentially, the 003 column is the state, and the 004 column is the city, but some cities are independent and do not belong in a state, which is why the 003 column is left blank for those two particular cities. My goal is to dissolve column 003 but also keep those two cities, to add population data later (so that those 2 cities and the other states are treated on an equal level). Would there be any way to do that?

This is what happens when I dissolve that column as is. I'm assuming the top row in the third column is those two cities, but they're not differentiated. The rest show the name of the states.


texts:
Can anybody tell me how to Download Qgis 2.2 new release version? In QGIS website they mentioned in download window QGIS 2.0 . if u checked released data it was 24 feb 2014.

So where i wll get QGIS 2.2 newly released version?


texts: Point sampling along a pole-wrapping coastline with PostGIS

texts:
I am working on a task that requires me to get sample points every 1000 km along coastlines, and have run into a problem with Antarctica. From what I can tell, it seems to be an issue with the use of geometry in the functions, when really I think geography should be used for this operation.

Using the function from this very similar question, I am able to produce a result that looks like this: .

As you can see, ST_AddMeasure() and ST_LocateAlong() don't seem to treat the geometry spherically, which results in many points that sit on the South Pole. A point was even added on the clip along the date line (left side). Per the documentation of these two functions, only geometry can be used.

The code used to generate the polygon and the points can be found here, but this is the SQL used to generate the points:

CREATE TABLE atest AS WITH line AS
  (SELECT
      id,
      ST_ExteriorRing((ST_Dump(geom)).geom) AS geom
    FROM line_sample_test),
linemeasure AS
    (SELECT
        ST_AddMeasure(line.geom, 0, (ST_Length(line.geom))::int) AS linem,
    generate_series(0, (ST_Length(line.geom))::int, 10) AS i
FROM line),

geometries AS (
    SELECT
        i,
        ST_LocateAlong(linem, i) AS geom
    FROM linemeasure)

SELECT
    * from geometries;

How can I generate points at every 1000 Km along this coastline?


texts: How to keep remain void islands in raster while mosaic?

texts:
I have raster images of DEM files

After running Mosaic To New Raster, I am getting images were all void spaces are getting filled(Highlighted with red lines in above image)

Below is the setting I have used

As a result, I am looking for a mosaic raster as it is.
I am using Arcgis Desktop 10.7.1


texts: Using a wildcard in an elif statement in pyscript...?

texts:
I’m trying label up a blank field with ‘yes’ or ‘no’ to the presence of names with either the word ‘Estuary’, ‘Loch’, or ‘Firth’ using * as a wildcard, in Field Calculator. I think I'm getting mixed up with SQL here... can anyone help?

The example of code I've used is:

arcpy.AddField_management(workingPath + "/GB_Catchment_merge.shp", "Offshore", "TEXT", "", "", 3)

expression = "LookUp( !MNCAT_NAME! )"
code_block = "def LookUp( a ):\n\

  if ( a == '*Estuary*' ):\n\

    return 'yes'\n\

  elif ( a == '*Loch*' ):\n\

    return 'yes'\n\

  elif ( a == '*Firth*.' ):\n\

    return 'yes'\n\

  else:\n\

    return 'no'"

arcpy.CalculateField_management(workingPath + "/GB_Catchment_merge.shp", "Offshore", expression, "PYTHON", code_block)


texts: Working with the FlowDir in the national hydography dataset

texts:
In the National Hydrography Dataset there's a FlowDir field.
According to https://www.townofnewipswich.org/sites/g/files/vyhlif4766/f/uploads/dataset_dictionary.pdf the two possible values for FlowDir are 1 for "with digitized" and 0 for "uninitialized".
According to https://nhd.usgs.gov/userGuide/Robohelpfiles/NHD_User_Guide/Interactive_Tutorials/Module_2/Navigating_the_Network/Navigating_the_Network.htm "Flow direction is from the ‘Start’ junction to the ‘End’ junction. If flow direction is incorrect,  NHDFlowline features may be ‘flipped’ to correct the flow direction.".
The shapefile has a sequence of coordinates. Maybe the first coordinate in this sequence is the "start"?
There's also an NHDFlowlineVAA attribute table. In this table there's fromnode and tonode columns. But...  when I try to open that file with DBF View I get the following error:

Exception executing: 'select * from NHDFlowlineVAA.dbf One or more errors occurred. - The provider could not determine the Decimal value. For example, the row was just created, the default for the Decimal column was not available, and the consumer had not yet set a new Decimal value.

OSGeo4W gave me 0.00000 for both of those for everything.
Any ideas?


texts: Using in-memory vector layer with QGIS processing / SEXTANTE

texts:
I'm trying to run the qgis:clip algorithm from console, but am getting an error when using an in-memory layer as the overlay parameter. Is this to be expected, or am I doing something wrong?

Code:

mem_layer = QgsVectorLayer("Polygon?crs=epsg:4326", "temp_layer", "memory")
if not mem_layer.isValid(): raise Exception("Failed to create memory layer")
mem_layer_provider = mem_layer.dataProvider()

clip_polygon = QgsFeature()
clip_polygon.setGeometry(QgsGeometry.fromRect(
    QgsRectangle(
        self.output_layer.extent().xMinimum() + 10,
        self.output_layer.extent().yMinimum() + 10,
        self.output_layer.extent().xMaximum() - 10,
        self.output_layer.extent().yMaximum() - 10
    )
))
mem_layer_provider.addFeatures([clip_polygon])
mem_layer.updateExtents()

output = self.output_layer_path + "2"
processing.runalg("qgis:clip", layer, mem_layer, output) # Fails

In the code above, self.output_layer and layer are vector layer objects (QgsVectorLayer - proper ones, loaded from shapefiles on disk), self.output_layer_path is a python string object with a path.

Here is the error I'm getting:

"C:/OSGEO4~1/apps/qgis/./python/plugins\processing\core\GeoAlgorithm.py", line 150, in
    execute self.processAlgorithm(progress)
File "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\algs\ftools\Clip.py", line 72,
    in processAlgorithm index = utils.createSpatialIndex(layerB)
File "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\algs\ftools\FToolsUtils.py",
    line 31, in createSpatialIndex features = QGisLayers.features(layer)
File "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\core\QGisLayers.py", line 211,
    in features return Features(layer)
File "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\core\QGisLayers.py", line 218,
    in __init__ self.iter = layer.getFeatures()
AttributeError: 'NoneType' object has no attribute 'getFeatures'

If I change my processing call to the following, it runs without error:

processing.runalg("qgis:clip", layer, self.output_layer, output) # Runs fine

Also, if of any help, this is the failing algorithm as it is logged in processing_qgis.log:

processing.runalg("qgis:clip","C:/path/to/shapefile.shp|layerid=0|subset=CONTINENT =
    'Europe'","Polygon?crs=epsg:4326","C:/path/to/output")


texts:
I'm looking for Poland zip-code data (polygons). I have found data in csv on the http://www.geonames.org/ but there are only point data.

Do you know any free source containing this data?


texts:
I have a route/path in my map that I need to draw a "corridor" on. Basically, the corridor are parallel lines on the route with a given distance in nautical miles.

The red line is the route, the blue line is the corridor.

I have the idea of making a radius with the given distance on each waypoint in the map and from there I would draw a line from the edges of the circles. So I was wondering if there is any official approach or calculation to this?


texts:
I'm using QGIS 1.8.0 and have some shapefiles with a "User Defined CRS" which I don't recognize. Is anyone able to recognize it?

Here is an image of the settings I see on QGIS.


texts:
I get XY coordinates of point features delivered to me in PDF format, this means that to bring them into ArcGIS I manually have to enter them first into an excel spreadsheet.

Is it possible to copy the coordinates straight from pdf?


texts: Exporting Images Linked with Survey data, from a GPS to GIS format

texts:
We were setting up a workflow system, wherein we attach Survey images to our collected point data (Leica GPS). The same can be viewed in Software like Leica Geo office, but a GIS export of the files (export shapefile option) doesn`t remember the Data - image link.

I have linked the photo to the survey points and later exported it to a shapefile from Leica Geo office 8.0, after importing the raw data.

But the shapefile does not recognise the point-image link.

Any expert advice on the best workflow for preserving the image link till GIS output?


texts: How to provide GDAL algorithms that require multiple layer input in the graphical modeler?

texts:
Some GDAL algorithms such as gdalbuildvrt or gdalmerge require the selection of input layers.
How to make an entry in a graphical model from a list provided by the processing flow?


texts: WFS getFeature request no longer working

texts:
I am trying to do a GetFeature with this parameters:
POST request:

Url: https://inspire.navarra.es/services/CP/wfs
Body:

<GetFeature
    xmlns="http://www.opengis.net/wfs" service="WFS" version="1.1.0" outputFormat="application/json"
    xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd">
    <Query typeName="CP:CadastralParcel" srsName="EPSG:25830">
        <Filter
            xmlns="http://www.opengis.net/ogc">
            <PropertyIsEqualTo matchCase="false">
                <PropertyName>nationalCadastralReference</PropertyName>
                <Literal>167090077</Literal>
            </PropertyIsEqualTo>
        </Filter>
    </Query>
</GetFeature>

But I'm getting this error:

Error occurred getting features
Duplicate mappingName or targetElement
across FeatureTypeMapping instances detected. They have to be unique,
or app-schema doesn't know which one to get.
Please check your mapping
file(s) with mappingName or targetElement of:
http://inspire.ec.europa.eu/schemas/cp/4.0:CadastralParcel

What is wrong?
It was working one week ago.


texts: R geographically weighted regression with negative binomial distribution

texts: I'm trying to find a R package that allows me to do geographically weighted negative binomial regression. I'm aware there's 'spgwr' and 'glm.nb' package. But is there a package that combines the two? If not, is it possible to change one of the package to meet my aim, for example adding negative binomial distribution to 'spgwr'?


texts: Extracting outer boundary of polygon (current polygon merged to one row but layer has holes creating small polygon needing removal) using ArcGIS Pro

texts:
I had to create boundary of raster layer, so I used raster to polygon to get the boundary. But the boundary has many holes based on the original raster which I want to get rid of as I only need the outer-boundary or extent for my analysis.
I have considered creating a new polygon and snapping the edge to this layer but I think there is a simpler way to do this?
Image below - referring to purple boundary


texts: One shapefile for African continent containing counties/regions?

texts: I use R to do my geo stuff with two limitations. First, I'm pretty new to all the geo stuff, and second I got lost while searching for a shapefile. I'm looking for a shapefile that contains all the gadm level 3 (possibly level 2 would suffice too) for all African countries. I could theoretically, download all the shapefiles from Global Administrative Areas and combine them into one masterfile (that would be another question). However, I was wondering if such a shapefile already exists... And if so, where?


texts: Where are the elements of an ArcGIS SpatialReference as a string representation defined?

texts:
This seems absurd, but I cannot find any documentation on what the various elements mean in the string representation of an ArcGIS SpatialReference object.  The first part is obviously the Well Known Text representation of the projection but there are 7 parameters, semi-colon separated, that follow it.  The last (7th) is typically "IsHighPrecision" indicating that the high precision grid should be used.

The string can be created using the object's exportToString() method and a SpatialReference can be created from a string using the object's loadFromString() method.  I see examples of using this string to create the reference or used in place of an explicit object (e.g. in CreateFeatureclass examples), but nowhere can I find a "dictionary" of the elements and their default values.

Can someone please point me towards the documentation on this?


texts:
I'm trying to copy a shapefile and give it a new name in QGIS python console, I am trying to use the QgsVectorFileWriter class' writeAsVectorFormat method

    writeAsVectorFormat(...)
    QgsVectorFileWriter.writeAsVectorFormat(QgsVectorLayer, QString,
    QString, QgsCoordinateReferenceSystem, QString driverName="ESRI
    Shapefile", bool onlySelected=False, QString errorMessage=None,
    QStringList datasourceOptions=QStringList(), QStringList
    layerOptions=QStringList(), bool skipAttributeCreation=False,
    QString newFilename=None, QgsVectorFileWriter.SymbologyExport
    symbologyExport=QgsVectorFileWriter.NoSymbology, float
    symbologyScale=1, QgsRectangle filterExtent=None) ->
    QgsVectorFileWriter.WriterError

I am unable to determine what the QString arguments (arguments 2 & 3) are supposed to be..I get that they are strings, but strings of what?? Is there a resource I could go to that would give more help than this? I am just trying to take a shapefile, copy it and give it a new name, is this the best way?


texts: Delete feature selection using Python (ArcGIS 10)

texts:
I'm having some issues with a Python script for ArcGIS 10. The idea is to delete points from a single shapefile according to the values of two columns ("Tracks" and "Tracks_1").
Therefore, I want to:

Select features who share a certain ID
Unselect some of the features depending on their values
Delete the remaining selection from the shapefile
Proceed with the next ID

.

import arcpy
env.overwriteOutput = True
centroids_thinned="O:\\centroids_thinned.shp"
arcpy.MakeFeatureLayer_management(centroids_thinned, "centroids_thinned_layer")
rows=arcpy.SearchCursor(centroids_thinned)
for row in rows:
    ID_centroids = row.getValue("FID_centro")
    ID_buffer = row.getValue("FID_cent_1")
    Tracks_centroids = row.getValue("Tracks")
    Tracks_buffer = row.getValue("Tracks_1")
    if ID_buffer == ID_centroids:
        arcpy.SelectLayerByAttribute_management ("centroids_thinned_layer", "NEW_SELECTION", "FID_centro = " + str(ID_buffer))
        arcpy.SelectLayerByAttribute_management ("centroids_thinned_layer", "REMOVE_FROM_SELECTION","FID_centro = FID_cent_1")
        arcpy.SelectLayerByAttribute_management ("centroids_thinned_layer", "REMOVE_FROM_SELECTION","Tracks_1 > " +  str(Tracks_centroids))
        count = int(arcpy.GetCount_management("centroids_thinned_layer").getOutput(0))
        if count > 0:
                arcpy.DeleteFeatures_management("centroids_thinned_layer")
                print str(count) + " Attributes Deleted."

The selection process works fine but for some reason those selected features are not deleted from the attribute table. How can i solve this problem?


texts:
I want to convert two double (double longitude_wgs84_radians, double altitude_wgs84_meters) which I get from a function.

The coordinates looks like

latitude_wgs84_radians: 0.924593
longitude_wgs84_radians: 0.176054

How can I convert them into UTM in C++? A code example or documentation would be nice.

I tried the example code from GeographicLib  1.49

// Example of using the GeographicLib::UTMUPS class
#include <iostream>
#include <iomanip>
#include <exception>
#include <string>
#include <GeographicLib/UTMUPS.hpp>
using namespace std;
using namespace GeographicLib;
int main() {
  try {
    // See also example-GeoCoords.cpp
    {
      // Sample forward calculation
      double lat = 33.3, lon = 44.4; // Baghdad
      int zone;
      bool northp;
      double x, y;
      UTMUPS::Forward(lat, lon, zone, northp, x, y);
      string zonestr = UTMUPS::EncodeZone(zone, northp);
      cout << fixed << setprecision(2)
           << zonestr << " " << x << " " << y << "\n";
    }
    {
      // Sample reverse calculation
      string zonestr = "38n";
      int zone;
      bool northp;
      UTMUPS::DecodeZone(zonestr, zone, northp);
      double x = 444e3, y = 3688e3;
      double lat, lon;
      UTMUPS::Reverse(zone, northp, x, y, lat, lon);
      cout << lat << " " << lon << "\n";
    }
  }
  catch (const exception& e) {
    cerr << "Caught exception: " << e.what() << "\n";
    return 1;
  }
}

but if Im trying to compile it, I get 3. errors.

undefined reference to `GeographicLib::UTMUPS::EncodeZone(int, bool, bool)'
undefined reference to `GeographicLib::UTMUPS::Forward(double, double, int&, bool&, double&, double&, double&, double&, int, bool)'
collect2: error: ld returned 1 exit status


texts: How to find the minimum pixel value within a region of interest?

texts: I am working on a Python script for dark object subtraction, a radiometric correction method for remotely sensed data.To accomplish this, I need to find the minimum band pixel value of a waterbody within a region of interest (ROI) (screenshot).  My initial thought was to create a polygon within a Landsat scene and use Zonal Statistics as Table to find the minimum value.  However, I cannot easily call this minimum value later on in the script.  Another thought was to convert the band to a numpy array and find the min value (excluding 0), although I am uncertain how to subset the array to my ROI.  How can I find the minimum value using a region of interest so that I can easily call this value later on in the script?


texts: Copying all attribute data within table and pasting into same table to create duplicates?

texts:
I have a spatial layer with an attribute table containing 55 rows.

How do I copy all attribute data within this table and paste it into the same table to create duplicates?

I have editor turned on and I can select the attribute table and click copy but I don't know how to paste into new rows in the same table.


texts: Local GeoJSON not showing on Leaflet map

texts:
I'm trying to load a local GeoJSON to a leaflet map. I'm new to leaflet so I was just following this example: https://www.igismap.com/add-load-geojson-file-point-polyline-polygon-map-leaflet-js/
So in the page html, I added this:

And in the overview.js file, I added this:
// loading GeoJSON file - Here my html and usa_adm.geojson file resides in same folder
$.getJSON("/maps/area.geojson",function(data){
// L.geoJson function is used to parse geojson file and load on to map
L.geoJSON(data,{style: myStyle}).addTo(map);
});

When I used the Brackets to see the live preview of the page, the features in the geoJSON file were shown. But when I accessed it from the web, the polygon was not shown. What could be the issue here?
Update:
Error found in browser console: