texts: Where is feature cache in ArcMap for ArcObjects to access?

texts:
It seems that ArcMap keeps a cache for every layer which contains features data and every time wants to refresh layer or apply a change, uses this cache. As I found, when ArcMap wants to draw a layer from the scratch, sends this query to geodatabase:

select ObjectId,Shape from "featureClassName" [where "layer definition query expression"]

In editing mode, when using snapping, it seems the features data in this cache are used to determine snapping state while mouse movement. I don't know how ArcMap stores spatial data in this cache, in ESRI shape format or as WKB, But I need to retrive these features data.

Is there any way to retrive feature data of a layer using its cache(layer cache) instead of using IFeatureClass::get_Feature method?
.....Maybe my idea isn't right and a layer doesn't have any cache (I don't mean cache which is optional for layers and boosts the speed of layer draw), but I'm almost sure spatial data of layers(map) is stored in ArcMap memory somewhere and ArcMap uses it for its operation like snapping during edit operation. Is there any way to access this part of ArcMap data(if really exists)?

I don't want to setup a cache, I want to retrieve data from ArcMap cache.ISpatialCacheManager doesn't offer any method for getting a specific data from it. Let me discuss what I mean in more details.

When I put this code inside my program

//C#

IFeature ipF=ipMyFeatureClassObject.GetFeature(10);

ArcMap issues a select command to geodatabase, let's suppose:

select * from MyFeatureClassTableName Where ObjectId=10

and gives back a record of data.

In edit mode, ArcMap can't use this method to get spatial data and calculate snapping state for features on the map because of time consuming this method has and snapping calculation needs a high speed data retrieve, also because mouse can moves every where on the map and snapping is enabled for any zoom, all these witnesses shows that spatial data must be available somewhere in memory and I'm looking for a method to access it.

If all these consumptions not come true, so there must be a perfect method that ArcMap uses to retrieve spatial data from geodatabase on the fly (for example 4 or 6 level spatial indexing based on the map extent and use the mouse position to retrieve only part of spatial data while mouse movement very fast, in this way only has access to geodatabase when the smallest indexing area of mouse changes).

I'm dealing with developing a tool like "Sketch Tool and "Edit tool" for a data source other than ArcObject data source and I don't want to setup a cache parallel with ArcMap. If possible I want to use same data which ArcMap stores in memory for this data source.

Now, do you guys think these are true assumptions and if yes is there any way to access ArcMap spatial data in its memory?


texts: Access QGIS menu functions copy paste style using python

texts:
Following question which was never answered:
https://stackoverflow.com/questions/27140520/how-to-access-qgis-menu-items-using-python

Is there any way to access functions in QGIS menu using python? Particularly in my case copy and paste style but best would be to get some general way / documentation where to find what function / tool is behind menu item and if / how is possible to get to it using python?

In my case because I frequently need to copy paste one or few layers from one project to other I made simple QGIS plugin for it. It can quickly copy layer path to clipboard (also raster/vector type and name from layers panel) and paste layer (using shortcut keys ctrl+shift+c / ctrl+shift+v). It works also between two QGIS windows which is exactly what I needed. Except one thing - it does not copy layer style from project. And because copy paste style is already as part of QGIS I would like to use that as step of copy paste layer without need to try to replicate the functionality in my own code.


texts: How can I include a map viewer into my web server from my ArcGIS Server?

texts:
I have a web server, what contains my website. I would like to put a map viewer into my site, what would connected to my ArcGIS Server. How can I do this in the easiest way?
The ArcGIS Server and my web server are different servers.


texts: PostGIS 2.3 - Union of lines if not starting/ending on a point table

texts:
I have a table with some lines (that are in fact multi-linestrings, even if it is not really needed, but that's another problem) that sometimes start and end on the location of a point table, sometimes not, like if what could be a single line what split in different lines. Here is how it could look like, I purposely drew the lines in different colors so we can see how each line is.

I would like to get the couple of points that are connected by the lines, and in order to do so I suppose the way to go is to first get a simpler table that union the lines in order to have only lines between the points, it would give something like that :

The way I see it, there are 3 different cases:

the line start and end on a point, so there is nothing to be done
the lines start or end on a point, and touches another line, in that case these two lines must be unioned
a line doesn't nor start nor end on a point, but touches two other lines that possibly do start or end on a point (even though I have some cases where there are more than one line with no starting or end on a point, but are on the same segment)

I try to conceptualize how to do so, I am guessing:

it would be possible to use plpgsql to iterate through the lines and create a new table with the proper geometry. But maybe there is way to do I don't know about with the PostGIS functions.
it would be needed, through some CTE, to first detect how many points a line intersects (or touches, but st_intersects is faster than st_touches isn't it?), so the three cases are segregated
it will be needed to consider that if the touching point of two lines are on a point, then this point should be excluded, if not all the lines would be merged regardless of the situation (in the second picture, the red and blue lines would just be one line since the middle point is a junction between two lines in the second case).

But from there I cannot figure the way to go, especially that it is possible that lines are intersecting without need to be considered (they intersects geometrically but not in the reality). Here is the piece of script I put together, but since I use a ST_Intersection, it also considers these not needed intersections, and anyway I am stuck with only possibility to merge a maximum of two lines:

WITH points AS (
SELECT id, geom, st_buffer(point.geom, 2) buff FROM point),
WITH lines AS (
SELECT id, st_linemerge(line.geom) geom, count(st_intersection(st_linemerge(line.geom),points.buff)) ngeom FROM line
LEFT JOIN points ON st_intersects(st_linemerge(line.geom), points.buff),
GROUP BY id, st_linemerge(line.geom))
SELECT a.id, b.id, st_union(a.geom, b.geom) FROM lines a
LEFT JOIN lines b ON st_touches(a.geom, b.geom) AND a.id != b.id
LEFT JOIN points ON st_intersects(a.geom, points.geom)
WHERE a.ngeom < 2 AND b.ngeom < 2 AND a.id < b.id
AND st_intersects(a.geom, b.geom)
AND st_intersects(st_intersection(a.geom, b.geom), points.buff) IS FALSE

Once again, my final need is to get the couples of points linked by lines, so in this example I would get 3 couples (I am saying this because it is always possible that my conceptualization of the problem is not the one to go).

EDIT: I looked into the approach proposed by GHRF which propose to merge all lines and split the result by the points.

The thing is that some of my lines are crossing each other, ending in two things:

The result see the lines cut not only on the points geometry, but also every time a line crosses a line
Merging the lines gives a network as a result, and cutting by the points then can only results in some kind of multilines, which nullyfies the process, given my need


texts: Getting error in finding spatial reference of shapefile

texts:
I am working on Ubuntu, when I tried to read spatial reference of shapefile it give me this error.

spatialRef = layer.GetSpatialRef().ExportToProj4()

  Traceback (most recent call last):

  File "< stdin >", line 1, in < module >

  AttributeError: 'NoneType' object has no attribute 'GetSpatialRef'

Help me to remove it


texts: AttributeError when I execute a python script

texts:
When I execute a python script  to assign my database in postgres, I get this Error !!

"AttributeError 'psycopg2._psycopg.cursor' object attribute 'execute' is read-only "

def ma_fonction(d):
    # connexion à la base de données
    conn = psycopg2.connect("dbname = 'test3' user = 'postgres' host = 'localhost' password = 'postgres'")
    cur = conn.cursor()
    for n in d :
        cur.execute = ("""
        PDATE douar SET distance_route = %s,  distance_piste = %s, distance_ecole = %s, distance_college = %s, distance_sante = %s, distance_souk=%s WHERE id =%s
        """,
        d[n]['route'] / 1000.0,
        d[n]['piste'] / 1000.0 ,
        d[n]['ecole'] / 1000.0,
        d[n]['college'] / 1000.0,
        d[n]['sante'] / 1000.0,
        d[n]['souk'] / 1000.0
        )
    conn.commit()
    print(commit)

Any solution plz !!


texts: Using Geotools or GeoserverManager to add a layer using just an URL

texts:
Lets say I have just an URL to a Layer in OpenLayers format:

http://www.geoservicos.ibge.gov.br/geoserver/CREN/wms?service=WMS&version=1.1.0&request=GetMap&layers=HidrogeologiaRegiaoNE

Is it possible to create a Layer in my Geoserver pointing to this external layer?

1) Using my Geoserver I've tried to use Web Feature Server (NG) but I don't want to list all layers before add it (GetCapabilities).

2) I need to do it using geotools. Again, GetCapabilities will be time spending because I already know the layer I want to add.

I can use Geotools or Geoserver Manager (REST API). I've searched everywhere for examples to add (publish) a WMS layer from a external geoserver layer but found nothing.

In other words, how can I do this using GeoserverManager API?
http://geoserver.geo-solutions.it/edu/en/adding_data/add_wfscascade.html


texts: How do I know the cause of my spatial error?

texts:
A (probably?) very basic question about spatial error it's about time I knew how to think about and answer correctly.

I have some data plotted from XY coordinates. I look at them and it turns out they appear in places where they shouldn't -- but only by a little bit. E.g. in the image below, the points appearing on land along the sides should be in the water.

The errors range from maybe 15 meters to 140 meters ... but usually about 30-50 m or so. How do I know why these are off and if there's any way to fix it? Is it:

1). Because I've incorrectly defined my coordinate system somewhere along the way? (I'm not sure if the original coordinate system of the points is known). But then, wouldn't they be astronomically off target, not by just 40 meters? Plus, it doesn't look like shifting them all in the same direction would fix the issue.

2). Because the resolution of the base layer doesn't capture the shoreline well enough and the points end up looking like they're on land? I've tried it with several sources for the background and none seem to match.

3). Error inherent to the GPS unit used to capture the points? I thought 100+ meters was a little high for anything above recreational-grade.

How do I know what's causing this? Have I done something wrong (#1) or is this sort of thing inevitable sometimes? How do I adjust for it to represent these points correctly, or verbally state an appropriate caveat about the accuracy?


texts: Export satellite images for an entire country from Google Earth Engine

texts:
How can I download all the satellite images covering an entire country(Kenya) for a period of years from Google Earth Engine? I drew a large geometry polygon over the entire country and then tried to get the images by using filterBounds() but doesnt seem to be working. What am I doing wrong?

var roi = /* color: #0b4a8b */ee.FeatureCollection(
    [ee.Feature(
        ee.Geometry.Polygon(
            [[[41.484375, -1.6477220517969224],
              [39.5068359375, 3.447624666646865],
              [35.2880859375, 4.981505049328245],
              [33.99169921875, 4.23685605976896],
              [34.98046875, 2.1088986592431382],
              [33.99169921875, -0.9447814006873896],
              [39.2431640625, -4.740675384778361]]]),
        {
          "system:index": "0"
        })]);

var image = ee.Image(ee.ImageCollection('MODIS/MOD09A1')
.filterBounds(roi)
.filterDate('2010-01-01', '2013-12-31')
.sort('CLOUD_COVER')
.first()
);
Map.addLayer(image, {bands: ['sur_refl_b01', 'sur_refl_b03',
'sur_refl_b04'], max: 0.5, gamma: 2}, "Modis Image");


texts: Adding variable number of layers with associated buffer distances with ArcPy for ArcGIS Pro

