texts: Has the function srid(geometry) been removed or has the name changed?

texts:
I am in the process of restoring a database that relied on PostGIS version 1.5 functions to a newer Postgresql engine that leverages PostGIS 2.2 functions.  In my restore process, I receive the following error:

pg_restore: [archiver (db)] could not execute query: ERROR:  function srid(geometry) does not exist.

I know this error is caused by a constraint on one of my tables that utilizes the srid(geometry) function.  I want to alter the constraint to use the new function name -- if it exists.  Has this function been removed or has it been renamed since PostGIS version 1.5?

I am unable to use the legacy PostGIS functions because I am restoring to an AWS RDS instance.  RDS instances do not allow the creation of functions that use C language libraries from the underlying OS for security reasons.  That is why I am in this situation.


texts: ArcScene does not update z value after edit

texts:
I measured a feature with x, y, z coordinates.
It's the green point cloud below

When I viewed the points in ArcScene, there was a known measurement error that I had to fix by adding 1.6 m to the Z field. (I selected all the too low points, then used field calculator to add 1.6 meters).

Elevation = [Elevation] + 1.6

The result was that the feature was updated, but the ArcScene display didn't.

You can see two selected points in the table and the workspace, with the visual top and bottom corresponding to top and bottom in table, but check the Z value after my editing, you can see the discrepancy.

I have tried reloading the file, exporting the data to a new feature, but the problem remains the same.

Any thought on why this happens?


texts: Artifacts when rendering TopoJSON geometries on Cesium

texts:
I am rendering a TopoJSON document, created via mapshaper.org from a shapefile, on Cesium. That works fine for all geometries except one.
I have tried re-creating the TopoJSON document and re-editing arcs manually but I always end up with different geometry artifacts.
To me, it looks like some of the arcs have an incorrect coordinate in between, which leads to the self-intersections. But I am not able to debug the coordinate in question.
Also, when I remove arcs from the polygon manually, I end up with similar artifacts at different locations, which seem to be completly random.
Using "Check validity" in QGIS does not give me errors or self-intersections for the initial geometry, neither in the shapefile nor in the TopoJSON.
EDIT:
I have displayed the outline of the polygon and it seems that the artifacts do not have an outline, which leeds me to the conclusion that this might be a Cesium rendering issue rather than an issue with the coordinates in the TopoJSON document.
Does anyone have an idea, how to debug or correct the polygon geometry?
Cesium rendered TopoJSON:

TopoJSON with some removed arcs:

TopoJSON with outline in black:

QGIS rendered TopoJSON:


texts: Associating Photo to Feature ID in ArcPad

texts:
I'm trying to associate a photo to a feature in ArcPad.

Currently I have three photo fields for my feature layer, titled "Poles".
Essentially, I have field technicians taking three different photos for each utility pole they come across, "COMM_PHOTO, POLE_PHOTO, and TAG_PHOTO" are the three different fields for the photos in the Poles feature layer. The problem I am facing is that when they take the photos with the tablet camera, it is saving all of them as simply the field name plus the photo number in succession, such as COMM_PHOTO_001, or COMM_PHOTO_002, which describes the field and the photo that was taken in succession with the camera.

What I need is the field photo to be COMM_PHOTO, POLE_PHOTO, AND TAG_PHOTO plus the unique pole_id which is populated in the Poles layer.

Anyone have a way of doing this?


texts: Are there project-specific svgPaths for QGIS?

texts:
I have a QGIS project that needs custom SVG symbols. These are all in a directory next to the myproject.qgs project file (and FWIW, everything is under version control):

myproject/
  .git/
  myproject.qgs
  symbology/
    icon1.svg
    icon2.svg

QGIS only finds the symbols icon1.svg, icon2.svg when I add the /absolute/path/to/myproject/symbology directory to the global SVG paths settings in QGIS.

However, this means that everyone else who checks out this repository will have to adjust their global QGIS configuration before they get correct map renderings.

Is there a more elegant way to do that? For example, is it possible to define custom SVG paths entries directly in the project file?


texts:
I have downloaded a file from OSM which contains geodata for around 50k addresses, much of which are incomplete.

+--------------------+-------+---------------+
| street             | house | city          |
+--------------------+-------+---------------+
| Forest Lane        | 75230 | Dallas        |
+--------------------+-------+---------------+
| Rendon Road        |       | Burleson      |
+--------------------+-------+---------------+
| Matlock Road       | 76018 |               |
+--------------------+-------+---------------+
| Hospital Boulevard | 75051 | Grand Prairie |
+--------------------+-------+---------------+

That table have two more columns for the coordinates, which I have not copied here.

Starting from there, how is the process to match an address? Has it to be an exact match? what it happens if additional info such as house number is provided? I was thinking on making a function for this, but it looks as a complex task.


texts:
I am new to GIS and I am trying to place the North and South Salt Lake City aeronautical charts in QGIS. The problem I'm having is that there is a black section where the two charts should touch. Am I missing something when I add the raster layers in QGIS? How do I get the maps to join without the black band in the middle? Sectional charts are supposed to have a small overlap between the North and South sides of the map but that doesn't seem to be the case. Can anyone help me fix this?

Thanks in advance


texts: QGIS graphical modeler: how to automatically enter parameters in r.drain.coordinate tool

texts:
i am carrying out a cost least path analysis and i am using QGIS. The processing tools i am using are: r.cost.full and then r.drain.coordinate
Considering i have more than 1000 layer to process (starting and ending points) i have built a model. Everything works except r.drain.coordinate which requires to input E,N coordinates of the starting point (which is the ending point for what i am doing). I have tried to add the coordinate as attribute values from a table (see images below) but although the algorithm runs i do not get what i want.
Anybody has got a clue? A script that can do the job is welcome too.


texts: Error exporting data from Oracle to shp using ogr2ogr

texts:
I have an oracle table which I would like to query/make a selection. I would then like to export the selection to a shp file.

I'm using gdal and ogr2ogr.

Following the guidelines here and here I have come up with the following:

ogr2ogr -f "ESRI Shapefile" c:/Testing/Scratch/GDAL/symbol.shp OCI:dbusername/db_password@tns_name -sql "SELECT* FROM tablename where VERSION_DATE >= CURRENT_DATE -8" -progress

I then get the following error: "Unable to open datasource 'OCI:dbusername/db_password@tns_name' with the following drivers"

There is then a long list of drivers including "ESRI Shapefile".

I have now added " after OCI as suggested.

ogr2ogr -f "ESRI Shapefile" c:/Testing/Scratch/GDAL/symbol.shp OCI:"dbusername/db_password@tns_name -sql "SELECT* FROM tablename where VERSION_DATE >= CURRENT_DATE -8" -progress

I no longer get the errors about it being unable to open the datasource with the drivers.  I now just get an "Access in denied".

I have a working connection to the database in ArcGIS using the same username & password.
Does gdal generate its own user?


texts:
I have two shapefiles with polygons (say S1 and S2). Some of the polygons are overlapping in the two shapefiles but not perfectly. In addition some polygons from S1 are touching polygons from S2. I need to join the two shapefile but if there is a polygon from S1 overlapping polygon from S2 keep the polygon as is from S1 and delete that same polygon from S2. if only touching delete the overlapping area.

Is there a way to do this in R?


texts: How can I combine together cells within the same raster in ArcGIS 10.2?

texts: I have a raster displaying Land Cover data for the UK at a cell size of 25m. The cells are currently categorised as one of 24 different types of land cover. Some of these are irrelevant to me, and others could be combined to have the same effect but make the map simpler. How can I join together 2 of these land cover categories within a single raster so they're effectively listed as a single land cover type?


texts: Exporting basemap in ArcMap, similiar to export dataset?

texts:
I am trying to export the basemap in ArcMap to use it with my mosaic dataset. I need these 2 to import and work with in Unity 3D.

When I am exporting my mosaic dataset consisting of DEM files, I can right click on the dataset and click export. There I have the option to "use renderer" and "square", together with choosing to use either entire mosaic or current view dataset.

However, from what I can notice, the only way to export the basemap is to click file - export.

What is the problem:

I need the exported DEM to match with the basemap when I export them as .png.
If the exported basemap is different (different size/pixels) the basemap which is used as a projection over the imported DEM will not match (for example the sea will be projected into the mountainous region of my DEM)

Could anyone help me match these 2?

I do use photoshop to edit my DEM/basemap for importing into Unity (16bit, grayscale for DEM, etc).

Maybe is it possible to photoshop the exported basemap to match with my DEM?


texts: Finding feature extents (not layer extent) in QGIS?

texts:
How can I find feature extents for each multipolygon in QGIS?

for feat in layer.getFeatures():

I want to calculate extents for each feature here.


texts: Why do I get missing data rows when I run this virtual layer area query?

