texts: Apply cloud mask to Landsat Imagery in Google Earth Engine Python API and output NAs instead of zeros

texts:
I am trying to mask clouds using the python API in Google Earth Engine, but the masked values are exporting as zeros instead of NAs. I am basically following the code here Apply cloud mask to Landsat Imagery in Google Earth Engine Python API
I have run the code in the GEE browser, and that actually exports the mask values as NAs, but not sure how to get the python API to do so.


texts: Is it possible to specify interpolation method for zoom functions in Styled Layer Descriptor (SLD)?

texts: To visualize images on Geoserver, I use xml file wrote in Styled Layer Descriptor (http://www.opengeospatial.org/standards/sld). I would like to specify what interpolation method would be used for zoom functions. Can I write it in xml SLD ?


texts: Defining fields and adding features to the output layer in a processing algorithm in QGIS

texts:
I have two input vector layers (one point and one line layer). I loop through every feature in the input-line-layer and calculate some things (with the help of the point layer).
What I am trying to do now is to generate an output line layer for which I want to have the following six fields:
1. QgsField('id', QVariant.String)
2. QgsField('vto', QVariant.Double)
3. QgsField('vfrom', QVariant.Double)
4. QgsField('priority', QVariant.Int)
5. QgsField('node1', QVariant.String)
6. QgsField('node2', QVariant.String)

Based on those defined fields, I want to add some features to the output layer,
but I do not know how to define the fields (of the output layer) and how to add features to the output layer in QGIS3.
Note: I am using the processing algorithm template for QGIS3 from Underdark: see https://anitagraser.com/2018/03/25/processing-script-template-for-qgis3/
Based on her template, she writes this:
    def processAlgorithm(self, parameters, context, feedback):
        inEdges = self.parameterAsSource(parameters,
                       self.INPUT_VECTOR_LAYER_EDGES, context)
        inNodes = self.parameterAsSource(parameters,
                       self.INPUT_VECTOR_LAYER_NODES, context)
        (sink, dest_id) = self.parameterAsSink(parameters,
                       self.OUTPUT_VECTOR_LAYER_MERGED, context,
                       inEdges.fields(), inEdges.wkbType(), inEdges.sourceCrs())

In the line which starts with (sink, dest_id), underdark has declared an outputlayer that takes the same fields as the inputlayer. But how can I declare an outputlayer with my customized fields (which are not ident with the fields from the inputlayer)?
EDIT:
As suggested from root676, I altered my code a little bit:
# 1. First of all, I defined a variable which holds the definition of my fields
outFields = QgsFields()

# 2. Then, I defined the fields
outFields.append(QgsField('id', QVariant.String))
outFields.append(QgsField('vto', QVariant.Double))
outFields.append(QgsField('vfrom', QVariant.Double))
outFields.append(QgsField('priority', QVariant.Int))
outFields.append(QgsField('node1', QVariant.String))
outFields.append(QgsField('node2', QVariant.String))

# 3. Then I created the output sink with the previously defined fields
(sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT_VECTOR_LAYER_MERGED, context,
                                           outFields, inEdges.wkbType(), inEdges.sourceCrs())

# 4. I created a new feature that will be added to the outputlayer after I set its attribute values.
newEdge = QgsFeature()

# 5. The ERROR happens when I try to set an attribute value like this:
newEdge['id'] = edgeIdTrimmed

The error I am getting is as followed:

Traceback (most recent call last):
File "", line 143, in processAlgorithm
KeyError: 'id'

EDIT: I solved the problem by initializing the new QgsFeature with my defined Fields variable. See below:
newEdge = QgsFeature(outFields)


texts: Assigning widths to string fields using writeOGR() from rgdal?

texts:
I am currently working with MapInfo files in the .tab format. I am using R (3.1.0 for Windows) to read and write files in this format, however, I am encountering an issue with writeOGR() from the rgdal package. When I write a MapInfo file in the .tab format, all the string fields are being assigned a width of 254 characters. This is very problematic for larger datasets with numerous string fields, as it dramatically increases the size of the output .dat file.

Below is an example:

require(rgdal)
shape <- readOGR(dsn = directory, layer = "Up")
writeOGR(obj = shape, dsn = directory, layer = "Up2", driver = "MapInfo File")

Here's an excerpt from the original Up.tab file:

!table
!version 300
!charset WindowsLatin1

Definition Table
  Type NATIVE Charset "WindowsLatin1"
  Fields 52
    ID Integer Index 1 ;
    unid_planejamento Char (4) ;
    codigo2 Char (4) ;
    descricao Char (100) ;

And one from my modified Up2.tab file:

!table
!version 300
!charset Neutral

Definition Table
  Type NATIVE Charset "Neutral"
  Fields 52
    ID Integer ;
    unid_planejamento Char (254) ;
    codigo2 Char (254) ;
    descricao Char (254) ;

Note that all the string fields have now been assigned a width of 254 characters. I should also note that the .dat file has gone from 43KB to 93KB.

I understand there are some creation issues when writing MapInfo files, however, it is not clear to me whether or not this can be resolved by somehow controlling the string field width.

I have searched for ways to do so using dataset_options and layer_options, but have not been successful.

Is there a way of controlling the string field width using writeOGR()? And if not, is there a way of using R to write a new .dat file (from my understanding, this is a dBase III DBF file) with the appropriate string widths?


texts: Why doesn't anything happen in using QgsRasterCalculator on a list?

texts:
I have a list of raster layers that I want to do a calculation on each layer with QgsRasterCalculator. My code is:

lstLyr = self.readRasters()

    for layer in lstLyr:
        registry = QgsMapLayerRegistry.instance()
        lyr = registry.mapLayersByName(layer)
        namlyr = lyr[0].name()

        output = "/LorDataPlugin/Data/Temps/Normlized_%s.tif" % str(namlyr)

        fileName = layer
        fileInfo = QFileInfo(fileName)
        baseName = fileInfo.baseName()
        lyr2 = QgsRasterLayer(fileName, baseName)

        entries = []

        ras1 = QgsRasterCalculatorEntry()
        ras1.ref = namlyr + '@1'
        ras1.raster = lyr2
        ras1.bandNumber = 1
        entries.append(ras1)

        formula = ras1.ref + '/' + str(self.sumValue(lyr2))

        calc = QgsRasterCalculator(formula, output, 'GTiff', lyr2.extent(), lyr2.width(), lyr2.height(), entries)

        calc.processCalculation()

But when run it without any error message, nothing happens. what is error in this code?


texts:
I am using TileStache to serve my own vector tiles using a postgis datasource. First I re-projected the shapefiles I had using ArcMap ( reprojected to Web Mercator (Auxilary Sphere)) and loaded them in Postgres using shp2pgsql (SRID:900913).

I am able to generate the Tiles using this config file:

{
  "cache":
  {
    "name": "test",
    "path":"/tmp/stache",
    "umask": "0000"

  },
  "layers":
  {
    "tpl":
    {
        "allowed origin":"*",
        "projection": "spherical mercator",
        "provider":
        {
            "class": "TileStache.Goodies.VecTiles:Provider",
            "kwargs":{
                "dbinfo":
                {
                    "host":"localhost",
                    "user":"postgres",
                    "password":"postgres",
                    "database":"tpl"
                },
                "queries":
                {
                    "7":"SELECT geom AS __geometry__,name, priority FROM isl_roads_sm WHERE priority IN (5,4) -- zoom 7+",
                    "8":"SELECT geom AS __geometry__,name, priority FROM isl_roads_sm WHERE priority IN (5,4)",
                    "9":"SELECT geom AS __geometry__,name, priority FROM isl_roads_sm WHERE priority IN (5,4)",
                    "10":"SELECT geom AS __geometry__,name, priority FROM isl_roads_sm WHERE priority IN (5,4)",
                    "11":"SELECT geom AS __geometry__,name, priority FROM isl_roads_sm WHERE priority IN (5,4,3,2)",
                    "12":"SELECT geom AS __geometry__,name, priority FROM isl_roads_sm WHERE priority IN (5,4,3,2)",
                    "13":"SELECT geom AS __geometry__,name, priority FROM isl_roads_sm WHERE priority IN (5,4,3,2)",
                    "14":"SELECT geom AS __geometry__,name, priority FROM isl_roads_sm  -- zoom 14+"
                }
            }
        },
        "preview":{"ext":"json"}
    }
  }
}

I am rendering the tiles generated using Leaflets TileLayer.GeoJSON. The problem I have is that the tiles dont overlay with the basemap, there is pretty large gap between my tiles and the features on the basemap.

When I use WGS84 as the projection in the config file (using data in WGS84) I get empty tiles, no other errors.

How can I get the tiles to overlay correctly? Why does WGS84 not work?


texts: Creating different diameter circles round points using QGIS

texts:
I want to map a tree survey. I am thinking that I will plot positional points of the trees and then add other data like tree type etc. to the attribute table later.
Is there a way I can have an attribute that describes a radius or diameter, whichever is easier, of a circle round the point?
In fact I need two circles, one to describe the canopy radius/diameter and another to do the same for the root protection area.
I am using QGIS 1.8.0


texts: Link QGIS to Anaconda Python Distribution: Spyder, iPython Notebook, pip, scikit-learn

texts:
I have Anaconda as my primary Python distro, and want to route Python related QGIS scripts through it.  With Anaconda, I  have Spyder, can easily install packages with pip, and can launch a iPython Notebook web-browser for both screen shares and storing code I intend to execute on my data before, during, and after manipulations in QGIS (which calls on packages such as scikit-learn).

I can not find any useful anywhere on how to link the two.

I'm using Windows 8.1, but if this can only be done in Linux, I can try to set it up on my desktop, which is running Mint 13. I just need to have a webcam and mic, found on my laptop but not desktop at the moment, for an interactive demo and interview on Tues.


texts: Using "Ruled Based Render" to display labels without the marker with python in QGIS

texts:
I'm writing a script to load and manipulate some shapefiles.
Unlike my previous question I now need to display labels only on certain zoom levels. My usual result is that the previous formatting gets over written after a rule has been set.

Outside of the rule set up I can display a point (without a marker) like this:

rendererP = layer.rendererV2()
symbolP = QgsFillSymbolV2.createSimple({'size':'0'})
rendererP.setSymbol(symbolP)
label = QgsPalLayerSettings()
label.readFromLayer(layer)
label.enabled = True
label.fieldName = 'Variable'
label.setDataDefinedProperty(QgsPalLayerSettings.Size,True,True,'6','') label.placement= QgsPalLayerSettings.OverPoint
label.writeToLayer(layer)
QgsMapLayerRegistry.instance().addMapLayers([layer])

I'm using the following to set up the display on certain zoom levels:

symbol = QgsSymbolV2.defaultSymbol(layer.geometryType())
renderer = QgsRuleBasedRendererV2(symbol)
root_rule = renderer.rootRule()
## copy the rule to edit
rule = root_rule.children()[0].clone()
rule.setScaleMinDenom(1000000)
rule.setScaleMaxDenom(4000000)
root_rule.appendChild(rule)
## delete the default rule
root_rule.removeChildAt(0)
## apply the renderer to the layer
layer.setRendererV2(renderer)

But the moment I attempt to apply the zoom level rules to this layer the markers returns.
Can anyone give me some guidance?


texts: Why OpenLayers v6 can not show layers when they secured by GeoServer 2.17.2 GeoFence extension?

texts:
I have two workspaces in GeoServer and two different users in GeoServer. For each user I created a rule in GeoServer by GeoFence Extension, each user must be see and access to only own layers are inside own workspace, that allowed by GeoFence rule, however if I set GeoFence like below images OpenLayers can show layers but each user can see content of another user's worksapce also.

I'd expect user 5f30ee36... can see only content of h0BUjvHBeM workspace, whereas it can see "Emeni" layers that is belong to j0BV2EjWBs workspace for user 5f5cd8...
However I can fix this problem by uncheck "Use GeoServer roles to geo authorization" and each user able to see only own layers are inside own workspaces, But my problem will happen after this step, "OpenLayers can not show layers!" but I can get link of layers! name of them and even see them in my browser (outside OpenLayers)"
My axios code is here
sysinfo: Os:Debian 10 x64, Openlayers V6, GeoServer 2.17.2, GeoFence 2.17.x,
my geoserver log


texts: How to normalize choropleth maps of census information?

texts:
I'm trying to improve on the census-based maps on this site, which shows census values at 5 scales from Collection Districts through to State.

Currently the maps show whole numbers (eg total population per polygon) which is not optimal. How should they be normalized?

Option 1: % of Total

I was advised to normalize by % of the total. This results in legend ranges such as "0.0004% to 0.0006% of the Total Population" which (IMO) make the legend hard to decipher - it's not an intuitive fraction, and is hard to put into perspective.

Option 2: Population Density

Another option could be to normalize based on area. Since the datasets are stored in meters, this also results in values like "0.000000019 to 0.000000035 persons per square meter".

I could flip this to show "area / person" - but this wouldn't make sense for datasets such as housing. And would it be offensive when applied to datasets such as religion (1km2 per X person)?

Is there a standard methodology to convert either of these values into something more meaningful? What do other people do on their census maps?

Thanks,
Steve


texts: Where does QGIS store info on a file's layers?

texts: I am using QGIS 2.16 and keep my GIS data (e.g. shapefiles) and my QGIS project files all on Dropbox because I work in two locations. When I open up a file made in the first location at the second location, new layers I was working on are missing, even though dropbox tells me my .QGS file and the shapefiles of the new layers have been updated.


texts: Geoserver envelope function not working in YSLD

texts:
I'm trying to use envelope function like this i my ysld:

- line:
    uom: metre
    geometry: ${envelope(geometry)}

My ysld won't validate with error:

java.lang.IllegalArgumentException: Invalid cql expression 'envelope(geometry)'

Functions: convexHull, boundary or even octagonalenvelope found on https://docs.geotools.org/latest/userguide/library/main/function_list.html works well.
envelope and octagonalenvelope definitions are very similar but still envelope doesn't work in ysld. I'm using Geoserver 2.16.2

envelope(geom): returns result
geom (Geometry) required

result (Geometry)

octagonalEnvelope(geom): returns result
geom (Geometry) required

result (Geometry)

Relevant log file part exceeded character body limit. Log file part available here

Caused by: org.geotools.filter.text.cql2.CQLException: Encountered "geometry" at line 1, column 10.
Was expecting one of:
    "-" ...
    <INTEGER_LITERAL> ...
    <FLOATING_LITERAL> ...
     Parsing : envelope(geometry).
    at org.geotools.filter.text.ecql.ECQLCompiler.compileExpression(ECQLCompiler.java:99)
    at org.geotools.filter.text.commons.CompilerUtil.parseExpression(CompilerUtil.java:172)
    at org.geotools.filter.text.ecql.ECQL.toExpression(ECQL.java:144)
    at org.geotools.filter.text.ecql.ECQL.toExpression(ECQL.java:125)
    at org.geotools.renderer.style.ExpressionExtractor.splitCqlExpressions(ExpressionExtractor.java:111)
    ... 140 more


texts:
Given a line (street, railroad, river...) I want to simulate a moviment inside this line.
Example:
See this picture:

Say I'm in "Nagqu" and want to move to "Kunlunqial".
In a time (t=0) I'm in Nagqu. By flowing the time, my point moves a little to northeast (lets say a velocity of 40MPH). In this time (t=1) I want a point exactly inside the railroad.
I know how to trace a straight line from a point to anohter, but I need to stay in the railroad (or line) no matter where it goes to. I need to find a point in that line until I reach my destination.
In other words: find a point inside a line given a direction (a line have only 2 directions).


texts: ArcPy: Selecting feature class by attribute then copy selection to new feature class

texts:
Using ArcPy I want to make to do a selection on a feature class by field to select everything 0.1 or greater, then copy the selection to a new feature class.
However, the output of my following code is exactly the same as the input and includes values under 0.1.
input_fc = r"C:\Geoprocessing\Model.gdb\input_fc"
output_fc = r"C:\Geoprocessing\Model.gdb\output_fc"

arcpy.SelectLayerByAttribute_management(in_layer_or_view=input_fc, selection_type="NEW_SELECTION", where_clause="Field1 >= 0.1")
arcpy.CopyFeatures_management(input_fc, output_fc)

is there a step I am missing?


texts: Merge several columns to one in attribute table, maybe coalesce with condition?

texts:
In my attribute table, I have 9 columns providing information on the same thing (group ID). One row does only have one group ID. I.e. one row has a group ID in one of these 9 columns, in the other columns it says 0. I need to have one column with the group ID information (an 8-10 digit number) for all rows. I have tried coalesce function, but it is (in it's pure form) not the right function for this, as each column has a number (i.e. 0 if not a group ID). Thus, the coalesce function provides me with the exact same column as the first stated column.

Is there any way to use the coalesce function but with a condition?
E.g. something like:

Coalesce ("Group_ID" if > 0, "Group_ID_1" if > 0, "Group_ID_2" if > 0, "Group_ID_3" if > 0 )

Or, is there any other way to do this?


texts: Sending requests including headers via PyQGIS by using QGIS Proxysettings

texts:
I am working on a plugin, where sending requests to an external server is part of the process. When using no proxy (like from my personal computer), this works just fine:
isochrone_headers = {"accept":"application/x-zip-compressed"}
isochrone_request = urllib.request.Request(isochrone_url, headers=isochrone_headers)
isochrone_response = urllib.request.urlopen(isochrone_request)

However, I have now tried this from my companies network where we are using a proxy. It is setup within QGIS settings and works just fine to browse plugin repository or use external basemaps. But my plugin is not working any longer.
So I was trying this (rather old) solution like
networkAccessManager = QgsNetworkAccessManager.instance()
req = QNetworkRequest(QUrl(isochrone_url))
r = networkAccessManager.get(req, headers={"accept":"application/x-zip-compressed"})

but I never get a valid response file from the OpenTripPlanner server. Unfortunately I really do need to send the headers to get my desired response file.
How can I send requests including headers by using the QGIS 3 proxy settings? I guess if thats not possible, using the proxysettings from operating system would be an acceptable workaround.


texts: i.maxlik cannot read i.class output signature file

texts:
I'm trying to perform a supervised classification using i.class. I'm then using i.maxlik, but this is returning the following error:

GRASS 6.4.3svn (India_01):~ > i.maxlik group=test_group@test subgroup=test_group sigfile=test_sigfiles1 class=test_results1
ERROR: Unable to read signature file <test_sigfiles1>

I'm not sure what other information I should provide. Can anyone suggest the possible reasons for this error? I've tried using i.cluster to generate an unsupervised signature file, and i.maxlik reads this perfectly well.


texts:
I just opened the QGIS and this appeared. What does it mean? Is it serious problems?

warning:/Applications/QGIS3.6.app/Contents/MacOS/../Resources/python/qgis/utils.py:672: DeprecationWarning: the imp module is deprecated in favour of importlib; see the module's documentation for alternative uses
          mod = _builtin_import(name, globals, locals, fromlist, level)
warning:/Applications/QGIS3.6.app/Contents/MacOS/../Resources/python/PyQt5/uic/objcreator.py:152: DeprecationWarning: 'U' mode is deprecated
          plugin = open(filename, 'rU')
warning:/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/psycopg2/init.py:144: UserWarning:

The psycopg2 wheel package will be renamed from release 2.8; in order to keep installing from binary please use "pip install psycopg2-binary" instead. For details see:


texts: Using C# and ArcObjects, what is the recommended way to convert a raster file into a polygon.  I need this for calling the Near tool so that I can move a station which may be inland to the nearest water feature.  I have googled and really not found anything useful yet.  There is a lot for taking a polygon to a Raster but not the other way around.


texts: How to Identify/Select features from WMS?

texts:
I'm trying to extract some features from the following WMS:

http://servicios.idee.es/wms-inspire/hidrografia

Namely, one layer called "Red hidrografica", which contains the traces of rivers and streams for all of Spain.

I've realized that by using the "Identify feature" option I'm able to select one particular river, copy the feature and paste it to a vector layer (even though the original is ostensible a WMS/raster layer).

This works, but it can be frustrating as the feature identification only seems to work when I click exactly on the river.

Any slight deviation results in nothing being identified, so I've had to spend a few minutes clicking at random around each feature trying to select them.

If I try to select a wider area, it fails to identify the only feature in the layer that falls within this area.

Is there a way to have this "Select by area" do what I need it to do?

Or some way to increase the size of the "Click-Identification-box-area"?


texts: Seeking open vector format for long term GIS data storage in off-line repositories?

texts:
Many field sciences require data management plans for their research data that specify how the data will be submitted for archival storage. These data come from various devices and are often in proprietary formats that require processing. For example, Trimble GPS devices produce SSF format files but these are not suitable for longterm storage because they are proprietary and require a Trimble software license.

If you worked for a library or the US Forest Service what vector format would you consider the most sustainable for people to use when they contribute data to your archive?

While not an open format Shapefile is a defacto standard and widely readable (in practice this is what people are using).
To be clear, these are not archives that are served online they're stored, perhaps burned onto archival DVD.

Relevant Links

GDAL - Geospatial Data Abstraction Library
Sustainability of
Digital Formats: Planning for Library of Congress Collections
OGC standards working group
FOSS Open Standards/Comparison of
File Formats


texts: Using geodatabase for processing 7200 polylines?

texts:
I have 7200 lines and I want to divide each of them every 100 meters and take the coordinates of points from them. For that task should I create a geodatabase and put them in there?

I am totally newbie with geodatabases and I don't have much time. Any good tutorial? Because then I will have to do that with python also.

Has the geodatabase the capacity for that task?

I have ArcGIS 10.2 for Desktop.


texts: QGIS Georeferencer adds padding to image after save

texts:
I am attempting to georeference an aerial image using satellite imagery. I'm successful in getting an output, however I notice there is an issue where QGIS's georeferencer is adding padding to all the edges around the image after the save. It looks like it's related to the gdal_warp command (as I separated out the commands from the georeferencer); I'm using a polynomial 2 with cubic spline resampling method with a total of 6 GCPs. The target SRS is EPSG:4326. Attached is an illustration of what I mean:

It's important to have this not require any human cropping as the human will not be perfect in generating the proper output.

Here's the zoomed out view to give the full picture with the georeferenced on on the left and the original on the right again:


texts: Merge overlapping rasters with ignoring no data values QGIS3

texts: I thought this would be extremely basic task but obviously it's not. In QGIS3 GDAL Merge (GUI), I am trying to merge two rasters which overlap and each of them has a different distribution of nodata values. I want the rasters to be mosaiced together such that nodata values are ignored if one the rasters does not have nodata value on that spot. In ArcGIS it is that simple, why not QGIS?


texts: County area calculated from NLCD rasters is too large (with R)

texts:
I'm trying to calculate landcover repartition for each US county.
I have obtained NLCD for the Apache county using the FedData package (devtools version) and I'm using county shapefiles from the Census bureau.
The problem is that I get an area that is much larger than the official one and the one indicated in my shapefile, namely 51,000km^2 instead of 29,000km^2 officially. There must be something I don't understand about the raster object.
The following describes the code used and the method used to calculate. The county data can be downloaded here:
https://www2.census.gov/geo/tiger/TIGER2016/COUNTY/
The following code assumes the county shapefile is saved and unzipped.

Get and read the data

#devtools::install_github("ropensci/FedData")
library(FedData)
library(rgdal)
library(dplyr)

#Get Apache polygone
counties<- readOGR('tl_2016_us_county/tl_2016_us_county.shp')
apache <- subset(counties,counties$GEOID=="04001")

# Get NCLD data
nlcd_data <- get_nlcd(apache,
                      year = 2011,
                      label = "Apache",
                      force.redo = TRUE)

nlcd_data #inspect the object, can see that number of cells is around 57 million

I have then extracted the values of the raster and put them into a frequency table. From there I calculate the resulting area. Since the NLCD data is 30m resolution, I multiply the number of cell of each category by 900 and divide the result by 1 million to obtain the area in Square Kilometer.

The calculated total area is too large.
# Calculating the landcover repartition in County
landcover<-data.frame(x2011 = values(nlcd_data)) #Number of rows corresponds to number of cells
landcover_freq<-table(landcover)

df_landcover <- as.data.frame(landcover_freq)

res <- df_landcover %>%
  mutate(area_type_sqm = Freq*900,
         area_type_km=area_type_sqm/1e6,
         area_sqkm = sum(area_type_km))%>%
  group_by(landcover)%>%
  mutate(pc_land =round(100*area_type_km/area_sqkm,1))

head(arrange(res,desc(pc_land)))

# A tibble: 6 x 6
# Groups:   landcover [6]
  landcover     Freq area_type_sqm area_type_km area_sqkm pc_land
  <fct>        <int>         <dbl>        <dbl>     <dbl>   <dbl>
1 52        33455938   30110344200      30110.     51107.    58.9
2 42        16073820   14466438000      14466.     51107.    28.3
3 71         5999412    5399470800       5399.     51107.    10.6
4 31          488652     439786800        440.     51107.     0.9
5 21          362722     326449800        326.     51107.     0.6
6 22           95545      85990500         86.0    51107.     0.2

## Total area calculated from raster is 51,107 square km

apache_area <- as.data.frame(apache) %>%
  mutate(AREA=(as.numeric(ALAND)+as.numeric(AWATER))/1e6) %>%
  select(AREA)

apache_area$AREA
29055.47 #Official area of apache county (in square km)

Visual inspection of the shapefile and the raster:

The difference doesn't seem large enough to justify the difference
apache <- spTransform(apache,proj4string(nlcd_data))
plot(nlcd_data)
plot(apache,add=TRUE)

[![Apache polygone on NCLD raster][1]][1]


texts: R when I reproject coastscoarse to Robinson using spTransform I get an unwanted horizontal extra line

texts:
When I run this code
library(rworldmap)
data("coastsCoarse")
plot(coastsCoarse)
RobinsonProj <-  "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
crsLoc <- CRS(RobinsonProj)
coastsCoarse <- sp::spTransform(coastsCoarse, crsLoc)
plot(coastsCoarse)

coastsCoarse is a SpatialLinesDataFrame. I get one plot that looks like the coasts of the continents. The second plot, after projecting into Robinson returns a horizontal line across the top of the plot of the continents. How do I get rid of that line in the projected version? Extra credit for an explanation of how to get rid of Antarctica.


texts: How to change the raster format using QGIS translate?

texts:
I tried to change a raster format from .tiff to .asc using Qgis's 'translate'.

As a result, I get 3 files: a PRJ, XML and a file without extension just a plain file...

What could it be happening? I used the same tool for a bunch of raster previously and it worked...


texts: Shortest path from points in one layer to points in another layer via roads in QGIS

texts:
I'm looking for a way to calculate the shortest path from addresses to subway stations via a network of roads.

I have the Road Graph plugin, however, the addresses and subway stations are point layers that do not overlap with vertices in the roads layer.

In QGIS, is there a way to calculate the shortest distance from the address/subway station points to the roads, create vertices at those locations, then calculate the shortest distance via the roads network?

black squares - addresses
white circles - subway stops
gray lines - roads


texts:
I'd like to find a function or plugin for QGIS that displays vector data in a dynamic report. For example, I have a point vector layer "sample plots" that consists of 4 fields: ID, plot-number, comments and grade.
I'd like to have a report or form that is linked to the dbf-file of that layer and displays the data dynamically in report. So if I add another point, it would show up on the report.

Does such a function or tool exist in QGIS?


texts: How to add a Mapinfo-style text halo in QGIS print composer?

texts:
Is there any way to imitate the Mapinfo halo (of something similar) in the Print Composer. I'm trying to put some text (copyright information) over the maps and I want to make it so that it is easier to read the text over the map.

I've tried using the "Background" option for the label, but this blocks out the entire label area, and not just what is printed (the text is being auto generated with a custom function).

I'm currently experimenting with text-shadow using the below HTML snippet, but the highlighting is too diffuse to make any clear difference

[% printCopyRights()%]


texts: Mapping data from coarse grid cells onto high resolution DEM (in QGIS)

texts:
I'm trying to figure out a way to map water depth data from a coarse triangular mesh cell onto a high resolution (HR) raster (in .asc format). My coarse model and the raster are 2D.

I created a very coarse (<-- intended) mesh for an urban flooding simulation by polygonizing a high resolution digital elevation model, smoothening the edges (using QGIS) and then meshing it (using gmsh).

Then I did the flooding simulation and now I have some results for the mesh in .vtu format. For each cell I have a water elevation and a bed elevation value
(--> substract the latter from the former to get water depth). The water elevation is of course higher than the bed elevation of a cell.

As mentioned the grid is coarse, I have big cells, for example let's say a cell is 100 m in characteristic length. The HR raster is very fine.

In a next step I'd like to get some water depth results (= water elevation - bed elevation). I already have results just for the coarse grid. But now I'd like to take my water elevation value from a coarse cell and the underlying bed elevation data from the HR model and compute the water depth this way.

The following sketch might elaborate my undertaking a little better

I have a couple of problems doing this:

Is this actually the right forum for this kind of task :D? Because I don't know if QGIS is capable of doing this, or if I only need something like ParaView for this task.
The raster calculator might come in handy I thought.
If QGIS fits for me, how do I get the .vtu data into QGIS? This would be nice because I could create the water depth for all cells but if it's not possible I'd just replicate a couple of cells as triangles in QGIS and give them the water elevation data by hand.

(I've seen a couple of questions concerning QGIS and VTK that go in the same direction and I will check them out now.)

If all of that should work out, how can I actually compute the water depth? My first thought would be to use the raster calculator to compute the difference between water elevation and bed elevation and create a new raster that has the resulting water depth as information on each raster point.

I've now managed to create a polygon with cell IDs and corresponding water elevations in QGIS. I used a gmsh plugin for QGIS and extracted the elevation data with ParaView into a CSV-File. I then created a table join within QGIs for the Cell IDs and their respective water elevation levels. Now I only need to use the raster calculator. Does anyone have an idea how I could compute the wanted water depth?


texts: Fill color not showing appropriately in folium?

texts:
The following code should have shown different colored circles in a map. But my result was circles with no color at all.

import folium
import pandas

data=pandas.read_csv("Volcanoes.txt")
lat = list(data["LAT"])
lon = list(data["LON"])
elev = list(data["ELEV"])

def color_producer(elevation):
    if elevation < 1000:
        return 'green'
    elif 1000 <= elevation < 3000:
        return 'orange'
    else:
        return 'red'
map = folium.Map(location=[38.58, -99.09], zoom_start=6, tiles="Mapbox Bright")
fg = folium.FeatureGroup(name="My Map")
for lt, ln, el in zip(lat, lon, elev):
    fg.add_child(folium.CircleMarker(location=[lt, ln], radius = 6, popup=str(el)+" m",
    fill_color=color_producer(el), color = 'grey', fill_opacity=0.7))
map.add_child(fg)
map.save("Map1.html")

Can you show me where did I did wrong?


texts: Is there a way to set rotation on inset maps to match data driven maps

texts:
I have been playing with the new map book functionality in ArcGIS 10 to re-create maps I have made in the past. The new Data Driven maps with page queries is a pretty powerful tool.

Is there a way to match the rotation of inset map, secondary data frame, with the main layer being used for the data driven maps? I found it is possible to match the extent in the 'Data Frame' tab of the 'Data Frame Properties', is this also the place to find the rotation option?

If this is something that you need to code, is it possible to build it into the mxd and to have it launch on refresh?

Not sure I understand how Python is supposed to work within a mxd.


texts:
How do you manage the following projects (i.e. what GIS software / platform would you use)?

 An extensive Footpath Construction & Renewal Program
 Improved drainage systems at x Street, y Street, z Avenue
 Traffic Management along A Avenue, B Road and C Road
 Reconstruction of the following roads: K Street, S Court, D Road


texts: How to synchronize a map (Google Map), Mapillary and Street View?

texts:
I've an HTML / Javascript page with three divs:

Mapillary images on the left
map (Google Map) at center
Google Street View images on the right

I'm using Google Maps API and Mapillary JS.

Here you are the code (NOTE: you need Mapillary key  and Google API Key to execute it ....):

    <meta charset="utf-8">
    <title>Street Image Compare View side-by-side</title>

    <link href='https://unpkg.com/mapillary-js@2.5.0/dist/mapillary.min.css' rel='stylesheet'/>

When I click on the arrows on the Street View image (on the right ...), I can see the pegman that moves on the map and this is what I want.

I'd like to see also the Mapillary images move in synchronize manner, but they don't move now.

Also I'd like to click on the arrows on Mapillary image (on the left ...) and see that pegman moves (at the moment I can move the Mapillary images but pegman doesn't move ...), and I'd like to see the Street View images move in synchronize manner.

Examples or suggestions?


texts: Scripting a conditional formula in field calculator QGIS

texts:
I'm trying to write down a script in the field calculator but I'm getting to nowhere.

In one field I have runoff data, and in two other ones I have the first quartile of the runoff data (Q1) and the third quartile (Q3) of the runoff data. Basically what I'm trying to achieve is to replicate my excel formula inside the field calculator:

=IF(N3<=Z3;"A";IF(N3>=AA3;"B";IF(AND(N3>Z3;N3<AA3);"C";0)))

with:

N3= runoff data

Z3= Q1

AA3= Q3

Inside the attribute table, I have the three fields (runoff, Q1 and Q3) labelled in the same way.
As I'm fairly new to the scripting method, I don't know the syntax I need to respect in my field calculator formula.


texts: ArcPy Script to update LYR data source not identifying all layers in folder?

texts:
I've put together a script that for the most part works as intended - to walk through a series of folders, find SDE layer files and update their data source to direct connects if the layers support workspacepath. However I've noticed that in some cases it is not finding specific layer files within folders. Here is the script:

import arcpy
import os
from arcpy import env

env.overwriteOutput = True
TheLayerDir = r"G:\GIS\Layers"
NewDataPath = r"G:\GIS\Data\vector_direct.sde"

for dirname, dirnames, filenames in os.walk(TheLayerDir):

    for filename in filenames:
        newFilename =  os.path.join(dirname, filename)
        desc = arcpy.Describe(os.path.join(dirname, filename))

    if newFilename[-4:] == '.lyr':
            lyrFile = arcpy.mapping.Layer(newFilename)

            for lyr in arcpy.mapping.ListLayers(lyrFile):
                if not lyr.isGroupLayer:

                    if lyr.supports("WORKSPACEPATH") == True:

                        if "sde" in desc.dataElement.catalogPath and lyr.workspacePath != NewDataPath:
                            print "------------------------------"
                            print "Layer: "        + lyr.name
                            print "Old Data Source: "   + lyr.workspacePath
                            lyr.replaceDataSource(NewDataPath,"SDE_WORKSPACE")
                            print "New Data Source: "   + lyr.workspacePath
                            lyrFile.save()

                        elif "sde" in desc.dataElement.catalogPath and lyr.workspacePath ==  NewDataPath:
                            print "------------------------------"
                            print "Layer: " + lyr.name + " has already been updated with path, skipping..."

                        else:
                            print "------------------------------"
                            print "Layer: " + lyr.name + " is not a SDE Layer"
                    else:

                        if "sde" in desc.dataElement.catalogPath:
                            print "------------------------------"
                            print "Layer: "        + lyr.name
                            print "Layer does not support workspacePath"

print "Script Complete! - That's all Folks!"

If I for example substitude the following code as follows:

 if newFilename[-4:] == '.lyr':

with

if extension.lower() == ".lyr":

then the script will successfully identify all layer files within a folder, however without declaring 'newFilename'.


texts: How to call Leaflet FullScreen plug-in from code

texts:
I'm currently using Leaflet along with leaflet.fullscreen to provide a full screen map.

I've got it working nicely with the button but I can't work out how to "call" the function manually from my JavaScript code.

I've tried: map.fullscreenControl.toogleFullScreen(); but get the following error and a "broken" fullscreen:

TypeError: this.invalidateSize is not a function
this.invalidateSize();


texts: Embed Esri map in MS Access Report or Form

texts:
I have forms and reports in an MS Access front-end that reference tables in an Oracle ESRI.ST_GEOMETRY back end. The forms allow users to tabularly navigate through construction project attributes in a more intuitive way than can be done in ArcGIS Desktop. The forms and reports have hyperlinks that open a popup form that contains project details. Among other things, I would like the popups to contain a map of the project.

Sure, I could embed or link to a Google map or Bing map, but this would only show the general location. I would like the map to contain several GIS layers, such as road centerline, road right-of-way, sewers, sidewalks, bridges, etc. This is where an embedded GIS map would work well.

The common link could be the project number in the record in the form, and a project number in a polygon layer in the database.

Is it possible to embed some sort of Esri-based map (desktop, web, or otherwise) into an MS Access report or form? Or as a Plan B, somehow hyperlink to an Esri map?


texts:
I am trying to use the following CQL filter in Geoserver for extracting all polygons (layer1) intersecting a line (layer2):

http://XXX.XX.X.XX:8080/geoserver/workspace/wfs?service=wfs&version=1.1.0&request=GetFeature&typeName=workspace%3Alayer1&srs=EPSG:4326&CQL_FILTER=INTERSECTS(geom, collectGeometries(querySingle('workspace%3Alayer2','the_geom','fid IN ('layer2.1')'))&outputFormat=application/json

Layer 2 only has 1 feature with fid = 1, Layer 1 has hundreds of polygons.

Possibly there are other ways of retrieving the same information, but why the above is not working?

It says Could not parse CQL filter list

I am using Geoserver 2.15.1


texts: ArcGIS Closest Facility with anisotropic round trips?

texts:
I'm working on a network analysis of an urban archaeological site where I'm running a closest facility analysis on some 2000 structures to specific features within the city. I've segmented the road layer and calculated travel times for each segment based on Tobler's Hiking Function, which I then enter as the time cost to calculate the quickest routes from each house to the closest feature (in terms of time). The problem I've run into is how to calculate these as round trips. I can't simply double the total time of each route because the time values are anisotropic based on slope and Tobler. I've tried running the analysis backwards (from Facility to Incident, instead of Incident to Facility), but this often results in different matches and routes, again because of anisotropy. I've also tried summing the From-To and To-From time values for the segments and using those as a kind of "round trip" cost, but this again results in different Incident-Facility matches from the original Incident-Facility analysis using the anisotropic times. I'm at a loss as to how to do this. I thought if there was a way to produce a list of network layer segments for each route, I could simply sum the opposite time values (TF vs FT), but I can't find a way to produce this, plus I'm sure that each route would consist of both TF and FT values for each segment. Does anyone have a solution or have I reach a limit of current ArcGIS capabilities?
UPDATE: I've found the "Copy Traversed Source Features)" tool that will export a table with all of the segment values, but the process of recalculating the accumulated inverse time values is daunting, to say the least. I have to link the segment IDs from each route to the original network layer, then somehow note if each segment's time value was TF or FT (via separate fields with 0 or 1 values) and replace it with its inverse. Even worse, sometimes segments are joined in the middle (not at a junction) and the time values have to be proportionally shortened.
I feel like there must be a way to write a script or something to automate this, rather than attempt to do it manually, but I have no idea how to do that.


texts:
I am working on a project for which I would like to plot some data on a map and I was told about QGIS by a friend that uses it. It seems pretty complex and loaded with features and I have been trying to figure it out through tutorials. But I really just need a simple plot, so I wonder if someone could help me out, or at least tell me where to look. Or maybe QGIS is not even the program I should use for this and there is something else out there.

I have biological data from about 2000 samples taken throughout Spain. Essentially what we have and want to display is each sample placed into one of 2 categories. I would like to have a map divided into the 17 autonomous communities and then I want to make a kind of heatmap of all the autonomous communities based on the percent of samples in each category. For example, if 80% of samples from Catalunya are in category 1 and 20% are in category 2, then on the map, Catalunya would be a lightish red (where 100% category 1 would be red, 100% category 2 would be blue and 50/50 would be white).

From one of the tutorials I got this data:

http://maps.cga.harvard.edu/qgis/data/older_data/workshop_data_2014-11.zip

which has this map of Spain:

workshop_data_2014-11/part_2/ESP_adm4.shp

The map is divided into cities and in the Attribute Table it has the autonomous community under NAME_1. So is there a way for me to change the map so it only divides into the different values of NAME_1? And then how can I introduce my data and color the communities accordingly?

And also, where do you obtain map files like this, in case I need to use a different one?

There are a number of other things I could plot on the map with our data, but I think if I can get this basic one down, I'll be able to figure out how to do some of the other things.


texts: left join returning one more row than expected

texts:
I have a table of single polygons (table rast_poly) derived from a raster and a table of polygons that represents buffered points (table pt_buff). The raster polygon geometry has been reduced by using an ST_Intersection() operation with the buffered points, so the geom is clipped to those boundaries. Each polygon from the raster contains an integer value that corresponds to a crop code. I want to sum the areas for each crop type within each buffered point polygon. I have to join the cropcodes table to the first query on 'value' and 'country' because the value-category pairs are unique to each country. I am almost there with the following code:

select t3.id, t3.country, t3.val, t3.km2::numeric(10,2), t4.category
from (select t2.id, t1.val, sum(st_area(t1.geom::geography)/1000000) as km2, t1.country
      from rast_poly t1, pt_buff t2
      where st_intersects(t1.geom, t2.geom) and t2.id = '17DYR07'
      group by t2.id, val, t1.country
      ) t3
left join cropcodes t4 on t3.val = t4.value and t3.country = t4.country;

Here is the output:

I've checked the results against both the dumped polygons from the raster, and the original raster itself. Everything is correct except the two duplicate rows at the bottom of the output.

How can I fix this?


texts: Cannot reuse output feature class in a custom geoprocess using ArcObjects

texts:
I'm developing an arcobjects 10.4.1 custom GeoProcess (fusion shift) which has to be able to run within a model builder.
The GeoProcess seems to run smoothly, but the problem is it needs to be outputting a resultant feature class that needs to be reused afterwards.
The resultant feature class path is being generated correctly and the feature class is there at the end, but the node doesn't turn green nor shades, so it cannot be reused.

This is the code i'm using to statically define this resultant feature class field:

            outputParameter = new GPParameterClass();
            outputParameter.DataType = new GPFeatureLayerTypeClass();
            outputParameter.Value = new GPFeatureLayerClass();
            outputParameter.ParameterType = esriGPParameterType.esriGPParameterTypeDerived;
            outputParameter.DisplayName = "Resultant feature class";
            outputParameter.Name = "ResultantFeatureClass";
            outputParameter.Direction = esriGPParameterDirection.esriGPParameterDirectionOutput;
            outputParameter.AddDependency("ShiftedWorkspace");
            outputParameter.AddDependency("ShiftedFeatureClassName");

            IGPFeatureSchema outSchema = new GPFeatureSchemaClass();
            outSchema.FieldsRule = esriGPSchemaFieldsType.esriGPSchemaFieldsNone;
            outSchema.FeatureTypeRule = esriGPSchemaFeatureType.esriGPSchemaFeatureFirstDependency;

            IGPSchema schema = outSchema as IGPSchema;
            schema.GenerateOutputCatalogPath = true;
            outputParameter.Schema = (IGPSchema)outSchema;

            parameters.Add(outputParameter);

And this piece of code is executed once in the update parameters method:

            // Output parameter feature class
            IGPParameter3 gpParamOutputMasterLayer = (IGPParameter3)m_Parameters.Element[++parameterCount];
            IGPParameterEdit3 gpParamOutputMasterLayerEdit = (IGPParameterEdit3)gpParamOutputMasterLayer;
            if (!gpValue.IsEmpty()) gpParamOutputMasterLayerEdit.Value = m_GPUtilities.UnpackGPValue(gpValue);

What am I missing? is there something else I need to define?


texts: How to update Shapefile's Properties on PostGIS?

texts:
Context: I'm building a WebApp with OpenLayers where users interact with some vector data (MultiLineStrings from Shapefiles) loaded from a GeoServer using WFS.

Basically they rate (like/dislike) some Features and choose the reasons why they like or dislike from a list. This information is easily stored in the Features using OpenLayers.

feature.set('Rating', 'Like');

For every user session I store all rated features as a GeoJSON into a PostGIS table. It works.

Problem: I'd like to also update a copy of the original Shapefile as a PostGIS table where each feature would contain a summary of all ratings. Let's say 87 users liked the feature with id = 42. I'd like my feat.42 to have a LIKES property with value 87. And I'd like to directly update those LIKES in the DATABASE. Something like:

UPDATE schema.shapefile_table SET Likes = newValue WHERE id = 42

But, when I use shp2pgsql & psql to insert my shapefile to PostGIS I get something like:

I suspect the original features' properties aren't in my table but just their geographical information. So,

Am I right?
Is there any way to import a Shapefile to PostGIS including its properties? (using shp2pgsql | psql or something as well non-interactive)

In any case, is there a better scenario to solve what I'm trying to do?

Edit: My SQL file generated by:

shp2pgsql -G -I myshapefile.shp

looks like:

SET CLIENT_ENCODING TO UTF8;
SET STANDARD_CONFORMING_STRINGS TO ON;
BEGIN;
CREATE TABLE "myshapefile" (gid serial,
"id" int4,
"length" float8,
"geog" geography(MULTILINESTRING,4326));
ALTER TABLE "myshapefile" ADD PRIMARY KEY (gid);
INSERT INTO "myshapefile" ("id","length",geog) VALUES ('0','184.68','0105000020E61000000100000001020000000C000000944DB982ED76184155BE671D0C8F514152499DE0CC76184172DEFFF60B8F5141DA5889C9B8761841EE0912E40B8F5141BDFDB9A8987618419641B5C60B8F51417638BA4A75761841AA622AA90B8F51414FB0FF6A4D761841245F098E0B8F5141A48B4D3B05761841C47B0E620B8F51412A8E03AFB175184149BDA7330B8F514182C64C525175184154AA44FD0A8F5141E083D77EED74184155DD23CB0A8F5141041C4215547418410DE02D840A8F5141A702EE690B741841F644D74C0A8F5141');
INSERT INTO "myshapefile" ("id","length",geog) VALUES ('1','83.37','0105000020E610000001000000010200000006000000D2FD9CB278721841BA8443AC3A8F51415C1E6B669D72184177853ED53A8F51417C28D142F57218410CCB9F023B8F5141D670914B4173184146ED7E233B8F5141651BB86381731841268FA7433B8F5141931804E6C573184174081C773B8F5141');
...
CREATE INDEX ON "myshapefile" USING GIST ("geog");
COMMIT;
ANALYZE "myshapefile";


texts: What is future of Silverlight with ArcGIS Platform?

texts:
At the September 2011 Build Conference, Microsoft said a lot about Windows 8 and the Future of JavaScript over Silverlight. No Plugins wanted!! JavaScript with html5 ftw!

What is the future of Silverlight with the ArcGIS Platform (ArcGIS for Server, ArcGIS Online, etc) in 2014 and beyond?


texts: Converting longitude and latitude coordinates to Indian Grid system?

texts:
How to convert lat/long coordinates to Indian Grid system?

I've googled for a couple of days but didnt find the exact algorithm for this.

and i got a link Indian Grid System.

but i didnt get a picture how the zones(there are 9 zones) are divided among india and how to convert the coordinates?


texts:
I have a thematic map layer Rekryterare in QGIS. The layer is defined as identifiable. The question is how can I make the thematic categories of this map layer unidentifiable when they are invisible on the map? (see the attached picture).

Right now I can click the invisible points on the map and read the feature attributes and this is not what I want, I would like to have the invisible points unselectable. Anyone know how I can solve this problem?


texts: QGIS atlas export as image with different page sizes / same dpi

texts:
I've made an atlas in QGIS where different pages have the same map extent, same dpi, but different page sizes / resolutions. When I export as image (png), however, QGIS overrides the different page sizes with one fixed image size, changing the resolution instead.

Simplified example:

Atlas page 1 is 4096x2048 pixels @ 96 dpi, scale 1:2257
Atlas page 2 is same extent, 2048x1024 @ 96 dpi, scale 1:4514

If I preview the atlas on screen, it's fine. And I can export page by page using Layout / Export as Image, as expected. But if I try Atlas / Export Atlas as Images, the Image Export options dialog that appears populates with the currently-being-previewed page's dpi and page width and page height, and then all pages are forced to that width and height in pixels. So if I was looking at atlas page 1 on screen, I get 2 exported images at 4096x2048 at different dpi. How do I instead get the dpi to be fixed and the image width and height to change?

The page width and height do change as expected if I export to PDF instead (where the image export options dialog doesn't muck with image size).

I'm doing this to generate pngs at different scale levels that I then convert to one .mbtiles with gdal_translate (plus sqlite3 to combine into a single file). This is to preserve resolution-dependent rendering that creating a single png and then using just gdaladdo would not allow, and to avoid the instability and cross-tile-oundary labeling issues of the QTiles3 plugin.

[QGis Windows 3.6.2]


texts: Extracting points bordering shape extent

texts:
I'm using ArcGIS 9.31.

I have a set of points as elevation and in order to get a soft clip area for the TIN, I would like to extract only all the points bordering the extent of this shape, or to get a shape deliminating the extent of the set of points. I don't have any other data than this point shape. Is there a way to select them?


texts: Arcpy.Buffer_Analysis Does Not Create Feature Class

texts:
import arcpy

arcpy.env.workspace = r'c:\data'

fcs = str(raw_input('Please enter the layers you would like to buffer: '))

dist = str(raw_input('What would you like your buffer distance to be? '))

output = str(raw_input('Name your output: '))

for fc in fcs:
    if arcpy.Exists(fcs):
    arcpy.Buffer_analysis(inputs, output, dist)
if not arcpy.Exists(fcs):
    print "Feature class does not exist"

The above script executes but does not create the feature class I want. Where am I going wrong?


texts: Flow Accumulation Python or MATLAB Script

texts:
does anyone have any resources for scripting the D8 flow accumulation in Python or MATLAB after obtaining the flow direction values and filling in 'sinks' and 'flats'?

I am trying to understand the logic so that I can apply a similar algorithm to a land use raster.


texts: How to validate high precision GPS systems location

texts:
I am currently working with a high precision RTK-GPS system which claims to be able to reach centimetre level accuracy, I am looking for a way to validate that it is actually reaching that precision and doesn't produce systematic errors etc.

My initial thought was to look for places nearby that had very precisely defined coordinates to compare the system against but I am struggling to find any (I live in Cambridge, UK). Any other ideas for how to go about validating it?


texts: Collection properties are lost when exporting in google earth engine

texts:
I have a general confusion over the 'id' field and how it is set and change in GEE, but I faced another strange thing: The collection properties are lost when exporting in GEE! I tried below sample code:

var features = ee.FeatureCollection([
  ee.Feature(ee.Geometry.Point(30.41, 59.933), {name: 'Voronoi'}),
  ee.Feature(ee.Geometry.Point(-73.96, 40.781), {name: 'Thiessen'}),
  ee.Feature(ee.Geometry.Point(6.4806, 50.8012), {name: 'Dirichlet'})
]);
features = features.set({'test_property':'Test'});
Export.table.toAsset({
  collection: features,
  description:'vectorsToAssetExample',
  assetId: 'exampleAssetId'
});

I can see 'test_property' and can not see any 'id' field by print(features), but  when I import 'exampleAssetId', I have an 'id' field and instead of 'test_property' in the properties, I have 'system:asset_size' property.
How can I have my properties retained when exporting?


texts: Most efficient method to get national-level data displayed with Leaflet

texts:
I'm looking for the most efficient way to get a national or U.S. and Canada level dataset working with a JS mapping library. My problem is there's a myriad of technologies that seem like they can work. I'm just not sure if they do or which of them is best.

My use case: my company does monthly crime mapping for several large clients and they use the data to understand risk around their businesses. The dataset itself is not at all large. We're talking maybe 800-1000 points among the U.S. and Canada but it is proprietary and sensitive and for that reason I'm hesitant to use geoJSON. Given they are using it to understand risk around a location, they need to be able to see a certain level of zoom - which is why MBTiles/TileMill kind of falls short. I'm looking to do zoom levels 0-12 or 0-14. That generates a 10GB tile file. I can split up the tiles into like two or more portions, but I'm customizing maps for several companies and that seems like it would create a bunch of additional steps. Additionally, I'm looking for the points to have teasers/tooltips so they need to be modifiable.

So in summary I need -

A decent amount of zoom (thinking 0-12 or 0-14 level) for about 1,000 points, which should have tooltips/teasers, dispersed among the United States and Canada. The data is proprietary and thus it should be in a secured format.

I'm looking for suggestions on what technologies can help me do this. I accept there's going to be a cost, but I want to avoid anything too expensive (> $10k).


texts: Can I define a style for one layer, taking into account a defined relationship with another? (QGIS 2.8)

texts:
I have two sets of data. A large set of lines with a line_ref field, and a small data-only (no geometry) table with a matching line_ref field (let's say for argument's sake that the data-only table marks recently surveyed sections of pipeline). I'm talking a about a one-to-one relationship on this occasion.

I'm using QGIS 2.8

I can specify a table join, and refer to the joined data from the data-only table in styling the lines.

I can do fancy things with input forms based instead on setting up a relation in the project properties.

What I was wondering was whether I could style the lines based on the relation (rather than specifying a join). I'm prepared to be told that this isn't a sensible thing to do if that's the case.... but I'm thinking ahead to when I start to want to style things where there's a one-to-many relationship. Surely at this point table joins start to get more complex?


texts:
I'm using GeoServer REST API to create a layer from a shapefile. That shapefile needs to be updated on a regular basis. Since there is no replace function on the REST API, I'm trying to delete the layer, replace the files and then recreate the layer with the new shapefile. pretty much like the answer at Uploading shapefile to GeoServer using cURL and replace file?

The problem I´m facing is that i´m receiving a 403 Forbidden message when I try to delete the layer. Since I have no problems publishing and editing layers I guess it is a security problem but have no clue of where I must change such security features.

The cURL command I´m sending is:

curl -v -u user:pass -XDELETE -H http://localhost:8080/geoserver/rest/layers/LayerName.xml

Sucth request gives me:

http/1.1 403 Forbidden


texts: How would I build the syntax inside a request in the case I wanted to multiply 2 variables?

texts:
I would like add something to my previous question
(here: How to build QgsExpressions in python?)

How would I build the syntax inside this request in the case I wanted field_name to be also a variable pointing to a field name?

And in the case I wanted to multiply two numerical variables inside a request?

expression = QgsExpression('variablepointingtofield * variablepointingtolistitem ')


texts: How to get the closest point to the map on the shapefile in Python

texts: The red dots represents the points in the shape file. I want to find point B (closest point on the shapefile to point A). The B might not be one from the given set of points in the shapefile (as we are finding the closest one.) Any solution to solve thus using Python? I want the coordinates of B.


texts: Interpolating specific area on point cloud in QGIS

texts:
I am making a heightmap and am trying to interpolate a specific area on a point cloud in QGIS. I have found some ways of using GRASS to interpolate using the attribute table. The problem is that the attribute table is empty.

Is there any way to interpolate with GRASS by just using the z-coordinates, or maybe importing the z-coordinate to the attribute table?


texts:
Long time qgis mac user.  I find with qgis 3.4, it will not change fonts within the Layout text box, nor will it show ESRI North arrows when I choose that font in a text box.
When I deleted 3.4 and reloaded 3.2 Bonn all the fonts are working as before.
Has this been noted and been fixed?  I cant use 3.4 at this point.

QGIS 3 is a vast improvement over 2.18.   Congratulations to the development team.

Tanoos


texts: Can Arcpy Reconcile Initiate the Interactive Conflict Review?

texts:
I'm pretty new to python but trying to hide a lot of the nitty-gritty of creating a grandfathered version and reconciling from my users. I've gotten pretty far, but can't figure out a way for a toolbox script in ArcMap to kick off the interactive conflict review window for my users - is this possible using python and ArcMap?

I can't seem to find any info online about this. I'm using AGS 10.2, python 2.7.
Here is the window I'm trying to initiate (conflict resolution window):  http://help.arcgis.com/en/arcgisserver/10.0/help/arcgis_server_dotnet_help/index.html#//009300001606000000


texts: Getting the segment of a line, crossed by another with QGIS expression

texts:
I have two line layers in QGIS 3.16: line 1 (blue), line 2 (red). The white points are the vertexes (nodes) of the two lines. They have no common vertex, are not connected. How can I get the segment from line 1 that is crossed by line 2 using QGIS expressions? I would like to generate a line by geometry generator of the highlighted yellow segment (see screenshot).

I tried this as a geometry generated symbol layer on line 1: I first segmentize line 1 to get the segments, than I cross it with the line_2 on the other layer. I get a boolean output (true/false) and I am also able to create a crossing point where the two lines intersect. However, I tried in vain to address directly the segment that is crossed.
crosses (
   segments_to_lines (
      $geometry),
   geometry (
      get_feature(
        'line_2',
        'fid',
        1
      )
   )
)

In this simple case, I could get the segment with the expression (that's how I generated the yellow line): geometry_n (segments_to_lines ( $geometry), 3). However, I have to set the segment-number (3 in this case) manually, but I want to create this number autormatically as I have a lot of different cases and can't do it manually. How can this be done?
For the background of this question, see: Re-locating points randomly at certain distance in QGIS (Step 5)


texts: Merging/Dissolving features with same ID using Transformer in FME?

texts: I was looking to merge some features to improve the performance of some of my layers on my ArcGIS Server. I have a few very large polyline layers for our fiber network that are exported out of a CAD based program (Mapcom M4) every week and into a GeoDatabase. I did notice that some of the features have the same attributes and fiber ID's, but are separated into up to between 4-15 different lines. When I tested this out and manually dissolved the featueres into one based on fiber ID it drastically improved performance. This would be a significant expense of time to do every week, though, and I was looking to make it part of my Weekly FME Workbench. I just wasn't sure what  transformer would make the most sense to use to merge all the features with the same Fiber ID before writing them to the Geodatabase.


texts: How to customize locator file to geocode unusual address?

texts:
I am very new to geocoding. I am trying to create address locator and geocode a bunch of address using ArcCatalog 10.1, and I get problem for the following addresses:

Road Type is in Street Name: eg: Glenview Heights Lane
There is direction in street name: eg: East Street
I have an address of 75 Romfield Cir (cir represents circuit), somehow it doesn't recongnize it even though I have added cir, circuit as a street type. I am thinking it is probably conflicting with the street type circle (cause they all start with cir).


texts: Checking if two WGS84 lines intersect in JavaScript

texts:
I am trying to work out given two GeoJSON lines (WGS84), if they intersect (i.e. cross one another). Let me give some pseudocode to illustrate:

// Providing Line 1 = lon1, lat1 and Line 2 = lon2, lat2
function linesIntersect(lon1, lat1, lon2, lat2) {

    var intersectionCheck = mathsToDoGeodesicIntersectionChecking()
    // This could be an external library or another code sample

    return intersectionCheck; // Return a boolean with true or false

}

The code samples I have tried are to long winded to post in the question, but you can see in the following GitHub repo (open line-crosses.html in a browser and open the dev tools to see some tests run).

It appears that the popular library Turf.js does not have a line intersect function. I am keen to make sure that the lines be created geodesically and checked in the same fashion (i.e. accounting for the curvature of the Earth rather).

The operation seems comparatively simple in a projected coordinate system scenario. I have a script (based off of code for in the comments of this post) that converts the lines coordinates to Web Mercator and checks for intersection, however I am uncertain as to the issues this may cause.

The code provided by moveable-type appears to get very close but unfortunately seems to check entire line paths and not line segments you provide (i.e. will these two lines ever intersect at their current bearings). Similarly GeographicLib looks like it may be able to help but I am struggling to find the solution using the library.

I am not specifically interested in where they intersect but more if they do. It would also be highly beneficial if they could handle discontinuity, as in to say  to be able to handle a line that starts at a latitude of -179 and ends in at latitude 179 (see here for more explanation of what I mean).


texts: Openlayers overview map interferes with loading

texts:
I'm using Openlayers 3.16.0 and this bug has really stumped me.

Whenever I turn on my overview map layers start loading more slowly and after a while just not anymore at all. Like this:

Is this an OL bug that got fixed in a later version, or am I doing something wrong here maybe?

EDIT:

Tested with 3.18 and it's still doing it...


texts: Getting wikimapia embeded in Leaflet map?

texts:
I have a leaflet map where I want to get wikimapia POI and data inside my map so I could identify some Points of interests and informations which are in wikimapia ,I have a wikimapia kml that loads in google maps, but How to have it in leaflet is this possible?

  <title>
    POI
  </title>
  <meta http-equiv='content-Type' content='text/html; charset=UTF-8' />
  <link rel="stylesheet" href="http://k2scc.onlinewebshop.net/plunkers/leaflet/leaflet.css" />
  <link rel="stylesheet" href="http://k2scc.onlinewebshop.net/plunkers/src/leaflet-search.css" />

here my pluncker code:
http://plnkr.co/edit/FNUDToR2XmSO2rIaK9pw?p=preview


texts:
I am currently working with the Union tool of ArcObjects but I do not know how to properly pass the input_features to the union object. I´ve tried several methods but all throw an Automation Exception (0x80004005 - Unspecified Error). I could not find a single sample for the Union tool in Java.

   GeoProcessor geoProcessor = null;
   try {
        geoProcessor = new GeoProcessor();
        geoProcessor.setOverwriteOutput(true);
        // Get featureclass objects of both shapes
        IWorkspaceFactory workspaceFactory = new ShapefileWorkspaceFactory();
        IWorkspace workspace = workspaceFactory.openFromFile("C:\\test", 0);
        IFeatureWorkspace featureWorkspace = (IFeatureWorkspace) workspace;
        IFeatureClass fc1 = featureWorkspace.openFeatureClass("a.shp");
        IFeatureClass fc2 = featureWorkspace.openFeatureClass("b.shp");
        // Create IGPValueTableObject and and add featureclasses
        IGpValueTableObject gpvtObject = new GpValueTableObject();
        gpvtObject.setColumns(2);
        gpvtObject.setRow(0, fc1);
        gpvtObject.setValue(0, 1, 1);
        gpvtObject.setRow(1, fc2);
        gpvtObject.setValue(1, 1, 2);
        /*
         * Also tried this:
         * gpvtObject.addRow(fc1);
         * gpvtObject.addRow(fc2);
        */
        Union union = new Union();
        union.setInFeatures(gpvtObject);
        union.setOutFeatureClass("C:\\test\\union.shp");
        geoProcessor.execute(union, null);
    } catch (IOException e) {
        e.printStackTrace();
    }


texts:
I have one country layer in a shapefile consisting of a lot of attributes, but the result needs to be in WKT file/format. I'm using QGIS 3.12 . I'd tried GetWKT Plugin but it limited to one selected feature only.
How can I convert the whole shapefile to a WKT file? Any suggestion of what plugins I can use to convert the shapefile?


texts: Points shapefile resulting in MultiPoints table in PostGIS

texts:
I have a Points shapefile. When opened with QGIS, the type is 'Point' (right-click on the layer and check properties). When imported in PostGIS using pgAdmin III plugin 'PostGIS shapefile and DBF loader', the created table is of Multi Point type... I can't change the type afterwards and this brings issues like updating other Points tables geometry column.

Any idea on this issue?


texts: Image Statistics Import in Google Earth Engine?

texts:
I've wrote this code for Nitrogen Index. but it returns error. how can I solved this?

code link: https://code.earthengine.google.com/ff6954984a933d11580d2717ba59670e

error massage:

Layer 2: Layer error: Image.subtract, argument 'image2': Invalid type. Expected: Image. Actual: Dictionary.

> Map.centerObject(geometry);
>
> var hyperion = ee.ImageCollection('EO1/HYPERION')
> .filterBounds(geometry) .filterDate('2000-01-01','2018-01-01')
> .toBands() .clip(geometry);
>
> print(hyperion);
>
> var NDRE = hyperion.expression(   '(B137 - B31)/(B137 + B31)' ,{
> 'B137':hyperion.select('EO1H1660372010258110PC_SGS_01_B137') ,   'B31'
> : hyperion.select('EO1H1660372010258110PC_SGS_01_B031')
> }).rename('NDRE');
>
> var min = NDRE.reduceRegion({   reducer: ee.Reducer.min(),   geometry:
> geometry,   scale: 30,   maxPixels: 1e9 }); print(min);
>
>
> var max = NDRE.reduceRegion({   reducer: ee.Reducer.max(),   geometry:
> geometry,   scale: 30,   maxPixels: 1e9 }); print(max);
>
> var CCCI = hyperion.expression(   '(NDRE - min)/(max - min)' ,{
> 'NDRE':NDRE ,   'min' :
> NDRE.reduceRegion(ee.Reducer.min(),geometry,30),   'max' :
> NDRE.reduceRegion(ee.Reducer.max(),geometry,30) });
>
> Map.addLayer(NDRE);
>
> Map.addLayer(CCCI)


texts: Add county attribute to gSSURGO soils raster

texts:
I have a high-resolution (10m cell size) raster of soils data for a US state. Pixel values are the map unit key (mukey) from gSSURGO. Every soil type has a unique mukey, but a given mukey can occur in one or more counties within a state.

I am trying to merge data into this raster attribute table, but the data I'm trying to merge has unique values for each mukey in each county that mukey occurs in. Consequently, I cannot make the merge until I first add a county field into the raster attribute table.

My current approach to doing this is to raster::rasterize a polygon .shp file of the state's counties to the same resolution as the gSSURGO data. Since the cell sizes are SO small, however, this is taking in incredible amount of time. Even with 8 cores on a fast machine...

Is there another good way to do this?!

I prefer to do this with R, but can use ArcMap instead if that can do it more quickly.

Here is the R code I am using for the raster::rasterize approach in parallel. I found this approach here. With this approach I was able to process one county on each of my 8 cores in 1.5 hours. And I've got a whole country to do!

IL_soils_raster <- raster("MapunitRaster_IL_10m.tif")
us_county_polygons <- readOGR("cb_2016_us_county_500k/cb_2016_us_county_500k.shp")

IL_county_polygons <- us_county_polygons[us_county_polygons$STATEFP == 17,]
IL_county_polygons  <- spTransform(IL_county_polygons, CRS = crs(state_soils_raster))

cl <- makeCluster(detectCores())
clusterEvalQ(cl, {library(raster)})
clusterExport(cl, list("ST_counties", "st_soils"))
rParts <- parLapply(cl = cl,
                    X = 1:nrow(IL_county_polygons[,]),
                    fun = function(x) rasterize(IL_county_polygons[parts[[x]],], IL_soils_raster, 'COUNTYFP'))
stopCluster(cl)
ST_counties_raster <- do.call(merge, rParts)


texts: How do I populate one attributes table column based on a value from another column? Qgis

texts: I am trying to populate all data points with an elevation based on the fid. Basically, all fid's of 0 have an elevation of 1700, all fid's of 1 have an elevation of 1740 (there are 48 fid's) ... There are over 5000 data points to populate so doing it manually is not an option. Is there a way to building an expression in orfer to populate the elevation column automatically?


texts: Can anybody help me with Spatial Join in ArcGIS 10?

texts:
I have a doubt, I do not find what's happening. Screenshot:

Using the Join (Join data from another layer based on spatial location) I have done an Spatial Join by location, and it is said that the FIELD "Count_" must to count the number of elements (in my case Points) within it. But when I see its metadata, I find that there more "number of elements" within the polygon! So, my question: Why do I have more countered elements?

Thanks for all.

----EDIT----

Done only with a simple example of 1 polygon with some points nearby and inside, it works. 5 elements and countered 5 elements. But still don't know why counts more elements with the raw data.  :(


texts: Selecting polyline between two point features using ModelBuilder?

texts:
I am looking to select a segment of a polyline that is between two GPS points in ModelBuilder. The line segment is created by enabling a user to define his/hers begin end point, and then using the split line at point tool to split out the appropriate line segment.

I am struggling to find a way to select this newly defined segment once I get to this point. I would like to select and export this segment out so I can run further analysis for the user.

Does anyone have a solution that can be shared?


texts: Creating FeatureLayer from JSON data using ArcGIS API for JavaScript 4?

texts:
I am trying to create a FeatureLayer from a JSON object that will be fetched from another service. I am not looking for something complicated, just something that works.

var layer = new FeatureLayer({
    fields: [{
        name: "ObjectID",
        alias: "ObjectID",
        type: "oid"
    }, {
        name: "fltid",
        alias: "FltId",
        type: "string"
    }, {
        name: "name",
        alias: "Name",
        type: "string"
    }],
    objectIdField: "ObjectID",
    geometryType: "polyline",
    spatialReference: { wkid: 4326 },
    source: [{
            geometry: new Polyline({
                paths: [
                    [-100, 38],
                    [-77, 35]
                ]
            }),
            attributes: {
                ObjectID: 1,
                name: "KATL",
                fltid: "UAL1"
            }
        },
        {
            geometry: new Polyline({
                paths: [
                    [-77, 35],
                    [-120, 40]
                ]
            }),
            attributes: {
                ObjectID: 2,
                name: "KZBW",
                fltid: "SW999"
            }
        },
        {
            geometry: new Polyline({
                paths: [
                    [-120, 40],
                    [-100, 38]
                ]
            }),
            attributes: {
                ObjectID: 3,
                name: "WKRP",
                fltid: "Fever1"
            }
        }
    ],
});

I was expecting three lines to be drawn, but nothing is shown and the error console is empty.

I am using the version 4 API and there are very few examples and most of them broken.


texts: postgis insert query to store geometry and other data

texts:
I'm using Postgis functions st_dump and st_split to split a polygon based on line. Now, after splitting; output geometry stores in a database table as well. Here is the query:

insert into aoo(geom)
SELECT ST_AsTEXT ((ST_Dump(ST_Split(circle, line))).geom) As wkt
FROM (SELECT
    ST_GeomFromText('LINESTRING(73.33393335342407 33.82005181729661,73.33685159683228 33.80960466650892)',4326) As line,
    (select geom As circle from aoo where gid=71) ) As foo;

But, the issue is I want to store other data with these geometries as well. Like

insert into aoo(geom,name)

but I don't what is the way to do it.


texts: How can I get one neighbourhood name using address and/or postal code?

texts:
this is my first time here.

I just need to get a neighbourhood/locality name, to get it I have the address and postal code (all over Argentina). -no polygons needed-

This is to change a background image acording to its neighbourhood/locality, this does not need to be a 100% accurate about neighbourhood/locality limits


texts: How do you scan/quantify vegetation using satellite images?

texts:
I am relatively new to GIS (about 4-5 months in) I did not receive any formal education about this so I am completely clueless about how to go about measuring/quantifying vegetation using satellite images. Can somebody please help? Is there any specific kind of image file I need to do this, if there is, how/where can I get them?

Thanks!

Enrique Roa
Undergrad biology student


texts: How to get correct coordinates system of drone images?

texts:
I have a folder containing images taken by a drone of a construction site. I have to convert those images into GeoTIFF.

For GeoTIFF conversion I know gdal_translate is used to do this task, but the problem is how I get or know the coordinate system of images taken by drone, so that I can convert images into proper GeoTIFF.

Is there any procedure to follow and solve this problem using Python GDAL bindings?

I used pix4dmapper to geo-referenced drone images, but now i want to make my own tool to perform this task so thats why i am asking all this.

EXIF Information of image:


texts:
I have a code to identify current date and its set to End_period. How can I set 60 days before the End_period to the Start_period?

var End_period = ee.Date(new Date().getTime())
var Start_period = ee.Date('2020-01-18')


texts: QGIS does not read ESRI Geodatabase correctly

texts:
I have a geodatabase from the EPA that does not seem to me reading correctly with Quantum GIS. QGIS Version 1.8.0 (Lisboa) and GDAL/OGR 1.9.2

When I open the attribute table the first two columns have the same data. The first column is CERCLA_ID and the second is SITE_NAME. Both columns have the SITE_NAME value. As far as I can tell all else is correct.

I opened various files in the .gdb folder using a hex editor and in the a00000004.gdbtable file the data is stored in an XML format internally and the correct CERCLA_ID value is indeed being stored along with the SITE_NAME.

Is there something wrong with the file or something wrong with the program?

I also heard back from my EPA contact. He says they are using "attribute domains" on their end.


texts: How to add a box with longitude/latitude tick marks in QGIS

texts: I am using QGIS version2.18.2 to generate a map. I need a box to embrace the map, and the box needs to show tick marks from the corresponding longitude and latitude of the mapped area. I've searched the way of adding this kind of box, but has no solid solution. Hopefully I can get some guidances here. It shall look like the box below:


texts: Using Planet Labs API tile service with leaflet map

texts:
I'm using Planet Labs API tile service with leaflet map. Need to add it as a tile layer/WMS layer.
But the map stays the same when zoomed in or out.

My Map: https://drive.google.com/open?id=1jodnQB9exuxxt6trxirGmbdF3FkPX-nx

Here's my code.

var mymap = L.map('mapid').setView([51.505, -0.09], 13);
var wmsLink='https://tiles0.planet.com/basemaps/v1/planettiles/global_monthly_2016_05_mosaic/gmap/0/0/0.png?api_key={your_key}'
var wmsLayer4 =  L.tileLayer.wms(wmsLink , {layers: 'layer_name',format: 'image/png',transparent: true  }).addTo(mymap);

Referred from https://developers.planet.com/docs/api/tile-services/

Providing a working example would be great.


texts: How to setup default toolbars and menu item?

texts:
I need to customize the UI for a group of users at another office.  I want to turn off a bunch of toolbars and menus that are not needed.  I know how to save my configurations (Settings -> Customization...) to a INI file.

How do I setup the other users computers to load the INI at startup?

I am using Windows 7 and Qgis 2.0.1.


texts: Is there a way to select only a certain polygon for editing?

texts:
I have a spatialite data base with 4300 polygons contained within. I wish to be able to edit the polygons but once i have highlighted the Polygon from the attributes table i find it hard to edit the vertices as when i try to select the cross hairs (vertex) unless i am extremely careful i end up selecting another polygon?

Is there a way of selecting the polygon i wish to edit from within the attributes table and making only that polygon editable? I know i could export it, edit it and then merge it back into the Spatialite database but this is additional steps that i would rather avoid if possible.

Any help or suggestions would be much appreciated.


texts: Converting raster data to have the attached results using ArcGIS Desktop?

texts:
I have a huge 30x30m reclassified raster dataset (10 classes) (in the size of Portugal) and I would like to find the easiest way to convert it to kml, having a similar result to the map we bought attached below.

If you look in the map, you will see that they poligonize each cell of the raster into 9 different shapefiles (the 9 small folders in the left side) according to a range of value, like "WSP < 3.4", "3.5 < WSP < 4.1" and etc.

Also they could keep a streched colour to each cell.

Normally when I transform a raster to a polygon it becomes one shapefile and if I want to transform it to a kml that have different folders I have to extract by attribute each value of my shapefile and convert each one of this shapes to kml. After I need to go in each one on Google Earth and set a color manually, paying attention if it is in the colour scale or not.

How can I make it easily and right to have a result like this one?

I have a ArcGIS 10.4 with spatial analyst license.

Desired Results


texts:
The write. function allows one to export a JSON file, but it does not seem to work for a GeoJSON file.

This is the code i ran:

write(map,"Map Layer2.GeoJSON")

and this was the error produced:

Error in cat(list(...), file, sep, fill, labels, append) :
  argument 1 (type 'S4') cannot be handled by 'cat'


texts: Converting Seismic Shot Points into polyline for linear distance for each Tract using ArcGIS for Desktop?

texts:
I have Seismic shot points (receivers and source points) that I am trying to convert from Points into a polyline for linear distance for each Tract in ArcView 10.1.

What is the best way to do that?


texts: Walking distance with pgrouting for several locations

texts:
I'm intending to calculate the walking distance (not the linear distance, but the "real travelled distance") to several locations (e.g. social infrastructure). Thus I will have to create a routable network.

Most implementations I found are based on a single point only. This would be hell of work, so I hoped there is a possible solution to do the calculation at once for a set of multiple points.

I think the way to go would be PostgreSQL and pgrouting with dijkstra, but I'm not sure. My input layer is a point layer with the given locations, then I have OSM data which has to be made routable. What is the best way to achieve getting isochrone or isodistance (travelled distance) maps?

I also took a look at grapphopper, which seems to be great and easy to configure for an isochrone map of a single point, but I didn't find an implementation for a multipoint isochrones case.

Any hints?


texts: Copy & Paste Feature in same Layer but not editable seperately?

texts:
First of all, I'm new here and hope this question is accurate enough and isn't answered already somewhere. Here's the problem:
I have an existing polygon and need to create another one of it (within the same layer), with a small part of the border edited.
Normaly, I just copy & paste the feature within the same layer. Then I can edit the copied entry and update the attribute table.
But when I do this, somehow QGIS also edits the borders of the original polygon. The entries in the attribute table remain updated, but the polygons are the same again.
The only way this worked was, when I copiy & pasted the polygon into a new temporary scratch layer, edit the borders there and then copied it back into the original layer. But with this method, I have to edit all entries in the attribute table for the copied polygon, altough I only would have to edit one row.
I use QGIS 3.10.8 longterm release.


texts: Sextante Modeler Error "global name os is not defined"

texts: I want to generate a model whic reclassify radter values. So i added Saga reclassify grid values tool to Modeler and saved the model. When i run the model i got this error both on Ubuntu OSgeolive 6.0 and Win7. But this tool runs succesfully directly in Qgis. the error occur when i want to run in Modeler. I use Qgis 1.8.0 and Sextante 1.0.9. I need your help...


texts: How can I plot GPS coordinates on Google Map

texts:
What is the best way to plot GPS coordinates on a Google Map?

Are there any tools that can do that (written in python specifically)?


texts: Problem with protocol Microsoft Azure Blob in QGIS

texts:
I'm trying to upload an image in TIFF format into QGIS using Azure Blob protocol with the name of the BLOB and the storage account key but I have this error:

Layer Invalid: provider GDAL You can not open GDAL dataset/vsiaz/cog2/ 3LiES........
Lack AZURE_STORAGE_ACCOUNT + AZURE_STORAGE_ACCESS_KEY
or AZURE_STORAGE_CONNECTION_STRING settings frame layer Provider is invalid
(Provider: gdal, URL:/vsiaz/cog2/3LiES.............

I write the name of the storage account and the key1 of that storage account but I get an error

What could be the problem?


texts: Does ST_GEOMETRY support Linear Referencing?

texts:
I am wondering if the ArcGIS ST_GEOMETRY provides Linear Referencing functions, such as creating a point based on a given measurement (like Oracle Spatial)?

I couldn't find any function like that in the ArcGIS documentation but there are pretty functions in ArcPy.


texts: Explanation of HAND (height above nearest drainage) model procedure

texts:
i tried to follow this procedure Spatial Analysis - Calculate a cells elevation (height) above the nearest stream cell but i have some problem to understand everything. I mean the 2th answer in the link.
I don't understand the 3th and 4th points. I have flowdirection, stream network and dem. I would to understand what about flow elevation difference and what i need to do to get the hand model.


texts: Criteria for accessing Workspace and its subfolders using ArcPy?

texts:
I am wondering when I can access the workspace and its subfolders with Arcpy tools (using ArcGis 10.2).

I have a raster in following directory:
E:/folder/sub1/sub2/sub3/raster.tif

Now I set the working directory to sub2:

import arcpy
arcpy.env.workspace = "E:/folder/sub1/sub2"

For some tools, I can access and write files in subfolder sub3 by referring to the subfolder of the workspace:

# load raster in subfolder
r = arcpy.sa.Raster("sub3/raster.tif")
# load and write rasters in subfolder
arcpy.gp.Reclassify_sa("sub3/raster.tif","Value","1 3;2 2;3 1;4 4;NODATA 0", "sub3/outRaster1.tif","DATA")

This works all very well. If I use the .save() method, however, I get an error:

outR = Reclassify("sub3/raster.tif", "Value","1 3;2 2;3 1;4 4;NODATA 0", "DATA")
outR.save("sub3/outRaster2.tif")

RuntimeError: ERROR 000875: Output raster: E:\sub3\outraster2.tif's workspace is an invalid output workspace.
From this error message, I can see that the .save() method is not using the working directory, but is using the E: root directory. I know that I can solve this problem by using

outR.save(arcpy.env.workspace + "/sub3/outRaster2.tif")

but I am wondering when the working directory is used and when not.  I looked for more information, but apparently this is not very well documented. Or is it just a bad practice to use subfolders in the working directory?


texts: point count from in Inverse Distance weighting procedure - QGIS

texts:
I have conducted an interpolation using IDW procedure. I have vectorized the interpolation in order to be able to specify a 3 class clasification.

I have additional point dataset indicating the families living in parts of the interpolated area. Each point repsresents a family. I would like to know, how many families live in each area, expressed as the proportion of total families.

How can I achieve this?

Is vectorization of the IDW appropriate here?


texts: Determining which QGIS layer has a spatial index?

texts:
I would like to loop over all layers within a QGIS project and determine which layer don't have a spatial index and cause QGIS to slow down.
I tried the method hasSpatialIndex() but this always returns SpatialIndexUnknown even if it's for example a shapefile with qix file present in the same folder (created from QGIS itself).
source=iface.activeLayer()
if source.hasSpatialIndex() == QgsFeatureSource.SpatialIndexNotPresent:
    print("not present")
elif source.hasSpatialIndex() == QgsFeatureSource.SpatialIndexUnknown:
    print("unknown")
elif source.hasSpatialIndex() == QgsFeatureSource.SpatialIndexPresent:
    print ("present")
else:
    print ("something else")

I guess it's a bug that QGIS does not recognize the qix files in the same folder where the shapefiles are stored?
Is there really no pyqgis-method to check layers if they are indexed? I am pretty sure that QGIS itselfs checks the existance of a spatial index because it uses the spatial index automatically if it exists.


texts: How to generate StreetFrom and StreetTo data from street segments?

texts:
Background:

I have a polyline feature class of about 4000 street centerline segments for which I would like to add data on what streets the segments intersect with.

I would like to add some new fields to my dataset that capture information on what street the polyline segment begins at and what street the polyline segment ends at (essentially the intersections)

What I have so far:

I wrote a quick Python script to calculate these values by doing Select By Location -> Boundary Touches:

# Import arcpy and set path to data
import arcpy
arcpy.env.workspace = "c:\Users\me\Documents\ArcGIS\Default.gdb"

#create a master search layer to encompass all streets.
searchLayer = arcpy.MakeFeatureLayer_management('StreetsCP', 'searchLayer')
#set the fields that we are going to be updating
fields = ('FULL_NAME', 'OBJECTID', 'StreetFrom', 'StreetTo')
#create a cursor to encompass all rows in streets
with arcpy.da.UpdateCursor('StreetsCP', fields) as uCursor:
        #loop over every row in streets
        for uRow in uCursor:
            #set where clause to create a feature layer from the row
            whereClause = '"OBJECTID" = ' + str(uRow[1])
            #create feature layer consisting of just the row we are on
            inputLayer = arcpy.MakeFeatureLayer_management('StreetsCP', 'inputLayer', whereClause)

            #select the streets that touch the boundary of the street we are on
            arcpy.SelectLayerByLocation_management(searchLayer, "BOUNDARY_TOUCHES", inputLayer, "", "NEW_SELECTION")

            touchingStreets = []
            uniqueStreets = []

            #populate the unique streets list
            with arcpy.da.SearchCursor(searchLayer, 'FULL_NAME') as sCursor:
                for sRow in sCursor:
                    if sRow[0] != uRow[0]:
                        touchingStreets.append(sRow[0])
                #eliminate duplicate street values
                uniqueStreets = list(set(touchingStreets))
            if len(uniqueStreets) <= 2:
                i = 2
                for street in uniqueStreets:
                     uRow[i] = street
                     i += 1
                uCursor.updateRow(uRow)

            arcpy.Delete_management("inputLayer")

The Problem:

The issue I am running into is the directionality of the street. In my street data I have address directionality using Left From/Left To/Right From/Right To, like in the below table:

The Street From and Street To should reflect this same directionality. What are some ideas for capturing the directionality of the street?

Some ideas I've been exploring:

Is there something exposed through ArcObjects that stores the endpoint and beginning point of a polyline using these metrics that I can use to select by location?
Can I somehow determine the directionality just in ArcPy?
As a last resort, Network Analyst might be able to offer me this functionality, but due to not having used Network Analyst before the ramp up time may be too much for my purposes.
If there is functionality that may help me in GRASS/QGIS I'd be willing to explore that as well.
I also have an "intersections" point file that could be used.


texts: create minimal spanning tree using table of distances between nodes rather than TSP

texts:
I am working on a transportation (linear programming) problem. Each potential solution also requires checking the distance of approx 20,000 routes with 5 to 15 nodes per route.
Calculating the shortest path for each route is likely to be time consuming using a TSP or similar call in pgrouting.

Would it be feasible to sacrifice some accuracy (for example < 10%) and instead build a minimal spanning tree for each route using a pre-calculated distance table of the shortest path from one node to another (say using OSRM).

Apologies if answer is blindingly obvious, I am new to GIS.


texts: Any issues in converting a AutoCAD 2010 DWG to FeatureClass?

texts:
I have a small .NET script that uses ArcObjects such as CADtoFeatureClass to convert DWG files to ESRI FeatureClasses.  This standalone program has been working great with DWGs created in previous versions of AutoCAD but I've been told that my code is now not working on DWGs from AutoCAD 2010.  As a temporary solution I can have the users save their DWG files in AutoCAD to previous versions but to permanently fix the problem I need to find out more details.  So does anyone know or has anyone had a similar issue in trying to convert DWG files to FeatureClasses using these types of ArcObjects on AutoDesk 2010 DWG files?

Thanks for any advice.

best regards,


texts:
Could you please help me with this error:

Traceback (most recent call last):
  File "C:\DDPeX\Export_GroupPDF.py", line 19, in
    finalPdf.appendPages(tempPdf)
NameError: name 'finalPdf' is not defined
Failed to execute (CFAsGroupDDP).

import arcpy
import os
mxd = arcpy.mapping.MapDocument("CURRENT")
output_folder = arcpy.GetParameterAsText(0)
ddp_lyr = mxd.dataDrivenPages.indexLayer
pdf_field = "CFA"

#Get a list of unique values to group by
pdf_list = set([row[0] for row in arcpy.da.SearchCursor(ddp_lyr,pdf_field)])

for pdfName in pdf_list:
    finalPdf = arcpy.mapping.PDFDocumentCreate(os.path.join(output_folder,pdfName + ".pdf"))
    for i in range(1, mxd.dataDrivenPages.pageCount + 1):
        mxd.dataDrivenPages.currentPageID = i
        cur = str(mxd.dataDrivenPages.pageRow.getValue(pdf_field))
        if cur == pdfName:
            tempPdf = os.path.join(output_folder,"temp_" + cur + str(i) + ".pdf")
            arcpy.mapping.ExportToPDF(mxd, tempPdf)
            finalPdf.appendPages(tempPdf)
            os.remove(tempPdf)
            finalPdf.saveAndClose()
            del finalPdf, tempPdf


texts: Creating perimeter around multiple polygons using ArcGIS for Desktop?

texts:
This may have a simple answer that I'm just missing.

I have multiple polygons in one shapefile I would like to make a larger Polygon around all of them.

Is there a tool that is able to do this?

Or even a Polyline would work if that is not possible, I can always do the Feature to Polygon tool after.

The images below better explain what I am looking to do (second is just a quick polygon I drew in a few seconds).

It is in order to show different sales areas for my company.