texts: Why TauDEM specialized tools (e.g. D-Infinity Distance Down) are only available with D-Infinity flow model?

texts:
I would like to know why most TauDEM specialized tools have been developed just for the D-Infinity flow model.
Is it because the D8 would not make sense for these algorithms, or they just have not been implemented?
In the first case, what is the reason?


texts: QGIS, Installation on Linux Mint 15 Mate

texts:
I've just tried to install the latest version of QGIS on Linux Mint 15 Mate using 'sudo apt-get update' and 'sudo apt-get install qgis'.
However, QGIS gives me the error message 'Couldn't load PyQGIS.
Python support will be disabled.'. Also, the installed QGIS version turns out to be 1.7.5 while I would prefer 1.8.0 or even 2.0.

I am pretty new to the world of Linux in general, so any advice how to get the latest QGIS version including a working Python plugin would be highly appreciated.

Best regards, T.


texts: I have deleted all the data in a SpatiaLite layer in QGIS. Can I get it back?

texts:
I am totally new to QGIS, and I am working my way through its tutorial manual.
After adding a spatialite layer in tutorial 2.2.4 (https://docs.qgis.org/3.16/en/docs/training_manual/basic_map/preparation.html), I somehow managed to delete all the data in the layer in module 3. The layer is still there but nothing shows on the map or in the features table. When I re-add the layer to the map, the data is still not there so I wonder if I have somehow permanently deleted the data from my file. I've been fiddling with it for a while.
Do you have a clue as to what I may have done or how to fix it?


texts:
When writing about GIS geometry I'm unsure whether to use the term "polygon" or "area".

In particular, does the word "polygon" suggest a specific type of geometry? For example, a colleague suggested to me that "polygon" is not an appropriate term for a donut feature!

But even if that's true, does anyone care? Is there any strong concensus either way on this?

Also, which word translates better into other languages such as German and French? I'm trying to consider users other than just English-speaking ones.


texts: How to count unit quantity in a certain polygon in QGIS

texts:
I have two shapefiles:

List item household units in each building in the whole city(green color building);
List item Polygons in the whole city(red color polygons);

There are 1 million buildings and 1000 polygons in the city.

How to count the unit in the certain polygons?


texts: Filtering out polygons that consists of less than 3 nodes

texts:
I have OSM data stored in PostGIS.

After running the following query:

select st_astext(st_forcepolygonccw(st_transform(geom,32633))), st_geometrytype(geom) from osm_versions.polygon_v1
limit 10

I get this result.

As you may see only 2 out of 10 spatial objects consist of more than 1 node yet all of them have GeometryType ST_Polygon. As a polygon is "a plane figure that is described by a finite number of straight line segments connected to form a closed polygonal chain" (after Wikipedia) I have no use for "polygons" with less than 3 nodes.

My questions are:
Is there any command that returns, given a set of geometries, a subset of geometries that consists of more than 2 nodes?
How should I proceed in order to filter out "polygons" with less than 3 nodes?

As for the question how come that points are counted as polygons, I have no idea as the data were provided by an external source.


texts: How to update attributes on spatialite table programmatically: SQLite Error

texts:
I am writing a plugin that updates 2 fields in an existing spatialite layer. The spatialite table contains several layers which are displayed on the current map. It executes to completion (very slowly), but the updates aren't appearing and I am getting an error that states:

  SQLite error: database is locked SQL:

I've tried putting the layer into edit mode before running the plugin and also calling startEditing() on the layer within the code.

def run(self):
   #bunch of code that works self.bRanges created here
   layer_name = self.dlg.cbStreets.currentText()
   self.layer = layer_list[layer_name]
   self.from2 = self.layer.fields().fieldNameIndex('from2')
   self.to2 = self.layer.fields().fieldNameIndex('to2')
   self.dp = self.layer.dataProvider()
   self.layer.startEditing()
   self.processEachVectorCount(self.updateRanges, 10)
   self.layer.commitChanges()

def processEachVectorCount(self, func, count):
    iter = self.layer.getFeatures()
    counter = 1
    for feature in iter:
        func(feature)
        counter = counter + 1
        if counter == count:
            break

def updateRanges(self, feature):
    # self.bRanges is a dict with a list as values
    bRange = self.bRanges.get(feature.attribute('road_id'))
    fid = feature.id()
    if bRange and len(bRange):
        bMin = min(bRange)
        bMax = max(bRange)
        print feature.attribute('full_name'), bMin, bMax
        #I'm seeing this print, so far so good
        attrs = {self.from2: bMin, self.to2:bMax}
        self.dp.changeAttributeValues({fid:attrs})

EDIT: The following code works.

def run(self):
   ... same setup as above, but no need to set the dataprovider
   self.layer.startEditing()
   self.processEachVectorCount(self.updateRanges, 10)
   self.layer.commitChanges()

def updateRanges(self, feature):
    ... again, same as above
    self.layer.changeAttributeValue(fid, self.from2, bMin)
    self.layer.changeAttributeValue(fid, self.to2, bMax)

The commitChanges call isn't strictly necessary. Without it the layer is left in edit mode.


texts: How do I edit a single label in QGIS 3.4 (Madeira)?

texts:
I'm making a map of North Dakota and trying to label the counties, but one of the names (Golden Valley) is extremely long in a very tall county. I don't want to turn the label vertically, but would rather put 'Golden Valley' on two lines so it will fit.

Is it possible to edit only that label so that it is on two lines? And how?


texts: How can I display points with points displacement renderer without circleoutline

texts:
I understand the point displacement, but I would like to display the circle without the circleoutline. Even when I change the "circle pen width" to 0,00 it displays a fine line. And there is no transparent color.
Do you have any ideas?


texts: Calculating Maximum Distance from Centre of a Polygon to the Boundary of a Polygon

texts:
I have a polygon ready in QGIS. What I want to do is to calculate the maximum possible distance in kilometres from the centroid of the polygon up to the boundary of that polygon.

I need two things:

Maximum possible distance covered from the centroid of the polygon up to their boundary
centre of the Polygon (Lat/Long)


texts: Explode/Decompose Polyline Network with MapInfo

texts:
I have a network of polylines which possibly include multi polylines as well. I wish to decompose these polylines into their fundemental line segments so I can produce the following table:

Link_ID | BegNodX | BegNodeY | EndNodeX | EndNodeY | Distance
--------------------------------------------------------------
0001    | 554612  |  456375  | 557864   |  455487  |  0.348
0002    | 554632  |  456355  | 557841   |  456522  |  3.454
....

I've noticed there is a Mapinfo tools such as Polybrk.mbx by Mr. Dressel but this tool only explodes one polyline selection. I need a whole network exploding and I do not wish to do each one manually.


texts: What is the country/region between Poland and Lithuania?

texts: What is the geographical area between Poland and Lithuania? It doesn't seem to have a name.


texts: I have a few shapefiles that I am analyzing in QGIS. They are in the WGS 84 coordinate reference system, which is what I want, however, the coordinates are in the seven digit range (they are all around 3829665, 5878362). How can I scale the coordinates all down by 10^5? e.g. I want them all around 38.29665, 58.78362


texts:
I am very new to Leaflet and coding and I am having some problems creating a map.

I have been able to create a map showing all my points but I would like to be to be able to distinguish between all different types of cheese via toggleable layers.

I have tried to use the solution shown in Leaflet: How to create toggleable overlays from a single GeoJSON FeatureCollection?, with no success.

This is my test JSFiddle, could you tel me what I am doing wrong?

Also I have all the Geodata saved in a cheese_map.geojson, that I would like to refer to in the javaScript rather than having 1000 lines of it.

I have managed to progress on this creating individual layers for each cheese type and using different icons as markers.


texts: finding the nearest point along the edge of a polygon

texts: I have a polygon shapefile which represents a study area (island). I also have a point shapefile representing some villages. What I am after is to locate for each village point the nearest point along the shoreline (the latter being represented by the edge of the polygon shapefile). Is there any tool to accomplish this in arcGIS 10.1 (advanced license)?


texts: How to use Leaflet Marker Cluster on markers from Leaflet-layerJSON (Overpass API)

texts:
I'm having troubles getting the Leaflet Marker Cluster working with Leaflet-layerJSON. I'm trying to make a Leaflet map with data from OpenStreetMap. The markers are showing on the map, but they are not clustering. It would be a great help, if you could point out my mistakes in the following code:

<meta charset="UTF-8"/>

<link href='https://api.tiles.mapbox.com/mapbox.js/v2.0.0/mapbox.css' rel='stylesheet' />

<link href='https://api.tiles.mapbox.com/mapbox.js/plugins/leaflet-markercluster/v0.4.0/MarkerCluster.css' rel='stylesheet' />
<link href='https://api.tiles.mapbox.com/mapbox.js/plugins/leaflet-markercluster/v0.4.0/MarkerCluster.Default.css' rel='stylesheet' />

I have tried to find answers to my question in other posts, but without any luck. I'm new to webmapping and stackexchange and apologize if my question is unclear.


texts: Performance Improvement of Postgres vs PostGIS

texts:
I am trying to run the following query thousands of times per second to determine which is any polygon contains the passed in point.

SELECT * from areas WHERE poly @> ?

Poly is the name of the polygon field and it has a GIST index on it.

I am wondering if I would see any performance benefit from using PostGIS over the regular postgress geospatial index for this simple query.  I am using Amazon RDS to host postgress and with 3 xl read replicas I can only get about 20,000 transactions per second and I need to triple that.  I have a webapp that makes the requests to postgres.  At 20k tps all the postgress replicas are at 100% cpu.


texts: Probabilities matrices and cellular automata in IDRISI

texts:
I apologize for asking about such elementary things, but I am utterly new to IDRISI yet am now having to learn how to perform certain tasks with it on very short notice.

Am I right that, when Markov chain analysis is performed in IDRISI, one of the outputs is a probabilities matrix? In other words, for a certain location L, for which land-cover data for years A and B are available, one can generate a matrix of possible changes for each land use from year A to year B? Is it possible, then, to use this matrix as an input to a cellular automaton that would predict land-cover changes in another location, M, based on (for example) current land-use data there? What tool or tools in IDRISI are needed for this type of analysis?


texts: Disjoint is very slow in Virtual Layer but fast on symbology

texts:
In QGIS, take a polygon grid (named 'grid') and a multipolygon layer (named 'surfaces').
For each feature in surfaces, I want to union all features of the grid that are disjoint (not intersected) with. I need to do this with a virtual layer (datasets regularly updated and exported with an atlas).
I have tried this two solutions, they work but are extremely slow (4-5 minutes):
select s.id, st_union(g.geometry)
from grid g, surfaces s
where not(st_intersects(g.geometry,s.geometry)) -- st_dissolve also works but is slower...
group by s.id

select s.id, st_union(g.geometry)
from grid g
join surfaces s on not(st_intersects(g.geometry,a.geometry))  -- st_dissolve also works but is slower...
group by s.id

The problem seems to be the intersects. The time is about the same with or without st_union.
For testing, I tried to add a rule-based symbology to my grid, with this rule:
not(intersects($geometry, @atlas_geometry ))

And I get the result instantly. Why intersects is slow on the Virtual Layer and fast in symbology? And how can I get a fast result in the Virtual layer?

EDIT
What I want to achieve (without union to see the grid):

And with union, to calculate the area of each part of the grid:


texts: Creating new GRASS locations while looping in Python script

texts:
I am making a Python script for processing my data.

import grass.script as grass
import numpy as np

def main():
     f = open("D:\\Glaciers_selected.csv", 'r')
     fl = f.readlines()
     for x in range(1, len(fl)):
           nf = fl[x].split(';')[1]

           #### Problem in next executable line
           ## I need to create a new location, using the projection from the nf file.
           grass.read_command('grass76' flags='c' nf ./UTM60S/     ## Here I am stuck

           ## The rest of the script works fine
           grass.read_command('r.import', input= 'D:\\GIS\Farinotti_glaciers\\' + nf + '.tif ', output='RGI60_11', overwrite=True)
           grass.read_command('g.region', raster='RGI60_11', flags='a', overwrite=True)
           ...
           ...

if __name__ == "__main__":
    main()

My script runs fine if I start GRASS, add a new location with the location wizard and then run each line, one by one from the Python tab in the layer manager. However, I do not know how to add a location once I am inside a loop of a Python script.

When I run grass76 -c  from the console, the location Wizard appears ...

I read Creating a GRASS location from a bash script and there seems to be a way of doing this, but I was not able to follow the explanation (which is for BASH, by the way, not for Python).


texts: Calculating field on if-else-condition checking input parameters for ?

texts:
I want to calculate a field depending on 3 input string parameters. Two of them can contain Null-values. My old VB expression was not able to work with Null-values, so the script stopped.
Now I want to develop a script which allows me to solve all my problems:

set an if-else-condition
calculate the field with Python expressions instead VB
eliminate Null-values by not computing/using these cells in the expression.

This is the input table to get my parameters

And here is my script:
import arcpy, os
from arcpy import env

#Get Parameters
importfeature = arcpy.GetParameter(0)
workspace = arcpy.GetParameter(1)
owner = arcpy.GetParameter(2)
criteria = arcpy.GetParameter(3)
keep_fieldname1 = arcpy.GetParameter(4)
keep_fieldname2 = arcpy.GetParameter(5)

#Set Workspace
arcpy.env.workspace = workspace

#Add Field
arcpy.AddField_management (importfeature, criteria, "text")

#Calculate field criteria
if keep_fieldname1 == None and keep_fieldname2 == None:
  expression = '"%owner%"'
elif keep_fieldname2 == None
  expression = '"%owner%" +" - "+ !%keep_fieldname1%!'
else
  expression = '"%owner%" +" - "+ !%keep_fieldname1%! +" - "+ !%keep_fieldname2%!'
arcpy.CalculateField_management (importfeature, criteria, expression , "Python")

I need to say that the parameters are automatically given by an iterator to the script via "get values". So don't worry about the input. This is working fine.
How do I solve my problems?


texts: Using the field calculator to copy one column's values to another's BLANK values

texts:
This is probably a basic and simple question which has been answered multiple times before, but my Google-fu has come up empty, so here we are.

Long story short, to simplify my analysis, I want to copy the data from the first zoning column field to the second into ONLY in cell values which are blank. I will attach a screenshot to show what I am talking about.

I have tried various [case when] and [if] examples, but nothing seems to get this thing to go.


texts: How to use search cursor in a given extent

texts:
I want to perform a selective search cursor for a feature class for a given extent.
I have tried setting env.extent in the below code but it is returning the whole feature class instead of a fixed extent.

I am using ArcGIS 10.3 and the featureclass is in PostGIS.

import arcpy
from arcpy import env

# Set the workspace
env.workspace = "C:\Users\ADevra\AppData\Roaming\ESRI\Desktop10.3\ArcCatalog\Connection to DAT-PGTEST-01.sde"
#Setting the Environment extent
arcpy.env.extent = arcpy.Extent(438859, 184490, 445034, 190295)

fc1 = "postgretest.postgres.boundary_line"

##Select the features in a given extent
with arcpy.da.SearchCursor(fc1,["fid"]) as innercursor:
    for searchrow in innercursor:
        print searchrow[0]

del(innercursor)

UPDATED ANSWER:

import arcpy

arcpy.env.workspace = r'DAT-PGTEST-01.sde'
extentPoly = "Test Shapefileszextent.shp"
fc = r'postgretest.postgres.boundary_line'
dataset = "postgretest.postgres_tmp.boundaryline"
desc = arcpy.Describe(dataset)
extent = desc.extent

array = arcpy.Array()
# Create the bounding box
array.add(extent.lowerLeft)

array.add(extent.lowerRight)
array.add(extent.upperRight)
array.add(extent.upperLeft)
# ensure the polygon is closed
array.add(extent.lowerLeft)
# Create the polygon object
polygon = arcpy.Polygon(array)
array.removeAll()
# save to disk
arcpy.CopyFeatures_management(polygon, extentPoly)
del polygon

arcpy.MakeFeatureLayer_management(in_features=fc, out_layer='fclyr')

arcpy.SelectLayerByLocation_management(in_layer='fclyr',    overlap_type='INTERSECT', select_features=extentPoly)

    for rows in arcpy.da.SearchCursor('fclyr', ["postgretest.postgres.boundary_line.fid"]):
    print rows[0]


texts: Customizing right click (context) menu of selected feature using ArcPy?

texts:
I want to confirm whether arcpy can customize the right click menu of the selected feature, which means I want to remove the default items and add in my own item. The following pictures show what I want to do:

(1) picture 1 - the default when right click a selected feature

(2) - picture 2 - what I want when right click a selected feature


texts: Export shapefile to KML for Google Earth - using QGIS

texts:
I am a real beginner with QGIS.
I have shape files of my study site.
I would like to have the nice google earth background.
I would like to be able to see my files in Google earth, then I can zoom and see details at a fine scale.
from what I understood, I need to export my files as kml files.
I clicked right on the layer and clicked on "saved as", I chose a place to save and a name for the file, I chose WGS 84, and saved.
When I open the kml layer in Google earth and/or QGIS it is completely wrong.
In Google earth, appears some segments every where on the globe; in QGIS, I have on the right a vertical line.
Can anyone explain to me what goes wrong?


texts: QGIS: how to combine a GeoTIFF with a country border shapefile?

texts:
I have two files. A GeoTIFF with different values and a shapefile with country borders of this area. Can I generate with these two files an image with the country borders and the data from the GeoTIFF?

I'm also interested in any other software/scripting solution for this.
I did the second image with matplotlib and basemap but the country borders are wrong (outdated).


texts: Converting multipart line into single parts using QGIS?

texts:
I have a shapefile in which a bunch of the lines are multipart. Sometimes this is due to joining ways that are not quite connected (gap less than a foot of real-world distance) and sometimes the endpoints of each part are at exactly the same place (thanks to snapping before joining).

This is really annoying because it means that splitting does not work properly, among other lesser problems.

I've tried every suggestion I've come across. GRASS v.clean.snap with tolerance 1 does nothing. Exporting to KML and removing the gaps in the KML file works but only if the parts are in the correct order and direction (otherwise I get garbage connecting lines). Exporting to GPX completely fails because GPX can only handle single-part lines.

Am I missing something here? Why does this problem even exist? I understand why it happens when there is a slight gap, but I can't see any situation in which joining two topologically connected ways in a multipart line is a feature rather than a bug.

Here (click "Click here to start download from sendspace") is a shapefile with two multipart lines. I-85 has a tiny gap (at the state line), while I-795 connects perfectly yet has three parts. Does anyone have a solution for this, one that can be batch-applied to a file with many lines?


texts: how to setup gdal_calc.py for command line use (windows)?

texts:
I've searched the web extensively (e.g. here)  but I am not able to find any documentation on how to setup gdal_cal.py.

When I try to run gdal_calc.py from my command line e.g.:

  gdal_calc.py -A input1.tif -B input2.tif --outfile=result.tif
  --calc="A+B"

i simply get prompted with which program i want to open the file.
how do i setup up gdal_calc.py so that I can run it from my command line? Or where can i find documentation about it.

INFO:

i have newest version of python
i have newest version of QGIS (3.0)
i have newest version of python (2.7.14)(using spyder IDE(3.2.6))
i have newest version of GDAL (2.2.4)
running on windows


texts: Customizing widgets in ArcGIS Viewer for Flex 3.4 to support French?

texts:
How can I customize my widgets in flex viewer such that the labels will be in French?

I work with flash builder 4.6.


texts: Anyone know the location of these Geo Vector Tiles from an iOS 13 image?

texts:
Last year I convert these files, but I cannot remember and I cannot seem to find any current information on them. You can find more here: https://cyberfenixtech.blogspot.com/2019/10/examining-maptiles-from-ios.html
Here is the header information for the file:

It doesn't look like a PBF format. I looked around at several tile sites with example data and I couldn't match the file header to it.
These are MapTiles located Data\mobile\Library\Caches\MapTiles\MapTiles.sqlitedb. SQL database in iOS.


texts:
I have a polyline dataset with lines of different lengths. I want to shorten all polylines the same distance on both ends (e.g. 10m each side). Is there a tool or script available that could help me do this? I want to automate this process, so manual editing is no option.

I´m using ArcGIS 10.0 with an ArcView license.


texts: Calculate probability of occurrence from multiple grids

texts:
I generate a lot of spatial grid files and I want to calculate the probability that a certain threshold isn't exceeded. In this example here I want to see what's the probability of each cell that the value is equal or smaller than 0. I do this by counting the occurrences and divide them by the number of rasters involved. This code works, but this looks terribly bloated. Is there a way to make this more efficient, especially with the aim of doing this for hundreds of really large raster maps? (and I need to work with SpatialGridDataFrame, not with raster files from raster package, thus the conversion)

library(raster)
library(sp)
seeds <- c(2013,2014,2016)
n <- 0
pfail <- matrix(0,5,5)
pfail <- raster(pfail)
pfail <- as(pfail, 'SpatialGridDataFrame')
for (i in seeds) {
  set.seed(i)
  n <- n+1
  xy <- matrix(rnorm(400),5,5)
  rast <- raster(xy)
  sgdf <- as(rast, 'SpatialGridDataFrame')
  rast <- paste("rast",n,sep="")
  assign(rast,sgdf)
  for (j in 1:length(sgdf)) {
    pfail[j] <- ifelse(sgdf$layer[j]<=0,pfail$layer[j]+1,pfail$layer[j])
  }
}
pfail$layer <- sapply(pfail$layer, function(x) x/length(seeds))


texts: JTX/Workflow Manager dependencies to control workflow

texts: When creating a new job in JTX, one of the required fields (under "Additional Information") tab we have is "Asset Type". The three choices for "Asset Type" are "Water", "Sewer" and "Water and Sewer". My question is, is it possible to setup JTX where depending on the Asset Type, certain step(s) in the workflow will be skipped. For example, if I choose Water as my asset type, I may want to skip a certain step and if I choose sewer as the asset type I may want to skip a different step.


texts: Leaflet multiple circle riseOnHover or bringToFront problem

texts:
I am drawing a multiple circle using Leaflet API. I need to apply hover effect so that I can hover on circle it will rise up or bring to front. I am using the following codes

for (var i = 0; i < count; i++) {
var circleLayer = L.circle([latitude[i], longitude[i]], bufferDistance[i], {
    color: 'white',
    weight: 1,
    fillColor: getColor(data),
    fillOpacity: 0.8,
    riseOnHover: true
}).addTo(map);
circleLayer.on('mouseover', function(circleLayer) {
    circleLayer.bringToFront();
});
circleLayer.addTo(map);
}

All the circles are drawing accordingly but there is no mouseover or hover action is taking place. What I am doing wrong?


texts: How can I see the coordinate transformation parameters in QGIS?

texts:
How can I see (and edit if possible) the transformation method and parameters used by QGIS for "on-the-fly" transformation, and then for saving data in another CRS? Like in ArcGIS?

Here is the similar question but answered differently, I don't want to create a new coordinate system, I want to check the transformation between existing in QGIS coordinate systems.

In ArcGIS I can select between several transformation methods for the pair of CRS, and what about QGIS?

I've changed the default WGS84 CRS to another local CRS. For example, if I set for the project EPSG 3763 (ETRS PT TM 06), described as +proj=tmerc +lat_0=39.66825833333333 +lon_0=-8.133108333333334 +k=1 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs, and add data in EPSG 20790 (Datum Lisboa) +proj=tmerc +lat_0=39.66666666666666 +lon_0=1 +k=1 +x_0=200000 +y_0=300000 +ellps=intl +towgs84=-304.046,-60.576,103.64,0,0,0,0 +pm=lisbon +units=m +no_defs, it is reprojected on-the-fly and looks more or less correctly. But if I cannot see is it really correct I want to know the transformation used. There are several recommended methods of different accuracy, and I want to know which one was used (if it was one of them).

(I am not fluent in QGIS, just learning).


texts: Overcoming spatial bias in ArcMap or ArcGIS Pro?

texts:
I have a shapefile showing invasive plant locations but they mostly collected from more accessible places such as near roads or hiking trails.

Is there a method to correct this bias in ArcGIS Desktop?


texts: How to get pixel coordinates in QGIS georeferenced output map?

texts:
using the QGIS Georeferencer GDAL plugin, I can select GCPs (ground control points) on a map (plain raster input map image in a format like PNG or TIF) in an unknown projection and boundaries, and specify the coordinates for those points in my desired output projection; then the Georeferencer plugin can reproject the map image from the unknown input projection into the specified output projection, epsg:4326 in this case:
(each GCP is the x,y coordinate pair of a pixel in the input image followed by the lon,lat coord of that pixel on the map in the desired target projection)

    gdal_translate -of GTiff -a_srs EPSG:4326 -gcp 738.445 571.363 -5.3664 39.8004 ...[lots of gcps]... -gcp 656.449 529.382 -5.2034 40.4102 inputmap.png  tmpmap.tif

    gdalwarp -r near -order 2 -co COMPRESS=NONE  tmpmap.tif  outmap.tif

When the input image with its ground control points is reprojected into the output image, the GCPs are reprojected too, and shown at their correct new locations.

Is there a programmatic way to obtain the pixel coordinates of the GCPs in the reprojected output image?
And more generally, is there a way to reproject a single point/coordinate pair from an unknown projection into a known specified projection, using GCPs to define the unknown projection like in the above commands?


texts: Getting differences between old and new versions of ArcSDE database?

texts:
Is there a tool or a process that can reveal changes made to the schema of an ArcSDE database when comparing two versions of the same database hosted on different servers?

I have an old database that has been running a site, and a newer version, and I'd like to be able to get a heads up about what has changed in terms of schema and contents.


texts: How to remove 'orphan' tables in SpatiaLite database

texts:
I have a SpatiaLite (4.3) database created in QGIS that I have modified numerous times using the SpatiaLite-GUI. I've realized today that when I add a table within QGIS from that DB, the selection dialog lists several tables that have been deleted long ago.
I can find references to the deleted tables with the following queries:
select * from vector_layers
select * from geometry_columns

I know that I can simply save the remaining tables into a new DB, but I'm wondering if there is some sort of bookkeeping query to clean up a database.


texts: Scale (or stretch) Leaflet map tiles to look bigger

texts:
I am using leaflet and I am using the same tile provider as the main openstreetmap's website (http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png). However, leaflet seems to fit many more tiles in 1 screen than openstreetmap. This results in the icons been super small, and I can barely read the letters in a 1080p monitor. Look at this screenshot:

That is the same exact tile in openstreetmap's website and my website.

This is the tile:

How can I scale (or stretch) the map tiles, so that my tiles look bigger ? It seems like Openstreetmap doubles scales their tiles by 200%. Which means that the 256x256 image ocupies a space of 512x512 in the actual website.


texts: Sudden increase in the download size of images from GEE

texts:
I am working on LandSat 8 data to build a classifier for built up and non built up areas in a particular district. I was using the following code for classification and then to export my classified image into my drive:

var cloud_masks = require('users/fitoprincipe/geetools:cloud_masks');
var maskClouds = cloud_masks.landsatTOA();

//Loading India image, the extracting data for Haryana (a state in India) and then subsequently Ambala (a district in Haryana)
var bands = ['B1','B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B10', 'B11'];
var india = ee.FeatureCollection('ft:1UDdgOCf8DoRJ9bVm-UVbR6CqxtkJToLQjTFd0r0Z','geometry')
.filter(ee.Filter.eq('Name','India'))
.geometry();
var india_image = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
.filterBounds(india)
.filterDate('2014-02-01','2014-09-01')
.sort('CLOUD_COVER')
.map(maskClouds)
.median();
var district = ee.FeatureCollection('ft:1PA2zwArj8EsplrX9eMxJ2H_TICyyx855KPnbJhC1','geometry')
.filter(ee.Filter.eq('name','Gurgaon'));
var haryana_image = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
.filterBounds(district)
.filterDate('2014-02-01','2014-04-01')
.sort('CLOUD_COVER')
.map(maskClouds)
.median();
var input = haryana_image;
input = addBands(input.select(bands));
india_image = addBands(india_image);

//Loading the points from the fusion table and training the classifier
var ft = ee.FeatureCollection('ft:1fWY4IyYiV-BA5HsAKi2V9LdoQgsbFtKK2BoQiHb0');
var ft_builtup = ft.filter(ee.Filter.eq('class',1)).limit(1200);
var ft_nonbuiltup = ft.filter(ee.Filter.eq('class',2)).limit(1800);
ft = ft_builtup.merge(ft_nonbuiltup);
var new_bands = ['B1','B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B10', 'B11','NDBI','NDVI'];

function addBands(image){
  var ndvi = image.normalizedDifference(['B4', 'B3']).rename('NDVI');
  var ndbi = image.normalizedDifference(['B5', 'B4']).rename('NDBI');
  var ndwi = image.normalizedDifference(['B6', 'B6']).rename('NDWI');
  return image.addBands(ndvi).addBands(ndbi).addBands(ndwi);
}

// Load a Landsat 8 image to be used for prediction.
var training = india_image.sampleRegions(ft,['class'],30);
var trained = ee.Classifier.cart().train(training, 'class', new_bands);

input = input.clip(district);
input = input.classify(trained);
input = input.expression('LC==1?1:0',{'LC':input.select('classification')});

Export.image.toDrive({
  image: input.clip(district),
  description: 'Test',
  scale: 30
});

I used this code for download in September and it gave me a .tif image of size 175.5 KB. Now I am running the same script again now in November but it gives me an error saying that the number of pixels limit is exceeded. So I modified it as:

Export.image.toDrive({
  image: input.clip(district),
  description: 'Test',
  maxPixels: 499295920080,
  scale: 30
});

This code now gives me 16 split .tif files, with their sizes summing up to a total of 298 MB.

My code remains untouched and I even checked my images of 3 months back, they are not corrupted or incomplete. I cannot figure out what happened in this span of 3 months that led to such differences. Could there be some changes in certain libraries I am using or the LandSat 8 images I am training and classifying on that are leading to such discrepancies?


texts: Unable to create field for PostGIS layer. PYQGIS

texts:
I'm having some troubles to create a new field in a postgis layer.This is the code I'm trying to use, but nothing happens, no field is created:

from PyQt4.QtCore import *

dbnam = "pruebapython"
uri = QgsDataSourceURI()
uri.setConnection("localhost", "5432", dbnam, "postgres", "xxxxxxx")
uri.setDataSource("", "masa", "geom","","gid")

layer = QgsVectorLayer(uri.uri(),"masa","postgres")
#layer = QgsVectorLayer('/xxx/xxx/xxxxx/xxx/xxx/DatosTest/MASA.SHP', 'MASA', "ogr")
layer.startEditing()

layer.dataProvider().addAttributes([QgsField("mytext", QVariant.LongLong )])
layer.updateFields()
layer.commitChanges()
idx = layer.fieldNameIndex( 'mytext' )
if idx == -1:
    print "No field"
else:
    QgsMapLayerRegistry.instance().addMapLayer(layer)
print 'Done!'

If I uncomment the layer which is a shapefile and I comment the posgres one, the field is created perfectly. What am I missing? Could anyone help me?

The QGIS log shows this:

Consulta errónea: ALTER TABLE "masa" ADD COLUMN "mytext"  devolvió 7 [ERROR:  syntax error at end of input
            LINE 1: ALTER TABLE "masa" ADD COLUMN "mytext"]

How can I solve this?


texts: Extracting the closest Point to a Line within a Polygon using ArcGIS for Desktop

texts:
I have the following shapefiles:

Polyline: Detailing a road network.
Polygon: Containing approx. 2000 polygons enclosing sections of that road network.
Point: Containing numerous (OS Addressbase) points falling within each of the polygons.

Using ArcGIS for Desktop, I would like to identify and extract to a new shapefile the point within each polygon that is closest to the road network, and also calculate the distance that each closest point is from the road network.


texts: Getting GRASS plugin working in QGIS for Windows?

texts:
I've installed QGIS 2.16.3-1 on Windows 10 via the OSGeo4W 64-bit installer and also selected the following packages:

grass (7.0.5-1)
qgis-grass-plugin7 (2.16.3-1)
qgis-grass-plugin-common (2.16.3-1)

I can't find the GRASS plugin in the list of plugins at Plugins > Manage and Install Plugins > Installed (panel)

And the QGIS Browser can see the grass location and mapset folder, but doesn't have the GRASS icon like it does in Linux. So layers can't be added to the map. I guess that means that the gdal-grass plugin is unavailable as well.

I tried launching QGIS Desktop 2.16.3 with GRASS 7.0.4 but get the same results as above. Maybe the GRASS version mismatch is a problem or just a shortcut name typo.


texts: Integrating QgsMapCanvas window in flask template

texts:
I have a code to show a vector layer in the canvas window using PyQGIS.
from qgis.core import *
from qgis.gui import *
from qgis.utils import *
from PyQt4.QtCore import *
from PyQt4.QtGui import *

app = QgsApplication([], True)
path = "C:/Program Files/QGIS2.18/apps/qgis"
app.setPrefixPath(path, True)
app.initQgis()
canvas = QgsMapCanvas()
title = "PyQGIS Standalone Application Example"
canvas.setWindowTitle(title)
canvas.setCanvasColor(Qt.white)
layer_info = 'LineString?crs=epsg:4326'
layer =  QgsVectorLayer(layer_info, 'MyLine' , "memory")
pr = layer.dataProvider()
linstr = QgsFeature()
wkt = "LINESTRING (1 1, 10 15, 40 35)"
geom = QgsGeometry.fromWkt(wkt)
linstr.setGeometry(geom)
pr.addFeatures([linstr])
layer.updateExtents()
QgsMapLayerRegistry.instance().addMapLayer(layer)
canvas.setExtent(layer.extent())
canvas.setLayerSet([QgsMapCanvasLayer(layer)])
canvas.zoomToFullExtent()
canvas.show()
exitcode = app.exec_()
QgsApplication.exitQgis()
sys.exit(exitcode)

I have a plan to develop a flask application. So I want this canvas window to be integrated into the flask Html template.


texts: Combining attributes of tiff layer file and shapefiles they fall into?

texts:
I am currently using ArcMap 10.5 and having trouble combining information of two layers. One is a tiff file with land classification data, and the other is a shapefile with broad scale regions.

Is there a way to combine the two layers so that the attribute table for the tiff file lists which of the shapefile regions it falls into or the percentage of cover for that specific shapefile?


texts:
I have a map showing the outline of an area which I'm using to clip other layers using the data frame > clip options.

In addition, I have other data frames to show some detail and want to show the extent indicators on the first map.

Is it possible to exclude the extent indicator from the data frame clip options as it isn't listed as a separate layer in the Exclude Layers... button?

So in the in the image below I want to keep the blue lines clipped to the shape, but have the entire red extent indicator visible.


texts:
I have a set of 32-bit TIFF files of DTM data (50 sm) and I need to save them as 16-bit.

In order to save some space, I compressed them using LZW algorith and saved as 16 bit tiffs using Arcpy.

But I see that the precision of height of compressed DTM is not preserved.

Is there a way to preserve it?

import arcpy
import os

arcpy.env.workspace = r"D:\Venera_Workflow\COMpression_TEST\TEST_uncompressed_input"
rasList = arcpy.ListRasters()
outWS = r"D:\Venera_Workflow\COMpression_TEST\TEST_Output_LZW"
arcpy.env.compression = "LZW"
for ras in rasList:
    arcpy.CopyRaster_management(ras, os.path.join(outWS, ras, ), "", "", "", "NONE", "NONE", "16_BIT_SIGNED")


texts:
In ArcMap, is it possible to display an image by using a path in the attribute table to show an image in MapTips?

Any Python or vbScript expression are welcome as long as it is MapTips which shows the image.

I do not want to use HTML Popus.


texts: Iso Cluster Unsupervised Classification error 000814 Invalid File Type in ModelBuilder of ArcMap

texts:
When trying to run the Iso Cluster Unsupervised Classification tool in ModelBuilder my model always fails at this step due to an invalid file type.
However, when I pull this process out and place it in a stand alone model everything works fine.
Has anyone experienced this problem before?


texts: Buffering multiple features with Geoserver WPS

texts:
I'm trying to create a buffer around more than one LineString feature at a time with geoserver 2.8.x. and openlayers 3. Currently, i can create a buffer around a single line feature with the following GML

<?xml version="1.0" encoding="UTF-8"?>
<wps:Execute version="1.0.0" service="WPS" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns="http://www.opengis.net/wps/1.0.0" xmlns:wfs="http://www.opengis.net/wfs" xmlns:wps="http://www.opengis.net/wps/1.0.0" xmlns:ows="http://www.opengis.net/ows/1.1" xmlns:gml="http://www.opengis.net/gml" xmlns:ogc="http://www.opengis.net/ogc" xmlns:wcs="http://www.opengis.net/wcs/1.1.1" xmlns:xlink="http://www.w3.org/1999/xlink" xsi:schemaLocation="http://www.opengis.net/wps/1.0.0 http://schemas.opengis.net/wps/1.0.0/wpsAll.xsd">
  <ows:Identifier>JTS:buffer</ows:Identifier>
  <wps:DataInputs>
    <wps:Input>
      <ows:Identifier>geom</ows:Identifier>
      <wps:Data>
        <wps:ComplexData mimeType="text/xml; subtype=gml/3.1.1">
            <gml:LineString>
                <gml:posList>0.0 0.0 10.0 0.0 10.0 10.0</gml:posList>
            </gml:LineString>
        </wps:ComplexData>
      </wps:Data>
    </wps:Input>
    <wps:Input>
      <ows:Identifier>distance</ows:Identifier>
      <wps:Data>
        <wps:LiteralData>2</wps:LiteralData>
      </wps:Data>
    </wps:Input>
  </wps:DataInputs>
  <wps:ResponseForm>
    <wps:RawDataOutput mimeType="text/xml; subtype=gml/3.1.1">
      <ows:Identifier>result</ows:Identifier>
    </wps:RawDataOutput>
  </wps:ResponseForm>
</wps:Execute>

How can i extend this GML to buffer more than one LineString feature at a time?

EDIT

Technically we could accomplish the same thing using a GeometryCollection, but this approach doesn't seem to work well when the vertices of a set of LineStrings in the data are spatially coincident, forming a closed sub component similar to a polygon. The result is that you get a buffer around those shapes, but for some reason, you also get the interior of that 'closed polygon' as part of the result, which would make it impossible to perform further geoprocessing like intersects or contains accurately on the result of that buffer.

In the image below, the green layer is the result of buffering the brown LineString features using a GeometryCollection. It shows in certain areas where the LineStrings are closed off that the interior of the closed area is added to the result when it really should have been ignored by the WPS/JTS buffer function


texts: GeoServer: CQL_Filter working with operator '=' but not with operator 'IN'

texts:
WMS 1:

http://localhost:8080/geoserver/mtm/wms?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&FORMAT=image%2Fpng8&TRANSPARENT=true&layers=copper_cable&tiled=true&tilesorigin=-20037508.34%2C-20037508.34&viewparams=import_job_code%3Anull&cql_filter=id%20%3D%274ABC-01-12-ABC-108%27&CRS=EPSG%3A3857&STYLES=&WIDTH=2858&HEIGHT=778&BBOX=17003975.290642627%2C-3209930.007007534%2C17005681.98616457%2C-3209465.4131950526

WMS 2:

http://localhost:8080/geoserver/mtm/wms?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&FORMAT=image%2Fpng8&TRANSPARENT=true&layers=copper_cable&tiled=true&tilesorigin=-20037508.34%2C-20037508.34&viewparams=import_job_code%3Anull&cql_filter=id%20IN(%274ABC-01-12-ABC-108%27)&CRS=EPSG%3A3857&STYLES=&WIDTH=2858&HEIGHT=778&BBOX=17003975.290642627%2C-3209930.007007534%2C17005681.98616457%2C-3209465.4131950526

In WMS 1 I am using CQL_Filer has operator '=' and is working fine

In WMS 2 I am using CQL_Filter has operator "IN" and is not workig.


texts: Use Field Calculator to change "Status" field in Annotations

texts:
I am trying to turn off a lot of annotations at once by using Field Calculator to set status to Unplaced. Just typing "Unplaced" into the expression, I get this error:

  The calculated value is invalid for the row with ObjectID = 408. For example the calculated value may be too large for the field or you may be trying to add a string to a number field. This row will not be updated.

I understand that I'm trying to drop a string value into the field and it won't accept that. What data type/value should I be using instead?


texts: Runtime Error Rerunning Arcpy Update Cursor Script

texts:
The following script works to update coordinates in my ArcGIS feature class, but if I run it twice in a row from ArcMap on the same geodatabase feature class, I get a "RuntimeError: Unexpected Operation[outputZoneTrack FeatureClass name]" error on the "updateCursor.updateRow(row)".  If I run it on a dataset, then run it on a different dataset, I can then re-run it on the first dataset again with no problem.  I just can't run it on the same dataset twice in a row. Why does this error occur when I run the script twice in a row on the same dataset?

outputZoneTrack = "C:\\Temp\\geodatabase.gdb\\ozt"
outputWorkspace = os.path.split(outputZoneTrack)[0]
outputZoneTrackName = os.path.split(outputZoneTrack)[1]

arcpy.env.overwriteOutput = True
if arcpy.Exists(outputZoneTrack):
    arcpy.Delete_management(outputZoneTrack)

spatial_reference = arcpy.SpatialReference(4326) #Use WGS84 so coordinates are calculated in Lat/Long
with arcpy.da.UpdateCursor(outputZoneTrack, ["ZSTARTLAT", "ZSTARTLON", "ZENDLAT", "ZENDLON", "SHAPE@"], "", spatial_reference) as updateCursor:
    for row in updateCursor:

        row[0] = row[4].firstPoint.Y
        row[1] = row[4].firstPoint.X
        row[2] = row[4].lastPoint.Y
        row[3] = row[4].lastPoint.X
        updateCursor.updateRow(row)

mxd = arcpy.mapping.MapDocument("CURRENT")
    for lyr in arcpy.mapping.ListLayers(mxd):
        if lyr.name == "Segment Oil Zones":
            if lyr.supports("DATASOURCE"):
                lyr.replaceDataSource(outputWorkspace, "FILEGDB_WORKSPACE", outputZoneTrackName, True)
                lyr.name = "Segment Oil Zones"
                lyr.visible = True
                arcpy.AddMessage("\nOiling Zone Snapped to ACP Segments: ON\n")
                arcpy.RefreshTOC()

Note - this problem only occurs if the script is run from ArcMap and the feature class is loaded within the ArcMap TOC - I can run it multiple times in a row on the same dataset from ArcCatalog.


texts: Sentinel-1 missing data gaps between two images GEE

texts:
I am trying to export a Sentinel-1 mosaic with the code below. I keep getting this missing data gap between the two images. Even when I export the images separately, and mosaic them myself in ArcGIS, I still have this problem.

https://code.earthengine.google.com/91fea33548352ade6dd961ed3d34d883

var start_date = ee.Date('2019-08-01');
var finish_date = ee.Date('2019-08-15');
var orbit = 'ASENDING';

// Get the VV collection.
var collectionS1 = ee.ImageCollection('COPERNICUS/S1_GRD')
    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
    .filter(ee.Filter.eq('instrumentMode', 'IW'))
    //.filter(ee.Filter.eq('orbitProperties_pass', orbit))
    .filterDate(start_date, finish_date)
    .filterBounds(polygons);

 // Get the VV collection.
var collectionVH = collectionS1.select('VH');

var VH = ee.Image(collectionVH.median());

Map.centerObject(polygons, 7);
Map.addLayer(VH, {min: -25, max: 5}, 'VH', true);

//Import Tool
var batch = require('users/fitoprincipe/geetools:batch');

// Change Folder Name

var folder = 'GEE_EXPORT';

// Download Collections
batch.Download.ImageCollection.toDrive(collectionVH, folder,
        {name: 'VH_'+orbit+'_'+'{id}',
               scale: 10,
              region: polygons,
                fileFormat: 'GeoTIFF',
                  maxPixels: 1e13,
              type: 'float'});


texts:
Has anyone seen any publically available benchmarks of MSDs vs MXDs?

If we ignore the fact that as part of the MSD publishing process, we run it through some validators with the map service publishing toolbar, in terms of rendering, how much quicker do MSDs render?

And can we break this down further.  I know MSDs are faster to render, but are there particular situations where they will render a lot faster?

Basic vector layers Vs Imagery
Simple symbols vs Multi-layer symbols
Lots of layers Vs Lots of features, etc

I thought there would be something in the Capacity Planning Tool wiki, but could not spot anything that goes into detail on what circumstances the MSD outshines the MXD.

I want to see some stats and graphs on MXDs vs MSDs.

As a sidenote, at 10.1 we only have the option of MSDs.


texts:
I'm trying to render tiles out of a KML file.
Normally, I'd render tiles out of an planet.osm extract. For example, I'd select a tag using osmosis and then process it with osm2pgsql. Then I'd generate the tiles using Mapnik's python script and I'd serve them with OpenLayers.

Maintaining this structure, i.e. PostgreSQL (with PostGIS) and OpenLayers, I'd like to know what's the best way to do this starting from a kml file, rather than an osm one. An option would be to transform the kml file into an osm file. Another option would be to find something like kml2pgsql. I haven't found any questions nor ways to do this (I'm specially surprised this is not a duplicated question).

Do you know how to do this? (might be using something else, as far as it's supported in Linux and uses Postgresql - Openlayers).

Answer:
I marked Andre answer as the final one because while I'm not using exactly what he said it put me in the right path. I ended up using TileMill rather than QGIS, but the process is practically the same...


texts: How to generate more pbf zoom levels with mbtiles file?

texts:
I have a 8GB canada.mbtiles vector file which contains all of Canada up to zoom level 14.
How do I generate zoom 15,16 from this canada.mbtiles file?  The goal is to serve smaller pbf on z15,16 for a mobile web application(mapbox-gl-js).
Tried mbutil from https://github.com/mapbox/mbutil, extracted all pbf files into folder 0~14, what to do next?  I assume it should be possible to generate z15,16 from z14 given nothing is changed, only getting smaller tile.
File: metadata.json
The mbtiles file is extracted into 0~14 folders +  metadata.json


texts: When opening more than two images in GIMP version 2.8.14 system gets hanged. Even if I change the preference folder path the problem persists. How can I keep the system from hanging?


texts: R in QGIS not working after change of computers

texts:
I have installed, uninstalled a couple of times, but still get the same issue.  Windows 7 and QGIS 3.6.1 and latest versions if R and R Studio.  It used to work perfect until I changed computers.
R execution console output
Warning in install.packages("sf", dependencies = TRUE) :
'lib = "C:/Program Files/R/R-4.0.3/library"' is not writable
Error in install.packages("sf", dependencies = TRUE) :
unable to install packages
Calls: tryCatch ... tryCatchList -> tryCatchOne -> <Anonymous> -> install.packages
Execution halted


texts: Extracting certain fields for selected feature from point shapefile file with multiple features using ModelBuilder?

texts:
I have a point shapefile that is filled with many features.

If I have 1 feature selected I need to be able to export certain fields of that selected feature (ID, lat, long) into a dbf file.

What tools in ModelBuilder can accomplish this?


texts: Comparing polygon attributes that share the same parcel?

texts:
I'm new to ArcGIS, or at least anything beyond the basics, especially ArcToolBox. so the answer might be extremely simple.

Basically I have one layer of building footprints, and another layer of parcel data.

Many residential parcels contain multiple buildings, one of which is the house (coded in field BTYPE as 401) and one or more accessory buildings (BTYPE 999). each building footprint has a height attribute.

What I'd like to do is highlight parcels that contain a 401 footprint and a 999 footprint, but only the parcels where the height of the 999 footprint is higher than the height of the 401 footprint.

I have no idea how to do this.. I have tried a spatial join between the two layers but all that gets me is a count of the footprints on the parcel. I don't need the exact answer but a rough outline of steps that I need to take to accomplish this would be welcome.


texts:
I have to find the ratio between the first raster layer in folder R:\prova\VH and the first raster layer in folder R:\prova\VV, then I have to repeat the operation for the second two layers belonging to the two folders and so on.

The layers are sorted alphabetically and have the same dimension and CRS.

Each folder contains the same number of raster files.

I wrote this script in Python for QGIS 2.18 but it doesn't sort any result.

import os, os.path
import processing
from qgis.core import *
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry

raster1_filepath = "R:\\prova\\VH"
raster2_filepath = "R:\\prova\\VV"

outputDir = "R:\\prova\\Out"

for i in os.listdir(raster1_filepath):

    layercount1 = os.path.join(raster1_filepath, i)
    layercount2 = os.path.join(raster2_filepath, i)

    raster_layer1 = QgsRasterLayer(layercount1, 'raster')
    raster_layer2 = QgsRasterLayer(layercount2, 'raster')

    entries = []
    # Define numerator
    ras1 = QgsRasterCalculatorEntry()
    ras1.ref = 'raster_layer1@1'
    ras1.raster = raster_layer1
    ras1.bandNumber = 1
    entries.append( ras1 )

    # Define denominator
    ras2 = QgsRasterCalculatorEntry()
    ras2.ref = 'raster_layer2@1'
    ras2.raster = raster_layer2
    ras2.bandNumber = 1
    entries.append( ras2 )

    # Process calculation with input extent and resolution
    calc = QgsRasterCalculator( 'raster_layer1@1 / raster_layer2@1', \
    outputDir, \
    'GTiff', \
    raster_layer2.extent(), \
    raster_layer2.width(), \
    raster_layer2.height(),
    entries )
    calc.processCalculation()


texts: Adding extra characters to elements in a list

texts:
I am trying to create a script that automates the generation of a MXD. It does this by the user copying and pasting several variables, running a definition query and exporting an MXD. I can get the script to work with hard coded values as in the code below:

#vari 1
A = "'A', 'B', 'C', 'D'"

#vari 2
B = "'1', '2', '3', '4'"

#Query 1
Query1 = "Field1 = '17' AND NAME in ("+A+")"
#Query 2
Query2 = "Field1 = '1' AND UNIQUE_ID in ("+B+")"

#get list of layers in the mxd
layers = arcpy.mapping.ListLayers(mxd1)

#for the specified list of layers in the mxd, update the definition queries with the current variables
for layer in layers:
    if layer.name == "layer1":
        layer.definitionQuery = Query1
    if layer.name == "layer2":
        layer.definitionQuery = Query2

#save a copy of the template mxd do the work request folder
mxd1.saveACopy(workFolder + "\\H_" + jobNumber + "_01.mxd")

#clear the lock on the mxd
del mxd1

The problem I have is that the user enters the data as A, B, C, D (copy and pastes from a report) and I need to transform that for the definition query into "'A', 'B', 'C'" (adding the apostrophes). I can do this for one value (A = "'" + A + "'") but struggling to work out how to format with multiple values. I have tried:

#add apostrophes to list
string1 = "'"
AA = [string1 + x + string1 for x in A]
BB = [string1 + x + string1 for x in B]

but that didn't work. Put simply I would like to add extra characters to the elements in the list (A, B, C, D.... into "'A', 'B', 'C',...") I would welcome any pointers on how I can resolve this?


texts:
How do I use python to find raster layers data type
The raster data type for a band like float32 or int16?
I know you can use rasterlayer.dataType(1) but how do I then lookup returned integer (eg 6)  and get the data type?
I'm currently using QGIS 3.10


texts: which Coordinate system should I use for 3D analysis?

texts: I am creating 3D surface below the earth's surface using ArcGIS and then I need to find distances to the plane from the surface. Which coordinate system shall I use? 3D interfaces like IProximity3D and IArea3D look like they work only with Projected coordinate systems which are 2-dimensional. I am new to this area and need direction to decide and use.


texts: Canging coordinate system units from mm to m using ArcMap

texts:
I've received a feature dataset, that has a wrong amount of digits/wrong unit as coordinates. The coordinate system is actually a variation of Gauß-Krüger, which coordinates normally look like in this example:

E 33570272.656, N 5652121.859 (in meter)

Unfortunately the feature dataset has coordinates with too many digits (millimeter instead of meter):

E 33570272656, N 5652121859 (here in millimeter!)

I've already tried changing the linear units of my coordinate system to millimeter, which had no effect. Maybe someone knows a solution for changing the coordinates to meter?

The parameters of my coordinate system look like this:

I'm using ArcGIS Basic 10.6 and I don't have any programming skills


texts: I have added a dxf file in QGIS and saved it as shp file. the problem it that the area of the the polygon in shp. file is 1% more than the area of polygon in dxf. just to double check i have exported the same shp. file into dxf. and check the area. the difference still exists. can someone please explain.


texts: Apply mask and compute statistics with GDAL

texts:
I've constructed a VRT that includes a band of data and a mask (0/1 representing water/land) with:

ds = gdal.Translate(masked_vrt, combined_vrt, bandList=[1], maskBand=2, format='VRT')

Where combined VRT has band 1 of data and band 2 of 0/1.

I would like to use the python GDAL bindings with this VRT to compute the statistics of the data after applying the mask but ds=gdal.Open('masked.vrt').GetRasterBand(1).GetStatistics(0,1)
gives only the stats of the original, maskless, data (band 1).

I can get the result I'm looking for like this:

ds = gdal.Open('masked.vrt')
data = ds.GetRasterBand(1).ReadAsArray()
mask = ds.GetRasterBand(1).GetMaskBand().ReadAsArray()
masked_data = data[mask>0].mean()

But reading in the arrays of course loses the efficiency of using the VRT.

Is there some way to use GDAL directly to apply the mask and compute the stats?

Here is a MWE:

import numpy as np
from osgeo import gdal, gdal_array
data1 = np.ones([5, 5]) * 2
data2 = np.ones([5, 5]) * 5
mask1 = np.zeros([5, 5])
mask2 = np.ones([5, 5])
data  = np.vstack([data1, data2])
mask  = np.vstack([mask1, mask2])

## dummy geotransform
gt = [-78.38, 2.7e-04, 0, 37.33, 0, -2.7e-04]

data_ds = gdal_array.OpenArray(data)
data_ds.SetGeoTransform(gt)

mask_ds = gdal_array.OpenArray(mask)
mask_ds.SetGeoTransform(gt)

## write datasets
gdal.Translate('data', data_ds, outputType=gdal.GDT_Float64)
gdal.BuildVRT('data.vrt', 'data')
gdal.Translate('mask', mask_ds, outputType=gdal.GDT_Float64)
gdal.BuildVRT('mask.vrt', 'mask')
gdal.BuildVRT('combined.vrt', ['data.vrt', 'mask.vrt'], separate=True)

ds = gdal.Translate('applied_mask.vrt', 'combined.vrt', bandList=[1], maskBand=2, format='VRT')
ds1 = gdal.Open('data.vrt', gdal.GA_ReadOnly)
print ('Original Statistics:', ds1.GetRasterBand(1).GetStatistics(0, 1))
print ('Masked Statistics:', ds.GetRasterBand(1).GetStatistics(0,1))
ds = None; ds1 = None

If working, original and masked would differ, of course.


texts: Visitors to city's website can digitize points?

texts:
I'm trying to collect input on a revision to my town's mainstreet.  Is it possible to create a map where a visitor can click a point, and answer a few questions about that point?  (OK I was told that this question is too broad, so I am going to try to refine it)

Essentially, I want to be able to embed this map on the city's website, then invite community members to visit it and place markers at points on Main Street that they either like/dislike.  Ideally they would be able to drag a red or a green marker onto the map, which prompts one or two questions about the site, then move on.

I expect this could be done many different ways, but I really don't even know where to start looking for a solution.


texts:
I have a DEM built from LiDAR data and the file extension is .tif. When I open this file in ArcMap it seems to be a shaded relief DEM (see the screenshot below).

It happens I want to visualize the DEM colored by elevation. How should I proceed?


texts: Postgis | Problems Reprojecting a Point from UTM ED50 to GEOGRAPHIC ETRS89

texts:
I've got a POINTS PostGIS table in the 23030 (ED50 UTM 30N) SRID.

I'm having a lot of problems trying to reproject those points from ED50 UTM coordinate system to ETRS89 Geographic (lat/long) coordinate system. The aim is to convert the points for Google Maps.

I take an example point for my tests:

x: 499462.1
y: 4799340.56

1.- If i just transform the point from UTM (23030) to Geographic lat/long (4230) using the ST_TRANSFORM("GEOMETRY",4230) PostGIS function. I get the following coordinates:

long: -3.00663639334493
lat: 43.3461345057409

You can see in the picture where it is located in Google Maps. It has sense to be located in a different place to the original point, because although being in Geographic projection system it is not in ETRS89/WGS84.

2.- So I reproject the original point first to ETRS89 (nearly the same as Google's WGS84) using ST_TRANSFORM("GEOMETRY",25830) and I get the following point:

x=499462.124763252
y=4799253.63249274

Ok, if I load that point in a GIS, I can see there is a difference between the previous point which is ok as I have changed the SRID.

3.-I thought that if I reprojected the ETRS89 UTM point to ETRS89 lat/long it would be enough for getting the coords valid for Google Maps. The surprise came when I got the same point of the step1!!! (just doing the 4230 transformation)

I used a ST_TRANSFORM(ST_TRANSFORM("GEOMETRY",25830),4258) UPDATE query applied to the original point.

I don't really know how to solve this mess, I'm really confused.


texts: How to set the coordinates at which Cesium.Viewer is centered?

texts:
In the official CesiumJS tutorial there is a sample HelloWorld application, which shows a map viewer.

When the user opens the page, I want the globe to be centered above a certain point and with a certain zoom:

How can I achieve this?

I assumed that it can be configured via Cesium.Viewer constructor parameters (like timeline below), but didn't find them in the reference.


texts: Adding a GeoTiff created in QGIS to Map Box

texts: I have created a heat map in QGIS and have successfully exported it as a Geotif file. I exported the file as a Rendered Image to include my classifications and colors. However when I pull the data into Tile Mill or Map Box, all of the pixels with "No Data" get colored as white or black. I have tired setting the "No Data" to be 100% transparent in QGIS but this did not change anything. I want the white/black space to be excluded.


texts:
I'm trying to process some OSM data extracted via Osmium. I've managed to process some of my smaller datasets, but as I've begun working on the bigger ones, I'm receiving some odd errors. Specifically, I've tried to batch-process SUMLINELENGTHS and to create an OD Matrix using QNEAT3. Oddly enough, the sumlinelengths tool works on individual layers, just not when applied as a batch.
Batch SUMLINELENGTHS Error:
Traceback (most recent call last):
File "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\algs\qgis\SumLines.py", line 119, in processAlgorithm
QgsFeatureRequest().setSubsetOfAttributes([]).setDestinationCrs(poly_source.sourceCrs(), context.transformContext())), feedback)
Exception: unknown

OD Matrix error:
Traceback (most recent call last):
File "C:/Users/dicki/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\QNEAT3\algs\OdMatrixFromPointsAsLines.py", line 207, in processAlgorithm
net = Qneat3Network(network, points, strategy, directionFieldName, forwardValue, backwardValue, bothValue, defaultDirection, analysisCrs, speedFieldName, defaultSpeed, tolerance, feedback)
File "C:/Users/dicki/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\QNEAT3\Qneat3Framework.py", line 130, in __init__
self.list_tiedPoints = self.director.makeGraph(self.builder, self.list_input_points, self.feedback)
Exception: unknown

Execution failed after 53.70 seconds

Looking at my data, I see nothing about this set that should make it any different than the other sets which were processed successfully. I've never run into errors that looked like this before. I'm not even sure if these two error sets are related, and the unhelpful exception: unknown tag doesn't help me pry them apart. Any advice?


texts: Positioning leaflet divIcon when adjusting size using css

texts:
I am using the L.divIcon when using iconSize:null which is causing iconAnchor:[0,0] to not be recognized,my code looks like this:

  var divIcon =L.divIcon({
    className:'current-location-icon',
    html:'<i class="fa fa-location-arrow fa-2x"></i>',
    iconAnchor:[0,0],
    iconSize:null,
    popupAnchor:[0,0]
});

How can I ensure that the icon is visible exactly at [0,0] when adjusting size using CSS?


texts:
My ultimate goal is to convert landcover raster (.tif) objects to an sf object representing the raster's grid and the original values of each cell within each geometry. I have been able to do this for smaller rasters doing the following:
library(sf)
library(stars)

# import raster using stars
landcover_stars <- read_stars(my_raster.tif)

# convert to sf object using st_as_sf
landcover_grid_sf <- st_as_sf(landcover_stars)

In larger rasters (e.g. my largest raster is currently 11482x12607 cells), however, the read_stars() function imports the raster as a "stars proxy", which is a step taken to handle large raster datasets by the package.  While stars proxy objects are not accepted by the st_as_sf function, it is possible to set "proxy = FALSE" in the function. If I do this in my largest dataset, however, running st_as_sf(landcover_stars) with the resulting object will crash my laptop {16 GB RAM, i7 2.70GHz processor}.
Is there a way I can proceed to ease the load on my machine when converting very large star objects to sf?
In addition - could it be that it is actually the newly generated sf object what is depleting my machine?
Here is a dummy raster in case youd like to test it, with integer values randomly generated ranging from 1 to 10:
raster(nrows=12000, ncols=12000, xmn=0, xmx=10, vals = floor(runif(12000*12000, min=0, max=11)))


texts: Setting Symbol Template for ArcGIS Pro GraduatedSymbolsRenderer from ArcPy

texts:
Is it possible to set a Symbol Template for a GraduatedSymbolsRenderer using ArcPy for ArcGIS Pro?

The reason I ask is that running the test code below gives an error of:

Traceback (most recent call last):
  File "C:\Temp\test.py", line 10, in <module>
    symTemp = sym.renderer.symbolTemplate
  File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\arcobjects\_base.py", line 90, in _get
    return convertArcObjectToPythonObject(getattr(self._arc_object, attr_name))
RuntimeError

To run the code I created a test project and added one layer into it which had the ne_10m_populated_places.shp from Natural Earth as its data source.

This is the test code, and it is very close to a code sample from the GraduatedSymbolsRenderer documentation:

import arcpy
aprx = arcpy.mp.ArcGISProject(
    r"C:\Temp\Projects\TestProject\TestProject.aprx")
mapx = aprx.listMaps("Map")[0]
lyr = mapx.listLayers("ne_10m_populated_places")[0]

sym = lyr.symbology
sym.updateRenderer('GraduatedSymbolsRenderer')

symTemp = sym.renderer.symbolTemplate         # Line giving error
symTemp.applySymbolFromGallery("Circle 1")
symTemp.color = {'RGB' : [38,115,0,0]}
sym.renderer.updateSymbolTemplate(symTemp)

sym.renderer.classificationField = "POP2020"
sym.renderer.breakCount = 10
sym.renderer.minimumSymbolSize = 4
sym.renderer.maximumSymbolSize = 16

lyr.symbology  = sym
aprx.save()

I ran my test using ArcGIS Pro 2.4.1 on Windows 10 Pro for Workstations.


texts: Password protect QGIS project file (.qgz)

texts:
I want to password protect my QGIS Project file (.qgz).

Goal: Without password no one can open or modify it. Maybe a read/write or just write protection.

I can't find an option in QGIS. Any idea? Did I missed the option in settings?


texts: Ghost images (composite images) in Google Earth?

texts:
I notice that Google Earth images often have ghost cars or ghost boats. For example:

So, obviously multiple images from different times are being combined. What is the reason for this?


texts: Error "generate features" with ArcScan vectorization

texts:
Windows 10, ArcGIS 10.4.1/ArcScene 10.4.1.5686

I have a geotiff that I pulled in from a pdf. I'm in an edit session of ArcScan and when I Show Preview, the output looks great. However, whenever I try to Generate Features, I receive the following error:

  The map does not contain an editable polyline layer. Please add one and try the command again.

Any suggestions to proceed with ArcScan?


texts: What can I use to calculate drive times /with/ traffic considered?

texts:
I am looking to obtain drive times for a large list of addresses and want the results to take traffic conditions into consideration. For example, drive times from a distribution centre to a store at various times throughout the day.

Google's free API doesn't take traffic into consideration and the business license which does is too costly an option.

Would ArcGIS's API supply this information and would it be easily accessible with ArcPy? I think its World Traffic Service covers drive times has the data I'm looking for.

Are there any other options I may have missed?


texts: View ArcSDE feature classes of published map service, when outside of corporate domain

texts:
I have a notorious problem, which keeps bugging me since days and I can´t find ANY proper answer on the net...

In my company we´re using ArcGIS Server 10.1 and publish map services throughout our organization. We set up a bunch of ArcSDE geoDB on SQL, which are used to build the published MXD´s. It has the distinct advantage that we only have to maintain the ArcSDE´s, no re-publishing necessary after feature updates. And this works all fine as long as the client application is connected to our domain.

The problem starts when external clients want to view map services, which contain feature classes from these ArcSDE. Even though the map service is shared with 'everyone' and security changed to 'Public', it doesn´t work.

Is it possible to configure ArcSDE DB in a way that they are publicly available (outside the domain) or is this a problem of security settings on the ArcGIS Server side? Read a bunch of confusing things on ArcGIS SOC permissions, but really didn´t get that...

And I refuse to believe that this is not possible.


texts: Is there a good GPS track analysis library?

texts:
I keep writing custom code to work with GPS tracks: GPX files and NMEA logs from my phone, dedicated GPS logger, etc. I feel like I'm reinventing the wheel and am wondering if there's some standard library / tool I can use to work with GPS tracks. Some of the things I like to do:

Calculate the distance travelled on the track
Calculate the average speed and instantaneous speeds along the track
Calculate time in motion vs. time stationary
Simplify the track, eliminate stationary segments
Segment a track into separate tracks, one per trip
Smooth out GPS sampling error
Create an elevation profile by looking up points in a DEM
Convert the tracks to GPX, KML, etc for visualization

There are a variety of applications that analyze GPS tracks. I'm looking for a software library or scriptable set of tools so I can write my own applications. I don't much care about input formats; any sort of timestamped lat/lon points is fine. I mostly work in Python on Unix but pretty much any language will do.

Existing tools I use include GPSBabel and GDAL/OGR for conversion and PostGIS for working with the track geometry as a line string. But mostly I end up writing my own Python code that works on arrays of lat/lon pairs and that feels too low level. Also some of the tasks above (like removing GPS error) are quite subtle. Is there some existing library out there I should be using?


texts: Are these 9.2 CASE tools instructions still valid using ArcGIS 10.1?

texts:
The Using CASE tools in ArcGIS 10 document links to some older documentation at Building geodatabases with case tools, which specifically mentions ArcGIS 9.2 throughout.

Are these instructions still valid for creating a geodatabase using CASE tools? Is there a more up-to-date set of instructions which I should follow instead?

thanks


texts: Extra file generated while writing exporting collection to csv

texts:
I am exporting csv file from simple feature collection with help of gt-csv the feature collection is coming from db. Export is working fine but i am seeing another temp file generated  in the same folder with my table name as prefix and secondly is there any option available in this library so that i can put double quotes on string values.

 CSVDataStore csvStore = null;

    try
    {
          SimpleFeatureCollection coll = getCollectionDB();
        String geometryPropertyName;
        geometryPropertyName = coll.getSchema().getGeometryDescriptor().getLocalName();
        File fileToSave = new File("Export.csv");

        Map<String, Object> cvsParams = new HashMap<>();
        cvsParams.put(CSVDataStoreFactory.FILE_PARAM.key, fileToSave.getAbsolutePath());
        cvsParams.put(CSVDataStoreFactory.STRATEGYP.key, CSVDataStoreFactory.WKT_STRATEGY);
        cvsParams.put(CSVDataStoreFactory.WKTP.key, geometryPropertyName);
        csvStore = (CSVDataStore) DataStoreFinder.getDataStore(cvsParams);

        csvStore.createSchema(coll.getSchema());

        CSVFeatureStore csvFeatureStore = (CSVFeatureStore) csvStore.getFeatureSource(csvStore.getTypeName().getLocalPart());

        csvFeatureStore.addFeatures(coll);

    }
    finally {
        if(csvStore != null)
            csvStore.dispose();
    }


texts: Is there a size limit when adding features to ArcGIS Server programatically using POST?

texts:
I'm inserting new points using the ArcGIS Server REST API addFeatures method - this is working correctly.
However, I'm finding that when the JSON reaches a certain (undefined) threshold the POST fails with a 404. If the JSON is reduced in size the (otherwise identical) POST succeeds.
I can't see any inherent limits in Python's request functionality. Is this a known limitation with the ArcGIS Server REST API, or are there any workarounds?
Below is some code to reproduce the problem, using a sample Esri server so you can test it.
The target field has a width of 60 so it's expected that the insertion will return an error message with a string 1000 characters long. But note that it errors with a 404 when the string is 3000 characters long.
import json, requests

for num in [10, 1000, 3000]:

    # The following works with 10, gives a legitimate error with 1000, and gives a 404 with 3000
    value = "a" * num

    data = [{"attributes": {'address': value}, "geometry": {"x": 0,"y": 0 }}]
    params = {"f": "json","features": json.dumps(data)}
    headers = {'content-type':'application/json','Accept':'application/json'}
    response = requests.post("http://sampleserver3.arcgisonline.com/ArcGIS/rest/services/SanFrancisco/311Incidents/FeatureServer/0/addFeatures", headers=headers, params=params)
    print("With number = " + str(num) + " the response is " + response.text)

Results:

With number = 10 the response is
{"addResults":[{"objectId":18209901,"globalId" :
null,"success":true}]}
With number = 1000 the response is
{"error":{"code":400,"message":"Unable to complete
operation.","details":["Unable to add features"]}}
With number = 3000 the response is  404 - File or directory not found.

Why does posting the long JSON string return a 404, and how can this be resolved?


texts: Convert ArcGIS Javascript api Geometry Object to GeoJson

texts:
Possible Duplicate:
  How to convert ArcGIS Server JSON to GeoJSON?

How to convert ArcGIS Javascript api Geometry Object to GeoJson using Javascript?


texts: Draw lines from points to points not by lon/lat in mapbasic

texts:
How to draw lines from points to points not by lon/lat?

What I need is a cross connection with cable like the image below:
I have N input and M output points


texts: How to show single title for double CSS stroke in map legend

texts:
I am using CSS for styling roads in my geoserver. Here is my css :

 /* @title HSAMB Roads*/
    * {
       stroke: #333333,#00CC00;
       stroke-width: 3px, 1px;
       z-index: 0, 1;
       label: [name_of_ro];
       font-family: Arial;
       font-weight: bold;
       font-fill: black;
       font-size: 10;
       halo-color: white;
       halo-radius: 1;
       -gt-label-follow-line: true;
       -gt-label-max-angle-delta: 90;
       -gt-label-max-displacement: 400;
       -gt-label-repeat: 150;
     }

Everything is working fine but the client want the legend to be single line legend not the double line legend i.e inner line and outer line :

Here is my legend
Please Tell me how to combine both black and green line to single in map legend so that they are not shown as two separate line for a single road type.


texts:
I was accessing USGS map services at http://isse.cr.usgs.gov/ArcGIS/services/Combined/SDDS_Imagery/MapServer/WMSServer, through a program.

This link no longer works. I found a redirect - that took me to http://viewer.nationalmap.gov/example/services.html.

What would be a replacement site for the one I was using ?


texts: How to convert Zipcode+4 into standard 5 digit zip code using Field Calculator

texts:
I have a ZCTA shapefile that I would like to join to a table by the zipcode. My problem is that some of the zip codes in the table have a 5+4 digit format, and the ZCTA shapefile has only the 5 digit format. I would like to remove the last 4 digits of the 9 digit zipcodes using field calculator. My current method was to use slicing in python, but the code will not run. This is what im currently using: I created a new field in the table ('Field9') of type text and made it equal to the zipcode field 'ZIP' and used field calculator with this code

Field9=
!ZIP!

for i in Field9:
    if len(i)>5:
        i = i[0:4]

Why wont this compile? I am fairly new to using python on arc.


texts: Programmatically setting CRS for processing algorithm with a memory output layer

texts:
On QGIS 2.18.1, Python 2.7.12, Ubuntu 64 bit, I am struggling with setting a CRS for a memory output layer within an iteration over several years of observations for a protected area (PA) polygon layer (UNEP-WCMC World PA database). The solution I have found so far (see, e.g. this) builds upon generating a multipart vector layer for a sub-selection of PA per year as a memory layer in order to intersect with GADM country polygons. This, however, requires a manual input of selecting the CRS for the memory layer. Is there a way of doing this automatically?

The code that requires solution:

from qgis.core import *
from PyQt4.QtCore import *

# singlepart to multipart algorithm
processing.runandload('qgis:singlepartstomultipart', my_PA_vl, "my_field_of_country_ISO3", "memory:my_output_layer_name") # for this step I would like to set CRS for mem layer automatically
my_output_layer_name = QgsMapLayerRegistry.instance().mapLayersByName("Multipart")[0]

Setting the output layer to None - which should set the memory layer CRS to input layer CRS (?) - somehow does not generate any output.

Potentially, there might even be a better solution / approach for what I am intending...


texts: How can I save a QGIS project with hyperlinks?

texts:
I have already created an hyperlink project using actions, it works, when I click on one point of my shapefile it opens a photo, but when I re-open the project, the action is gone and I have to create it again :( I tried to open the project from qgis and with double clicking direct on the project.

What can I do to mantain the action? Am I missing an step?


texts: How to solve the error "ERROR 160196: Invalid Topology" when running the split tool

texts:
I am trying to split a multipart features using Split (Analysis) tool. When launching the process I come across the error message "ERROR 160196: Invalid Topology. Failed to execute split".
After some reading (this one for example: https://support.esri.com/en/technical-article/000006932) , I run the "check geometry" tool. I successfully fixed the identified error (incorrect ring ordering).
I tried again the Split (Analysis) tool but the module still failed. I still get the error "ERROR 160196: Invalid Topology. Failed to execute split". Is there another approach to solve it?
I am using ArcGIS Pro 2.7


texts: How to query a specific race/age profile from Census data?

texts:
I want to count Asian males, 15-30 years old by Census Tract. I can get the data sorted, but I don't know how to make a Boolean type of query (i.e., race is Asian AND sex is male AND age is 15-30). I currently have a table of all Asian, all female, all male aged 15-30, but this is not useful.

I am using ESRI Arcmap Business Analyst to pull ACS data to tracts.

Data shapefile


texts: Processing on files output by qgis:reprojectlayer have the original CRS from before reproject tool was used

texts:
I am writing a QGIS standalone Python command line tool. I check that all input shapefiles have a valid CRS with layer.crs().isValid() and then I reproject any files that are not already in ESPG:4326. One input file is originally in EPSG:3338. So, I call:

processing.runalg('qgis:reprojectlayer', input.shp, 'ESPG:4326', reprojectOut.shp)

the output 'reprojectOut.shp' has the correct CRS (4326) and correct .prj file and also displays properly even with 'On The Fly' disabled. So all seems well... until I try to extract a few features from the reprojected shapefile with:

processing.runalg('qgis:extractbyexpression', reprojectOut.shp, expression, extractOut.shp)

and somehow this output 'extractOut.shp' contains the correct features but has gone back to the original 'EPSG:3338' CRS. Is there something I am missing? Is there a reference to the original input file hiding somewhere? I have verified the filepaths are correct and tried passing both direct filepaths and also QgsVectorLayer() objects as input parameters. I am running QGIS 2.18.13 on Mac OSX High Sierra.


texts: Finding and snapping nearest point to polyline/point

texts:
I'm new to Python programming in ArcGIS, and I'm currently trying to write a script which I think is beyond me right now.

Script should take each point from first layer and:

Snap this point to nearest point from second layer (with given
distance)
If in given radius there is no point, only line, it should snap to
this line
If there are two points in radius it should snap to the closest one
If distance between points from first and second layer is 0, no
snapping is needed

It's hard for me to come up with the right solution. I thought of doing for loop which will iterate through points in first layer and for every point it should check the distance to points/lines in second layer and then snap to it. But the problem may be a huge amount of points in both layers, and maybe there is a faster and easier way of doing it.

Situation where point should snap to nearest point

Situation where point should snap to nearest line if there is no point in distance

EDIT:

I know that I could do it with tools from ArcToolbox but the problem is that it has to be done with a script which can work with hundreds of points in both layers.

So far in my code there is nearly nothing. Just some thoughts of doing it using cursor and for loop.

import arcpy

arcpy.env.workspace =
layer1 = path to layer1
layer2 = path to layer1

cursor = arcpy.SearchCursor(layer1)
cursor2 = arcpy.SearchCursor(layer2)

tol = 0,3

for row in cursor:
##getting 1 point from first layer

    for row in cursor2:
    ##  getting all points from layer 2

## measure distance between these two points
dd =
    if dd <= tol

## the same again but with other point from layer1


texts: Getting notified when ArcGIS REST map layer is updated?

texts:
I am trying to figure out how to get notified when a map layer is updated on this county government map server.

I am not a developer, so I am in way over my head, but I figured out the various map layers here. I want to know when layer 14 is updated but can't figure out how to do that without downloading the whole map and using a md5 checksum to see if it changed.

Is there a command I can send the API that will tell me when the layer has last been updated?

Ideally, I would like to create a script that checks every 30 minutes or so and can notify/email me when it has changed.


texts:
can anybody help me with setting a default value for a python toolbox parameter:

param1 = arcpy.Parameter(
            displayName="Exportgebiet",
            name="in_area_set",
            datatype="GPExtent",
            parameterType="Required",
            direction="Input")

What I want to do ist set the default value to the extent of the currently loaded map document.

I tried

    def updateParameters(self, parameters):
    ...
      mxd = arcpy.mapping.MapDocument("CURRENT")
      parameters[1].value = mxd.activeDataFrame.extent
    ...

but it did not work.

The dropdown box for the parameter should be "As Specified Below" or "Same as Display".

What am I doing wrong?


texts:
I'm using Leaflet API and geojson. I defined a function which zooms to the feature when clicked:

function zoomToFeature(e) {
            map.fitBounds(e.target.getBounds());
        }

 function onEachFeature(feature, layer) {
                layer.on({
                    mouseover: highlightFeature,
                    mouseout: resetHighlight,
                    click: zoomToFeature
                });
            }

But, I don't need to zoom on the feature, but to go to the URL. So, when user click on polygon on the map, it takes him to some url.

Is that possible and how? I thought I can maybe define some function which will lead to the url, and then pass it to layer.on click event, but not sure how.

Thank you for any help.


texts: Count the total number of trees from DSM drone data

texts:
We have a DSM, DTM and ortho data processed from a drone survey with a resolution of 5cm.

Are there any methods to "automatically calculate" the total number of estimated trees inside a dense vegetation forest area using Global mapper, QGIS, or GRASS software and mark it as a vector point?


texts: Finding Pseudo Nodes using ArcGIS for Desktop?

texts:
Trying to write a python scrip to find and dissolve pseudo nodes for a large dataset.

Anyone have any code to perform such a task?

For example, I have a road line with several end/start nodes (pseudo nodes) created by the cartographers.

Easy enough to run a topology but I have upwards of 4,000,000 +/- pseudo node to get rid of.

Software: ArcGIS 10.2.2
Python 2.7


texts: Selecting multiple vertices by dragging box doesn't work

texts:
I'm doing a QGIS course at the moment using QGIS 3.10 on a Windows 10 Pro 64bit system at the moment.

Settings are default, except for the "filter invalid geometries" option is set to "ignore".

With a shape selected and editing+vertex tool active, I can select a single vertex (shift+click) and edit or delete it. I can also select all vertices on a line by pressing shift+r and then selecting the first and last vertex.

I cannot, however, select multiple vertices by dragging a box. I have tried for 2 hours now, and I just don't get why it doesn't work. Has anyone else seen this problem with QGIS 3.10, and has anyone solved it?


texts: SpatiaLite / QspatiaLite: get subgroup row ids

texts:
I'm using QspatiaLite plugin to join polygon ids to points that are within those polygons. Now I' d like to get subgroup rownumbers for each subgroup of points with same polygon id! Like so..

grp_ID | sbgr_rowID | geom
--------------------------
1        1            ...
1        2
--------------------------
2        1
2        2
2        3
--------------------------
3        1
..

Here's my query:

   SELECT
     plgns.ID AS grp_ID,
     -- ?subgroup row number? AS sbgr_rowID --,
     pts.Geometry
   FROM pts JOIN plgns
   ON CONTAINS (plgns.geometry, pts.geometry)
   ORDER BY grp_ID, X(pts.geometry)


texts: Where to execute sdetable command in Windows?

texts: I am trying to fire some SDE Command in ArcGIS Desktop 10.1 but I am not aware where to execute the command in windows operating system.


texts: How to extract 1 band image (Boundaries) in Google Earth Engine Code Editor

texts:
I am quite new in remote sensing and trying to get boundaries of fields to build a prediction model in R. I tried Canny and Hough algorithms but the results are not satisfying. Lastly, I tried the model in 'Segmentation Session EEUS18' by Noel Gorelick. Here is my problem:

Gorelick's code: code.earthengine.google.com/?accept_repo=users/gorelick/examples
Example images are in https://drive.google.com/open?id=1bFVgedyVfzwXCyp6ZVXkCp_KkRl3Uq66 -

Suppose that I got perimeter pixels after the snic algorithm as [PerimeterZoomOut]. I want to extract the lines as vectors in a shapefile.

I tried reduceToVector but it needs at least 2 bands whereas perimeter image has only 1 band.

I also want to add a constraint so that it would not extract unnecessary lines as [PerimeterZoomIn].

The difference is also as between [NDVIZoomIn] and [NDVIZoomOut]. I don't want the noise in the fields.

Are there any way of extraction NDVI and perimeter lines images as vectors? They both have only 1 band.


texts: Iterating in a points geodataframe and then check if point is inside polygon

texts:
I'm trying to make a spatial verification.
I have a geodataframe with 245 points, and I'm having a bad time trying to check if they're inside a polygon or not.
I've been attempting this code, but I just won't work. As far as I understood, geodataframes are slightly different from dataframes when in a for loop.
for x in event_log_c.geometry:
    i = event_log_c.iloc[x].within(hexa)
    print(i)

After running this, it brings me this error:
TypeError: Cannot index by location index with a non-integer key

hexa is the polygon geodataframe, and event_log_c is the 245 points geodataframe.
How can I make a spatial verification iterating a geodataframe? Are there any way of doing this?


texts: Finding nearest point A to Point B and then excluding all point Bs without nearest point A?

texts:
I'm using GIS for the first time for a university project, ArcGIS specifically. As far as I'm aware I have access to the full suite of Arc software but I'm typing this out from home so can't check right now.

I've got two datasets; one is all of the kingfisher sightings uploaded to a biodiversity database in 2015 (~18000 data points) and the other is a list of all the water bodies in England surveyed by the Environment Agency in 2015 and a rating for each (poor, moderate, excellent) based on various factors (~800 data points).

Firstly, I'm not bothered about how far the kingfisher sightings are from the water bodies in terms of m or km etc. What I want to do is find the nearest water body to each kingfisher sighting and then separate out all of the water bodies that aren't the nearest water body for any kingfishers. Basically I want to know which of the water bodies on my list aren't used by kingfishers.

I've looked at some other forum posts and from what I've seen I think what I want to do is possible, but I haven't seen any that want to extract data about the points that aren't linked.


texts: I have 3D matrix file where XYZ are lon,lat and height. I plot these XYZ value using scatter3. Now I want to plot profile of height (along a line in km) of the data where xaxis:distance along line in km and yaxis:Z value . I tried mapprofile and improfile but could not able to get what I want.


texts: Creating PZ-90.02 and PZ-90.11 geodetic datum for QGIS

texts:
I would like to plot maps based on PZ-90.02 and PZ-90.11 in Russia, but in QGIS there is only PZ-90 to choose.

Can anyone explain how to create a proj4 record for these datums?
For PZ-90 proj4 looks like this:

  +proj=longlat +a=6378136 +b=6356751.361745712 +towgs84=0,0,1.5,-0,-0,0.076,0 +no_defs

but I have no knowledge about these values
towgs84=0,0,1.5,-0,-0,0.076,0 +no_defs


texts: Can't find plugins button on Ubuntu 16.04 - QGIS Las Palmas

texts:
I recently installed QGIS 2.18 'Las Palmas' on my Ubuntu system. I have located the "plugin toolbar" but can't seem to find the actual plugins button.

In this plugins toolbar, the only button present seems to be the python console.

I believe I installed everything properly but I suppose I may have missed some download.


texts: How do I cycle through the available overlays and base layers using OpenLayers

texts:
I would like to have my map cycle through the display of loaded overlays and base layers using JavaScript or jQuery.  I would like to do this with a delay of 5 seconds between cycling.

Does anyone have any idea on how to do this?


texts: Importing map information with osm2pgrouting?

texts:
I'm importing the cartography of Spain with osm2pgrouting.

I want to import address information, that is, the tag addr = *.
I modified the file mapconfig.xml follows:

  <type name="junction" id="4">
    <class name="roundabout" id="401" />
  </type>
  <type name="addr:housenumber" id="5"></type>
  <type name="addr:street" id="6"></type>
  <type name="addr:flats" id="7"></type>
  <type name="addr:city" id="8"></type>
  <type name="addr:country" id="9"></type>
</configuration>

I did not get any result.


texts: Value Error with Float/Integer/double Fields in Python?

texts:
When running this Code:

import arcpy
from arcpy import env
import os

#Specify the map document and the data frame
mxd = arcpy.mapping.MapDocument(r"G:\\PythonCode\\mxd_test.mxd")
env.workspace = r"G:\\LoopTest.gdb"
#output file folder
outPath = r"G:\Mxd_MapAutomation\\Test4\\"

FCs = arcpy.ListFeatureClasses()

# loop through all Featurecalsses in env.workspace
for FCLASS in FCs:

    #Skip all Text. OID and Geometry Fields - unfortunately WID, Shape Area etv is still included
    AllFields = [f.name for f in arcpy.ListFields(FCLASS)]
    LeaveOutField = [f.name for f in arcpy.ListFields(FCLASS,field_type='String')]
    LeaveOutField.extend([f.name for f in arcpy.ListFields(FCLASS,field_type='Geometry')])
    LeaveOutField.extend([f.name for f in arcpy.ListFields(FCLASS,field_type='OID')])
    attrlist = [x for x in AllFields if x not in LeaveOutField]

    # layer which has to change the attribute field
    flayer = arcpy.mapping.ListLayers(mxd, "*")
    featurelayer = flayer[7]

    # text element that has to change (e.g. title...)
    title = arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT", "map_title")[0]
    for attr in attrlist:

        # change value field
        featurelayer.symbology.valueField = attr

        # change tilte
        title.text = attr

        # update view
        arcpy.RefreshActiveView()
        arcpy.RefreshTOC()

        # export as jpeg
        outJPEG = os.path.join(outPath, FCLASS + "_" + attr + ".jpg")
        arcpy.mapping.ExportToJPEG(mxd, outJPEG, "PAGE_LAYOUT", 200, 200, 300, False, "24-BIT_TRUE_COLOR", 90, False)

del mxd

I get following error:

    Traceback (most recent call last):
  File "G:/DDP_loop_Code.py", line 35, in <module>
    featurelayer.symbology.valueField = attr
  File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\arcobjects\_base.py", line 87, in _set
    return setattr(self._arc_object, attr_name, cval(val))
ValueError: hhfishing

and I have no idea why.. the code basically runs through all FC in the GDB. It does it correctly with the first FeatureClass and stops with the error in the second Feature Class. Even if I would put in other Feature Classes the Error would appear in the same way. All the fields are either float or integer and the field where the error appears is not different to any other field where the code works on..


texts: Geoserver WPS - Tiled Barnes Surface with geospatial index

texts:
I'm using Barnes Surface to compute a set of points from my PostgreSQL db. The problem is that my whole table is processed for each tile.

I added a geospatial index on my points, how can I tell geoserver to compute only the points present in the bounds of the tile?


texts: Fixed distance buffer version 2.16.1 problems

texts: I am trying to buffer a simple line and the program simply won't complete the process after 10 attempts. I am buffering to 20,000 (I assume it's metres). Anyone had this issue before? My internet is fine and I restarted QGIS too. I am new to the programme.


texts:
I'm tasked with creating a .pyt to distribute across our office and field computers (state geological survey) for converting .gpx to features but constraining the choice of projections/datums to seven (the point being that such functionality will be more efficient). The tool must be able to write locally to a specified gdb. I'm not sure what best practices would be:  create a file in the specified .gdb with the .prj files and point to that or make the .prj variables global. Have not had success with portability either way (new to .pyt, I can fashion a script tool that works fine but is more trouble for the user than its worth). I seek insight or sample code.

Or is there a better way?