texts:
I have a shapefile of circles called circlesdixie. I have a polygon layer called FCC (the shp download link at https://www.fcc.gov/reports-research/maps/auction-904-updated-jun20-eligible-areas/. I'm creating a new virtual layer which calculated what percent of each polygon each circle covers: If a polygon is covered equally by two circles this new layer will have a row saying: Circle 1 Polygon 1 50 Pct and a row saying Circle 2 Polygon 1 50 Pct. That's what this does:
with circle_land_intersection as (
select c.name circlename, l.name landname, st_area(l.geometry) total_land_area,
    st_intersection(c.geometry,l.geometry) geometry
from circlesdixie c, FCC l
where st_intersects(c.geometry, l.geometry)
),
area_calc as (
select *, st_area(i.geometry) intersect_area
from circle_land_intersection i
)

select *, round(intersect_area/total_land_area*100) pct from area_calc;

This works on other polygon layers but when I run it on the FCC layer, I get a shapefile with a DBF of only one feature, even though 98 features correctly display on a map. The table says

Features Total 98 Filtered 1.

It's strange.
Is there any way to fix this? Is it bad geometry or corrupt data?


texts: Add X, Y column to a ESRI shapefile using ogr2ogr

texts:
I've a point ESRI shapefile in ESPG32632 spatial reference.

I need to add X and Y column to it and I'd like to use something from command line like ogr2ogr but I don't know the syntax.

Any example?


texts: Is it possible to rollback an ArcGIS Pro Update instead of uninstalling/reinstalling?

texts:
I had a moment of weakness today agreed to update my ArcGIS Professional with the latest version: 1.0.0.1808.  But unfortunately, it is crashing so bad that its unusable and I need to rollback.  Is this possible or is my only option to uninstall completely and reinstall from an earlier version?

If I have to reinstall, can I install and then update to a specific version but not the latest version?  I have the 1.0.0.1472 installer.

Is there a list of ArcGIS Pro releases versions and when they were released?


texts: What is the best way to test if a CPLString is not empty?

texts:
If I define a CPLString this way:

CPLString osString;

And I want to check if osString is not empty, what would be the best way?

My bet:

if (osString.c_str() == NULL)

But I find it a little rusty... Is there any other better way?


texts: Determining validity of odd polygon using ArcObjects?

texts:
Using ArcGIS 10.4.1 a user in our department created a polygon and saved it to a feature class as shown below. The area covered is correct as at times it could overlap itself. I have calculated the number of parts of this self-intersecting feature and it only has 1 part. I also ran the Check Geometry tool and that resulted in no issues. Through ArcObjects I used the ITopologicalOperator4 on the polygon and that results in:

IsSimple=True
IsSimpleEx=True with reason=esriNonSimpleOK
IsKnownSimple=True

A couple of my questions:

Is this a valid polygon?
If so how can I programatically determine that a polygon is of this type that self-intersects?
Can I run the Intersect tool on itself to determine I have a polygon of this type?


texts: QGIS x ArcGIS - Geoprocessing Vectors in Geographic Coordinate Systems

texts:
Is there a possibility of in future versions of QGIS, the software be able to geoprocess vectors (clip, buffer, intersection, etc.) in Geographic Coordinate System, and not only in UTM?

In ArcGIS, the software calculates automatically the geographic vectors and do the maths behind the scenes and geoprocess your data, even if all of them are geographic.

Until today, the QGIS software is only capable to geoprocess vectors if all of them are UTM, and in the same datum.


texts: C# - Geocode pharmacy by name and location

texts:
How can I find the coordinates of a pharmacy by it name and location?

I'm trying to search with the Google geocode API like this:

var pharmacyName = "Farmácia Ereirense";
var address = "Cartaxo, Santarém";

var url = "http://maps.googleapis.com/maps/api/geocode/xml?address=" + pharmacyName + ", " + address + "&sensor=false";

but I only got ZERO_RESULTS on the GeoResponse Status, and if I google "Farmácia Ereirense, Cartaxo, Santarém" it found the right location...

I already tried to do:

var pharmacyName = "Farmácia Ereirense";
var address = "Cartaxo, Santarém";

var url = "https://maps.googleapis.com/maps/api/place/nearbysearch/xml?name=" + pharmacyName + ", " + address + "&key=" + apiKey;

but I got the INVALID_REQUEST result.

Documentation I based on


texts:
I have two point feature classes (PointsA and PointsB) and one completely contiguous polygon feature class that cover all the points. For each point in PointsA I need to find the closest point from PointsB that is within the same polygon.

I feel like I need to use an iterator in model builder to iterate through each feature in PointsA - select the polygon the point is in, then select the points from PointsB within that selected polygon, and then some how add the value of the closest point from that selected group of points to the feature from PointsA, and then iterate to the next.

However, I am using 10.2 standard edition and do not have the near tool.

A Spatial Join using "closest" seems like it would work, but I don't believe I can apply that feature by feature in an iterator, can I?

Any suggestions are welcome.


texts: Generate equal distanced sections on polygon

texts:
I generated a polygon by buffering a reference line. What I want to do now is generating multiple equal-distanced points along the reference lines, upon which to generate sections on the polygon. A sketch map is here:

EDIT:

polygon data: polygon data link
reference line data: reference line link


texts:
I add a new column, naming it "reverse_cost", and give it the same data in the shape_leng column.

But as I run the shortest path algorithm in pgadmin, I got an error message.

here's my sql code

SELECT * FROM shortest_path('
            SELECT gid as id,
                     start_id::int4 AS source,
                     end_id::int4 AS target,
                     shape_leng::double precision as cost
                    FROM network',
            10, 17, true, true);

and here's the error

Error:  Error, reverse_cost is used, but query did't return 'reverse_cost' column

********** Error **********

Error: Error, reverse_cost is used, but query did't return 'reverse_cost' column
SQL status: XX000

As I change the code to this one, everything went well,

SELECT * FROM shortest_path('
            SELECT gid as id,
                     start_id::int4 AS source,
                     end_id::int4 AS target,
                     shape_leng::double precision as cost
                    FROM network',
            10, 17, false, false);

so I want to know how to use reverse_cost in PostGIS？


texts: Simplify PostGIS authentication process from QGIS

texts:
In QGIS I have a project with data from various (~15) PostGIS databases, with no login credentials (authentication user/pwd) stored for each of the connections as shown below.

Now when another user opens the project, he (or she) has to enter credentials for EACH database used:

I strongly fear that legitimatly none of my users is going to accept entering his (or her) credentials 15 times (but in an enterprise environment with potentially 100s of users, acceptance is key for success).

So two questions arise for me:

is there a way to suggest the Username according to the global project variable user_account_name and
is there a way to tell QGIS to use the credentials entered for the first db connection for any other connection as well?

I can't estimate if it would be a great deal to modify the 'Enter Credentials' dialogue as proposed below, but is this possible?


texts: Does leaflet.js offer an gutter option just like Openlayers

texts:
Currently we are switching from Openlayers to leaflet.js and there is a problem using a WMS layer with leaflet. The WMS being used is Deegree. When a layer is visualised within leaflet there are artifacts at the tile edges as seen in the added picture.

With openlayers we fixed it with a gutter option. Is there a similar solution for leaflet?


texts: Linking search box to dataset column in Carto?

texts: The search box in a published map of CartoDB isn't linked to the dataset which is shown in the CartoDB map. How can I link the search box to a certain column of a dataset?


texts: downscaling of resolution of surface roughness

texts:
I am using NCEP/NCAR reanalysis data for mesoscale modeling at a resolution of 5 Km.

I have calculated the surface roughness for the model and using this surface roughness, I will calculate the potential zones for wind farms.

Now, I wish to use these results in microscale model which requires to downscale the data.

What methodology should I adopt?


texts: So I have spend quite a bit of time reassigning F1 to F12 to mapinfo tools that suits our work line using Command Editor. Now I need to do same for couple of other computers. I was wondering if I can easily export this shortcut settings to another computer.


texts: How to convert lines from dxf to polygons?

texts:
I have several dxf files. Each file, contains the 0 layer. The 0 layer contains lines (which make a closed polygon). In each dxf the lines have their color and line width.
I do not have autocad or else like.
I have Qgis 1.7.4.
I can put these dxf files into vector layers (as line vector) into Qgis.

I do not know the way (using Qgis 1.7.4 ):
1. To turn each vector layer (with lines) to vector layer with polygon.
To “copy” some of these polygons (vector layers) into a new vector layer.
2. Can each polygon into the new vector layer to have different border color?
thank you


texts: How can I get QGIS to release/unlock a geopackage with 1 layer after use in python? Geopackage/Permission Error (WinError 32)

texts:
I have a python script to update 4 authoritative geopackage files with additions and changes made to copies of the geopackages in the field (field geopackages).  Once the new records and edited records within the field DB are copied to the authoritative geopackages my script closes the geopackages from QGIS, followed by cutting and pasting the processed field geopackages to an archive folder.  This process works fine for 3/4 of the geopackages, but fails on the 4th.  The geopackage I have this issue with is the the only 1 of 4 with a single layer/table.
The function I use to move the field geopackages to an archive folder is this:
def moveUpdatedInputs(incomingDataPath, processedDataPath):
inP = incomingDataPath
outP = processedDataPath
##get all geopackage files
gpkgIn = findGPKG(inP)
##get time stamp
t = time.localtime()
tryAgainList =[]
tStamp = time.strftime("_%Y_%m_%d_%H_%M_%S", t)
##create new save as names
saveOutNames = []
for item in gpkgIn:
    insertSpot = item.find(".")
    interimName = item[:insertSpot] + tStamp + item[insertSpot:]
    newName = interimName.replace(inP,outP)
    shutil.copyfile(item,newName)
    try:
        os.remove(item)
    except:
        print ("locked")
        tryAgainList.append(item)
for item in tryAgainList:
    insertSpot = item.find(".")
    interimName = item[:insertSpot] + tStamp + item[insertSpot:]
    newName = interimName.replace(inP,outP)
    shutil.copyfile(item,newName)
    try:
        os.remove(item)
    except:
        print ("locked")
        #tryAgainList.append(item)

I access this script after creating the data syncing object, then deleting the data syncing object to unlock the files from QGIS.  This works when entered line my line in the command line:
v1 = compareFMO_TableSettings(incomingDataPath,outgoingDataPath)
del(v1)
moveUpdatedInputs(incomingDataPath, processedDataPath)

But if I try to automate the commands with a function, it fails because the geopackage with one layer remains locked:
 def runner():
       v1 = compareFMO_TableSettings(incomingDataPath,outgoingDataPath)
       del(v1)
       moveUpdatedInputs(incomingDataPath, processedDataPath)

After adding a blank second table to the geopackage of concern everything works, so I have found a work around, but I cannot figure out why QGIS maintains a lock on a geopackage with a single layer when the geopackage is used within a python script but does not behave this way with a geopackage containing 2+ layers.
What I am looking for is a solution that can cut and paste geopackages containing only 1 layer after processing, whereas currently my solution only works for geopackages containing 2+ layers.
I am using QGIS 3.10.5, if that is a useful detail.


texts: Finding shortest path while avoiding polygon features

texts:
I'm developing an extension to ArcMap where I have a set of static, concave polygon features and I'm trying to find the shortest path between 2 points without crossing any of the polygons.

Is there an algorithm to do this?


texts:
Does anyone know a source of districts in Spain?

I can find Madrid and Barcelona each at a different site. If I try to collect the districts for the whole country it would be a massive undertaking.

Is there any central information site with these? or does anybody has this shapefile?


texts: How to clear cache in ArcPy so visual appearance updated after processing

texts:
I have written a script which updates lines and points, the script runs fine through ArcGIS PRO but the issue is after running the script the visual appearance of the updated lines and points are cached in memory and remain in their original positions (unless you scroll to a different zoom level, I am running the script as I am looking at it on the screen)

I have noticed that scripts within ArcGIS pro automatically update the appearance after running the script, how do I implement this into my custom script?

So far I have tried calling arcpy.ClearWorkspaceCache_management() before ending the script but it does nothing, I am not convinced this is the correct thing I should be using because it did not crash.

Thanks.


texts: Creating mutliple connector lines from point to nearby lines

texts:
I need to create multiple connector lines from a suburb centroid point to closest two or three road lines in arcmap.

How can I go about this?


texts: align tool on Arcview 3.2 Window Seven missing

texts:
after following the necessary procedure to install  arcview 3.2 on my windows seven...
I activated all extensions on including the spatial tool, make it as default, then applied. But i still got an error saying that some spatial... script is missing. hence i dont see the align tool on the toolbar.
i guess this tool is used to enter co-ordinate or georeferencing based on the online tutorial have been watching.
Pls how can i overcome this error


texts: Creating new empty vector layer that has same set of atributes as another layer and then copying subset of features into it

texts:
I am new to both Python and QGIS.

Here's a sample of my code:

unique = layer.uniqueValues(index)
attributes = layer..dataProvider().fields().toList()
assistant_layer = QgsVectorLayer("Point", "temporary_points", "memory")
assistant_pr = assistant_layer.dataProvider()
assistant_layer.startEditing()
assistant_pr.addAttributes(attributes)
for i in unique:
    features = vector.features(layer)
    for f in features:
        idVar = f[fieldName]
        if unicode(idVar).strip() == unicode(i).strip():
            assistant_pr.addFeatures([f])
    assistant_layer.commitChanges()
    delone_triangles = processing.runalg("qgis:delaunaytriangulation", assistant_layer, None)['OUTPUT']
#rest of the code is irrelevant to the question

Basically I'm trying to create a layer that has the same set of properties as my input layer. Much in the same way as here How to create a new empty vector layer programmatically? except there the Attributes have been defined in code, I want mine to load them from another input layer. Then I want to copy a specific subset of features onto that layer. Something however seems to be wrong with my method.

The broader issue that I'm trying to solve here is that I need to feed a subset of features from a layer into an algorithm that takes a whole layer as its input (specifically "qgis:delaunaytriangulation"), so my intended solution is to create another layer and then copy specific features that fit the criteria onto it. However, whereas normally qgis:delaunaytriangulation outputs a path to a file containing the triangulation, when I do the above it returns a path to a file that doesn't exist. I figure it's because either my assistant_layer or the method of copying features onto it is faulty.


texts:
I want to retile big TLS files into quite small chunks of maximum 20 m x 20 m size. I have to normalize the point cloud, thus I need to classify the points and create a DTM for each tile. For this process, especially smaller tiles seem to cause me trouble (I guess it's hard to fit  a plane in only few points). But I noticed when splitting up the data into chunks, there are quite many chunks with only little information.
ctg <- readTLSLAScatalog(path_points)
opt_chunk_buffer(ctg) <- 0
opt_chunk_size(ctg) <- 20
plot(ctg, chunk=TRUE)

Can I somehow reposition the tiles so the bottom left corner of the data is also the bottom left corner of the first chunk? Or (even better) a possibility to center the tiles and the chunks so all chunks on the borders have similar sizes? I saw there is an option opt_chunk_alignment(), but I don't understand how it works exactly. Also, I guess I could play around with the chunk size because it might help, but I want to automize the process for several LAS files.


texts: Set Z value to a linestring from a postgis raster

texts:
I have a table with multilinestring geometries that can be updated or inserted by external users, and a raster table with elevation values.

I want to write a trigger that sets the z value of each point in the line from the DEM.

Probably I can write a function that iterates over the linestring points with ST_NPoint or ST_DumpPoints. Gets the elevation from the raster with ST_Value. And recreates the point with ST_SetPoint.

But, is there any easier way to to this?


texts: How to install OSGeo4W libraries in older version of QGIS (2.16)?

texts: I am trying to install older version of QGIS (2.16) with OSGeo4W. But from the latest OSGeo4W Setup i can choose only 2.18 or 2.14 version. I want to install few libraries that are provided by OSGeo4W setup file to older version of QGIS. How can i do it?


texts: gdal missing jp2 driver files how to install them

texts:
Am operating on OSX(latest) have gdal installed and rasterio (both from conda versions 2.1.3 and 1.0xxx.  I can read tiff perfectly but jp2 files (aka from Sentinel-2 sat).

When I search for installed gdal drivers the ones for jp2 are on the list.  I have attempted to remove gdal via conda and reinstall with no luck.  Am not looking to use brew as it has a pre 2.0 version and i need 2.1 functionality.

How can i install this driver?


texts: Positional Accuracy of linear features using buffers in ArcPy

texts:
I am kind of lost on how to go about this, but I have two vector shape files of a road and want see how accurate they are by estimating the positional accuracy. The files are actualRoad as the true road and testRoad is a generated road. I want to create a buffer around the actualRoad with a width of x amount on both sides so I can calculate the proportion of the buffer intersecting to the testRoad file. I will have to loop through this process in order to increase the x width of the buffer on both sides until I get close to 99% of the testRoad file intersected with the buffer. To do this, I want to create a table in arcpy with columns of increasing buffer width, length of road within the buffer, and % of testRoad in the buffer. The table should show a new row for each iteration the actualRoad buffer gets larger until 99% or prop .99 intersects the testRoad.

My issue is trying to get the loops in working order where the script can go through both files, compare the two shapefiles and properly execute the buffer and calculate the % intersected.

import arcpy
from arcpy import env

import os

path = 'C:/filepath'
env.workspace = path

#import files into list
#actual and test, both end in .shp
linelist = [f for f in os.listdir(path) if f.endswith(".shp")]
print(linelist)

#create an empty table for info
table = 'outputTable.dbf'
#add the appopriate table attributes: buffer width, length of road,
#buffer width, % of road shapefile in buffer
arcpy.CreateTable_management(path, table)

arcpy.AddField_management(table, "Buffer_Width", "FLOAT", 10, 1)
arcpy.AddField_management(table, "Length_Road_Buffer", "FLOAT", 8, 1)
arcpy.AddField_management(table, "%_Road_in_Buffer", "FLOAT", 8, 1)

tableCursor = arcpy.da.InsertCursor(table, ["Buffer_Width",
"Length_Road_Buffer", "%_Road_in_Buffer"])

for r in linelist:
  bufferWidth = 0 #reset buffer between files
  lenRoad = 0 #reset road length between files

   with arcpy.da.SearchCursor(r, ["SHAPE@LENGTH"])as cursor:
      #proportion of intersect between bufferWidth and lenRoad need to be
      #<.99
      with a while loop
        #buffer analysis(road in question(acutal or test), bufferName, r)
         ##bufferList.append(b)

       #calculate intersection
       for r in cursor:
         #twidth
         bufferWidth += r[0]
         lenRoad +=1
   outrow = [bufferWidth, lenRoad]
   tableCursor.insertRow(outrow)

del tableCursor


texts: How to make maps in QGIS or Openstreetmap, etc. to wrap around a physical sphere?

texts: How do I go about creating a map of the earth in QGIS or Openstreetmap, etc. and then export art in the appropriate projection so that it can be wrapped around a sphere?


texts: can I add additional data to an existing pgrouting postgis database using osm2pgrouting?

texts:
This is just a simple question but I could not find the information in the documentation. If I have an existing pgrouting database table named routing, can I add additional OpenStreetMap data to this table through the osm2pgrouting software? I have some data in the database for Syria roads and I wanted to add additional road information for Turkey, Jordan, and Lebanon. I was not sure if I can just add these osm files to the existing database table, or if I need to drop and load some large osm file that has all of the countries I want.

Note that I did try using the osm2pgrouting tool on a subsequent osm file. The tool seems to have completed okay, though there were some errors on creating some of the relationships and keys, etc. (because they already existed). However, I was not sure if this approach was the right way to go. With so much data added, it would be nice to have a diagnostic to check whether the routing is correct. According to the osm2pgrouting workshop, the osm2pgrouting tool can produce certain types of errors if you are not careful.


texts: Getting the area of a polygon from Google Earth or .xyz digital elevation model INCLUDING area due to topography

texts:
I'm doing some analysis for a study on geomorphology. For this I'm comparing the 'bumpiness' of lava fields to the 'bumpiness' of the surroundings. What I'm trying to do for this is measure the ratio of the actual area inside a polygon (so including topography) to the area inside the polygon disregarding topography. So a a value >> 1 would be quite bumpy etc.

I have all the cut outs in .xyz so Lat Long, Elev and all the cut outs as .kml but I'm wondering if there is a straightforward way to do this before I launch into trying to code it myself. I have and can use MatLab!! I can take the area disregarding elevation straight from Google Earth so this bit is not a problem!


texts: GRASS "Projection of dataset does not appear to match current location" Error

texts:
I tried to import a virtual raster (that's possible right?) into GRASS v. 7.4.0 loaded as a plugin of QGIS v. 2.18.18. It gave me this error however. I know the rasters in the VR are NAD83 (I checked the metadata to be certain), and my location/ mapset is also NAD83 (I double checked this as well), so why is there a projection error? I didn't get this error when I loaded the rasters individually either

g.proj gave me:

C:/PROGRA~1/QGIS2~1.18/apps/grass/GRASS-~1.0/bin/g.proj.exe -p
-PROJ_INFO-------------------------------------------------

name       : GRS 1980(IUGG, 1980)

ellps      : grs80

proj       : ll

towgs84    : 0,0,0,0,0,0,0

no_defs    : defined

-PROJ_UNITS------------------------------------------------

unit       : degree

units      : degrees

meters     : 1.0

Successfully finished


texts:
I'm trying to draw a circle with OpenLayers with a Google Maps Layer. I'm able to transform a simple marker and have it draw at the correct location, but I can't figure out the circle. I'm trying to draw a 50 meter radius circle:

function InitializeMap() {

// create map
map = new OpenLayers.Map({
    div: "map_canvas",
    theme: null,
    units:'m',
    projection: "EPSG:900913",
    displayProjection: new OpenLayers.Projection("EPSG:4326"),
    controls: [
        new OpenLayers.Control.Attribution(),
        new OpenLayers.Control.TouchNavigation({
            dragPanOptions: {
                enableKinetic: true
            }
        })
    ]
});

// load layers
var gmap = new OpenLayers.Layer.Google(
    "Road", // the default
    {numZoomLevels: 20}
);
map.addLayer(gmap);
if (map.layers.length > 0) {
    var gphy = new OpenLayers.Layer.Google(
        "Terrain",
        {type: google.maps.MapTypeId.TERRAIN}
    );
    map.addLayer(gphy);
    var ghyb = new OpenLayers.Layer.Google(
        "Hybrid",
        {type: google.maps.MapTypeId.HYBRID, numZoomLevels: 20}
    );
    map.addLayer(ghyb);
    var gsat = new OpenLayers.Layer.Google(
        "Satellite",
        {type: google.maps.MapTypeId.SATELLITE, numZoomLevels: 22}
    );
    map.addLayer(gsat);
    map.addControl(new OpenLayers.Control.LayerSwitcher);
}
var blankmap = new OpenLayers.Layer("Blank", {
    isBaseLayer: true,
    maxExtent: new OpenLayers.Bounds(
        -128 * 156543.03390625,
        -128 * 156543.03390625,
        128 * 156543.03390625,
        128 * 156543.03390625
    ),
    sphericalMercator: true,
    maxResolution: 156543.03390625,
    units: "m",
    projection: "EPSG:900913"});
map.addLayer(blankmap);

// Set default layer if not offline
if (map.layers.length > 1) {
  switch (default_map_type) {
    case 1: map.setBaseLayer(blankmap); break;
    case 2: map.setBaseLayer(gmap); break;
    case 3: map.setBaseLayer(gsat); break;
    case 4: map.setBaseLayer(ghyb); break;
    case 5: map.setBaseLayer(gphy); break;
  }
}

// Add Layers to map
    var currentLocationLayer = new OpenLayers.Layer.Markers("Current Location");
    var currentAccuracyLayer = new OpenLayers.Layer.Vector("Current Accuracy");
    map.addLayer(currentLocationLayer);
map.addLayer(currentAccuracyLayer);

    // Add marker
    var current_point = new OpenLayers.LonLat(-81,37.transform(map.displayProjection,map.getProjectionObject());
    map.setCenter(current_point);
    map.zoomTo(17);
    current_loc_marker = new OpenLayers.Marker(current_point,image_GPS);
    currentLocationLayer.addMarker(current_loc_marker);

    // Add circle
    var hA_circle_geometry = new OpenLayers.Geometry.Polygon.createRegularPolygon(current_point, 50, 50, 0).transform(map.displayProjection,map.getProjectionObject());
    hA_circle = new OpenLayers.Feature.Vector(hA_circle_geometry);
currentAccuracyLayer.addFeatures([hA_circle]);
}

I tried the following for the circle and it still does nto show. Does it need to be transformed due to the Google Layer?:

var hA_circle = new OpenLayers.Control.DrawFeature(currentAccuracyLayer,
                                        OpenLayers.Handler.RegularPolygon,
                                        {sides: 40, radius: 50, angle: 0});
map.addControl(hA_circle);


texts: Enabling colleagues access to tabular data in QGIS?

texts:
This question may sound like it's a bit vague, but I think that lots of people seem to face similar issues, and on this occasion a general answer seems better than a specific focused one:

I use QGIS. I need to be able to collaborate with colleagues who use typical office software - working with QGIS to visualise their data. I've struggled for months to find a robust and simple way to do so. Last time I tried I ended up reverting to using Mapinfo, despite the other problems/limitations this caused me with symbology.

I'm in a typical organisation - most colleagues are limited to using typical office tools. This includes MS Office (Excel, Access). They are familiar with these tools and unlikely to use anything new. Anything I do has to work within their boundaries or it'll fail before it starts. (I've converted several to using QGIS themselves - but for the majority this isn't an option).

I use QGIS. I'm competent with working with shapefiles, TAB files, CSV files, linking data layers, relationships, etc - but pretty much a newbie when it comes to coding, python and other database solutions.

I (and others I think) need a simple and robust method to allow colleagues to work with tables of data - data which includes a field which I can use to link it to GIS records. I don't even need to be able to cope with geometry fields in their data. It's my job to manage geometry... all I'm working toward is relationships/links in their data with an existing GIS layer.

I've completely failed to find a simple way to do this... repeatedly getting bogged down in trying (incompetently) to connect MS Access to text files, or QGIS (64bit, 2.8, Windows 7) to MS Access in one way or another (no I can't change to 32bit). I'm currently working with asking the Excel user to save as a csv file so that QGIS can recognise this, but this is messy for the Excel user who has to handle odd (for them) messages about csv saves. I had QGIS talking directly to MS Excel for a while: I'm not sure how I achieved it but I think it knew the column headings should be field names and I could even edit the Excel file using QGIS at one point. But things fail and fall apart regularly. After one save to the Excel table (using Excel I think) QGIS stopped recognising the column headings (taking them to be data). I can't find how to re-set this, so am faced with columns labelled "field 1" etc.

Is there any simple solution to this?

I can see from a range of more specific questions that many other people face the same issue (e.g. for example see questions about connecting Access to QGIS). On this occasion I'm free to create a new solution... so long as it's simple for me to set up and doesn't rely on colleagues learning anything new.

This simple task seems to me to be a key requirement of a good GIS system for use in the real world, which is why it seems worth such a wordy/general question

Previously I asked a vaguely similar question about enabling the editing of text data, and although not tested much I think I found a way to work with some office software using a dbf file (see my own answer) but I can't make dbf work with MS Office software at the moment.


texts:
I have a Baroque tombstone for which I have to map the damage phenoma. Normally one can use a vector drawing programm but one has always two boundaries between the different polygons. Questions:
Is it possible to add a layer of a photograph of the tombstone and to scale it afterwards in Quantum GIS? If one adds a layer, QGIS asks for the CRS and, of course, tombstones has no CRS.
Besides the „boundary question“ it should then be possible to calculate the percentage of the different damage phenomena which would be very nice.

Thanks in advance,

Jörn


texts: QGIS Georeferencer, unexpected transformation, and GCP placement

texts:
I'm trying to get the sample map which I believe to be EPSG:3112 into QGIS and georeferenced. In my screenshot, on the right is the raster in georeferencer with the points I have added. On the left is how that same raster appears now those Ground Control Points (GCPs) have been added.

Two things are happening that I'm not sure why or what I can do about:

The raster is being stretched
And the GCP's are being arranged in a parallel type grid

Both my project, and raster layer are set to EPSG:3112 and when I do the referencing I have EPSG:3112 selected as well in the transformation settings.

Ive since tried setting these settings but no difference:

This is what i get with "projective" transformation (still deformed)

This is the result after making the grid and changing on the fly to 3112, there is also pictured the settings i try for "densify geometries", when i click ok it asks for crs, and after that it i dont see anything new and the except in the layer window, except the dialog stays open and progress bar continually slides back and forth...


texts: Error in .bboxCoords(coords) : nrow(coords) > 0 is not TRUE using raster::sampleRandom()

texts:
I am trying to sample a raster (pixel values= 1 or NA, where 1 means opportunity as defined for my research).
The specs of my raster are
class      : RasterLayer
dimensions : 15838, 3808, 60311104  (nrow, ncol, ncell)
resolution : 60, 60  (x, y)
extent     : 2402625, 2631105, 788651, 1738931  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=43 +ellps=WGS84 +units=m +no_defs
source     : <folder> not including in this post
names      : <name> not including in this post
values     : 1, 1  (min, max)

I want to sample 30 random points in all pixels with value=1 for further analyses and research.
sampled_points <-raster::sampleRandom(above raster, size = 30 ,na.rm=TRUE, xy=TRUE, sp= TRUE)

But the Error in .bboxCoords(coords) : nrow(coords) > 0 is not TRUE is present. I have used sampleRandom() multiple times before in other analyses an dhave not encountered this problem. How do I solve this problem?


texts: Deploying Custom dll Commands etc in ArcGIS 10

texts:
Having followed the tutorial for deploying custom commands etc at the following link:

How to deploy a custom component using a setup project

It's a little different to how things were at version 9.  Well I've tried many things and I guess I am missing something stupid - but I cannot seem to get the custom command to show up in the Customise dialog.
I look in the Categories.exe app - and there is nothing under MxCommands - so I definately know that it has not installed propely.  Jeez I hate these installers - they are so fiddly to setup properly!!!!  I have a suspicious feeling that the custom installer class is not being fired properly to trigger the regasm utility but I can't be sure because I haven't been able to properly debug the installation process.

Are there any glaring omissions or gotchas I should be aware of?


texts: How to use CARTO map visualization in Odyssey.js sandbox

texts: I have a story map in Odyssey's sandbox GUI. I want to use a visualization I have in CARTO as a base map for this map. How can I find my viz key? All the documentation seems to refer to an option that is no longer available. The documentation says to find it in the share options and click on CARTO.js, which is no longer an option inside CARTO builder.


texts:
I am working with QGIS and in particular with river basins.

If you see my picture below, each label shows a particular basin but actually only the biggest is shown (because I think its area covers the areas of the smaller ones).

e.g. 4103600 is the largest basin and its boundaries are shown correctly, but within it there are 7 more nested basins (i.e. 4103631, 4103650, 4103615, 4103610, 4103750, 4103751 and 4103570).

How can I show ALL the basins' boundaries (i.e. nested and non-nested) with QGIS?

I believe there should be a trick for doing so but I don't know which one!


texts: Converting geographic coordinate system in R


texts: white image while Trying to mosiac two hyperspectral data in ENVI

texts:
I am trying to mosaic two aviris-ng georefrenced hyoerspectral data but one of them become white in preview. so I am not able to mosaic the image.enter link description here

error ink


texts: How to select "incidents" in network analysis based on their routes?

texts:
I am using Closest Facility tool in Network analysis to locate closest routes between houses and bus stops.

I am using Length as accumulation attribute to know distance between each house and bus stop. Houses are incidents and bus stops facilities.

When I am done with analysis I can see each route from house to bus stop but I can't find way how to select houses according to their calculated routes. For example I would like to do selection of houses which can reach bus stops with distance no more then 300m.

Since the points of houses are not attaching to road network there is no way how I can do selection based on location. If I would like to do selection based on attributes I can't see any match between attribute tables of "routes" and "incidents".

I guess there must be some way!


texts:
I want to convert DXF from AutoCAD to KML. It works but the coordinates don't come out right when I compare them to a KML produced with global_mapper and consequently it is shifted some when uploading the KML to Google Earth. Maybe someone can illustrate me how to set/calculate the parameters right to adjust the shift?

ogr2ogr -f KML -s_srs "+proj=utm +zone=30 +ellps=WGS84 +towgs84=-36.7193107889714,-4.418216852011938,0,0,0 +units=m +datum=WGS84 +no_defs " -t_srs epsg:25830 -dialect SQLite -sql "select * from entities where layer='SUPERFICIES  ZV'" SUPERFICIES_ZV.kml D02-32.dxf

Southern Spain ETRS89 30 EPSG:25830


texts: How to avoid background zero values in erdas imagine model maker?

texts: Does anyone know how to avoid background black color/ zero value in background pixels in processed images using erdas imagine? Particularly in model maker. It's confusing in statistical calculation.


texts: functionality of Surface volume tool in ArcMap?

texts:
I'm using the surface volume tool to calculate multiple volume below a plane height.
Starting from the lowest value of my DEM, I perform several calculations in 0.1m intervals.
However I'm not satisfied with my output, as I don't understand how its possible that my volumes and 2D area outputs do not match with the simple:
V = A*h
As my h = 0.1;0.2 ... I'm confused my 2D area is not equal V/h.
What is wrong here? My linear unit of the DEM is meters (EPSG 25832).
For illustration i include the following example (xyz data) which I converted into a Raster (1m resolution):
647724  5812415 45;
647724  5812416 45;
647724  5812417 45;
647725  5812415 45;
647725  5812416 44;
647725  5812417 45;
647726  5812415 45;
647726  5812416 45;
647726  5812417 45

I ran the surface tool for four plane heights (44 m, 44.5m, 45m & 46m) and received the following volumes:
plane height, volume & 2d Area
44, 0, 0
44.5, 0.1, 0.66666666666667
45, 1, 4
46, 5, 4

for 44m respectively 45m I receive correct results, while the remaining ones for 44.5m and 46m make no sense at all.
Is there anyone who could reproduce that and is willing to explain how those results come up?


texts: Install latest GDAL on docker image of CentOS

texts:
I am trying to install gdal latest version on a docker image of CentOS. Below is the minimal code.

# Start and connect .
docker run --rm -it centos:latest bash
yum -y upgrade
yum -y install wget
yum groupinstall "Development Tools"

cd /tmp
wget http://download.osgeo.org/gdal/3.0.4/gdal-3.0.4.tar.gz
tar zxvf gdal-3.0.4.tar.gz
cd gdal-3.0.4/
./configure

At ./configure step, I am getting error. configure: error: PROJ 6 symbols not found. Ideally I would like to install gdal latest version just by using yum, i.e. without building it. I am looking for a solution that would either fix the mentioned issue of allow me to install the latest version using yum.


texts: pyscripter : add a legend to my vector layer

texts:
I work on PyScripter, i managed to display a vector layer from Postgresql with a symbology.
But I have a difficulty to display a legend.

def main():

def init():

# appel de l'application et initialisation
   a = QgsApplication(sys.argv, True)
   QgsApplication.setPrefixPath("C:/OSGeo4W/apps/qgis", True)
   QgsApplication.initQgis()
   return a

def show_shape(app):

# definition a connection
uri = QgsDataSourceURI()
uri.setConnection("localhost", "5432", "stat", "postgres", "admin")
uri.setDataSource("public", "com", "the_geom", ' ', "gid")
uri.uri()

# crÃ©ation d'une fenÃªtre et chargement d'une couche
carte = QgsMapCanvas()
layer=QgsVectorLayer(uri.uri(), 'com', 'postgres')
if not layer.isValid():
   print "Layer failed to load!"

# add layer to the registry
QgsMapLayerRegistry.instance().addMapLayer(layer)
carte.setExtent(layer.extent())

# affichage de la couche dans une fenÃªtre
carte.setLayerSet( [ QgsMapCanvasLayer(layer) ] )
carte.show()

#personnaliser la symbologie
myTargetField = 'code'
myRangeList = []
myOpacity = 1

# Make our first symbol and range...
myMin = 0.0
myMax = 100.0
myLabel = '0.0-200'
myColour = QtGui.QColor('white')#ffee00
mySymbol1 = QgsSymbolV2.defaultSymbol(layer.geometryType())
mySymbol1.setColor(myColour)
mySymbol1.setAlpha(myOpacity)
myRange1 = QgsRendererRangeV2(myMin,myMax,mySymbol1,myLabel)
myRangeList.append(myRange1)

#now make two symbol and range...
myMin = 100.1
myMax = 102
myLabel = '200.1-500'
myColour = QtGui.QColor('yellow')#00eeff
mySymbol2 = QgsSymbolV2.defaultSymbol(layer.geometryType())
mySymbol2.setColor(myColour)
mySymbol2.setAlpha(myOpacity)
myRange2 = QgsRendererRangeV2(myMin,myMax,mySymbol2,myLabel)
myRangeList.append(myRange2)

myRenderer = QgsGraduatedSymbolRendererV2('', myRangeList)
myRenderer.setMode(QgsGraduatedSymbolRendererV2.EqualInterval)
myRenderer.setClassAttribute(myTargetField)

layer.setRendererV2(myRenderer)
QgsMapLayerRegistry.instance().addMapLayer(layer)
print "fine"

app.exec_()

app = init()
 show_shape(app)

if name == "main":
    main()


texts: How to make combined symbols in the legend in ArcGIS Desktop 10.7

texts:
I am looking for a method to make more compact legend entries. Suppose 3 point layers "inner", "middle", and "outer" in ArcGIS 10.7. Assume each layer has a field called measurement with a value range from 1 to 1000 symbolized by different colors. The coloration is the same for all three layers. The size of the symbol for each layer is set to a fixed amount. Layer "inner" has a small circle etc., c.f. image below.

What I would like to do, is make a merged/joined/combined symbol in the legend so that it becomes more compact, exempli gratia like so:

If you know of a better design, I am open for suggestions.

Explanation and a little bit more background

The problem arises because I want to have legend entries to explain the different colors, but as the colors appear in all three sizes, I want to also show the sizes. The automatic legend by Arc gets very very long when showing all colors for all sizes, so I'd like to first show all the colors and what they represent and then say, look, there's three sizes and the colors apply to all three sizes respectively.


texts: QGIS crashes when toggle editing for child layer in QgsDualView

texts:
The crash occurs for:

Ubuntu 18.04 bionic
QGIS 3.4.2

For reproduce the crash, I use the data from Natural Earth Data:

https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/cultural/ne_110m_admin_0_countries.zip
https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/cultural/ne_110m_populated_places_simple.zip

Here the code:

canvas = qgis.utils.iface.mapCanvas()
for ml in QgsProject.instance().mapLayers().values():
    if ml.name() == 'ne_110m_admin_0_countries':
        layer = ml
        par = ml.id()
    elif ml.name() == 'ne_110m_populated_places_simple':
        chi = ml.id()

rel = QgsRelation()
rel.setReferencingLayer(chi)
rel.setReferencedLayer(par)
rel.addFieldPair('adm0_a3', 'ADM0_A3')
rel.setId('countries_places')
rel.setName('Relation beetween places & countries')
if rel.isValid():
    QgsProject.instance().relationManager().addRelation(rel)

dialog = QgsDualView()
dialog.init(layer, canvas)
dialog.setView(QgsDualView.AttributeEditor)
dialog.show()

In a QGIS Project with the two vector layers loaded, the code above generates a QgsDualView dialogue type and sets up a QgsRelation between vector layers.

When I click on the toggle Editing button on the QgsRelationEditorWidget, QGIS crashes, even if the relation is set manually. What could be a problem?


texts: How can get correct elevation values from Digital elevation models

texts: I am a bit new to spatial data analysis. I am trying to analyse elevation data for a coastal (almost flat) city with the highest elevation of about 30m. I have DEMs from ASTER, SRTM, and a DSM from ALOS W3D-30. I have realized that the values are too high: for ASTER, 0-250, for srtm, 0-101, for ALOS, 0-92. Are the values not in metres (m) or is there something I need to do to have the correct (approximately) elevation values?


texts:
With MapServer, is it possible to use an expression to round a number for a label?

For example, my dataset has the float value 174.000000000 which I want to show in a label as 174.

I have tried a few mathematical expressions for a label, but, now think that MapServer may not accomodate mathematical expressions in labelling.


texts:
After hours and days of trying to build TileMill from source on Windows, I'm more than frustrated. Is there anyone willing to share his/her up-to-date TileMill for Windows with me?

I tried to build a version a few weeks ago. But basically, I'm only interested in a version which supports setting a transparent color for raster data (which current 0.10.1 for Windows does not). It's "image-filters:colorize-alpha", that I'm interested in.


texts: Set default values for value table in Python Toolbox tool

texts:
I created a Python Toolbox tool to reorder fields and create a new feature class with the reordered fields.  The tool works nicely and I am able to use a value table to let the user arrange the fields in the order they choose or they can fill in a rank value for each field.  However, the annoying part of this tool is that all fields need to be added to the value table one at a time before reordering.

I am trying to set this up to bring in all the fields to the value table by default and any unwanted fields can be removed before reordering them.  Has anyone had success doing something like this before?  I am trying to achieve this in the UpdateParameters method.  Here is the code I am trying:

import arcpy
import os

class Toolbox(object):
    def __init__(self):
        """Define the toolbox (the name of the toolbox is the name of the
        .pyt file)."""
        self.label = "Reorder Fields"
        self.alias = "Reorder Fields"

        # List of tool classes associated with this toolbox
        self.tools = [ReorderFields]

class ReorderFields(object):
    def __init__(self):
        """Define the tool (tool name is the name of the class)."""
        self.label = "Reorder Fields"
        self.description = ""
        self.canRunInBackground = False

    def getParameterInfo(self):
        """Define parameter definitions"""
        fc = arcpy.Parameter(displayName='Features',
            name='features',
            datatype='Feature Layer',
            parameterType='Required',
            direction='Input')

        vt = arcpy.Parameter(
            displayName='Fields',
            name='Fields',
            datatype='Value Table',
            parameterType='Required',
            direction='Input')

        output = arcpy.Parameter(
            displayName='Output Features',
            name='output_features',
            datatype='Feature Class',
            parameterType='Required',
            direction='Output')

        vt.columns = [['Field', 'Fields'], ['Long', 'Ranks']]
        vt.parameterDependencies = [fc.name]
        params = [fc, vt, output]
        return params

    def isLicensed(self):
        """Set whether tool is licensed to execute."""
        return True

    def updateParameters(self, parameters):
        """Modify the values and properties of parameters before internal
        validation is performed.  This method is called whenever a parameter
        has been changed."""
        if parameters[0].value:
            if not parameters[1].altered:
                fields = [f for f in arcpy.Describe(str(parameters[0].value)).fields
                          if f.type not in ('OID', 'Geometry')]
                vtab = arcpy.ValueTable(2)
                for field in fields:
                    vtab.addRow("{0} {1}".format(field.name, ''))
                parameters[1].value = vtab
        return

    def updateMessages(self, parameters):
        """Modify the messages created by internal validation for each tool
        parameter.  This method is called after internal validation."""
        return

    def execute(self, parameters, messages):
        """The source code of the tool."""
        fc = parameters[0].valueAsText
        vt = parameters[1].valueAsText
        output = parameters[2].valueAsText

I want to bring all the fields as shown in the value table above by default.  I also tried using the parameters[1].value to add rows the the specific value table from the GUI, but that gave me errors.  I'm using ArcGIS 10.2.2.


texts: Extracting index/position of maximum value in annual time series in Google Earth Engine?

texts:
I am trying to extract the integer index/position of the maximum value in an annual MODIS NDVI time series in Earth Engine but I am having a difficult time to figure out a plan to do this. It is quite easy to calculate the maximum value of the image collection using the following javascript code:

var bbox = ee.Geometry.Polygon([
  [-8.44299,42.20410],
  [-8.43200,41.65239],
  [-7.75634,41.65239],
  [-7.76733,42.20817],
  [-8.44299,42.20410]]);

var getNDVI = function(img){
  return img.normalizedDifference(['sur_refl_b02','sur_refl_b01']);
};

var MOD = ee.ImageCollection('MODIS/006/MOD09A1').
          filterDate('2015-01-01','2015-12-31').
          filterBounds(bbox);

var NDVI = MOD.map(getNDVI);

var NDVImax = NDVI.max();

For example, in R raster package this would correspond to an operation with calc() or stackApply() with which.max() as a function to retrieve the position of the maximum value in the raster stack time series. However, I would prefer to do this directly in GEE without having to download the data.

Suggestions?


texts: PHP/MySQL, PostgreSQL spatial queries. Find a polygon that contains a coordinate?

texts:
I am building a non mapping javascript app and I am wanting to take a lat/long coordinate and query a MySQL table of polygons for the one that contains the point.

I have looked into MySQL spatial datatypes and how to query it, however, it seems to me that you can only return a true/false value for whether or not a specific point is within a specific polygon. What I want to do is look in a massive table for the correct polygon. My goal is to find the polygon.

I have also looked into maybe using PostgreSQL and PostGIS, but haven't found a solution yet.

Can anyone here point me in the right direction?


texts: PostGIS erase/difference vector polygons from raster

texts:
I have a raster and a vector table in Postgres 12 PostGIS 3.0

The vector polygon(orange) overlaps the slope raster(grid in black, whiteish)

All I want to do is cut out/erase the orange polygon from the black/white raster.

**My question is just focused on this specific raster and vector but in general this going to be a big process of hundreds of thousands of polygons erasing slope rasters

My process so far:

Convert the polygon to a raster
use MapAlgebra and subtract the layers which is essentially acting as a clip

query

drop table if exists corners.rast_test;
create table corners.rast_test as
with a as(
SELECT ST_AsRaster(m.shape, (SELECT rast FROM corners.slope_clip LIMIT 1), '32BF', -9999) rast,id
FROM (select m.shape,m.id
from corners.mappluto m
join corners.corners_final c on st_intersects(m.shape,c.n9x9) and m.id=2019) m
)
SELECT ST_MapAlgebra(b.rast,a.rast, '([rast1]-[rast2])') rast,cornerid
from a, corners.slope_clip  b where cornerid=1085547

the result

So this is the piece of the raster I do not want. I also want to keep the original slope value during this whole process

I found the expression I need to use from this question I posted about the ArcGIS Raster Calculator

Con(IsNull(rast2),rast1)

how do I turn this above expression into the correct st_mapalgebra statement


texts: How to copy values from field to another using arcpy?

texts:
Can some one please let me know how I can do simple field calculation like Copying From "A" to "B" in ArcPy?

I have found lots os example on web which they all used extra expression and they were complicated.what I need is simply doing what we can do in ArcMap GUI like

B = !A!

update:

Here is the code I have so far but I am getting error on running this as

  NameError: name 'A' is not defined

point_shp = "G:\\Temp\\All_Provinces.shp"
arcpy.AddField_management(point_shp, "B", "TEXT", "", "", "25", "", "NON_NULLABLE", "NON_REQUIRED", "")
arcpy.CalculateField_management(point_shp, "B", "A", "PYTHON_9.3")


texts: Leaflet GeoJSON bindtooltip from feature attribute

texts:
Leaflet allows you to use .bindtooltip for geoJSON. Is it possible to access the tooltip options like {permanent: true, opacity: 0.5} etc when using this method?

// My JS
var layers_grid = L.geoJSON(false, {
      style: function(feature) { return feature.properties}
})
.bindTooltip(function (layer) {
    return layer.feature.properties.name;
 }
).addTo(map);

// Later in the script
$.getJSON(url, function( data ) {
  layers_grid.addData(data);
});

Edit: JSFiddle


texts: How to view and change the linear units of a feature class in ArcGIS 10.0

texts:
According to the documentation for Near Analysis the output units for the tool are in the linear units of the feature.

When accessing the properties of the feature class I can see where angular units are defined (as a result of the coordinate system), however I do not see where I can view or change the linear units of a feature.

The Data Source textbox from the Source tab for the Layer Properties of my point feature class does not contain a Linear Unit specification.

How do I access the linear units of a feature class?


texts: single sided buffer in pyqgis 3 gives "_core.QgsProcessingException: Error calculating single sided buffer" error

texts:
I am trying to create a single sided buffer around drains (polygon shapefile). The Distance parameter must be taken from a field in drains shp. Here is my code:

buf_out = self.dlg.lineEdit_3.text()
b1 = drainLayer.fields().indexFromName('b1')
processing.run('qgis:singlesidedbuffer',{'INPUT': drainLayer, 'DISTANCE': b1,'SIDE':1, 'SEGMENTS':25, 'JOIN_STYLE':0, 'OUTPUT':buf_out})

I tried giving an int value for distance but even that didnt work. An empty shapefile is created and qgis processing exception is raised. Normal buffer is working but one sided buffer is not working.


texts: Uploading large shapefiles to ArcGIS Online?

texts:
I am trying to upload a large shapefile to ArcOnline, however I found out ArcGIS Online won't upload any shapefile that has more than 1000 features.  My 30 shapefiles have over 1,000,000 features each.

Is there a way for me to condense these files into a smaller number of features that won't require me to go through manually and break the file down into individual shapefiles of 1000 features?


texts: Changing opacity for KML layer fetched in Leaflet map

texts:
Is it possible to change the opacity for the .kml file, that has been attached via leaflet-plugin?

Basically the plugin is here:

https://github.com/windycom/leaflet-kml

but after implementation into the Leaflet we have a full opaque image. It's fine for placemarks, but how about the icons/images included, which cover some box on the map?

My code looks like this:

var kml = fetch('Coventry.kml')
      .then( res => res.text() )
      .then( kmltext => {

            // Create new kml overlay
            parser = new DOMParser();
            kml = parser.parseFromString(kmltext,"text/xml");

            console.log(kml)

            const track = new L.KML(kml)
            map.addLayer(track)
            // Adjust map to show the kml
            const bounds = track.getBounds()
            map.fitBounds( bounds )
      }).setOpacity(0.6).addTo(map)

      $('#sldOpacity').on('change', function(){
$('#image-opacity').html(this.value);
kml.setStyle({opacity: this.value, fillOpacity: this.value});
 });

whereas to the pure code straight from the leaflet-kml example I added up sth like:
      .setOpacity(0.6).addTo(map)

but isn't working.

The console states, that:

Uncaught TypeError: fetch(...).then(...).then(...).setOpacity is not a function
    at (index):143

I was also trying this example:

Change opacity using leaflet without plugin

but it didn't help me unfortunately

In the L.KML.js file I found a snip:

 var value = e.childNodes[0].nodeValue;
                if (key === 'color') {
                    options.opacity = parseInt(value.substring(50, 50), 16) / 255.0;
                    options.color = '#' + value.substring(6, 8) +
value.substring(4, 6) + value.substring(2, 4);
                } else if (key === 'width') {,

where in value.substring I put 50 instead of 0. No result at all.

Is it some way to make this .kml layer more transparent?


texts:
It is possible to create a style with default color but with variable quantity with GeoServer?
I've a lot of grids that reproduce surface temperature. I don't know previously the minimum and maximun temperature but I would like to use Magma colormap from QGIS.
My aim is divide in ten classes the variability of temperature like this:

(MaxTemperature - MinTemperature) / 10

Then, with this classes I need to obtain a dynamic coloramp.
Something like this but I've no clear what is ${env
Now I'm be able to use a JSON from GeoSever with which I can create a legend from the styles with OpenLayers and I hope that is possible to do the same thing also whit the dynamic coloramp.


texts: Efficiently referencing/connecting data without overlap

texts:
I am doing an analysis that places a buffer around X property, then identifies all the parcels within that buffer to store in a database and perform analysis on later.

Currently, I have two separate databases. One for all the properties I run the buffer on, then a second which includes all of the parcels that are within the buffer. Each record in the second database has a reference field to the first so I know where they all belong. The issue is there is massive duplication of properties in the second database since many of the buffers would end up overlapping each other. In turn I'm adding the same parcels to the second database, they just reference different records to the first. I feel like this is hugely inefficient.

Is there a better way to do this?

I thought I could just add a field to the Parcel data and use Calculate Field to just append any reference fields together. Is my "Better Way" suggested below the best way to do this, if so is there a reason why the Model Builder does not properly select the parcels so I can then calculate the field?


texts: Replicate elevation profile extraction (3D Analyst toolbar) in ArcPy

texts:
How can I replicate in ArcPy the 3D Analyst Toolbar operations which extract an elevation profile from a transect line drawn atop an elevation layer? This will be incorporated into an automated workflow that will loop ~500K times. Hence, streamlining is critical.

Here's the procedure using the toolbar:

Point toolbar to the elevation layer
Click 'Interpolate Line' to draw the cross-section
Click 'Profile Graph' to generate the profile
Right-click graph, 'Export', 'Data'

(Inputs: elevation layer, transect. Outputs: graph, XY table.)

Rather than manually drawing a transect, as in the toolbar functions, I need the ability to specify the transect start & end by entering coordinates. Similar to the tool, I want the output to be XY pairs (ideally a numpy array), with x=0 being the transect origin, X-values being distance from origin, and Y-values being elevation.


texts: How can i publish polygon layer on Geoserver with polygon geometry

texts:
I can not publish polygon layer (ESRI Shapefile) on Geoserver with Polygon geometry.

Every time is published as MultiPolygon. Shapefile has really Polygon geometry.


texts: Extracting values from raster stack and aggregating results using R

texts:
I'd like to extract values from a rasterstack, cell by cell, and create a list of the values:

ext <- extent(0,100,0,100)
r1 <- raster(nrows=100, ncols=100,ext)
r1[] <- sample(seq(from = 1, to = 6, by = 1), size = 10000, replace = TRUE)
r2 <- raster(nrows=100, ncols=100,ext)
r2[] <- sample(seq(from = 1, to = 6, by = 1), size = 10000, replace = TRUE)
r3 <- raster(nrows=100, ncols=100,ext)
r3[] <- sample(seq(from = 1, to = 6, by = 1), size = 10000, replace = TRUE)
r4 <- raster(nrows=100, ncols=100,ext)
r4[] <- sample(seq(from = 1, to = 6, by = 1), size = 10000, replace = TRUE)
s <- stack(r1,r2,r3,r4)

using extract I can create a matrix telling me the values in each layer of the stack, cell by cell:

y <- extract(s,c(1:ncell(s)))

now what I want to do is return only the unique rows of the matrix (unique(y)) with a count of how many times that row occurs.

In other words, stick a pin through every cell in a rasterstack and count how many times each 'vector' of values occurs, and also what those vectors are.


texts:
I would like to use cartopy to display a map of the cities in the contiguous United States. In order to do so, I have written the following code (the data comes from here):

import geopandas as gpd
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

cities = gpd.read_file("./data/citiesx010g.shp")
cities_m = cities[cities['STATE'].map(lambda s: s not in ['PR', 'AK', 'HI', 'VI'])]

fig = plt.figure()
ax = plt.subplot(111, projection=ccrs.AlbersEqualArea())
ax.scatter([p.x for p in cities_m.geometry], [p.y for p in cities_m.geometry], transform=ccrs.PlateCarree())
plt.show()

This is the output:

This leads to the following three questions:

In this case, what would I have to do to "flip" this plot right side up?
In the general case of having some kind of data in a AlbersEqualArea projection, what would I have to do to "flip" that data right side up?
In the most general case of having just about any projection, how would I go about "flipping" the data right side up? (only the naive PlateCarree always comes out "right").


texts: get and change style property openlayers

texts:
i need to change the property of some style after its created, and i need to change it with knowing only the string of the property name.
i am getting string like "fillColor" and value like "blue", and i need to change the style object. how can i do it? i didnt find anything in the api.
the style code:

            var selectStyle = new OpenLayers.Style({
            externalGraphic: 'styles/tree_selected.png',
            graphicWidth: 16, graphicHeight: 26, graphicYOffset: -24,
            title: 'tree'
        });


texts: Using expression to number unique string values consecutively in QGIS?

texts:
I have a set of string values in a shapefile's attribute column and want to assign a number to them in a way that each unique string gets an integer number, like this:

strings        numbers
-------        -------
Cloud          1
Water          2
Water          2
Vegetation     3
Water          2
Cloud          1
Cloud          1
Vegetation     3

The order does not matter (but of course, if it's possible, alphabetically would be perfect).

Is there a way to achieve this without knowing the actual strings or the count of unique strings (assuming I have 20-30 different strings)?

I am using Python in QGIS 3.

Using the suggestion of BERA, I was able to run the Python script from the console. Then I tried to implement it as an expression:

from qgis.core import *
from qgis.gui import *
from qgis.utils import iface

@qgsfunction(args='auto', group='Custom')
def numberUniqueStrings(layer_name, src_column, dst_column, feature, parent):
    layer_list = QgsProject.instance().layerTreeRoot().children()
    layer = [lyr.layer() for lyr in layer_list if lyr.name() == layer_name][0]
    idx = layer.getFeature(1).fieldNameIndex(dst_column)
    unique = []
    with edit(layer):
        for feature in layer.getFeatures():
            value = str(feature[src_column]).strip()
            if value not in unique:
                unique.append(value)
            counter = unique.index(value) + 1
            feature.setAttribute(idx, counter)
            layer.updateFeature(feature)

and called it using the "Field Calculator" (created the field "numbers" as integer before):

numberUniqueStrings(@layer_name, "strings", "numbers")

Unfortunately nothing happens, all the fields stay empty.


texts:
I tried to sent request to server Overpass API:
req.body()= "(node(51.249, 7.148, 51.251, 7.152);< ;);out meta;";

I got this error:

The data included in this document is from www.openstreetmap.org. The data is made available under ODbL. Error: encoding error: Your input contains only whitespace.


texts: Request Multiple WMS for Raster Calculations in Leaflet JS

texts:
I need to display in leafletJS a Layer that is the difference of two rasters received from a regular WMS. With the help of this repo I figured out how to manipulate the values of a raster tile in Leaflet using L.TileLayer.extend:
import L from 'leaflet';

L.TileLayer.Manipulated = L.TileLayer.extend({

  initialize(url, options) {
    L.TileLayer.prototype.initialize.call(this, url, options);
    this.on('tileload', (e) =>  {
      if (!e.tile.getAttribute('manipulated-tile')) { // If tile is not manipulated yet
        this._manipulate(e.tile);
      }
    });
  },
  _manipulate(img) {
    const canvas = document.createElement("canvas");
    canvas.width = img.width;
    canvas.height = img.height;
    const ctx = canvas.getContext("2d");
    ctx.drawImage(img, 0, 0);

    const imageData = ctx.getImageData(0, 0, canvas.width, canvas.height);
    const pix = imageData.data;
    for (var i = 0, n = pix.length; i < n; i += 4) {
      pix[i] = pix[i] * 0.8; // red
      pix[i + 1] = pix[i + 1] * 0.8; // green
      pix[i + 2] = pix[i + 2] * 0.6; // blue
    }
    ctx.putImageData(imageData, 0, 0);
    img.setAttribute('manipulated-tile', true);
    img.src = canvas.toDataURL();
  },
});

L.tileLayer.manipulated = function (url, options) {
  return new L.TileLayer.Manipulated(url, options);
};

Usage:
L.tileLayer.manipulated('https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png');

Does anybody know a simple approach how I can get two geographically matching tiles of different layers and then calculate the difference between those? The calculations need to be done in the frontend.


texts: Loop through a range of rasters and use tool output as input for the next tool run

texts:
I'm using the particle track tool and limiting it to a specific trackingtime. For each tool run, I want the rasters input (magnitude and direction) to change to the next rasters in the lists, and the starting point of the tracking (x,y) to be derived from the last polyline point of the previous run.
This is what I have so far but I get an error. It works fine if I remove the section: for row in cursor: but then it won't loop.

import os
import arcpy
from arcpy import*
from arcpy.sa import*
arcpy.CheckOutExtension ('Spatial')
arcpy.env.overwriteOutput= True

arcpy.env.workspace= r'C:/Time_slices/January_newrun/Magnitude2/Rectified/Stripped'
mag= arcpy.ListRasters('*.tif')
arcpy.env.workspace= r'C:/Time_slices/January_newrun/Direction360'
direction= arcpy.ListRasters('*.tif')

sourcePoint = arcpy.Point(418350.56, 84472.21)
outTrackFile = "C:/Time_slices/January_newrun/trackparticle/trackfileMod10.txt"
trackingTime = 86400
OutputFolder = r'C:/Time_slices/January_newrun/trackparticle'

for idx in range(0, len(mag)):
    inMagnitudeRaster= mag [idx]
    inDirectionRaster= direction [idx]
    basename= inDirectionRaster.split(".tif")[0]
    outTrackPolylineFeatures = os.path.join (OutputFolder, basename + '.shp')
    ParticleTrack (Raster(inDirectionRaster), os.path.join(r'C:/Time_slices/January_newrun/Magnitude2/Rectified/Stripped',inMagnitudeRaster), sourcePoint,outTrackFile, "#", trackingTime,outTrackPolylineFeatures)
    arcpy.AddField_management (outTrackPolylineFeatures, "LastX", "DOUBLE")
    arcpy.AddField_management (outTrackPolylineFeatures, "LastY", "DOUBLE")
    arcpy.CalculateField_management(outTrackPolylineFeatures, "LastX", "(!SHAPE.lastpoint.X!)", "PYTHON_9.3")
    arcpy.CalculateField_management(outTrackPolylineFeatures, "LastY", "(!SHAPE.lastpoint.Y!)", "PYTHON_9.3")
    cursor= arcpy.SearchCursor(outTrackPolylineFeatures)
    for row in cursor:
           endx= row.getValue("LastX")
           endy= row.getValue("LastY")
    sourcePoint= arcpy.Point(endx, endy)
    inMagnitudeRaster= mag[idx+1]
    inDirectionRaster= direction[idx+1]
    basename= inDirectionRaster.split(".tif")[0]
    outTrackPolylineFeatures = os.path.join (OutputFolder, basename + '.shp')
    ParticleTrack (Raster(inDirectionRaster), os.path.join(r'C:/Time_slices/January_newrun/Magnitude2/Rectified/Stripped',inMagnitudeRaster), sourcePoint,outTrackFile, "#", trackingTime,outTrackPolylineFeatures)

The error message I get is:

Traceback (most recent call last):
  File "C:/particle_cursor_15_02.py", line 40, in <module>
    ParticleTrack (Raster(inDirectionRaster), os.path.join(r'C:/Time_slices/January_newrun/Magnitude2/Rectified/Stripped',inMagnitudeRaster), sourcePoint,outTrackFile, "#", trackingTime,outTrackPolylineFeatures)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.6\ArcPy\arcpy\sa\Functions.py", line 2149, in ParticleTrack
    out_track_polyline_features)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.6\ArcPy\arcpy\sa\Utils.py", line 53, in swapper
    result = wrapper(*args, **kwargs)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.6\ArcPy\arcpy\sa\Functions.py", line 2140, in Wrapper
    out_track_polyline_features)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.6\ArcPy\arcpy\geoprocessing\_base.py", line 510, in <lambda>
    return lambda *args: val(*gp_fixargs(args, True))
ExecuteError: ERROR 999999: Error executing function.
Failed to execute (ParticleTrack).


texts: Adding classes to raster visualization using QGIS

texts:
This is my first time working with rasters in QGIS. I'm working with a single band gray raster.
It has 5 classes. When I open it in QGIS, I can see the gray scale for the 5 classes. However, in the layers panel I can only see 2 classes, which are the 2 extremes (white and black).
How can I visualize all 5 classes in the layers panel?


texts: Presenting geospatial data in Powerpoint

texts: I am looking for a more elegant way to display geospatial data in Powerpoint.  KMLs seem to be becoming the norm for quickly showing interactive geospatial data to colleagues and the general public. Is there a better way to integrate these incredibly useful KMLs into Powerpoint presentations than the example I've shown below?  Maybe there is a better solution for presenting geospatial data altogether in Powerpoint--I look forward to your ideas.


texts: Is there a typescript interface for "GeoXML" formats (such as WFS and WMS) apart from GeoJSON?

texts:
I want to work with JSON behind TypeScript interfaces but I'm not sure what the spec should be outside of GeoJSON.  Is there a spec which is more encompassing then GeoJSON and includes things like a GetCapabilities response?

For example, it is very obvious with ESRI because ArcGIS Server services have built-in json support via adding a query string.  For example, http://sampleserver3.arcgisonline.com/ArcGIS/rest/services?f=json&callback=foo

But GeoServer limits json support to GetFeature/GetFeatureInfo requests so it's not obvious what a GetCapabilies response might look like if a tool were to actually produce JSON.

So my question is this:

Is GeoJSON format limited to spatial data?
If so, is there any other JSON standards/schemas related to WFS/WMS services that is more encompassing?

NOTE: OpenLayers has WFS/WMS readers for various services.  Presumably the resulting JSON conforms to some standard?

https://github.com/openlayers/ol3/blob/master/test/spec/ol/parser/ogc/wfs_v1_1_0.test.js
https://github.com/openlayers/ol3/blob/master/test/spec/ol/parser/ogc/wmscapabilities_v1_0_0.test.js


texts: Layers Labeled "Override" in Printed PDF Legend

texts:
I have a Javascript function that creates a GraphicsLayer, adds some features to it, and then adds the layer to my ArcGIS Web Map. When I use the Print Widget to print my map I see the layer show up in the Legend labeled "Overrride 1". Is there a way to add a title or name to a GraphicsLayer so that it will be properly rendered by the print task?

Here's a plunk to demonstrate: http://plnkr.co/edit/DM7U7H


texts: QGIS Modeler Selecting based on Geometry Shape: Line, Point and Polygon

texts: I am in the process of building a model, where three layers (in a batch progress) there is currently one. Go through the same Refactor Field Alogrithm. However to make it more compeltly and to kind of prevent new empty layers being generated.I need to have a selection based on the Type of Geometry, for example the Ducts and Cabinets are supposed to be points so I need to do a selection based on point. I think I am supposed to select by attribute and then Geometry but is that really the sollution?


texts: Using GDAL in Windows with command longer than 8191 characters

texts: I have a large number of files, over 7000. I need to run GDAL commands on those files, but I am limited by 8191 characters in Windows. Is there any way to avoid this limitation?


texts: Why does fasterize gives different output for multipolygons with and without internally resolved boundaries?

texts:
I was thrilled to be clued into the package fasterize, which speeds up the rasterization process in R tremendously.

However, I cannot figure out why the fasterize() function outputs two different rasters when generated from a multipolygon with internally resolved boundaries and a multipolygon without internally resolved boundaries. I have tried specifying fun = "any" in the fasterize parameters, as well as fun = "sum". The former should (I believe) take any area where there's a polygon; the latter should take the sum of any overlapping polygons. In both cases, the coverage should be the same. However, the raster from the multipolygon generated from st_combine() drops overlapping areas whereas the raster from the multipolygon generated from st_union() keeps any/all areas.

Am I missing something? Why are there holes in the middle plot (rasterized from st_combine()?? There's reproducible code below the plot, with sample shapefile here (of Yellowstone and Grand Teton NPs -- provided as separate shapefile components and as zipped folder).

# Read in shapefile.
YNPGT <- st_read("YNPGT.shp")
bbox <- st_bbox(YNPGT) # Keep bounding box for plotting b/c global extent
# Rasterize. "Sum" so any spot with poly is retained (e.g., not just 'first' poly).
YNPGT.r <- YNPGT %>% fasterize(raster = ref.ras, fun = "sum")

# Combine that to flatten all polygons into one -- resolves internal boundaries
combo <- st_combine(YNPGT)
# Rasterize
combo <- st_sf(combo) # combine outputs sfc; need sf for fasterize.
combo.r <- combo %>% fasterize(raster = ref.ras, fun = "sum")

# Union that to flatten all polygons into one -- drops internal boundaries.
union <- st_union(YNPGT)
# Rasterize
union <- st_sf(union) # combine outputs sfc; need sf for fasterize.
union.r <- union %>% fasterize(raster = ref.ras, fun = "sum")

par(mfrow=c(2,3))
plot(st_geometry(YNPGT),
     xlim=c(bbox[1], bbox[3]), ylim=c(bbox[2], bbox[4]),
     col = c("green", "pink"))
title("YNP-GT orig")
plot(st_geometry(combo),
     xlim=c(bbox[1], bbox[3]), ylim=c(bbox[2], bbox[4]),
     col = c("green", "pink"))
title("YNP-GT  st_combine")
plot(st_geometry(union),
     xlim=c(bbox[1], bbox[3]), ylim=c(bbox[2], bbox[4]),
     col = c("green", "pink"))
title("YNP-GT  st_union")
plot(YNPGT.r,
     xlim=c(bbox[1], bbox[3]), ylim=c(bbox[2], bbox[4]),
     col = c("green", "pink"))
title("YNP-GT raster from orig")
plot(combo.r,
     xlim=c(bbox[1], bbox[3]), ylim=c(bbox[2], bbox[4]),
     col = c("green", "pink"))
title("YNP-GT raster from st_combine")
plot(union.r,
     xlim=c(bbox[1], bbox[3]), ylim=c(bbox[2], bbox[4]),
     col = c("green", "pink"))
title("YNP-GT raster from st_union")


texts: Copy multiple rasterfiles into geopackage

texts:
I have a folder with ~100+ raster images which covers my study area (the data comes in .ecw or .tiff, i have both at hand). How do i get them into one geopackage without having to go the export -> geopackage -> select the geopackage -> select a layer name etc. way for every single file?

So far i tried to create a .vrt with all the rasters and export this into the gpkg, which suceeded, but when i want to open it, it is loading for ages. Is there a better way to achieve this?


texts: Using Arcade in ArcGIS Online to customize labeling?

texts:
I am trying to clean up the labeling for a city zoning map using ArcGIS Online. To do this I need to trim the property address to only the numeric values, which has been accomplished, but due to incomplete data some values in the Property_Address field do not have street address number so the the expression will display the first four characters of the name of the street address. I need this second expression to be able to determine if numerical or text characters will be output and if there are no numerical characters and text is present, then those characters need to be hidden.

Left($feature.PropertyAd,4)

IIF(Left($feature.PropertyAd,4) !=


texts: Why did my Visualization tools disappear from the CartoDB Sidebar?

texts:
I was working with the "Intensity" visualization tool yesterday.

Using the same data, I created a new map, and now I can only access: Simple, Choropleth, Category, and Bubble.

Where did Intensity, Torque, and the others go?


texts: Advanced android app for mapping and GPS logging?

texts:
I am looking for an android app able to log GPS data like coordinates, altitude, current speed and would be nice if also the GPS accuracy/number of satellites. Many GPS apps are able to display these data, but not to log them.

Reason is that while doing a monitoring (no matter if helicopter, car or on foot), there is a variable accuracy of the coordinates and so this also influences the measurement results and the output maps. So it would be helpful to know these data to take this into consideration.

I am currently using Oruxmaps, but it logs only date/time, position and elevation and not the other data. Just found Antimap and Geopaparazzi which I will surely try, but I am after other recommendations.

I also need to work in offline mode.


texts: Legend Graphic for a multi polygon layer styled the same with unique (different) labels

texts:
When creating an SLD for a single polygon layer in QGIS and using that in Geoserver, by default the legend description is "single symbol". When you categorize the single polygon layer symbology and classify that based on a single attribute, you then can change the symbol name to whatever name or description you would like. This solves the single symbol description issue for your GetLegendGraphic displaying in either OpenLayers/Leaflet. This works for a single polygon. Multiple polygons layer categorized based on unique values and different styling also work fine. The challenge I have is creating an SLD for a multi polygon layer that shares the same styling with unique labels so that on the legendgraphic/png, there is a single symbology and description next to it. Is there a solution that I can possibly do in QGIS without having a polygon SLD and creating a point layer for the labels?

This is the single polygon layer in Leaflet:

This is what I would like to have in leaflet, basically the same as the single polygon but with different labels, tried categorizing but did not work:


texts: I have a city map with digital boundaries. Each polygon has some numbers between 1 to 7. For example one of the polygon has 4,4,6,1,4,7,3,4,5,4,2 (the mode is 4 for this polygon). I wonder is there a way in Arcmap to find mode value for each polygon and make it colored by these number's value. (i.e. 1-yellow, 2- red, 3-green, etc..)?


texts: Opening multiple instances of ArcGIS Pro?

texts:
Is there a way to run concurrent ArcGIS Pro instances under the same window account?

It appears that ArcGIS Pro 1.x doesn't allow this with the default configuration.


texts:
I would like to display my GML file but I don't understand why nothing appear in the map:

code:

import Map from 'ol/Map'
import Layer_Tile from 'ol/layer/Tile'
import source_OSM from 'ol/source/OSM'
import View from 'ol/View'
import VectorS from 'ol/source/Vector'
import VectorL from 'ol/layer/Vector'
import WFS from 'ol/format/WFS'
import GML from 'ol/format/GML32'

var map = new Map({
  layers: [
    new Layer_Tile({
      source: new source_OSM()
    })
  ],
  target: 'map',
  view: new View({
    center: [0, 0],
    zoom: 2
  })
});

const vectorSource = new VectorS({
    url: 'JeuxTestv2.gml',
    format: new WFS({
      gmlFormat: new GML({
        srsName: 'EPSG:3948'
      })
    }),
});
var vector = new VectorL({
  source: vectorSource
});
map.addLayer(vector);
vectorSource.on('change', e => {
  if (vectorSource.getState() === 'ready') {
    map.getView().fit(vectorSource.getExtent());
  }
})

My GML file:

<?xml version="1.0" encoding="UTF-8"?>
<!---JEUX TESTS DU PCRS V2 : Données fictives-->
<PlanCorpsRueSimplifie xmlns="http://cnig.gouv.fr/pcrs" xmlns:gml="http://www.opengis.net/gml/3.2" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:xlink="http://www.w3.org/1999/xlink" gml:id="SchemaPCRSv2" version="2.0" xsi:schemaLocation="http://cnig.gouv.fr/pcrs ../CNIG_PCRS_v2.0.xsd">
    <gml:boundedBy>
        <gml:Envelope srsName="EPSG:3948">
            <gml:lowerCorner>2050999.431282464 7273998.947281252</gml:lowerCorner>
            <gml:upperCorner>2051701.173832032 7274399.972199726</gml:upperCorner>
        </gml:Envelope>
    </gml:boundedBy>


texts: Controlling XY Resolution on file geodatabase feature class output from Append?

texts:
This question arises from research into Appending polygon feature classes with already repaired geometry => WARNING 000442: self intersections from Check Geometry?

Before running the code snippet below the properties of feature class rwc_aurukun show:

inputFileGDB = <file_geodatabase_where_all_work_is_done>
arcpy.env.workspace = inputFileGDB
fcList = <list_of_feature_classes_including_"rwc_aurukun"_as_first_member>
arcpy.CreateFeatureclass_management(inputFileGDB,"RWC","POLYGON",fcList[0])
coord_sys = arcpy.Describe(fcList[0]).spatialReference
arcpy.DefineProjection_management("RWC", coord_sys)
arcpy.env.XYDomain ="-180 -90 180 90"
arcpy.env.XYTolerance = "0.000000001 DecimalDegrees"
arcpy.Append_management(fcList,"RWC","NO_TEST")

After running the code snippet the properties of feature class RWC (which is 88 feature classes appended - including rwc_aurukun) show:

How can I control the XY Resolution on the output feature class (RWC)?

My experimentation with arcpy.env.XYDomain and arcpy.env.XYTolerance are not achieving the result that I am after.


texts: Get intersects of indexed grouped polygons between 2 SpatialPolygonDataframes in R

texts:
I am testing 2 different methods (m2b & hmm) for classifying the movement behaviour of penguins to ultimately classify feeding behaviour from GPS points and from these create kernel polygons to identify foraging areas. I want to check the overlap (area) between the foraging kernels from the 2 different methods. I have 1 SpatialPolygonsDataframe for each method which has a number of foraging kernels for each animal ID. Ultimately for each ID (which can be > 1 polygon) I want to calculate the overlap area between the 2 SpatialPolygonsDataframes.
I have tried using rgeos::gIntersection:

o1 <- gIntersection(ud50_m2b,ud50_hmm,byid = T, id = 'ID')

but it returns the following error:

Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, unaryUnion_if_byid_false,  :   ids vector has incorrect length of 1 expected length of 35910

My data looks like this:

> str(ud50_m2b)
Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
  ..

@ data       :'data.frame': 2 obs. of  4 variables:
  .. ..$ ID    : Factor w/ 2 levels "CB1201","CB1203": 1 2
  .. ..$ area  : num [1:2] 1971 1178
  .. ..$ survey: Factor w/ 15 levels "","April2012BI",..: 2 2
  .. ..$ n     : int [1:2] 3024 4326
  ..

@ polygons   :List of 2
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 3
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 435288 6250135
  .. .. .. .. .. .. ..@ area   : num 5857355
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:30, 1:2] 434056 433968 433882 433922 434056 ...
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 442476 6251628
  .. .. .. .. .. .. ..@ area   : num 13697898
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:64, 1:2] 439588 439576 439588 439740 439983 ...
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 442305 6247942
  .. .. .. .. .. .. ..@ area   : num 156406
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:8, 1:2] 442353 442112 442185 442353 442474 ...
  .. .. .. ..@ plotOrder: int [1:3] 2 1 3
  .. .. .. ..@ labpt    : num [1:2] 442476 6251628
  .. .. .. ..@ ID       : chr "CB1201"
  .. .. .. ..@ area     : num 19711659

I think there is an indexing issue here in that the gIntersection function does not register the same IDs for multiple polygons in each SpatialPolygonsDataframe. Is there a trick to getting the overall intersection area for each ID between the 2 SpatialPolygonsDataframe?


texts: Creating checkbox parameter in Python script tool for ArcGIS Desktop?

texts:
I am working on creating an ArcGIS tool from a Python script I am writing.  I am wondering if it is possible to have a checkbox parameter.

I want to have a parameter where the user selects a feature class, then from the feature class the user will choose the field for the upper most layer in their model, then I want the user to be able to choose what layers they want the script to run on with a checkbox structure derived from the upper most layer field.

Is this possible with python, and ArcGIS Desktop?


texts:
Directions Magazine, Pitney Bowes (MapInfo), and many blogs are talking about 'Cloud' in relation to GIS. I have read about Cloud-Computing but still don't understand the concept.

What is the GeoCloud that everyone is talking about?


texts: Openstreetmap Nominatim web service API - Result language

texts:
I notices something strange:

When I execute simple request:

https://nominatim.openstreetmap.org/reverse?format=xml&lat=45&lon=10&zoom=18&addressdetails=0

From the PC web browser for location located in my country (not 45/10) I get the result in English, although we have different language here.

But when I execute the very same request, without any other data like headers from application running on a mobile device - I do get the result in our local language.

Why this is happen? I familiar with the accept-language param of Nominatim, but as seen - I'm not specify him nor touching headers in my request.


texts: Counting number of points in polygon using R?

texts:
I have two classes sharing the same CRS (Latitute and Longitude):

bolognaQuartieriMap: a SpatialPolygonDataFrame containing data of a city boroughts.
crashPoints: a SpatialPointsDataFrame containing data of accidents.

They are well plotted using:

plot(bolognaQuartieriMap)
title("Crash per quartiere")
plot(crashPoints, col="red",add=TRUE)

What I need is to get the number of points (crashPoints) in each polygon that constitute bolognaQuartieriMap. I was suggested to use over() but I did not succeed.


texts:
I have a grid layer and a shapefile representing a closed boundary. I need to make a new grid as big as my current grid but with values equal 1 inside the closed boundary shapefile and -999 for grids outside the closed boundary. Would any body help me with this using arcmap 10?

Thanks a lot.