texts:
I am using ArcGIS Pro 2.7.2
I am not an experienced python programmer. My goal is that the user can select a various number of layers in the python tool and specify different buffer distances for each. I imagine it visually so that there are two columns:
Column A: Upload layer
Column B: Individual buffer distance in m
This process can be carried out until the user has specified all layers. This is important as it is currently not known how many layers are needed. Already thought about python dictionary but I do not know if this works for this case.
Finally, the buffered layers have to be loaded into a GDB.
import arcpy
from arcpy import env
env.workspace = "file_path"
env.overwriteOutput = True
fcs = arcpy.ListFeatureClasses()
bufdist = arcpy.ListFeatureClasses()
for fc in fcs:
    arcpy.Buffer_analysis(fc, OutputName, bufdist,"FULL", "ROUND", "LIST", "Distance")


texts: Finding Raster Pixel Neighbors using GDAL?

texts:
I'm new to GDAL Python, and I'm working on a little recursive algorithm which requires finding the neighbors of a a pixel on a DEM raster. I've been able to do this by creating a little PixelLite class which is essentially just a tuple representing the raster coordinates of a pixel. Here's a bit of its implementation (feel free to skip):

    # a lighweight version of a pixel for quick computation
    class PixelLite:

        def __init__(self, x, y, raster_meta):
            # set x and y
            self.x = x
            self.y = y

            # set max_x and max_y
            self.raster_meta = raster_meta
            self.max_x = raster_meta.x_max
            self.max_y = raster_meta.y_max

            # set error flag
            self.flag = False

            # the set of all neighbors
            self.neighbors = set()

            if self.x > self.max_x | self.x < 0 | self.y > self.max_y | self.y < 0: # something has gone terribly wrong
                self.flag = True
                return

            # !!! Below, I'm going to set out of bounds values to the value of the pixel itself, which will later be removed from the set !!!

            # neighbor columns
            if self.x > 0: # if left value will not be negative
                left = self.x - 1 # add legit left
            else:
                left = self.x # add left value that will be removed

            center = self.x

            if self.x < self.max_x: # if right value will not be larger than raster dimensions
                right = self.x + 1 # add legit right value
            else:
                right = self.x # add right value that will be removed

            # Note: bottom and top are arbitrary names, what matters is that values do no exceed bounds; it may be the case that, because of raster indexing 'bottom' actually represents the pixel above the current pixel
            # neighbor columns
            if self.y > 0: # if bottom value will not be negative
                bottom = self.y - 1 # add legit bottom value
            else:
                bottom = self.y # add bottom value that will be removed

            middle = self.y

            if self.y < self.max_y: # if top value will not exceed raster extent
                top = self.y + 1 # add legit top value
            else:
                top = self.y # add top value that will be removed

            # top row of neighbors
            self.neighbors.add((left, top))
            self.neighbors.add((center, top))
            self.neighbors.add((right, top))

            # middle row of neighbors
            self.neighbors.add((left, middle))
            self.neighbors.add((right, middle))

            # bottom row of neighbors
            self.neighbors.add((left, bottom))
            self.neighbors.add((center, bottom))
            self.neighbors.add((right, bottom))

            # the method I'm using will sometimes add the center, middle element as I've used that as a default to set values to if there is an out of bounds error. So let's just pop that off...; a pixel cannot be its own neighbor
            self.neighbors.discard((center, middle))

        ...

Basically, what this does is create the tuples representing each neighbor for a given pixel whenever the PixelLite class is instantiated. Although this is pretty lightweight, because it is implemented in very basic Python without any kind of binding, I would imagine it lacks the optimization of a more sophisticated solution such as I would expect to find in GDAL.

Is there a faster getPixelNeighbors function provided by GDAL or another raster library i.e. one which would allow me to do the same but which might have some neat little optimization tricks up its sleeve?


texts: Does ArcGIS Pro have an alternative to ArcMap's "Copy Selected" option for attribute tables?

texts:
I'm in the process of trying to switch fully from ArcMap to ArcGIS Pro (version 2.0). One small piece of functionality in ArcMap that I used very frequently was the "Copy Selected" option to copy selected rows from an attribute table into Excel. (See this answer if you don't know what I'm referring to.)

In ArcGIS Pro, I cannot find an analogous method to copy tabular data. The only method I can see is using the "Copy Rows" geoprocessing tool, which would majorly slow down my workflows.

Am I missing something?


texts: Find Closest Index in line string to insert new vertex using PostGIS

texts:
I want to find the index from Line String to insert new vertex on it.
I am using ST_Add Point function but it requires specified index for inserting vertex.How to get those index based on given input points from line string.
My Input Line String:
"LINE STRING Z (496328.772173402 4622876.01448448 208.4005,496327.883130287 4622814.67050952 208.4005,496494.578714408 4622814.67050952 208.4005,496497.898850613 4622875.89518039 208.4005)"
My Input Point:
'POINT Z(496458.12794668 4622814.67050952 208.4005)'
Expected Output:
"LINE STRING Z (496328.772173402 4622876.01448448 208.4005,496327.883130287 4622814.67050952 208.4005,496458.12794668 4622814.67050952 208.4005,496494.578714408 4622814.67050952 208.4005,496497.898850613 4622875.89518039 208.4005)"
I need to insert on 3rd position on given line string. But i am not able to
find the exact segment index for inserting new point.
I am using following query with manual index 2.But i need to get index automatically
update public.linetype
set geom=ST_AddPoint(geom,  ST_Make Point(496458.12794668 4622814.67050952 208.4005),2) where id=4
from public.linetype

public.linetype is my table name and geom is linestring geometry column.
PostGIS Version 2.4.4 and PostgreSQL 10.4


texts: Osm2po : Querying generated table to get the closet line from a point

texts:
I generate the OSM data of Africa postgis and I want to launch  queries for a list of  gps points to project them into the closest line.
I made this request to a single point but took around 3 seconds which is not acceptable because i'll repeat it for dozens of points.

SELECT osm_name,ST_Distance(st_geomfromtext('Point('|| 10.426925 || ' ' || 36.611049 ||')',4326),geom_way) AS dist
FROM af_2po_4pgr
WHERE geom_way && expand(st_geomfromtext('Point('|| 10.426925 || ' ' || 36.611049 ||')',4326), 0.01)
ORDER BY dist limit 1;

1) Is there a better approach to achieve what I wanted to do.
2) How to optimize this query.

I want the result illustrated in the first image.

But i have this one.


texts: Adding raster layer to map in ArcGIS Pro project using ArcPy not working?

texts:
I have a raster dataset in a geodatabase and I would like to add it to my map using arcpy. I have tried with the following lines of code but it doesn't work (no error thrown):

raster_data = r'D:\ArcGIS_Pro_Projects\Optimity_Alf\Optimity_Viewsheds.gdb\Accelerator_Viewshed_AGL_raster'
aprx = arcpy.mp.ArcGISProject(r'D:\ArcGIS_Pro_Projects\Optimity_Alf\Optimity_Alf.aprx')
map = aprx.listMaps('Optimity Map')[0]
map.addDataFromPath(raster_data)

Any ideas?


texts: Calibration annotation data for Sentinel-1 images downloaded from Sentinel Hub?

texts: I received some Sentinel-1 tiff images that were downloaded via the WCS service, so they're just the tiff files, not the full SAFE-format products.  They're for the IW GRD product, and my understanding from the documentation of that product is that to use the data, you're first supposed to radiometrically calibrate the digital numbers from the tiff files, to get eg backscattering coefficient sigma-naught.  The documentation says that the calibration vectors for doing this conversion are in the Annotation section of the product.  Is there a way to get that information from Sentinel Hub, or maybe to specify when doing the download that we want it to apply the calibration before giving us the file?


texts: Building custom form for QGIS in Python plugin?

texts:
I want a custom "QGIS 2.14.6 Essen" layer edit form in my own python plugin.

I use layer.setEditForm( NODE_ui_path ) and layer.setEditFormInit('NODE.formOpen')

Function NODE.formOpen does not work!

MyPlugin is in ~\.qgis2\python\plugins\MyPlugin

Ui Form is in ~\.qgis2\python\plugins\MyPlugin\NODE.ui

formOpen is in ~\.qgis2\python\plugins\MyPlugin\NODE.py


texts: Matching Trimble Data Dictionary File (.DDF) to File Geodatabase?

texts:
I'm implementing a new database model in our company's GIS (Utility). We collect field data with Trimble units. I'm having trouble figuring out how to create a Trimble data dictionary that matches the schema (fields and domains) of the database.

I want to easily be able to load data from the field into GIS with out much editing or modification.

Is there a way to export the schema from Esri's file geodatabase to a Trimble .ddf file?


texts: Creating points in interval 0.125 decimal degree using ArcPy?

texts:
I have a polygon shapefile. I want to create points at an interval of 0.125 decimal degrees in this polygon.

How can I create a script with python in ArcGIS 10 to automate the process of create points?

The start point is the south-west corner and the end point is north-east of the polygon.


texts:
I want to find features that are inside a buffer around a line. The idea is to find how far along the line the intersection with the housepolygons happens. I have tried the intersect tool with using a distance in XY Tolerance for 10 m (the buffer distance), but it does not return the desired output.

In the image above I want all the blue areas to get a point and to join it along the line using the "Locate features along routes tool". The projection is utm 33N wgs 84.


texts: How to update spatial point locations in ArcGIS

texts:
In ArcGIS, is there a way to quickly reflect changes made to XY coordinate values in my table?

In edit mode I can change the values in my X and Y columns but I'm not able to refresh the map. Currently I export the table, make changes in an editor (e.g. Excel) and re-"create a feature class from XY table" using my new source file.


texts: Export Web Map Task: Configure HTTP response headers

texts:
Is there any way to configure the http headers (response headers) of a geoprossessing task hosted by ArcGIS for Server?

I would like to make the export web map task, return an image which would be directly downloaded by my browser.


texts:
I am making several maps from a template with a PyQGIS script and I am having problems to align the scale bar on them. My code is:

# create scale bar
item = QgsComposerScaleBar(comp)
item.setStyle('Single Box') # Other possibilities are: 'Single Box', 'Double Box', 'Line Ticks Middle', 'Line Ticks Down', 'Line Ticks Up', 'Numeric'
item.setComposerMap(composerMap)
item.applyDefaultSize()
item.setAlignment('Right')
comp.addItem(item)

I've found the QgsComposerScaleBar.setAlignment() QGIS API Documentation, but when I introduce 'Right' as a parameter I get a type error message. I am using QGIS 2.6.0


texts: Carto trying to set syncing with a PostGIS database connection

texts:
I'm trying to set up a syncing connection between Carto and a PostGIS Table. However, I'm not seeing updates in Carto, and the table does not show up in my list of syncs when I make a Get request.

This is the code I'm using to set it up

curl.exe -v -H "Content-Type: application/json" -d '{\"connector\": {
 \"provider\": \"postgres\",
 \"connection\": {
   \"server\":\"[SERVERURL]\",
   \"database\":\"[DBNAME]\",
   \"username\":\"[USERNAME]\",
   \"password\":\"[PWD]\"
 },
 \"schema\": \"[SCHEMA]\",
 \"table\": \"[TABLE]\"},   \"interval\": 3600}' "https://[DOMAIN].carto.com/api/v1/imports/?api_key=[APIKEY]"

The response I get is

  "item_queue_id":"[QUEUEID]","success":true}* Connection #0 to host [DOMAIN].carto.com left intact

What am I doing wrong here?


texts:
I got a postGIS error like this,

ERROR:  transform: couldn't project point (100.496 13.7118 0): latitude or longitude exceeded limits (-14)
SELECT ST_Length(ST_Transform(ST_GeomFromText(E'LINESTRING(100.495995129 13.7117836894,100.495962221169 13.7117761471941)',4326),27700));

This lat/lng comes from Bangok, Thailand and the 27700 is a british national grid. Do you think this casued the problem? May I mention that problems like this did not occur when I got lat/lng from Spain or USA.

Any input would be much appreciated.


texts: Null vs. blank in click vs. hover info window in CartoDB?

texts:
I've noticed that in my CartoDB info windows, if a given column does not have a value for a point, that lack of value shows up as a blank in the hover window, but it shows up as a null in the click window. Is CartoDB set up so that this intentionally happens? Is there a way to change this setting?

Edit: I also am working on a project with over a dozen columns, any of which could potentially be blank at a given time, so to go in and change the HTML for each one of them, although doable, is rather tedious.


texts: SQL Error: Latitude values must be between -90 and 90 degrees

texts:
I'm more of an ArcGIS user than a SQL Server user but I have polygon data in a 'geography' type in SQL 2008 R2.  For a couple of these poly layers, I get this persistent drawing error in ArcMap:

  One or more layers failed to draw:

  Cultural Resource Site Polygon:  Underlying DBMS error
  [[Microsoft][SQL Server Native Client 10.0][SQL Server]A .NET
  Framework error occurred during execution of user-defined routine or
  aggregate "geography":  System.FormatException: 24201: Latitude values
  must be between -90 and 90 degrees. System.FormatException:     at
  Microsoft.SqlServer.Types.GeographyValidator.ValidatePoint(Double x,
  Double y, Nullable1 z, Nullable1 m)    at
  Microsoft.SqlServer.Types.Validator.BeginFigure(Double x, Double y,
  Nullable1 z, Nullable1 m)    at Microsoft.SqlServer.Typ]
  [db.layer]

I thought at first this might be related to other posts I've seen on incorrect ring orientation in polygons but those errors tend to say that the data exceeds a hemisphere.  This error seems to suggest that I have vertices inside some of my polys that are perhaps backwards.

Unsure how to find them in SQL (any help would be great!), I turned to ArcMap and used the Feature Vertices to Points tool to create a set of points at every vertex of my polys.  I then calculated X and Y on those points.  I can't find any instance where my Y exceeds either 90 or -90.

Unsure at this point...


texts: arcpy fieldmappings: can't set output field name for fieldmap

texts:
I am testing out some code in the Python window of ArcGIS so that I can insert it into a standalone script later. I want to create a copy of a table and add a new field that is
the result of concatenating two of the existing fields using a space as a delimiter.

I am sort of following the suggestion below to use field mapping which is from this page http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//001700000065000000 because I ultimately want to use the JoinField command to join two tables based on two attributes:
"To merge two or more fields in the join table before the join, first export the table or feature class using the Table To Table tool and merge using the tool's field map."

My code does the concatenation or merging just fine, but I cannot get it to rename the output field. In my resulting table, the name of the new field is the same as the first field I add to the fieldmap (ie, "COC_NAME" in the code below.)

fieldMappings = arcpy.FieldMappings()
fieldmap_join = arcpy.FieldMap()
fieldmap_join.joinDelimiter = " "
fieldmap_join.mergeRule = "Join"
fieldmap_join.addInputField("COC_DATA", "COC_NAME")
fieldmap_join.addInputField("COC_DATA", "GRID_ID")
fieldjoin = fieldmap_join.outputField
fieldjoin.name = "JOINID"
fieldmap_join.outputField = fieldjoin
fieldMappings.addTable("COC_DATA")
fieldMappings.addFieldMap(fieldmap_join)
arcpy.TableToTable_conversion("COC_DATA", "C:\\HEA_GIS.mdb", "COC_MOD", "", fieldMappings)

Can someone point out what am I doing incorrectly or perhaps provide another solution?


texts:
I am trying to include an epsg_operations.properties definition into my Java /Geotools Project. I thought it would be the best solution to include the gsb-File for a GRID-Transformation. So I inserted the following into my epsg_operations.properties File:

25832,31467 = PARAM_MT [ "NTv2" , PARAMETER [ "Transformation XXX" ,  "XXX.gsb" ]]

I added the epsg_operations.properties File as well as the gsb-File to my resources root, als well as the org.geotools.referencing.epsg package within the resources.

Nothing helped so far. Any other suggestions are very welcome.


texts: How to summarize GIS skills on a resume?

texts:
I'm a recent graduate and I've used ArcGIS 10.0-10.2 (I just say 10.X) for a variety of tasks:

Analysis
Data Conversion
Georeferencing
Editing
Data Management (i.e., creating and organizing geodatabases in ArcCatalog)
General map creation

Does anyone have any suggestions on how best to summarize my GIS skills on my resume?


texts:
How can PostGIS shorten the end segments of polylines?

I already figured out how to calculate the length of the shortest end segment of any end of any polyline in my_table:

select min(least) from
    (select least
         (ST_Length
              (ST_MakeLine
                   (ST_PointN(geom, 1),
                    ST_PointN(geom, 2))),
          ST_Length
              (ST_MakeLine
                   (ST_PointN(geom, ST_NPoints(geom) - 1),
                    ST_PointN(geom, ST_NPoints(geom)))))
     from (select (ST_Dump(geom)).geom as geom from my_table) as polylines)
    as a;

Now I just want to shorten both end segments of every polyline with 1/3 of that value. (So that a polyline with a single segment will still consist of the middle 1/3 of the original.)


texts: After creating tiles in QGIS by enabling VRT, I'm getting black color tiles

texts: I have a big TIFF image and I created tiles of size 512 X 512 but all the tiles colour is black I don't know why.


texts: Select form the attribute table in with OGR

texts:
I want to select the (Sumava) from the attribute table and then make a shapefile. My problem is I don't know which function should I use to select this area.

If any one can help me?

import ogr
import os
direction_parks = 'ArcCR500_v32/ChranenaUzemi.shp'
direction_result = 'ArcCR500_v32/result.shp'

open_parks = ogr.Open(direction_parks)
layer_park = open_parks.GetLayer(0)

layer_park.SetAttributeFilter("???????")

number_park = layer_park.GetFeatureCount()

driver = ogr.GetDriverByName("ESRI shapefile")

if os.path.exists(direction_result):
    driver.DeleteDataSource(direction_result)

dataset = driver.CreateDataSource(direction_result)

new_layer = dataset.CreateLayer(direction_result, geom_type= ogr.wkbPolygon)

number = 0
for j in range(0, number_park):
    h = layer_park.GetNextFeature()
    h_geometry = h.GetGeometryRef()

    new_layer.CreateFeature(h)

dataset.Destroy()


texts: Can a Geopackage style be deleted through SQL window in DB Manager?

texts:
I am trying to clean up a Geopackage and as part of this I would like to delete some styles from the "layer_styles" table. I have tried doing this with the below query in the SQL window of the DB Manager in QGIS 3.4.3:

Delete from "layer_styles" where id = 48

When I run this I get the a database error popup window with the error message

  SQLite objects created in a thread can only be used in that same thread. The object was created in thread id 29468 and this is thread id 31144.

Does anyone know if my query is wrong or whether this is a limitation with the Geopackage format?


texts: Counting the number of pixels for each distinct class

texts:
I created a land cover transition map that has over 200 distinct transition classes. I would like to count the number of pixels for each distinct transition class within the region and I was wondering if there is an efficient way of deriving the number of pixel counts for each distinct class.
//------------------ Extracting MYR Boundary -----------------------------

//  Call ee.FeatureCollection (i.e. vector data) - MYR boundaries
var ssea = ee.FeatureCollection('users/jjohanness1992/ssea/ssea-adm0');

//  Filter ee.FeatureCollection only to Southeast Asia
var countrylist = ee.List(['KHM']);
var CAM = ssea.filter(ee.Filter.inList('GID_0',countrylist)); // Filters to Myanmar

Map.centerObject(CAM,6);

// 2001 land cover
var lc2001 = ee.Image('users/earthorbgis/rlcms_lc/rlcms_lc_2001').select('Mekonglandcover2001_lc').rename('lc').clip(CAM);

// 2018 land cover
var lc2018 = ee.Image('users/earthorbgis/rlcms_lc/rlcms_lc_2018').select('Mekonglandcover2018_lc').rename('lc').clip(CAM);

//  Confirm projection and nominal pixel resolution
print('Projection: ',lc2001.projection()); // EPSG:4326 // confirmed as correct
print('Scale: ',lc2001.projection().nominalScale()); // 30.000000000000004 m // confirmed as correct

//  Define colour palette for RLCMS landcover maps
var palette = [
//colorcode   // pixel value // landcover
'6f6f6f',     //      0      // Unknown
'aec3d4',     //      1      // Surface Water
'b1f9ff',     //      2      // Snow and Ice
'111149',     //      3      // Mangroves
'287463',     //      4      // Flooded forests
'152106',     //      5      // Deciduous forests
'c3aa69',     //      6      // Orchard or Plantation Forests
'7db087',     //      7      // Evergreen broadleaved forests
'387242',     //      8      // Mixed forests
'cc0013',     //      9      // Urban and Built Up
'8dc33b',     //      10     // Cropland
'ffff00',     //      11     // Rice
'cec2a5',     //      12     // Mining
'674c06',     //      13     // Barren
'3bc3b2',     //      14     // Wetlands
'f4a460',     //      15     // Grassland
'800080',     //      16     // Shrubland
'51768e'      //      17     // Aquaculture
];

Map.addLayer(lc2001,{palette:palette, min:0, max:17}, '2001 ori',false);
Map.addLayer(lc2018,{palette:palette, min:0, max:17}, '2018 ori', false);

// Making a transition map

var transition = lc2001.multiply(17).add(lc2018);
print(transition)

// see sample transition
//Map.addLayer(transition.eq(91),{},"DCD_OPF")

var stat = transition.eq(91).reduceRegion({
  reducer:ee.Reducer.count(),
  geometry:CAM,
  scale:30.000000000000004,
  maxPixels: 1e13
});

print(stat)

In my mind, I would have to run an ee.Reducer.count() for each class manually. However, I was hoping if there could be a method that tabulates all the counts for each class.


texts: Cannot get the temporal plugin in the plugin manager

texts: I am using QGIS 3.4.8 and I want to animate time series data using temporal controller but I couldn't get the temporal plugin in the plugin manager list. I couldn't install in the software. What might be my problem?


texts: How to get distinct values from multiple columns in GeoTools?

texts:
Is there a way to get distinct values from multiple columns in geotools?

I know I can get unique values from a column like this:

SimpleFeatureSource featureSource = dataStore.getFeatureSource(layer);
Query query = new Query(layer, Filter.INCLUDE, new String[] { "PROPERTY_A",
"PROPERTY_B", "PROPERTY_C" });
UniqueVisitor visitor = new UniqueVisitor("PROPERTY_A");
SimpleFeatureCollection collection = featureSource.getFeatures( query );
collection.accepts(visitor, null);
Set set = visitor.getResult()

But what I need is an equivalent to select distinct PROPERTY_A, PROPERTY_B, PROPERTY_C.

Edit:
I've tried to get unique values from the first column, and then create cql filters with the unique values from the first column to query unique values from the second column and so on, an then join it all in the end. But the performance, as expected, is very poor.


texts: Adding ArcGIS Feature Service connections to QGIS 3 using Python API?

texts:
I am trying to deploy multiple ArcGIS Feature Service connections to all of the QGIS 3 users in my workplace.

Is there a way to do it using the Python API?

At the moment I'm manually adding the connections into the QGIS3.ini, like this :

[qgis]
connections-arcgisfeatureserver\My%20Service%20Name\url=https://services.somewhere.com/arcgis/rest/services/My_Service_Name/MapServer
ARCGISFEATURESERVER\My%20Service%20Name\username=MyUsername
ARCGISFEATURESERVER\My%20Service%20Name\password=MyPassword
ARCGISFEATURESERVER\My%20Service%20Name\authcfg=

but I would prefer to add them using Python if possible.


texts:
When buffer analysis is performed on a polygon what determines the order of assignment for FID?
This came up with the following:

arcpy.Buffer_analysis(inFC, outFC_50, "50 Meters")
arcpy.MultipleRingBuffer_analysis(outFC_50, outFC_100, [100], "Meters", "Distance", "NONE", "OUTSIDE_ONLY")

During the first buffer a polygon goes from FID 5 to 9, not a big deal since there is the added ORIG_FID field. During the second buffer the ORIG_FID is of course the FID from outFC_50. I can join the FID from inFC to outFC_100 and work around it, but I'm curious why the FIDs are changed during a buffer where no dissolving occurs?


texts: How to compare polylines attribute from the same GRID cell in ArcGIS 10?

texts:
I need some help because I do not know how to model this issue:

I have two polyline shapefiles. They have an attribute that I need to be compared. These are roads from different sources. So, I have their names and attributes.

My intention is to make this comparison within each GRID cell, so, I created a GRID with the fishnet tool, and then "Intersected" each polyline with the GRID, I have something like this:

BUT, I think that it is important to use the "NAME" attribute, for a better comparison between elements with same name. I should need to distinguish polylines, like this example:

For example, If it compares the SE attribute with the one that crosses (it has SW), it gives me something wrong.

Maybe, another possibility is taking into account their bearing, or something related with its angle, because it happens to be, in most cases parallels.

Sorry, but I'm a bit lost. I think it is something tricky, but I hope there is a possible methodology for this purpose. Any ideas?

---- Edit ----

I'm thinking about the bearing question. This is my idea:

Select all the GRID cell that intersect with my desired geo (the black one).
Export to a new layer the selected. So I only have grid cells spatially containing my geo and possibly the geo to contrast (red one).
Compute the bearing for each line Black and Red (previously intersected with the GRID cells).
Start a WHILE-clause loop with:

WHILE counter > 0

Count the number of Black features contained into the cell. "var counter"

Select a black feature.

Take value of the QUESTIONED attribute (for example NE);

IF the red geo has a bearing value +-5% from the blacks value THEN

take (in another var) the red questioned attribute (for example NE);
IF var1 == var2 THEN
write "SAME" into a String field from Red geo.
ELSE IF
Write "ERROR" into the same string field.

END IF

END IF

Counter = Counter - 1
Next feature

END WHILE

I have never implemented such a pseudocode, any help for this?
maybe is a good starting...

Anyone would adventure in that?


texts: ogr2ogr: Convert geojson to shapefile - additional features?

texts:
I have a geojson-file (geojsonfile.geojson) with points, lines and polygons.

I would like to convert the geojson-file to 3 shapefiles (points, lines and polygons) and keep all the attributes.

Therefore i used ogr2ogr, see this topic:

GeoJson to ESRI Shapefile using ogr2ogr

with the following commands:

ogr2ogr -nlt POINT -skipfailures points.shp geojsonfile.geojson OGRGeoJSON
ogr2ogr -nlt LINESTRING -skipfailures linestrings.shp geojsonfile.geojson OGRGeoJSON
ogr2ogr -nlt POLYGON -skipfailures polygons.shp geojsonfile.geojson OGRGeoJSON

the exported point shapefiles show the exact same amount of features (2142 in shapefile, 2142 in geojson), but the linestring and polygon feature count is different (linestring: 511 in shapefile, 132 in geojson line, polygon: 426 in shapefile, 379 in geojson polygon). The additional features are also visible.

How can a conversion from geojson to shapefile lead to additional features? Also some errors have been reported while converting, but I guess this is because I forced the SHP geometry type (explained in other topic)?

Any ideas?


texts: Load domain values from .xml when creating PostGIS database with ili2db from .ili

texts:
I'm trying to create a PostGIS database with ili2pg. This creates a proper database schema. There is an additional .xml file with domain-values. I think these values aren't loaded in the process. How can I load them to the database schema?

I'm using PostgreSQL 11 with PostGIS 2.5 and ili2pg 4.3.0 on Windows 10
ili-Model is "Minimales Geodatenmodell - Langsamverkehr" --> Hpm_Network_V1.ili with the additional Hpm_Catalogues_V1.xml
from https://www.astra.admin.ch/astra/de/home/fachleute/weitere-bereiche/geoinformation/geobasisdaten/langsamverkehr.html or Model Repository: https://models.geo.admin.ch/ASTRA/


texts:
Is there anyone who has a experience with relational operation in eCognition. I'd like to ask 2 questions related to eCognition.
1. How can I classify the pixel values between 0.23 and 0.46 as a class 1, the pixel values between 0.46 to 0.64 as a class 2 etc. How can I achieve this purpose in eCognition.
2. Is it possible to use two variables for classification of one class. For example, can I use given ranges of temperature and precipitation to classify class 1. How should I build rule set for this purpose?

Thanks a lot


texts: OSError: [Errno 24] Too many open files (gdalbuildvrt)

texts:
I tried to build a virtual raster from 2600 GeoTIFF's and got this message:

  OSError: [Errno 24] Too many open files: [...]

gdalbuildvrt -resolution average -r nearest -input_file_list /tmp/processing/buildvrtInputFiles.txt /home/map/Ar/virt.vrt

I tryed with 1000 files and got the same error. Using 250 files is working fine!
The files are in average 1.5MB. I got 16GB of RAM and running GDAL 2.2.3 on Linux Mint.


texts:
I made a map using tilemill and osm data. This map is indented for print, so it will be exported to pdf, so no zoom levels used in mss files.

I need to prioritise labels used on this map but can't figure out how. i need street names and subway stations to be always on map, and fill all unused space with amenities markers. Metro icons are made as shields

#spb-point [railway='subway_entrance']
  {
  shield-face-name:"Proxima Nova Regular";
  shield-unlock-image: true;
  shield-file: url(maki/svgs/subway-spb-15-3.svg);
  shield-name:[name];
  shield-size:10;
  shield-wrap-width:20;
  shield-halo-fill:#fff;
  shield-halo-radius:2;
  shield-text-dy:20;
  }

Amenities are made as markers

#spb-point [amenity='school']
  {
  marker-file: url(maki/fill/svgs/college-11.svg);
  }

Street names are made as lines with labels

#spb-line [highway='primary_link'],
#spb-line [highway='primary']
{

line-width:3;
line-color:#fff;
    ::label{
  text-name: [name];
  text-margin:100;
  text-size: 10;
  text-face-name: "Proxima Nova Regular";
  text-fill: #000;
  text-placement: line;
  text-halo-fill: #fff;
  text-halo-radius: 1;
  text-transform: uppercase;
  text-character-spacing:0.8;
  }

What I want to achieve is that street names and subway shields won't disappear as I export map to pdf.


texts: How to convert .tif to .pix or .gpr or open a .tif file in Focus

texts:
I was wondering if anyone knows how to convert a raster .tif file into some thing that PCI Geomatica (Focus) can use, such as a .pix or .gpr file?

I have tried to run it through FME, but I cannot seem to come up with a transformer that lets me do it.

Alternatively, does anyone know of a way to open a .tif file in Focus, as it already is in RGB?


texts: Spatial join of polygon to point not transferring attributes in ArcGIS for Desktop?

texts:
What I have:

1 point layer with adresses.
1 polygon layer with zones.

What I want: 1 point layer with the polygon info added.

I just did the spatial join, the columns are added but there's no info joined. Everything is NULL.
I Googled and there seemed to be a bug, but that was for ArcMap 10 and I have 10.3. I can't find any patches nor servicepacks.

What am I missing here?


texts: Remove or delete selected feature in openlayers

texts:
Using OpenLayers 3.9.0.
I have two vector layers which are both desplayed on the map and can select any feature with a singleclick var select = new ol.interaction.Select();

Now I need to make it possible to remove selected feature with Del key.

I ended up with creating the folloing function and listener

var deleteFeature = function(e){
    if(e.keyCode == 46 && (select !== null)){
        layer.getSource().removeFeature(select.getFeatures());
    }
};
document.addEventListener('keydown', deleteFeature, false);

but the problem is that I don't know the layer of the selected feature to call the removeFeture method on.

Any suggestions how to find out the layer of the selected feature or another way of removing selected features?


texts: Create a heatmap with R for printed media

texts:
I am looking for a way to represent density points using heatmaps in R for printed media (e.g. using Rmarkdown to generate a PDF).
The map should have a basemap (OpenStreetMap or Stammen) and instead of presenting all points I want them to be presented as a density heatmap overlayed.
I know there are several ways to create heatmaps in R, like leaflet.extras's addHeatmap() for HTML maps. Unfortunately, they can't be exported to PDF (or other formats different to html).
ggmap, on the contrary, is aimed for printed maps, but unfortunately, it has two main problems: 1) it requires google API key (I do not want to do do it, neither I want to use google services) 2) it is quite cumbersome, as it requires knowing bounding box beforehand and, since I have to do plenty of maps, I need an automated way to get that bounding box (ideally without requiring other packages, such as osmdata's getbb()).
tmap would be a great solution, as it allows re-using same code for printed and interactive maps, but I haven't found any documentation about generating heatmaps.
Is there any way to achieve what I am looking for?


texts: QGIS -Split vector layer tool not generating an output

texts:
I'm trying to split a vector layer into new layers based on one field. The field is string type but when I run the tool, a message saying 'Creation of field fid failed (OGR error: Wrong field type for fid)' comes up. Before running I created a folder to store the output layers.
Any idea about what I have to change or what the right type field is, please?
The output should be many layers as unique values the field have.


texts: Setting spatial reference to OSGB_1936_British_National_Grid in ArcGIS Desktop?

texts:
I have two datasets in ArcGIS Desktop. One is from the Ordnance Survey (UK National Mapping Agency) and, when I view the properties of the dataset it says that the Spatial Reference is

OSGB_1936_British_National_Grid

I want to set another file to the same spatial reference. I know I can use the Import function to import it from the first file, but how would I select this coordinate system from the built-in ArcGIS ones? I can find British National Grid under Projected/National Grids/Europe and I can find OSGB 1936 under Geographic/Europe, but I can't seem to find the combination of the two anywhere?

I'm not hugely experienced with coordinate systems, but I think that BNG is the co-ordinate system, and OSGB 1936 is the datum. Is that correct?

I know that when I select spatial references in ENVI I can choose a coordinate system of United Kingdom and a Datum of OSGB1936.  Is there a way to do this in ArcGIS Desktop?


texts: I'm running into problems running the r.neighbors GRASS command in QGIS. Although it will run for small neighborhoods as soon as I increase my neighborhood above a certain number it becomes unbearably slow. I don't understand why Focal Statistics on Arcmap will do this in the blink of an eye while R.Neighbors takes ages?


texts: Splitting data into geographic zones for multiprocessing

texts:
I have been looking at multiprocessing and tools like Find Duplicates in Big Data and they seem to simply split the data into sections based on rows/number of records. I am trying to do an intersect on datasets ranging from 30 million to 100 million records see this link for details on the project)  and for that the data (say for aspect, slope and vegetation) needs to be split into zones which are geographically the same. Any idea on how to do this?

Ideally in ArcGIS 10.


texts: How to Calculate length of polyline geometry for several tables in PostGIS?

texts:
PostGIS database contains several geometry (polyline) tables (with names "D1_r", "D2_r", "D3_r"). I calculate length for one "postgis layer" use request:

SELECT
sum(ST_Length_Spheroid(the_geom,'SPHEROID["WGS 84",6378137,298.257223563]'))/1000 AS km_roads
FROM
"D1_r";

How to make a request to get a result for each table ("D1_r", "D2_r", "D3_r") and the total result for all tables ("D1_r"+"D2_r"+"D3_r")?
Thanks!


texts: ArcGIS No font options for Character symbol

texts:
I can't get the font options to display when trying to use character symbol in ArcMap. This problem is local to my machine only. See screenshot below.

ArcMap 10.2.2 & Windows 7 Enterprise 64bit


texts: Features becomes bigger on zooming in Openlayers 2.13

texts:
i'm new in developing Openlayers applications.I want to add Points and LinearRing objects to the same layer. My Problem is, when I add a LinearRing and I am zomming in, the geometry gets bigger and bigger on every zoom level. But the LinearRing should have the same behaviour as a Point. When I am zooming, the point is automatically scaling so it becomes a little bit smaller and not bigger!

Here are some pictures to understand my problem:
First picutre show how it looks like without zoom:

Second picutre show how it looks like with zoom:
So as you can see, the LinerRing beomces very big.

I hope you can understand my problem and somebody can give me a good solution.


texts:
When I try to draw polygons why am I losing vertices? I'm delineating stream terraces on some large DEMs but have clipped the data frame to an area of around 3 X 8 km. I want to be able to create vertices that are as close as 5 meters on the ground.

I adjusted the snapping tolerance but that didn't work. Must be the x,y tolerance of my layers. Should have mentioned I'm using both 10m and 1m DEMs. Guessing then that the 10m DEM is simplifying my small polygons. Do I join/relate the 10m and 1m DEMs?


texts: Importing OpenStreetMap data in ArcGIS for Desktop?

texts: It is possible to import data in QGIS but ArcGIS has more features and easy to use, so is there any way to import OSM data and we can use ArcGIS tools for analysis and database creation?


texts: GEE – Smoothing curve of detected water area

texts:
I am trying to plot the time series chart of the detected water area by otsu's method in Sentinel-1 images as a smoothed average curve.

I found two approaches to this for the Google Earth Engine:

Using the Savitzky-Golay filter or
Smoothing via moving time window

Actually the first variant with Savitzky-Golay filter is very brilliant. Unfortunately the time scale from a certain point in time is considered here and I would use the original 'system:time_start' for the chart.

So I used the second variant with moving time window. The problem here is that the calculated average is far to high and the bad performance. When viewing a longer period of several months I get the following error 'Error generating chart: Computation timed out'.

Here is the code in Google Earth Engine:  https://code.earthengine.google.com/94d17969df4c88abc7b59dfd9ecda3a9

Can someone tell me where the problem is or does anyone know a better solution?

I would be very happy about a clue.


texts: Converting selected multipart features to single parts in QGIS

texts:
I'm looking for a way to convert one or many selected multipart feature to different single part features without recreating another layer, the way Arcgis allow a user to select a multipart feature and convert it to single parts while editing.

Please note I'm not looking for the existing processing QGIS tool doing this for the whole layer, but a tool who could do it on choosen entities while I'm editing the layer.


texts:
public static void SelectBy(this IActiveView view, ILayer2 targetLayer,
    string filterClause, esriSelectionResultEnum selectionResult ) {

    IFeatureSelection featureSelection = targetLayer as IFeatureSelection;
    if ( featureSelection == null ) return;

    IQueryFilter queryFilter = new QueryFilterClass();
    queryFilter.WhereClause = filterClause;

    IColor selectColor = featureSelection.SelectionColor;
    selectColor = Local.ColorToIColor( System.Drawing.Color.OrangeRed );
    featureSelection.SelectionColor = selectColor;
    view.PartialRefresh( esriViewDrawPhase.esriViewGeoSelection, view.FocusMap, null );

    featureSelection.SelectFeatures( queryFilter, selectionResult, false );
    ISelectionSet selectionSet = featureSelection.SelectionSet;
    view.PartialRefresh( esriViewDrawPhase.esriViewGeoSelection, view.FocusMap, null );
}

If I use this code to make a selection calling it with:

view.SelectBy( layer, "[KeyField] = 10000.0",
    esriSelectionResultEnum.esriSelectionResultNew);

My feature polygons are always FILLED. If I set the featureSelection.Symbol field to a SimpleLineSymbol (and set SetSelectionSymbol=true) it shows nothing. How do I do a feature selection that is NOT filled (as happens when I use IMap.SelectByShape()?)


texts: Passing list to viewparams variables in GeoServer SQL views?

texts:
I have a sql view which take only one parameter, but this parameter is a list, and a need to know how can i pass this parameter in the url

Here is my sql view

SELECT * FROM shop WHERE id IN %list%

When I put the default value for %list% as (1, 2 ,3) as example, the layer preview works and show the shopa with the default ids,
but when I try this in the url

viewparams=list:(1,2,3)

does not work.

I already tried other things like

viewparams=list:"1,2,3"

but nothing works.

Does anyone know how to do this?


texts: Loading CSV with X&Y and Save as Shapefile using PyQGIS?

texts:
I am attempting to transform a CSV file with joined fields of Eastings and Northings into a shapefile with PyQGIS (QGIS V3).

Here is where I've got to:

1: Loaded the postcode shapefile which contains postcode and X Y fields

2: Loaded the XLS file holding address/contact details which I will be converting to a shapefile

3: Run algorithm join attributes table which saves output as CSV

My next step(s):

4: Load text delimited layer (CSV) with X & Y and save as shapefile with CRS EPSG:27700

5: Load the resultant shapefile

My code for steps 1-3 is below.

# 1: loads the postcodes SHP file from V drive
from qgis.core import QgsVectorLayer, QgsProject
layer = QgsVectorLayer('V:/GIS - Files/3. Data/OS Data/UK Postcodes/UK Postcodes (2017).shp', 'Postcodes', "ogr")
QgsProject.instance().addMapLayer(layer)

# 2: loads the XLS spreadsheet from the server
from qgis.core import QgsVectorLayer, QgsProject
layer =     QgsVectorLayer('//xxxxsql/xxxxx/xxxxx/xxxxx/WeeklyContacts-2019-07-15.xls', 'Contacts', "ogr")
QgsProject.instance().addMapLayer(layer)

# 3. performs join field attribute algorithm, joining X and Y of 'Postcodes.shp' to 'Contacts.xls' and saves to 'Contacts.csv'
params = { 'DISCARD_NONMATCHING' : False, 'FIELD' : 'Postcode', 'FIELDS_TO_COPY' : ['Easting','Northing'], 'FIELD_2' : 'postcode', 'INPUT' : '//xxxxsql/xxxx/WeeklyAuditExtracts/xxxx/WeeklyContacts-2019-07-15.xls', 'INPUT_2' : 'V:/GIS - Files/3. Data/OS Data/UK Postcodes/UK Postcodes (2017).shp', 'METHOD' : 1, 'OUTPUT' : 'C:/Users/xxxx/OneDrivexxxx/Desktop/Temporary/WeeklyContacts-2019-07-15.csv', 'PREFIX' : '' }
processing.run("native:joinattributestable", params)

The steps of code where I speculate things to go wrong:

# 4. converts the joined CSV to SHP
uri = "C:/Users/xxxx/OneDrivexxxx/Desktop/Temporary/WeeklyContacts-2019-07-15.csv?delimiter=csv&xField=Easting&yField=Northing"
layer_csv = QgsVectorLayer(uri, "WeeklyContacts-2019-07-15", "delimitedtext")
QgsVectorFileWriter.writeAsVectorFormat(layer_csv, 'C:/Users/xxxx/OneDrivexxxx/Desktop/Temporary/WeeklyContacts-2019-07-15.shp', "UTF-8", layer.crs(), "ESRI Shapefile", layerOptions=['SHPT=POINT'])

# 5. loads the SHP file derived from the joined CSV
from qgis.core import QgsVectorLayer, QgsProject
layer = QgsVectorLayer('C:/Users/xxxx/OneDrivexxxx/Desktop/Temporary/WeeklyContacts-2019-07-15.shp', 'WeeklyContacts-2019-07-15', "ogr")
QgsProject.instance().addMapLayer(layer)

The Python console reads this successfully but generates an empty points shapefile without headers.


texts: ILayer to IFeatureLayer casting problem in VBA

texts:
I am facing problem with ILayer to IFeaturelayer casting in ArcObjects(VBA) after casting debugging time it is showing error message. Run-time error(Type Mismatch). error getting at ****

 For i = 0 To pMap.LayerCount - 1
     Set pLayer = pMap.Layer(i)
     MsgBox (pLayer.name)
     Dim pFlayer As IFeatureLayer
     **** Set pFlayer = pLayer ****
     If (Not pLayer.name = "ParcelNumbers") Then
         Set OtherFClass = pFlayer.FeatureClass
         InterCount = OtherFClass.FeatureCount(pslt)
         If (Not InterCount = 0) Then
             objTextStream.WriteLine pFeature.Class.AliasName & "," & pFeature.oid & "," & "Fuse Overlapping with other feature."
         End If
     End If
 Next i


texts: GeoServer on access antivirus scanning exclusions

texts:
The powers that be have decided to implement on-access virus scanning on all servers - i.e. a file is scanned each time it's accessed, unless it's on an exclusion list, and I've been asked to come up with the exclusion list - just ask the supplier, they say.
One of the servers is running GeoServer, and there is no supplier.
Files accessed at GeoServer start aren't really a problem. It's those that will be accessed during normal map serving, and which will cause response times to suffer that I'm worried about (I've got all the mapping data covered)
Has anyone any clue about where I could find the info?


texts:
I have a dataset with points in NYC. I extract the bounding box (max and min lat/lon) and plot it on a map. It works for the basemap stamen watercolor, but not for google's. Why?

Here is the code to create the bounding box polygon:

library(ggmap)

bbx <- c(left=-74.06077,bottom=40.63317,right=-73.76324,top=40.84902)
x <- c(bbx["left"], bbx["left"], bbx["right"], bbx["right"])
y <- c(bbx["bottom"], bbx["top"], bbx["top"], bbx["bottom"])
df <- data.frame(x, y)

Here is the code to plot the bbox on the stamen+watercolor map, and it works:

myMap <-get_map(location=bbx,source="stamen",maptype="watercolor", crop=FALSE)
ggmap(myMap)+geom_polygon(aes(x=x, y=y), data=df, fill="red", alpha=.5)

Here is the code to plot the bbox on the google+satellite map, and it doesn't work:

myMap2 <-get_map(location=bbx,source="google",maptype="satellite")
ggmap(myMap2)+geom_polygon(aes(x=x, y=y), data=df, fill="red", alpha=.5)


texts: Reducing height through ArcGIS raster calculator?

texts:
I am using ArcMap 10.4.1

I have a DTM (10x10)m, and a raster (rasterized from polygon) that I want to subtract the DTM with.

In the picture below, four areas have different values I want to reduce the height of the dtm with, and the beige area is NoData. I tried simply just subtracting the rasters, but ended up with a different result than desired.

All data is lost except where data overlays. So I wonder how to set up the formula to get the original DTM to have its original values, except from the areas where the four values from the new raster will reduce the height.


texts: QGIS Dissolve: Limit output file format to .shp only

texts:
I'm new to QGIS and recently I've been trying to dissolve two large shapefiles together after they have been merged. It's taking an awful lot of time for each operation and I'm trying to speed up the process.

I've noticed that, once the dissolve operation completes, the destination file has a mix of different data types, such as .shx, .dbf, .qpj, and of course .shp. I am wondering if I can make QGIS to generate .shp files only in a dissolve operation and, if I can, would doing so make the dissolve operation complete faster?


texts: Calculating angle between two points (according to reference point)?

texts: I need to calculate the angles between two points on ArcGIS. I have ArcInfo License, too. I have 5 points on a diagram. One of them is my reference point. The other four creates angles between each other. I need to calculate each angle automatically. Because I have 7500 of these points. I'm also adding a picture of my points. I need to calculate the angles of point number 1-2, 2-3, 3-4, 4-1.


texts: Need to put points on overlapping boundaries

texts:
I have two layers: a route (Line) and a county boundary (polygon).

I need to intersect them so that the output layer is a point layer. The ArcGIS Intersect tool places a point where the route passes the county boundary.

But there are some locations where the route goes along the county boundary and then enters the boundary polygon. In such a case, ArcGIS fails to place any intersection point.

Is it possible to place two intersection points in such cases, one where the common segment starts and the other where the common section ends?

Thanks a lot for your help.


texts: Putting different id in different paths on my Leaflet map

texts:
I have many polygons to display on map like this:
      L.geoJSON(polygon, {
        style: function (feature) {
          return {
            stroke: true,
            fillColor: feature.color,
            fillOpacity: opacity,
            color: feature.color,
            opacity: opacity,
            weight: 2,
          };
        },
      })
        .bindPopup(function (layer) {
          return layer.feature.properties.name;
        })
        .addTo(map);

and each polygon generate a <path> like this
<path class="leaflet-interactive" stroke="red" stroke-opacity="0" stroke-width="2" stroke-linecap="round" stroke-linejoin="round" fill="red" fill-opacity="0" fill-rule="evenodd" d="M477 171L469 171L468 176L455 175L455 171L459 171L461 167L465 167L467 163L471 163L472 160L479 153L488 148L493 147L494 158L489 159L487 163L477 167L477 171z"></path>

Is there a way that I can put a different id for each path?


texts: Add attachment to feature service feature issue

texts:
I have a feature service that contains a number of features for which I want to allow the attaching of files/images etc. The layer is editable and it has attachments enabled, both as a property on the layer and as an attachment table in arcmap.
The code to attach the file/image is as follows:

             layer.addAttachment(featureOBJECTID, document.getElementById("attachmentUploader"), callback, function (err) { });

The element being passed is the form:

<form id="attachmentUploader" style="display: inline-flex; margin-top:10px;">Add:&nbsp;&nbsp;<input type="file" name="attachment'" onchange="uploadAttachmentFeature()"><input type="hidden" name="f" value="json"></form>

This used to work up until around last October, however it had suddenly stopped working and I had found a thread that said that the issue was with arcgis server 10.2. The server was just upgraded to 10.6, however the issue is persisting with the following error:

Error: Error performing add attachment operation
at Object.e.load (init.js:1005)
at init.js:87
at c (init.js:103)
at d (init.js:103)
at b.Deferred.resolve.callback (init.js:105)
at c (init.js:104)
at d (init.js:103)
at b.Deferred.resolve.callback (init.js:105)
at init.js:92
at l (init.js:107)

An error is also being shown in the arcgis server logs saying that the attachment process was not successful.

Is there any solution for this issue?

EDIT - When accessing the add attachment service from within the url the file attaches fine from the add attachment window in the service itself.


texts: Repointing Datasource of layers using ArcObjects

texts:
I'm currently working on a ArcObjects C# project that will repoint datasources for layers of various types. Basically a database is maintained that holds the old datasource connections and the new connections. When the user presses a button, the application goes through all the current layers in ArcDesktop and compares the layers current datasource with the new one. If there is a difference, it resets the datasource.

However, I'm not quite sure how to access the datasource, and in my perusings of this stackexchange, I noticed someone who was working on something similar, but they were using raster layers. I'm not too familiar with ArcGIS and didn't know if I needed to handle feature layers differently as he mentions he already found a solution for feature layers.

Here's the section I'm trying to change in ArcDesktop:

I have everything working just fine except the changing of those datasources. I found this resource on arcgis, but it doesn't appear to change those strings, though I could be mistaken.

This is what I currently have:

//Loop through all layers and get their information
pEnumLayer.Reset();
pLayer = pEnumLayer.Next();
while (pLayer != null)
        {
            Log.Debug("Manipulating the layer: " + pLayer.Name);
            IDataLayer dataLayer = (IDataLayer)pLayer;
            var datasource = dataLayer.DataSourceName; // this is what I believe to be the current datasource

            Log.Debug("Repaired layer: " + pLayer.Name);
            pLayer = pEnumLayer.Next();
        }

Is that the correct datasource I'm looking for? If so, would all I need to do to change the datasource is to set that variable to the new datasource string and then redraw the map?


texts: Draw a polyline around a route to act as geofence with leaflet routing machine

texts:
I am using leaflet routine machine to get a route between 2 points. I want to be able to draw an invisible polyline around the route and use this as a geofence and if the user goes out of the bounds then reload another route.

Not quite sure how to go about this and the closest thing I have found is this

https://coderwall.com/p/zb_zdw/buffered-polyline

Where do I start?

Most examples use googles routine machine.


texts: Add date to an image collection to export to CSV

texts:
I calculate the number of fires for a region and export the result to Google Drive as a CSV file, but I can't add the date to the image collection to export, it gives the following error:

  ImageCollection (Error) Error in map(ID=2019001): Image.addBands,
  argument 'srcImg': Invalid type. Expected: Image.
  Actual: Date.

Here is the code:

// Load country shapefile
var lsib = ee.FeatureCollection("USDOS/LSIB_SIMPLE/2017");
var vn_shape = lsib.filterMetadata('country_na', 'equals', 'Argentina');
print (vn_shape);

// Load fire counts image
var fire = ee.ImageCollection('FIRMS')
             .filterBounds(vn_shape)
             .filterDate('2019', '2020')
print(fire)

var scale = fire.first().projection().nominalScale();

// Filter fire with more than 50% confidence and add a new band representing areas where confidence of fire > 50%
var filterConfidence = function(image) {
  var line_number = image.select('line_number');
  var confidence = image.select('confidence');
  var conf_50 = confidence.gt(50).rename('confidence_50');
  var count_band = line_number.updateMask(conf_50).rename('count');
  // HERE IS THE PROBLEM!
  var date = image.date();
  return image.addBands(count_band)
              .addBands(date);
};

var fire_conf = fire.map(filterConfidence);
print('fire_conf', fire_conf);

// Count for individual image.
var countIndividualImg = function(image) {
  var countObject = image.reduceRegion({
    reducer: ee.Reducer.countDistinct(),
    scale: scale,
    geometry: vn_shape,
    maxPixels: 1e9
  });
  return image.set(countObject);
};
var fire_ind_count = fire_conf.map(countIndividualImg);
print('fire_ind_count', fire_ind_count);
Export.table.toDrive({
  collection: fire_ind_count,
  description: 'FireCounts_gt50_conf',
  selectors:['count','date'],
  fileFormat: 'CSV'
});
print('Total fire count', fire_ind_count.aggregate_sum('count'));


texts: Labeling edges/borders of polygons using ArcMap?

texts:
How would I label polygon features, such as Counties, by its edge in ArcMap 10.1?

This is similar to what one might find in a street atlas.

For instance, have there are two counties (Lamb and Castro), at the same spot I want the labels to be inside along the boundary in the respective counties.


texts: How to exclude bounding boxes of plots in lascatalog

texts:
I have 96 1 ha (200m by 50 m) plots that I collected TLS data (each file separate). Most of the plots are close to each other and therefore they have overlapping bounding boxes (point clouds don't overlap themselves). I would like to calculate some grid metrics for all 96 plots using lidR package.
ctg <- readLAScatalog("infolder", filter = "-drop_z_above 1.3 -drop_class 2", chunk_size = 0,
                      chunk_buffer = 0)
opt_chunk_size(ctg) <- 0
opt_chunk_buffer(ctg) <- 0
opt_output_files(ctg) <- "outfolder/{ORIGINALFILENAME}"
summary(ctg)

hmean <- grid_metrics(ctg, mean(Z), 0.1)

Since I am giving chunk size of 0, catalog class considering bounding boxes as tiles and calculating metrics based on the bounding boxes. For this reason, I am getting raster image of a file which also includes part of neighboring image. Also, following warning message from R also referred to this problem.
Be careful, some tiles seem to overlap each other. lidR may return incorrect outputs with edge artifacts when processing this catalog.

If I want to make small chunks, size of 10m, lidR gives warning:
Be careful, a chunk size smaller than 250 is likely to be irrelevant.
Be careful, some tiles seem to overlap each other. lidR may return incorrect outputs with edge artifacts when processing this catalog.

How to calculate grid metrics for each file separately in my case?


texts: Saving Command History/Outputs to a file - GRASS

texts:
Is there any way to save the commands typed and the outputs obtained in terminal under GRASS GIS to a text file? I have read that GRASS saves history for each mapset but I have no much idea about how to obtain that.

Although there is "Save Output" in "Command Console" in GRASS GIS Layer Manager, many times I have faced a situation when the saved output has only one line!


texts: Is there a way to use image sprite for different marker icons?

texts: I just made the switch from Google Maps V3 to OpenLayers. I realise that OpenLayers API uses a lot of inline styling. If an image file consisting of several images can be used, it would cut down server request. Has anyone implemented this on marker icons of varying sizes?


texts: Error transforming geometry into lat-long coordinates

texts:
I'm trying to transform some geometries from the postgis database into lat long coordinates:
This is the query:
SELECT st_x (geometry) as lng, st_y (geometry) as lat FROM final

and this is the error I receive:
SQL error:

ERROR: Argument to X () must be a point


texts: What should I set vertical exaggeration to in gdaldem?

texts:
In the comments to another question I suggested that I would use a value of between 5 and 7 for the vertical exaggeration in gdaldem. The questioner pointed to the ESRI docs where they suggest a value related the pixel size of your data (and recommend a value of 0.00000956 if your data is geographic and the heights are in metres).

I had always assumed that value went into the scale parameter and that vertical exaggeration was used to overcome the fact that even the Himalayas are "small" in height when compared to the extent of the map and have traditionally used a value of 7 (or so) when creating hill shades.

So who is right?


texts: Incomplete data - lines missing from osm2pgrouting

texts:
I am importing OSM data into PostGIS using osm2pgrouting but some data is missing/did not import. I am using the default mapconfig.xml and this appears to work for most roads.

For example in below screencap, blue is the ways created by osm2pgrouting, and red is the line when the data is imported using osm2pgsql (using default.style and red filtered on "highway" = 'motorway' for this questions purposes). The missing section is a bridge, however other bridges have imported correctly, and some other roads (like 'service') are also missing.

The osm id of one of the missing lines in the pic is 158307203 if that helps, data sourced from Mapzen. The XML of that line is:

<way id="158307203" version="11" timestamp="2015-03-28T10:51:02Z" changeset="29802641" uid="1242214" user="jc86035">
        <nd ref="291167188"/>
        <nd ref="2355888876"/>
        <nd ref="1834210085"/>
        <nd ref="291167189"/>
        <nd ref="1274007359"/>
        <nd ref="1274007366"/>
        <nd ref="2355888862"/>
        <nd ref="291167190"/>
        <tag k="lit" v="yes"/>
        <tag k="ref" v="8"/>
        <tag k="foot" v="no"/>
        <tag k="name" v="xxxx Tsing Ma Bridge"/>
        <tag k="lanes" v="3"/>
        <tag k="layer" v="4"/>
        <tag k="bridge" v="yes"/>
        <tag k="oneway" v="yes"/>
        <tag k="bicycle" v="no"/>
        <tag k="highway" v="motorway"/>
        <tag k="name:de" v="Tsing-Ma-Brücke"/>
        <tag k="name:en" v="Tsing Ma Bridge"/>
        <tag k="name:es" v="Puente Tsing Ma"/>
        <tag k="name:fr" v="Pont Tsing Ma"/>
        <tag k="name:sv" v="Tsing Ma-bron"/>
        <tag k="name:zh" v="xxxx"/>
        <tag k="surface" v="asphalt"/>
        <tag k="cycleway" v="no"/>
        <tag k="maxspeed" v="80"/>
        <tag k="wikipedia" v="zh:xxxx"/>
        <tag k="start_date" v="1997-05-22"/>
        <tag k="name:zh-yue" v="xxxx"/>
    </way>

(Chinese characters replaced with x to post on StackExchange).

Below is the log when I imported the data:

$osm2pgrouting -file ~/osm2pgrouting/mapconfig.xml -dbname gis -user gisuser -passwd gisuser -clean
host=127.0.0.1 user=gisuser dbname=gis port=5432 password=gisuser
connection success
Trying to load config file /path/to/mapconfig.xml
Trying to parse config
Trying to load data
Trying to parse data
unknown maxspeed value: 56; 72; 72
unknown maxspeed value: 56;108
unknown maxspeed value: 56;108
Split ways
Dropping tables...
Creating tables...
Nodes table created
2create ways failed:
Types table created
Way_tag table created
Relations table created
Relation_ways table created
Classes table created
Adding tag types and classes to database...
Adding relations to database...
Adding nodes to database...
Adding ways to database...
Creating topology...
NOTICE:  PROCESSING:
NOTICE:  pgr_createTopology('ways',1e-05,'the_geom','gid','source','target','tru
NOTICE:  Performing checks, pelase wait .....
NOTICE:  Creating Topology, Please wait...
NOTICE:  1000 edges processed
...
NOTICE:  282000 edges processed
NOTICE:  -------------> TOPOLOGY CREATED FOR  282284 edges
NOTICE:  Rows with NULL geometry or NULL id: 0
NOTICE:  Vertices table for table osm2pgr.ways is: osm2pgr.ways_vertices_pgr
NOTICE:  ----------------------------------------------
Create Topology success
#########################
size of streets: 128763
size of splitted ways : 282284
finished

I noticed the "2create ways failed" message but not sure why no error message was generated. Hoping someone can explain why these roads failed to import.


texts: Changing name output with ModelBuilder Iterator?

texts:
I have several DEMs where I want to calculate their slope degrees.
I used iterate rasters with slope tool and I put %Name%_slope as Output Name.
The problem is the Name: it creates a file name too long ("dtm000628_slope")...
I would have slope rasters with name 335_%n%, but "%n%" doesn't work.
Do you know how I could change the name in an rasters iterators?

**The solution is found:

In the output raster, %n% makes a sequential number in the file name of output rasters, but not if you start the file name with a number:

-Example 1:

  in output raster = "C:\Output.gdb\335_slope%n%"

it will not work

-Example 2:

  in output raster = "C:\Output.gdb\slope335_%n%"

it will work

-Example 3:

  in output raster = "C:\Output.gdb\s335_%n%"

it will work


texts: Automatically create point based on location of polygon using ArcGIS Desktop?

texts:
I have a task that makes me so confused to figure out.
I have a layer Houses (Polygon), and i have to create points at the locations of the house. Like the picture i attached here, left or right of the house. Because there are thousands of house so i wonder if there is any ways to automatically create those points. Of course i have a field of location of the house. So i want to create points based on those locations.
I think the first thing to do is to identify the location of the house and set them as x,y coordinate then we create points based on x,y coordinate. But i still cannot figure out how to do it.

Actually, i couldn't determine those point which one is on the north or south. I have attached one more picture. The task is create the points which are opposite to the line and at the left corner or right corner of each polygon. In fact, i think when we cut the left or right corner of polygon we can solve the problem but i don't know how to do it. I also use Create Fishnet to make small polygon inside The House Polygon but still don't figure out how to choose the exact left corner or right corner of each house. The idea of using ArcObject, i researched it too and it seems like the IHit Test_interface could do something but i don't know how to use it. I some kind of newbie of programming with arcobject (T T).


texts: QGIS offset dash line style - dash spacing does not vary with rotation

texts:
I'd like to make a USGS-style dirt road marker but with smaller dashes and bigger spacings (i.e. less duty cycle) to indicate lower quality road.

First attempt = two simple lines, custom dashes of the same value, offset same distance in opposite directions.

Problem is, when the line goes around a curve, the spacings don't automatically adjust (stretch on outside of curves, compress on inside of curves) to account for the curve, so you get something that looks like phase shifting or 'pulling' of one side relative to the other

So, some attempts to fix that:

solid black line, centered, say width = 1; solid white line on top, centered, say width=0.6; so, this works to keep the dashes across from each other as the line curves, but, the center line that covers/hides the black is solid white - not transparent as in a USGS.  Anyone know of a way to do that?
change it to a marker layer - I used a font marker, picked a colon, set the y offset so it's basically centered, and good to go.  The 'rotate' checkbox, on by default, seems to keep the dashes 'across' from each other as desired.  Problem is, it's pixellated and a bit chunky looking.

So, is the cleanest solution to make a svg that looks nice and use that?  Not sure how it will handle compressing on inside of curves vs stretching on outside of curves.  What would really win is a checkbox in the dashed line type that says 'keep offset dashes across from each other around curves' or such, akin to the 'rotate marker' checkbox for marker lines.


texts: Defining invalid regions of VRT file with rasterio in a way that allows windowing

texts:
I'm using a VRT file to represent a bunch of individual map tiles as one big image. I noticed that some of the data is invalid but not marked as such (via the invalid data value). I'd like to define these regions as invalid using a file (e.g. a shapefile) in a way that would allow me to call dataset.read(..., masked=True) and get a masked NumPy array that masks out both data that has the dataset's invalid data value, and data that I manually define to be invalid in my file.

I found an example in the rasterio documentation using Fiona to load a shapefile and use it to mask some data. I figured that would be useful, so I built a shapefile with my invalid regions as polygons. I'm now having trouble reading from my VRT file with these shapes masked out. I've got the following code:

import rasterio
import rasterio.mask
import fiona

VRT_FILE = "/path/to/my/vrt/file.vrt"
SHAPE_FILE = "/path/to/my/shape/file.shp"

with rasterio.open(VRT_FILE) as dataset, fiona.open(SHAPE_FILE, "r") as shapefile:
    shapes = [feature["geometry"] for feature in shapefile]
    masked_data = rasterio.mask.mask(dataset, shapes, invert=True)
    # Error thrown on line above.

but this throws the error MemoryError: Unable to allocate array with shape (216640, 282250) and data type bool because the rasterio.mask.mask call is trying to load the whole (huge) VRT file, not a small section of it as I would like.

It seems there is apparently no way to define a Window object with the mask method. Until now (before starting to use shapefiles to define invalid data) I've been reading the dataset with dataset.read(... window=Window(...)), which lets me define a window to avoid reading the whole (huge) dataset at once. I'd like to do this for my shape-masked data too.

I have two questions:

Is my approach to defining these invalid regions - i.e. with a shapefile - reasonable? Or is there a canonical way? Ideally I'd like to avoid having to manually write the original data back to disk with these regions set to the dataset's real invalid data value, since it's 50 GB compressed. Perhaps there's a way to define this invalid data while generating the VRT file in the first place?
If my approach is indeed reasonable, then is there a way to create a DatasetReader-like object that supports .read(...) with a window parameter, using data masked from my shapefile?


texts: Geoserver integrated GWC not 'seeing' tiled layers

texts:
I am having trouble using Geoserver 2.4.1 integrated GWC, Ubuntu 13.04 and tomcat7

My data directory is set by setenv.sh with this string:

export CATALINA_OPTS="$CATALINA_OPTS -DGEOSERVER_DATA_DIR=/opt/gisdata/data"

When I check the /gwc folder in that directory, there are seeded folders and images. IE:

However, geoserver's tiled layers screen does not show any disk usage.

How can it be writing to the correct directory but not see those layer that it just wrote? Chmod 777 does not seem to help. Any advice, need more information?


texts: Removing return/newline (\n) character from Field using Python and Field Calculator?

texts:
I have a geodatabase table with a field that I'm trying to remove return characters (new line).  I found this post (How can I remove (chomp) a newline in Python?) of how to do it, however it is not working within field calculator.  Here are the code snippets that I tried:
Note return character is not at the end of the string.

!myField!.rstrip()

OR

!myField!.rstrip('\n')

OR

!myField!.rstrip('\r\n')

OR

!myField!.replace('\n', '')

000539 error given for this option:

  Description The calculation used by the Calculate Field or Calculate
  Value tool is invalid. This error message provided will list the
  specific Python error.

  Solution This error code covers a number of Python errors:

  Example error 1: exceptions.TypeError: cannot concatenate 'str' and
  'int' objects.  The above is a Python-specific error. The calculation
  is attempting to add or concatenate a string and a number.

  Example error 2: Invalid field shape@distance  The above is an error
  using the geometry object. The distance method is not a valid method
  of the geometry object.

  For specific Python issues, consult the external Python help for more
  information, or consult the Calculate Field or Calculate Value help
  for more information on these tools.

OR

import os
def removeReturn(myField):
  s = myField.rstrip(os.linesep)
  return s

Any ideas of how to remove return characters using field calculator?


texts: How is Landsat NDWI calculated in Google Earth Engine?

texts:
I am using Google Earth Engine to compute NDWI time series, in particular by plotting "Landsat 7 Collection 1 Tier 1 8-Day NDWI Composite" data values

The product description says:
These Landsat 7 Collection 1 Tier 1 composites are made from Tier 1 orthorectified scenes, using the computed top-of-atmosphere (TOA) reflectance.
It is derived from the Near-IR band and a second IR band, ≈1.24μm when available and the nearest available IR band otherwise.

This seems strange to me, since Landsat 7 has no 1.24μm band, but only a SWIR band at 1.55-1.75 μm
https://landsat.usgs.gov/what-are-band-designations-landsat-satellites

How does the Google algorithm calculates it?


texts: How to add point attributes to containing polygon?

texts:
I have three layers in a QGIS project:

layer 1 is a boundary layer of polygons each with and id I will call "p_id"

layer 2 is a point layer each with an id called "p_id"

layer 3 is a .csv layer with count data that I am joining to both layer 1 and layer 2 on p_id (both the points and the polygons are the same type of thing)

What I would like to do is for every polygon add the count data from each point contained within it to the polygons count data

Is there a tool in QGIS that will allow me to do that?


texts:
I am a beginner in eCognition and am working on trying to extract volume for certain classified objects. When I extract volume with the ruleset given in the program it is the same value as the area (ie. 200m2 for area and 200m3 for volume). From what I understand in the handbook there needs to be voxels available to create a 3D volume. I am currently using a handful of DEM derivatives (slope, roughness, aspect, a change map, curvature) and the DEM itself.

Can anyone suggest a way to determine volume in eCognition (or any alternative methods)?


texts: Getting TypeError:  from arcpy.mapping.ListLayers()?

texts:
I want to use ArcPy to turn on the labels for a layer.

I use this code but it fails and returns an error.

import arcpy
mxd = arcpy.mapping.MapDocument(r"C:\Users\P51\Desktop\testingLABEL.mxd")
df = arcpy.mapping.ListDataFrames(mxd,"Layers")[0]
layer = arcpy.mapping.ListLayers(mxd, "point10", "Layers")[0]
layer.showLabels = True

This is the error :

Runtime error  Traceback (most recent call last):
File "<string>", line 4, in <module>
File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\utils.py", line 182, in fn_
    return fn(*args, **kw)
File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\mapping.py", line 1517, in ListLayers
    result = mixins.MapDocumentMixin(map_document_or_layer).listLayers(wildcard, data_frame)   File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\arcobjects\mixins.py", line 843, in listLayers
    raise TypeError(str(type(dataframe))) TypeError: <type 'str'>

Does anyone know how to fix this?

I have many layers (100 layers/feature class/ file geodatabase) and I want to label those layers then I convert those layers into annotation and after that I convert them to DWG CAD (point and annotation), but I want to do that automatically (especially when do labeling and turn on those label), because it will be wasting time if labeling those 100 layers 1 by 1.


texts: "QGIS quit unexpectedly" at startup on Mac

texts:
I am having for a few months an issue with QGIS on my mac. Every time I start the program, QGIS crash with “QGIS quit unexpectedly”. The program starts only once I select “Reopen” in the message box. This is not dramatic but a bit annoying…

I have been using QGIS for a while on this laptop, and I suspect that the problem comes from the fact that I still have some traces of previous versions that are conflicting.

I already tried several things:

Reinstall QGIS
Remove all the plugins
Remove the content of .qgis folder

These steps did not solve the problem and I am now out of ideas.

My config: MAC OS 10.11.5 QGIS 2.18


texts: Can the ArcGIS Viewshed Analysis Sample take height as a parameter?

texts:
The ArcGIS JavaScript API has a viewshed analysis sample that calculates a viewshed after clicking the scene view. The geoprocessor takes an Input_Observation_Point and a Viewshed_Distance, executes its task, and draws the viewshed:
var params = {
  Input_Observation_Point: featureSet,
  Viewshed_Distance: vsDistance
};
gp.execute(params).then(drawResultData);

Here is the geoprocessor: https://sampleserver6.arcgisonline.com/arcgis/rest/services/Elevation/ESRI_Elevation_World/GPServer/Viewshed
Does the Input_Observation_Point parameter have any variables for height? If so, how can I set it in the sample program?


texts: How to export Landsat Simple Composite output?

texts:
I want to extend the "Landsat Simple Composite" example by exporting the stuff inside a rectangle.

var geometry =
/* color: #d63000 */
/* displayProperties: [
  {
    "type": "rectangle"
  }
] */
ee.Geometry.Polygon(
    [[[-47.8334884033203, -0.5160872232277051],
      [-47.8334884033203, -0.6455132239753004],
      [-47.62577813720702, -0.6455132239753004],
      [-47.62577813720702, -0.5160872232277051]]], null, false);

// Composite 6 months of Landsat 8.

// Note that the input to simpleComposite is raw data.
var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1');

// The asFloat parameter gives floating-point TOA output instead of
// the UINT8 outputs of the default simpleComposite().
var composite = ee.Algorithms.Landsat.simpleComposite({
  collection: l8.filterDate('2015-1-1', '2015-7-1'),
  asFloat: true
});

// Pick a spot with lots of clouds.
Map.setCenter(-47.6735, -0.6344, 12);
// Display a composite with a band combination chosen from:
// https://landsat.usgs.gov/how-do-landsat-8-band-combinations-differ-landsat-7-or-landsat-5-satellite-data

Map.addLayer(composite, {bands: ['B6', 'B5', 'B4'], max: [0.3, 0.4, 0.3]});

// Export the image, specifying scale and region.
Export.image.toDrive({
  image: composite,
  description: 'imageToDriveExample',
  scale: 30,
  region: geometry
});

The lines for exporting the image are from https://developers.google.com/earth-engine/exporting.  However, the .tif file, when opened, displays:

In the browser is seen:

where what is inside the red rectangle is the stuff I want to export.  What am I doing wrong?


texts: How to update GeoWebCache as part of GeoServer for maptiles

texts: I read the documentation for GeoServer but I can't seem to find a way by means of which one can update GeoWebCache dynamically by removing/deleting old (dirty map tiles) and then replace those dirty tiles with new maptiles and also how to dynamically create new maptiles in place of maptiles which are already removed being dirty and update the GeoWebCache dynamically?


texts: Layer from server as input to ArcGIS Pro Intersect

texts:
I am working in ArcGIS Pro 2.6.1 and I have managed to pull in data from a server found here: https://spartagis.ncem.org/arcgis/rest/services. The layer in particular is the NC_AllBuildings layer. I am trying to use this layer as an input to an intersect tool but when I try to use the layer, I receive an error ERROR 000840 The value is not a Feature Layer. I have tried to directly download the data from the website using a spatial subset but have not had any luck. I have even let it run overnight but when I come back in the morning it has not downloaded.
Is there a way I can use the data from the server as an input to the intersect tool that I'm overlooking?


texts: PostGIS - Calculate progressive distance between points on a GPX path

texts:
I am attempting to build an elevation profile using PostGIS-PostgreSQL only.
I have imported GPX data in my PostGis database, so now I have a table full of the GPS coordinates of the track and the timestamp of every single GPS point of the track. Now I would like to build a second table that contains, in a single row:

the track id
an array of floats representing the incremental distance from the starting point
an array of floats representing the incremental time from the starting point

Is it possible to do it in a single SQL stored procedure?

Thanks,

Andrea


texts: Is there a command to duplicate a map in Carto?

texts:
I'm hoping to make a quick duplicate of a map in Carto, the same way you can 'make a copy' of a spreadsheet or a word doc.

My plan is to add a new map layer (with new data) to the duplicate and delete the old data to create a new map with all the same attributes as before, but new geography.

I can't modify the original because it's embedded in a live webpage.


texts: Using pyproj inverse transform to get azimuth, row iteration and lambda apply (I'm close to the answer!)

texts:
I am using pyproj inverse transform to add azimuth and distance "info" to an ordered geodataframe.
I am trying to iterate over all the rows in the geodataframe to add the "info" as a tuple to a new column in the geodataframe, but there is something amiss with my lambda function.
Running the following:
Windows 10
conda 4.8.2
Python 3.8.3
shapely 1.7.0 py38hbf43935_3 conda-forge
pyproj 2.6.1.post1 py38h1dd9442_0 conda-forge
Example geodataframe:
%matplotlib inline
import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
from shapely.geometry import LineString
import pyproj

myid = [1, 1, 1, 1, 1]
myorder = [1, 2, 3, 4, 5]
lat = [36.42, 36.4, 36.32, 36.28,36.08]
long = [-118.11, -118.12, -118.07, -117.95, -117.95]
df = pd.DataFrame(list(zip(myid, myorder, lat, long)), columns =['myid', 'myorder', 'lat', 'long'])
gdf_pt = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df['long'], df['lat']))
gdf_pt = gdf_pt.set_crs(epsg=4326)
print(gdf_pt.crs)
display(gdf_pt)
ax = gdf_pt.plot();
ax.set_aspect('equal')
ax.set_xticklabels(ax.get_xticklabels(), rotation=90);

Expected output:
I am iterating over i row number, starting at the second row, and then using that row & the previous to calc the azimuths & distances. The results are as expected.
g = pyproj.Geod(ellps='WGS84')
for i, r in gdf_pt.iloc[1:].iterrows():
    myinfo = g.inv(gdf_pt.long[i], gdf_pt.lat[i], gdf_pt.long[i-1], gdf_pt.lat[i-1])
    gdf_pt.loc[i, 'az_fwd'] = myinfo[0]
    gdf_pt.loc[i, 'az_back'] = myinfo[1]
    gdf_pt.loc[i, 'dist'] = myinfo[2]

display(gdf_pt)

Wrong output, need help here!
I am trying to iterate over i row number, starting at the second row, and then using that row & the previous to calc the azimuths & distances. However, as seen in the table display, it is not passing through i as I expected, when compared to the expected results above. It is just reading the last effort and/or over-writing
g = pyproj.Geod(ellps='WGS84')
for i, r in gdf_pt.iloc[1:].iterrows():
    gdf_pt['mytuple'] = gdf_pt.apply(lambda x:
                                     g.inv(gdf_pt.long[i], gdf_pt.lat[i],
                                           gdf_pt.long[i-1], gdf_pt.lat[i-1]), axis=1)
display(gdf_pt)

I also tried the following, which threw the error shown:
g = pyproj.Geod(ellps='WGS84')
for i, r in gdf_pt.iloc[1:].iterrows():
    gdf_pt.loc[i, 'mytuple'] = gdf_pt.apply(lambda x:
                                     g.inv(gdf_pt.long[i], gdf_pt.lat[i],
                                           gdf_pt.long[i-1], gdf_pt.lat[i-1]), axis=1)
display(gdf_pt)

ValueError: Incompatible indexer with Series


texts: Segmentation Fault (segfault) when using OGR CreateField() in Python

texts:
Receiving a segfault when running this very short script in Ubuntu.

from osgeo import ogr, osr

shpfile = 'Census_County_TIGER00_IN.shp'

def cust_field(field):
    '''cust_field(shpfile, field) creates a field definition, which, by calling cust_field(), can be used to create a field using the CreateField() function.
    cust_field() DOES NOT create a field -- it simply creates a "model" for a field, that can then be called later. It's weird, but that's GDAL/OGR, as far as I can tell.'''
    fieldDefn = ogr.FieldDefn(field, ogr.OFTInteger)
    fieldDefn.SetWidth(14)
    fieldDefn.SetPrecision(6)
    return fieldDefn

ds = ogr.Open(shpfile, 1)
lyr = ds.GetLayerByIndex(0)
field = cust_field("Test")
lyr.CreateField(field)

Everything runs smoothly until that last line, when iPython, normal shell Python and the IDLE command line all dump to a segmentation fault. Is this an error on my end or an issue with the underlying C that I'm not addressing properly?