texts: How to get route plan with certain distance with pgRouting?

texts:
I am working for a project. I have a road network which has been imported in PostgreSQL/PostGIS database. I can query shortest path between two points. I want to get route plan with certain distance. Suppose, I want to go from point A to point B. There may have different way to go from point A to B. But, I want to fix the length, suppose, 10 km. That means, I want to know the route plan to go from point A to point B but in 10 km range. Is it possible to get this kind of query using pgRouting? Or I should use any other tools? If it is possible how should I write the query?

Note: in my database table 'route', i have 'gid', 'source', 'target', 'cost' etc columns.

Thanks, in advance.


texts: How to create an attractive choropleth map in R?

texts:
I've been experimenting with choropleth mapping techniques in R, having reaped the benefits of ggplot2 for creating beautiful graphs within a powerful data analysis package.

When it comes to mapping I've failed to produce results comparable to those I routinely get from QGIS. Using data of flow in Sheffield (fully replicable data and .qgs file here), QGIS easily produced the following:

The best I could produce using R (using code described here) were these two maps, the first using GISTools, the second using ggplot2:

To me, the ggplot2 option seems far more attractive, if only I could solve the problem of the faulty lines (probably a problem with the fortify() command, or not reading in shapefiles using readOGR() described here.)

So the question is 2-fold: is the ggplot2 option the best choropleth mapping solution in R and, if so, how can I solve the problem of the faulty white lines?

Replicable code to find out what I've done is here.

Edit - Since made choropleth() output option more attractive:

The ugliness of this option can be reduced by exporting at higher resolution and removing the legend (other R legends are available). Still not managed to change line thickness but it's getting better! Red lines represent flows to zones which employ > 5000 people.


texts: Developing QGIS Python plugins using Jupyter?

texts: Is it possible to run PYQGIS commands using Jupyter notebook and, if so, is there any documentation to develop QGIS plugins using Jupyter?


texts: WAB's AGOL search bar dropdown list No Results error

texts:
After following all the steps in Web Map > Settings > Application Settings > By Layer > Save.... the search bar's dropdown list returns the error message No Results
Have attempted multiple combinations is the Application Settings and it is still not working...
I found a thread on the topic, but still looking for a solution.
Is this still a known issue?
We are using a Feature Service (not a Map Service) due to selection tolerance issues if that matters.


texts: Getting multiple shortest paths with PgRouting in one query?

texts:
I want to run shortest path algorithm on multiple source&target pairs at once and get a result as a table and process it then.

How do I do this? This query does not work:

SELECT a.source, a.target, paths.*
FROM all_to_all a, shortest_path_astar('...', a.source, a.target, false, false) paths;

ERROR:  function expression in FROM cannot refer to other relations of same query level

(btw, all_to_all does not literally mean all to all, :) it's just a number of random pairs)

This does not work either:

SELECT *
FROM all_to_all a, (
   SELECT * FROM shortest_path_astar('...', a.source, a.target, false, false) yyy
) AS t2;


texts: Problem with LoadLibrary and encrypted DLL

texts: I'm using QGIS and I try to implement a new function into my plugin. For that I need to use a DLL which is encrypted. In my C++ test application I don't have any problems to use that DLL with LoadLibrary. In my plugin I don't have any problems too until my dongle is connected. If no dongle is connected an error message is displayed because DLL can't be decrypted. This is correct too, but LoadLibrary gets no return and hangs. I don't have this behavior in my C++ test application. QGIS does not react too. Is this a general problem of Python? Is there a possibility to fix this behavior?


texts: Using select by attributes then calculate field on selection in ArcPy?

texts:
ArcGIS Pro 2.3

I want to create a standalone script which runs outside of ArcGIS Pro (for automation purposes), which selects features (with blank attributes) in a feature class and runs calculate field on those selections.

The code below doesn't work. It updates all fields, not the selected ones.

import arcpy
arcpy.management.SelectLayerByAttribute(r"C:\Automation\Categories\Coal\Coal_automation.gdb\CLAIMS", "NEW_SELECTION", "Claim_Spend_in_£ = ''", None)

arcpy.management.CalculateField(r"C:\Automation\Categories\Coal\Coal_automation.gdb\CLAIMS", "Claim_Spend_in_£", '"Unrecorded"', "PYTHON3", None)
print ("Done")

How can I write the code to run calculate field on a selection?


texts: How to access the terminal in Grass 7.0.4?

texts:
Just installed Grass 7.0.4 via "
Download GRASS GIS for MS-Windows: Stand Alone Installer (64 bit)".

It's been many many years since I last used Grass, I recall on my old computer there a start menu icon and it would start in terminal/shell/command mode (whatever). With this latest install, there appears to be no option to start grass in command mode? There doesn't appear to be an option to type commands in the GUI, either.

The GUI is handy, if all you want to do is...make a simple map...but I'm trying to run some addons (successfully installed), that can only be run via command. Not an experienced Grass user here. What research I have done on this issue resulted in mind-numbing, pages long, multiple link discussion threads that went way over my head.

All I want to do is access the Grass Command Prompt: How do I enable this feature in Grass 7.0.4 for Windows?


texts: How can I generate a series of maps from time series data in QGIS3?

texts:
I've got a polygon layer containing some time series data in cells and I'm trying to create a series of maps from it.

Data structure: every polygon has a timestamp and a value that I'm displaying on a map. There are usually multiple polygons "stacked" in one place with different timestamps.

What I'd like to do is to create one map for each unique timestamp using the same print layout, extent, layer styling, etc.

I tried using the Atlas feature of the layout manager by using the filter function with the expression "date" = attribute(@atlas_feature, 'date') and sort by "date", but that didn't work.
The resulting pages all look exactly the same when I do this, but I don't understand why.

The only way I found that does exactly what I need is to use the filter function from the layer context menu in the layer overview. In there, I can use "date" = '2012-01-01T00:00:00.000' for example and then export the resulting map. I can't do that for all my layers by hand, though.

Any ideas how I could automate this approach?


texts: How to bring feature selection into Printing Composer?

texts:
I'd like to make maps of Finland's municipalities. Our National Land Survey provides data I can make into a Shapefile, and when a municipality is selected in QGIS/Identify Features, three different levels of associated administrative boundaries are helpfully highlighted, even the exclaves and enclaves of the selected municipality.

My question then is, how can I export this exact map view? When I open a new Print Composer, the selection/color scheme is gone. The exclaves are an important feature I'd like to keep.

You can probably see that there are hundreds of those tiny exclaves, whose names aren't easily readable in the GML, so I really wouldn't bother checking up manually which of them belong where, and then color them in Inkscape.

I tried exporting with the plugin SimpleSVG, but it just compresses the shapes into uglyness. Compare original and after:


texts:
I want to use QGIS to extract lat/lng data from the vertices of polygons.  There are about 400 polygons, each with maybe 200 vertices.  I am having trouble with step one, getting the data into QGIS.

The polygons are viewable in ArcGIS Explorer using a zip file that contains files with the following extenstions:

.shp, .dbf, .sbn, .sbx, and .prj

I can load the .shp file into QGIS, but without the other files there is nothing to view.  How do I load the .dbf, .sbn, .sbx and .prj files into QGIS?


texts: Can you store a Classifier in Google Earth Engine for later use?

texts:
Training a classifier in Google Earth Engine is slow and can sometimes show temporary errors like "Capacity exceeded". To combat this, I would like to train the classifier, store it, and then load it and use it in other projects.
I haven't seen an API call to do this, but maybe I'm missing something.
(Or maybe it's not possible now, but one day this question will get a positive answer)


texts: How to fix distortions at the edges of a UTM zone?

texts:
I have a grid of ~48000 points, however this grid has this break which coincides with two different UTM zones (36S and 37S):

I don't know how this grid was obtained but it's pretty obvious to me that the problem is related with the junction zone of the two UTMs (and eventually was originally obtained in one UTM projection since the points of one side seems to be well-positioned - see next figure).

Is there any way I can fix this (I'm more into R but any suggestion would be great)?


texts:
I am attempting to georeference this parcel map from a local assessor's office.
In doing so, I am using the georeferencer in QGIS 3.4.
I am using the DMS coordinates on the image to reference the corners of the parcels. The format that QGIS requires is as follows: (dd mm ss.ss).
I think the corners on the parcel map are as follows:
Top left corner:  N0° 13' 08"E, N89° 52' 06"E
Top right corner:  S0° 15' 08"W, N89° 52' 06"E
Bottom left corner: N0° 13' 08"E, S89° 52' 52"
Bottom right corner: S0° 15' 08"W, S89° 52' 52"
I reference the parcel as follows:
Top left corner: X/East = 00 13 08.00, Y/North = 89 52 06.00
Top right corner: X/East = 00 15 08.00, Y/North = 89 52 06.00
Bottom left corner: X/East = 00 13 08.00, Y/North = 89 52 52.00
Bottom right corner: X/East = 00 15 08.00, Y/North = 89 52 52.00

I set the target SRS to EPSG: 4326

When I finish georeferencing the image, I expect it to be in Cochise County, AZ, but it is way off.
I asked the assessor what the image's CRS is and he said, "The mapping system used by the county assessor is based on the U. S. Public Land Survey System (PLSS). This system uses a section, township and range which is tied to the Gila & Salt River Base & Meridian (GSBM)."
How can I georeference this image? Where am I going wrong, or what more information do I need?


texts: Is there a plugin or a software where I can visualize my spatial queries from PostGIS? I remember I have used something similar to this in the past for the Oracle database but I can not recall the name of it.


texts:
We know what SRID 4326 means and what 3857 means but I do not know what these numbers mean. What is the minimum and maximum number for these SRIDs?

I know it is a standard but how did they select 3857 for Mercator, or 4326, or 900913...?

In fact I think this is not an ordinary list of numbers, it could be some kind of binary number or things like that. I have searched but could not find any answers.


texts: Road data missing in shapefile downloaded from OpenStreetMap

texts:
I am trying to prepare a base map of road network of Cooch Behar District of West Bengal, India. I have been trying to download the OpenStreetMap (OSM) for the district but found that many roads are missing in the downloaded shapefile. I have also tried to download the same area from geofabrik and the issue remains the same. I saw when I downloaded the same from a small area all roads came but I need the complete road network shape file for the entire district.
I have decided to purchase the road network shapefile for the area but the cost is huge.


texts:
I'm following Building a Python Plugin (QGIS3) tutorial to set up the development environment for building QGIS3 python plugins in Win10.
I have a problem compiling "pyrcc5".
As explained in the tutorial I wrote the compile.bat:

@echo off
call "C:\Program Files\QGIS 3.4\bin\o4w_env.bat"
call "C:\Program Files\QGIS 3.4\bin\qt5_env.bat"
call "C:\Program Files\QGIS 3.4\bin\py3_env.bat"

@echo on
pyrcc5.bat -o resources.py resources.qrce

I copied compile.bat in the plugin folder and run from cmd but I had the follow error:

"C:\PROGRA~1\QGIS" is not recognized as an internal or external command,  an executable program or batch file.

I also follow this suggestion and substituted the content of compile file with this:

@ECHO OFF

set OSGEO4W_ROOT="C:\Program Files\QGIS 3.4"

set PATH=%OSGEO4W_ROOT%\bin;%PATH%
set PATH=%PATH%;%OSGEO4W_ROOT%\apps\qgis\bin

@echo off
call "%OSGEO4W_ROOT%\bin\o4w_env.bat"
call "%OSGEO4W_ROOT%\bin\qt5_env.bat"
call "%OSGEO4W_ROOT%\bin\py3_env.bat"
@echo off
path %OSGEO4W_ROOT%\apps\qgis-dev\bin;%OSGEO4W_ROOT%\apps\grass\grass-7.2.2\lib;%OSGEO4W_ROOT%\apps\grass\grass-7.2.2\bin;%PATH%

cd /d %~dp0

@ECHO ON
::Ui Compilation
call pyuic5 save_attribute_dialog.ui --from-imports -o save_attribute_dialog.py

::Resources
call pyrcc5 ui\resources.qrc -o gui\generated\resources_rc.py

@ECHO OFF
GOTO END

:ERROR
   echo "Failed!"
   set ERRORLEVEL=%ERRORLEVEL%
   pause

:END
@ECHO ON

But I have the same error. How do I solve it?


texts: Filling Field2 to icorresponding values in Field1 of shapefile using Field Calculator of ArcGIS for Desktop?

texts:
I have to give explanation in a new field of the existing field which contains abbreviation.

Example: Field "Geology" has values Q, Qs,N....(abbreviation), I have to give explanation of those in a new field "Geologic Age", where Q is Quaternary Sediments and Qs is Quaternary sand & dunes.. and it need to be done for more than 100 entries, so manual is practically time taking,

something like "IF[GLG]=Q;THEN[Geologic_Age]=Quaternary sediments"


texts:
I used "Prepare Network" function before conducting "Integral Analysis". However, the analysis still failed because of the tolerance problem. What's wrong?

The error message:
sDNA processing

sDNA is running in 64-bit mode

Building network and checking for tolerance errors...

ERROR: Edges closer than Arc's XYTolerance/Cluster Tolerance detected
these are not necessarily viewable in ArcMap
and will lead to incorrect results

Please inspect and fix using sDNA Prepare

TIME 0.0 for start_gs=;end_gs=;metric=ANGULAR;custommetric=;disable=;intermediates=;oneway=;;origweight=;destweight=;weight_type=Link;radii=400;radmetric=EUCLIDEAN;;nojunctions;nohull;bidir;arcxytol=0
.001000;arcztol=0.001000;

Calculation failed.

External command completed

Process did not complete successfully

Failed to execute (sDNAIntegral).


texts:
Is it possible to work on QGIS with QTCreator since most of the build documentation is assuming Linux environment and I have to use windows 7?

all of the old posts in this forum didn't solve this issue it is still not clear how can I integrate Qgis with QTcreator.


texts:
I am trying to retrieve all points p that are within polygons c. While building the query, I came across this strange behavior.

Question: Why does query 2 takes so long compared to query 1? Is it also doing a full tablescan of table points even though the WHERE clause is using c.osm_id and not table p? If so how can we prevent this? My final query probably requires both tables in the FROM clause.

Query 1

SELECT *
FROM counties As c
WHERE c.osm_id IN (-1842139,-1933719,-2362362,-1933718,-2362263,-2306361,-1933745,-1942017,-1931997,-1865772,-1933055,-1933746,-2315704)

Query 2

SELECT *
FROM counties As c, points As p
WHERE c.osm_id IN (-1842139,-1933719,-2362362,-1933718,-2362263,-2306361,-1933745,-1942017,-1931997,-1865772,-1933055,-1933746,-2315704)

This is my final query, it seems to take forever (207 sec) to execute... I'm wondering if the slowness of this query is related to the slow speed of the above query 2.

Table p countains 1 million points and table c contains 300 polygons.

SELECT *
FROM counties As c, points As p
WHERE c.osm_id IN (-1842139,-1933719,-2362362,-1933718,-2362263,-2306361,-1933745,-1942017,-1931997,-1865772,-1933055,-1933746,-2315704)
AND ST_Within(
        ST_Transform(ST_SetSRID(ST_Point(p.lng, p.lat), 4326), 2163)
    ,c.geom)

Update

Query using indexed geometry column p.geom and JOIN still takes 144 seconds. Am I expecting an unreasonable increase in speed?

SELECT *
FROM counties As c
JOIN points As p
ON ST_Within(p.geom,c.geom)
WHERE c.osm_id IN (-1842139,-1933719,-2362362,-1933718,-2362263,-2306361,-1933745,-1942017,-1931997,-1865772,-1933055,-1933746,-2315704)

I figured the long query time is due to transferring the geometry data for 17k rows from the Pgsql database to my computer running pgAdmin3. Transferring the ids alone took 480ms.


texts:
By following the ZOO Client Documentation, I have set up the zoo-client-demo folder to be accessible to localhost. I have also set up the assets folder and its subfolders to contain the required js and mustache files.

The problem is, I also need an HTML file to be present somewhere so that the client actually displays something when the appropriate URL is entered into the browser.

The only bit of HTML code in the documentation is the following:

  Before using the ZOO-Client, you will first have to include the javascript files from your web page. With the use of requirejs you will need only one line in your HTML page to include everything at once. This line will look like the following:

<script data-main="assets/js/first" src="assets/js/lib/require.js"></script>

1) Can someone please explain what this line means?
2) Is the above line enough, or do I need some more code (both the HTML and the .js files)?


texts: Transpose (copy) Google Earth geometry to different location?

texts: I just tried to copy a folder of concentric radius circles to a different location (different centerpoint) in Google Earth Pro & quickly discovered this does not work. Is there anyway to transpose (I guess) a folder of coordinates/shapes easily? ...


texts: What tools in QGIS can be used for calculating habitat structural diversity?

texts: I have a canopy height model derived from a LiDAR DTM and DSM and I would like to calculate or illustrate habitat structural diversity from the canopy height model file. I have tried using roughness but would like to explore other tools. Are there other tools available in QGIS/GRASS etc to do this?


texts:
I have a line feature class that I am trying to symbolize with a Marker Line Symbol type.  I have my symbology set up but I am having a problem with the rotation angle.

When I enter an angle in the box and click ok, the symbol does not rotate.  If I go back into the Symbol Selector dialog box, the angle is still shown as 0.  In this example, I am trying to enter 180 in the Angle box to make the arrow turn 180 degrees.  Is this a bug, or is there a workaround?


texts: Filling Polygon with Lines in ArcGIS Desktop?

texts:
i am trying to fill a polygon with lines and calculate total distance of these lines.

like this example. but i have different distances between these lines like 1 meter, 2 meter, 3 meter.. so i need to take this as an array from client side, calculate length for each of these parameters and return a json response, should i use arcgis modal as a geoprocessing service? or Server object extensions?

Actually i want to do this operation
Filling a polygon with lines using Postgis
in ArcGIS but not sure how


texts: PostGIS Nearest Neighbor query performance - can it be quicker?

texts:
I want to canvas the community to see if my nearest neighbor query can be optimized at all.  I am happy with my current results, but I'm always up for improving.  I have 2 tables, lakes and towns.  I want to get the nearest town to each lake within 0.5 degrees.  I use a calculation to get a cardinal direction output as an integer (ceil function).  Below is the query.  Any other suggestions are also welcome.

SELECT DISTINCT ON (lakes.wkb_geometry) lakes.wkb_geometry, lakes.ogc_fid AS lakes_id, lakes.name AS lakes_name,
towns.ogc_fid AS towns_id, upper(towns.name) AS towns_name,
ceil(degrees(ST_Azimuth(lakes.wkb_geometry, towns.wkb_geometry))/90) AS card_dir
FROM lakes, towns
WHERE ST_DWithin(lakes.wkb_geometry, towns.wkb_geometry, 0.5)
ORDER BY lakes.wkb_geometry, ST_Distance(lakes.wkb_geometry, towns.wkb_geometry);


texts:
I am trying to optimize my python code using pypy. While I am able to import most libraries and run succesfully, I am not able to import GDAL libraries. The same libraries when imported in Python 2.7 run fine.

I have tried adding GDAL location to PYTHONPATH, and even tried downloading GDAL Python bindings to PyPy directory, somehow pypy does not want to recognize GDAL, and comes with the following error.

Traceback (most recent call last):

  File "app_main.py", line 51, in run_toplevel

  File "ImportENC-P.py", line 50, in <module>

    from osgeo import gdal, ogr, osr

  File "B:\PYPY\site-packages\osgeo\__init__.py", line 21, in <module>

    _gdal = swig_import_helper()

  File "B:\PYPY\site-packages\osgeo\__init__.py", line 13, in swig_import_helper

    import _gdal

ImportError: No module named _gdal


texts: Geographic coordinates in PostGIS to QGIS and shapefile layers

texts:
I've got in PostgreSQL/PostGIS the table "points" with latitude and longitude fields. I've done

update  points
set geom = st_point (longitude, latitude);

and also

select updategeometrysrid ('points', 'geom', 3395);

In QGIS I've added that table with the Data Source Manager to a project in EPSG::3395 with a shapefile layer already loaded. However, points from PostGIS table do not appear correctly projected. Any suggestions to fix it?


texts: When I create a map book, using data driven pages, when I place teh actual pages side by side there is some overlap, so the pages do not line up exactly. Is there a way to clip the overlap so that the pages are fairly seamless in the map book?


texts: Display points labels starting at certain zoom level

texts: I use WMS to display my points on the map. I would like to display points labels(name_ID) starting at certain zoom level. From What I have read I understand I can do it with GeoServer rules, but I need a clear example to clarify this.


texts:
I'm fairly new to QGIS and was hoping someone could help me load some MrSID photos into QGIS.

I have read some posts and have installed the OSGEO4W Installer and selected the GDAL-MRSID Lib.

I also have the GDAL plugin installed.

However, after it installed I rebooted but have the same result.

I know the .sid file is good because I can load it into MapWindow without any issues.


texts: GRASS GIS 7, Windows 8.1, Graphical Modeler, Export to Python

texts: I need help. I´m trying to export my model made in Graphical Modeler in Grass GIS 7.0.0 beta3 to Python Script. The exported file contains only header.


texts: I need to know the exact date of image acquisition for NDVI purpose, the date on image name is " LC80920852014222LGN00" is in Julian date, but in meta data it says"acquisition_date>2014-08-10" , I assume that I can use the date in meta data, any ideas?


texts: How to count polygons in multipart features?

texts:
I'm trying to count the number of polygons in a multipart feature. My aim is to have a 'count' attribute to quantify fragmentation of multipart land parcels.

I've looked at the field calculator and for suitable plugins, but can't see anything that would help. There are lots of ways to count points, lines and vertices, but not polygons.


texts:
I'm trying to figure out a way to draw archeological sections in QGIS (example). We currently only draw the layers and features digitally (like this). When we make a section of features or a profile we draw a line from point a to b with the gps and make an analog drawing of the actual section or profile.

The idea is to draw section in the field on a rugged tablet within QGIS and link them to a section line (point a to b).

How do I set this up?


texts: Highlighting point using ArcGIS API for JavaScript?

texts:
How can I highlight a point using ArcGIS API for JavaScript?

When I use the function of centerAndZoom, I also want to highlight this point.


texts: How to map Traffic Count Points to Streets?

texts:
I've got point data representing traffic count locations and I am trying to figure out the best way to represent them on a map.I thought about kriging to provide a continuous surface but I don't think that is a good option.

I could do a spatial join to join the points to the streets but I'm not sure how I could represent the actual counts. Any ideas?


texts:
Recently I've started playing with OSM2PO and pgRouting but I've stumbled upon a weird problem which I'm not really sure where it comes from…

I've installed OSM2PO on Ubuntu, with no errors, and I've also imported latest Portuguese data from Geofabrik, also with no errors at all. When Osm2poService started running, I've tried a very simple routing from 40.64013,-8.61054 to 40.64021,-8.61112 and get a “no routing found”-like alert box. If you check http://goo.gl/maps/k6qyr and http://osrm.at/1KY you’ll see what the result should be.

I've also seen this problem while trying routings between other points, I found some routes that are really strange (not very optimized).

I think that for some reason there are some ways missing in my OSM data.

I'm using OSM2PO v4.5.25.


texts: I am trying to print out a network map using black/white hatching patterns, similar to Black/white hatching pattern fills in QGis but for Lines. Is there some trick I can use?


texts: Generate ArcGIS Token without Username/Password Each Time

texts:
I have an ASP.NET C# web application that authenticates a user via their username and password.  Once the user is authenticated, the app generates an ArcGIS token using Esri's Generate Token operation (http://resources.arcgis.com/en/help/rest/apiref/generatetoken.html).  Once the ArcGIS token is received, the web application can access secured map services from the ArcGIS server.

The web application is authenticating against the same LDAP directory as the ArcGIS server.

My problem is that if I need to generate a new ArcGIS token from the web application, I cannot do so without forcing the user to login again because I don't have the password (and I don't want to store the password).  If the ArcGIS token expires, then the web application loses access to the secured map service.  If in the future, I implement an authentication system that doesn't require a username or password (Oauth via Facebook for example), then I would never be able to generate the ArcGIS token.

Note: I have to use ArcGIS Server - I cannot use Portal or ArcGIS Online.

Is there another way to generate the ArcGIS token?  Is there a way for me to specify, upfront, what the token value will be?  Has anyone experience this before and come up with a solution?


texts: Shapefile extent units in decimal degrees but values in meters

texts:
I have a shapefile that had its geographic coordinate system defined as the projection NAD 83 UTM Zone 18. The scale of the layer was correct but it was at the wrong point in space. I redefined the layer's coordinate system as GCS_North_American_1983 using Define Projection tool. When I did this, the layer's extent values stayed in meters but the units changed to decimal degrees, resulting in the scale to change to 1:167,078,124,575 in order to bring the layer into view.

Data Type:  Shapefile Feature Class
Shapefile:  F:\ADKCouncil\RawData\Natural\Hydrography\NoSource\ADKlakes.shp
Geometry Type:  Polygon
Coordinates have Z values:  No
Coordinates have measures:  No

Geographic Coordinate System:   GCS_North_American_1983
Datum:  D_North_American_1983
Prime Meridian:     Greenwich
Angular Unit:   Degree

The layer has the correct values for the projected coordinate system in meters however it is in the geographic coordinate system using decimal degrees.

Extent
Top: 4969303.000856 dd
Left: 474011.435231 dd
Right: 4766549.500856 dd
Bottom: 637414.935231 dd

I tried to project the layer into NAD 83 UTM Zone 18 using the Project tool however an error results:

  invalid extent for output coordinate system


texts: ArcGIS API for JavaScript QueryTask for user's own map

texts:
I am using ArcGIS online to create a web map and then I want to query this web map using ArcGIS API for JavaScript. In the tutorials I read I see

var queryTask = new QueryTask("http://sampleserver1.arcgisonline.com/ArcGIS/rest/services/Demographics/ESRI_Census_USA/MapServer/5");

But I want to use my own web map which is "http://pawebgate.maps.arcgis.com/home/webmap/viewer.html?webmap=a0c647159b46456b963d8b180e062b57"

I have tried using the whole http://... it did not work.
I have tried using only a0c647159b46456b963d8b180e062b57 it did not work.

Has anyone got any idea about how to use my web map instead of the http://samplserver1 ...


texts: Export data from image collection Google Earth Engine

texts:
I have a function to get an average value per month for a band from an image collection. Full code is here.

var dat = ee.ImageCollection('NOAA/CDR/SST_WHOI/V2').select('sea_surface_temperature')
    .filterDate(ee.Date('1990-01-01'), ee.Date('2000-01-01'))
    .filterBounds(ROI);

// month data for dataset

var startyear = 1990;
var endyear = 1999;

var startmonth = 1;
var endmonth = 12;

var startdate = ee.Date.fromYMD(startyear, startmonth, 1);
var enddate = ee.Date.fromYMD(endyear , endmonth, 1);

var years = ee.List.sequence(startyear, endyear);

var months = ee.List.sequence(1, 12);
var monthlySST =  ee.ImageCollection.fromImages(
  years.map(function (y) {
   return months.map(function(m) {
      var w = dat.filter(ee.Filter.calendarRange(y, y, 'year'))
                .filter(ee.Filter.calendarRange(m, m, 'month'))
                .mean();
      return w.set('year', y)
          .set('month', m)
          .set('system:time_start', ee.Date.fromYMD(y, m, 1));

    });
  }).flatten()
);

In order to extract the data I've been plotting a graph then exporting as a csv.

print(ui.Chart.image.series(monthlySST, ROI, ee.Reducer.mean(), 4000)
     .setOptions({
     title: " Sea surface temperature mean by month of the year ",
     hAxis: {title: 'Time'},
     vAxis: {title: ' Sea surface temperature (°C)'},
     colors: ['#EF851C'],
     pointSize: 3}));

However, I'm trying to export a very large time period. And the graph times out or exceeds user memory. I tried exporting the image collection as a table but it doesn't include the the mean monthly band value in the table. Is there anyway of doing this?


texts: Error printing pdf in OpenGeo Suite client SDK app

texts:
I am using OpenGeo Suite 3.0.2. I created an app using Client SDK and stored the raster and vector data on postgreSQL/PostGIS.

The print tool is not working, giving error the following error:

  Java.io.IOException: the byte array is not recognized imageformat.

These are the dependencies, I included:

@require plugins/Print.js
@require GeoExt.ux/PrintPreview.js
@require GeoExt/widgets/PrintMapPanel.js
@require GeoExt/plugins/PrintProviderField.js
@require GeoExt/plugins/PrintPageField.js
@require OpenLayers/Control/ScaleLine.js

This is the code I included for print service:

{
   ptype: "gxp_print",
   customParams: {outputFilename: 'DAA Map-PDF'},
   printService: "/geoserver/pdf/",
   actionTarget: "map.tbar",
   showButtonText: true
}


texts: Difference in the elevations for the same building on two DSMs saved in *.asc and *.laz formats

texts:
The following question concern the same type of data as the one I've asked in Differences in the longitude and latitude between a LiDAR DSM and an optical image.

Basically, I've downloaded a tile of the DSM for the UK from UK Survey Open Data. The DSM is available under two extensions: *.laz or *.asc files.

The thing I cannot understand is the difference in the elevation values between these two types of data. Below are the two DSM files together with the elevation of the rectangular building seen in white. Note that the *.asc file was opened with QGIS and saved in GeoTiff format while the *.laz files was opened using ENVI LiDAR and processed to produce the DSM from it. The elevation of the building in each image was retrieved roughly from its pixel intensities in the produced GeoTIFF images.

*.asc 2014 DSM, 100cm resolution (file: LIDAR-DSM-1M-2014-SO83/dsm_d0177253_20141212_20141212_mm_units.asc)

White building elevation: around 25m.

*.laz 2014 point cloud, 100cm resolution (file: LIDAR-LAZ-2014-SO83/SO8832_P_9938_20141212_20141212.laz)

White building elevation: around 75m.

Update1:

The *.asc files was downloaded directly from UK Survey Open Data, and used as is with no further processing on my part.
The only processing I carried out on the *.laz files was to open it with ENVI LiDAR and to produce a DSM in the GeoTIFF format (menu Process > Process Data..., tick only the dsm checkbox) from the given *.laz point cloud (The Export Coordinate System was changed to Geographic Lat/Lon).
To follow the terminology in DSM vs DTM, both *.asc and *.laz correspond to a DSM and not a DTM as in both cases the buildings were not flatten. a DTM is available on UK Survey Open Data where you can clearly see that the location where the building is supposed to be is almost the same elevation as the ground near it.
The tile in the screenshots is SO83. It's in the city of Tewkesbury (Gloucestershire, England).
The (latitude, longitude) of the building are (51.999409, -2.162770), in case one needs to visualize it on Google Earth.

Update2:

In order to investigate what had caused the differences in the heights between the two files (.las and .asc), I tried normalizing the .las file (the uncompressed version of the .laz file) using a software called Fusion, as described in Normalizing point cloud data:

Since the DTM is available in .asc format, it needs to be converted first to .dtm format so that it could be given as an input to the second Fusion command. Below is the command I typed on the prompt (The *.asc is the input file, and the *.dtm is the output file):

ASCII2DTM <my-folder>\lidar.dtm m m 1 10 0 0 <my-folder>\dtm_f0177253_20141212_20141212_mm_units.asc

The previous command worked well, but not the one that actually normalized the .las file (i.e. flatten the ground):

ClipData /height /dtm:<my-folder>\lidar.dtm <my-folder>\SO8832_P_9938_20141212_20141212.las <my-folder>\normalized-lidar.las

Update3:

This update concerns a difference in the georeferencing between both .asc and .laz files.

If I convert both .asc and .laz files to the GeoTIFF format as explained above, and I open them with the Sentinel Application Platform (SNAP) to put a Ground Control Point in one of them, it doesn't match the same location in the other. Something might be responsible for the differences in the georeferencing!
If I add the minX, minY, maxX, maxY arguments to the previous command clipData by taking these coordinates from the provided .las file (in ENVI LiDAR) I get the following error: No input tiles overlap the sample area: (556547.00,5759837.00) - (558568.00,5761855.00).


texts: How to center the map on a point while their crs differ?

texts:
Just starting with Python, plugins and actions in QGIS.

From an postal code api I get a position in WGS84: lat and lon.

I want my map (in WGS 84 pseudo-Mercator) to center on this point.

Since several hours I'm googling everywhere I have tried everything but to no avail.

My code now is (based on several answers on different forums since I've not seen the straight answer I need):

map_pos = QgsPoint(lon, lat)

map_pos.crs = QgsCoordinateReferenceSystem(4326)

rect = QgsRectangle(map_pos,map_pos)

self.iface.mapCanvas().setExtent(rect)

self.iface.mapCanvas().refresh()

I get this message:

  rect = QgsRectangle(map_pos,map_pos) TypeError: QgsRectangle():
  arguments did not match any overloaded call: overload 1: too many
  arguments overload 2: argument 1 has unexpected type 'QgsPoint'
  overload 3: argument 1 has unexpected type 'QgsPoint' overload 4:
  argument 1 has unexpected type 'QgsPoint' overload 5: argument 1 has
  unexpected type 'QgsPoint'

Question: how to center the map on a point while their crs differ? Should be very simple!


texts: How to make selection tool with 3 variables which depend on each other (ModelBuilder)?

texts:
I try to create model to selection feature by attribute with 3 variables which depend on each other but i don't know how.
for example:
variables: kat_1, kat_2, kat_3
Kat_1{ gastronomy, education, sport, religion}
Kat_2{bar,restaurant, kindergarten, elementary school, middle school, high school, church, cemetery, playground, pool, tennis}
 Kat_3{names}

how to create some filtr which works that when I choose gastronomy from kat_1 then it only allow me to choose bar or restaurant from kat_2?


texts: Extracting multiple values of raster to polygons using ArcGIS for Desktop?

texts:
I was hoping I could get some help on how to extract information from a land use raster using polygons. I have a layer with 157 polygons and a land use raster with multiple land use categories for the whole region. I need to calculate the area for each land use type in each of the polygons (either the area or the number of pixels). I already extracted the land use information using my polygons as a mask, but I have been unable to get the counts/area for each land use category PER polygon.

I have tried all the zonal and extraction tools, but I either get pixel counts of each land use in the whole layer (not counts per polygon) or I get the overall area for each polygon without any detail on the area for each land use category per polygon.

Could someone point me in the right direction please?


texts: How to place line labels underneath point symbols in QGIS 2.16.2?

texts: I would like to make a map, where the street names are labeled but point symbols shouldn't be covered (see attached image). I learned from another post, that labels are always rendered last and that there is no way to change the rendering order. Is this still true?


texts: Accessing feature class from file geodatabase for intersect in arcpy

texts:
I want to intersect one polygon file in a file geodatabase with multiple point files in a folder.

To access the file in the file geodtabase I do this:

arcpy.env.workspace=r'F:\Sheyenne\Sheyenne_lidar\for_histos\lidar.gdb'
shapes=arcpy.ListFeatureClasses()

when I do this:

for shape in shapes:
    print shape

my file is returned, named lidar_final_poly.

To access the files in the folder I want to intersect with I do this:

arcpy.env.workspace=r'C:\path_to_intersects'
intersects=arcpy.ListFeatureClasses(*.shp)

and then to intersect I do this:

for i in intersects:
    arcpy.Intersect_analysis([shapes,i]), outpath)

but this returns:

ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Input Features: Dataset F:\Sheyenne\Pixel_Regression\Viking\GDD_regress\NDII\shapefiles\NDII_ts.shp #;lidar_final_poly # does not exist or is not supported
Failed to execute (Intersect).


texts: QGIS NNJoin - Add Point to Single Nearest Polygon

texts:
I know there have been tons of posts on this, but none of them seem to answer my question/situation. I have a layer of many polygons (building footprints) and I'm trying to join a of points (less points than polygons). I am trying to join the attributes from each point to the closest polygon. When I use NNJoin, I get a resulting layer where each polygon has assumed the attributes of the closest point to it. Similar posts are here and here. In many instances, my points are not within a defined polygon. Maybe I'm missing some settings in NNJoin?

NNJoin interface here; 'Building Outlines' is my polygon layer and 'Building Properties' is my point layer. Do I have the settings correct?

Am I right in understanding that if I check 'Approximate geometries by centroid' it will find distance based on centroid, but if I leave it unchecked it will base it on closest line of the polygon?


texts: pgr_drivingdistance seems to be excluding nodes that should fall within the cost parameters set

texts:
I am getting a result that I don't understand and I am not sure if its my lack of understanding of OSM.

I am using the pgrouting functions pgr_drivingdistance to get the list nodes reachable within a certain time frame. This is my query:

  create temp table prouting as select * from pgr_drivingdistance('select gid as id, source, target, traveltime_min as cost from ways', 3605, 5, false);

Travel time was calculated as follows:

  sudo -u postgres psql -d planet -c "update ways set traveltime_min = length_m / 5000 * 60"

The contents of this table look as follows:

I wanted to now see a list of the nodes that attach to all source nodes of this graph but are not part of the result set of the pgr_drivingdistance function. I then paired these up with the costs of the result set in the prouting table. The idea was that I would now be able to have a cost traversing from the origin point up until the matching node.

select wi.gid,
       name,
       wi.source,
       pr.agg_cost,
       traveltime_min,
       5 - agg_cost as rcost,
       length_m,
       traveltime_min - (5 - agg_cost) as acost
from (
   select gid from ways wi join prouting pr on (wi.source = pr.node)
except
   select edge as gid from prouting) as sexclude
join ways wi on (sexclude.gid = wi.gid)
join prouting pr on (wi.source = pr.node);

rcost is the remaining cost that is left over when subtracting the agg_cost from the prescribed max cost (5) that was fed into pgr_drivingdistance. This gave me back the following result:

As you can see from the highlighted/underlined elements, those nodes should have been included in the original query as their cost when added to the value of the connected node falls within the cost limit I prescribed, 5.

Can anyone allude or give me a hint as to why these nodes are not included in my pgr_drivingdistance query.


texts: Showing only current index polygon extent in Data Driven Pages

texts:
I am currently working on putting together a few map series in Data Driven Pages. I can't figure how to only show the index layer outline for the current page. I would like to hide the outline for any index polygon that is not the current focus of the data frame.

In the example picture below, the first page of my series (index layer is the yellow polygon labelled 1), I would like to hide the index polygon for page 2 (the yellow box in the upper right of the image). As I move on to page two, I would like to hide the index polygon for pages 1 and 3, as all of index polygons adjoin and are sequential.

Is there something simple I am missing or a way to do this programmatically?


texts: How to calculate weighted mean coordinates with negative weights in QGIS?

texts:
I'm pretty new to prof. GIS, but have OSM experience. So if I just miss something obvious, please don't punish me :)

I have a .CSV file with measurements of wifis along a track. I like to guess the wifi accesspoint positions by using the signal levels during the tracking.
Adding the locations works fine (coloured circle markers), but applying the mean coordinates feature to the point layer results only in the usual centeroid (blue star) of the convex shape polygon(soft yellow) and seems to take not the signal level into account:

But I expect that the new position should be usual more close to the dark blue circles in the south.

So my CSV structure looks like

trackpoint_id, lat, lon, ..., signal
1, 54... , 12... ,..., -90.0
2  54... , 12... ,..., -74.0

The signal col is recognized as double. As signal is in dbm so negative (smallest negative number is best) and logarithmic scale, I might need to create a new column containing a cleared signal_strength?
Sadly the log doesn't tell me if the mean coordinates tool has problems with computation :(


texts:
I am trying to work in QGIS and I am very inexperienced with it. I recently was able to import several shapefiles into a SpatiaLite database and I need metadata to go with these shapefiles. 4 out of 5 my shapefiles were created in QGIS and I put metadata in under their Properties, but I do not see this metadata when I open them from SpatiaLite. The other shapefile was download from a county GIS data portal but had no metadata entered in.

I have already tried the Metatools plugin, but the buttons where I can try and edit the metadata appear to be grayed out and I'm not sure how to fix it and use it.

How can I put in metadata for the shapefiles in the SpatiaLite database in QGIS?


texts: Do ArcGIS Online feature services support a DISTINCT or returnOnlyUniqueValues option?

texts:
ArcGIS Server REST API has:

returnDistinctValues
Description: If true, then returns distinct values based on the fields specified in outFields. This parameter applies only if the supportsAdvancedQueries property of the layer is true.
Values: <true | false>
Example: returnDistinctValues=true

http://resources.arcgis.com/en/help/arcgis-rest-api/index.html#//02r3000000r1000000

Is there an equivalent when querying agol feature service?


texts: GeoDjango and PostGIS distance query efficiency concerns

texts:
For example, if I want a query that says "show me all houses within 50km from here and order them from nearest to furthest", these are the questions and concerns I have:

Assuming we are using EC2 (application server with GeoDjango) and RDS (database server with PostGIS), it seems that this query would execute on RDS. This is a problem, isn't it? Not only is RDS optimized for IOPS, but we only have one RDS instance serving all our EC2 instances. It would only be worth it if this distance query is only slightly more expensive than simply reading all the latitudes and longitudes.
Will GeoDjango or PostGIS be smart enough to run precise distance calculations only on houses that are within a reasonable distance, skipping all the houses that are too far away to matter? For example, it can calculate the maximum bounds given a radius, and only compute distances for houses within those bounds.
Does there exist a less precise method to calculate distance in GeoDjango/PostGIS, so we can sacrifice accuracy for speed? If not, would it be significantly better to roll our own solution?


texts:
I want to georefrence streets in Autocad.
I used _Align command for a single object, but it is not working with multiple objects.
How can we align multiple objects together in autocad?


texts: Inserting accumulative attribute in dynamic text element based on Date field using ArcGIS for Desktop?

texts:
I've got a line feature class that contains Date field. I've made a layout view in order to display time animation of this feature class. I've managed to enable the time option in the layer and to create a Data Frame Time text in correspondence to the elapsed time.

I'm looking to enrich the layout by creating a dynamic text element of the feature class'  accumulated 'Shape_Length' field, according to the elapsed date\time - time elapsing > accumulated 'Shape_Length' gets bigger.

I've tried several ways to do it, but didn't succeed so far:

I've tried inspecting the Data Frame Time element Text syntax --- <dyn type="dataFrame" name="Layers" property="time" emptyStr="[off]"/> --- but couldn't figure how to change it to show the 'Shape_Length' field.

I've try an alternative approach: to create a date based graph. I've followed this link, and search quite along time in the Advanced Properties of the graph how to make it a time enabled and also how to make it accumulative. Didn't work as well.

I think the answer would be in connection to the LayerTime function.


texts: Multiply Rasters with Scalar Numbers as Input in QGIS Modeler

texts:
I built a little model in the QGIS modelbuilder which uses the raster calculator with the command

  gt(a,threshold)

I need to run the model with varying thresholds from 2 to 30 so it would be great to create a Number input in the model for this threshold so I'm asked for its value when starting the model, but the raster calculator only accepts rasters as inputs.

Maybe someone has an idea how to get around this problem? For example by generating a raster containing only the threshold value in every cell in a previous step if there is an easy way to do this?


texts: How to fix buffer selection issue in ESRI ArcGIS Javascript API and LGIM applications

texts:
I had strange LGIM application behavior where the map spatial search did not highlighted the parcels within the Geometry Buffer. What could be the possible reasons?

LGIM stands for the Local Government Information Model. It is a geodatabase schema ESRI uses to support the ArcGIS JavaScript applications they have.


texts:
I already asked in previous question QGIS table join issue, how to join shape file data and CSV data which contain multiple SAME ID for one object?

Shape file contains columns: object id, area, height etc

CSV contains columns: object id, materials as brick, metal, concrete etc

So first I was joining them by JOIN option, but the numbers did not match as the software only recognises one value of the field and overwrites it to all other same IDs.

I was offered solution Project>Properties>Relations and I used instructions on: https://docs.qgis.org/testing/en/docs/user_manual/working_with_vector/attribute_table.html#id40. But I have over 84000 rows and don't have time to go through all of them and make connections. Further this option does not seem to have all the information I need for example in the attribute table. I would ideally like to be able to do this one-to-many join such that I could have overview of each material (concrete, steel, brick etc) for each field (I did not figured that out with relation).

Join table would be ideal but as stated the numbers do not match so it is issue with representation. Also it is important to have all materials for each field as they will be intersected with different layers and the correct amounts are needed for further calculations.


texts: How to create a reference grid in QGIS print composer?

texts: Is there a way to create what is called a 'reference grid' in ArcMap using QGIS print composer:


texts: Using Python codeblocks in ArcGIS to evaluate if all values in many columns is 0

texts:
I'm working on zoning data and I've been evaluating if the percentage for various zoning is more than 90% of the area, so I can define it as that zone in a new field, like so:

def ifBlock(RS, RM):
if RS >= 90:
    return "RS"
elif RM >= 90:
    return "RM"
else:
    return "Mixed"

There are approximately ~40-50 lines of the above code for the different zones. However, for some of the parcels, the values for all of the zones come out to 0. I want to make a new field to identify these areas as "no data". I thought  that I could define all of the zoning columns in a list, and I wanted to write in Python something like "if all of the values in the list = 0, return no data".


texts: Getting Unique Intersections of Multiple Polygons Efficiently in Python

texts:
I would like to get a list of all polygons that are intersections of other polygons.  However, that list of intersections cannot self-overlap.  For example, if I have a triple Venn diagram, I would want four polygon intersections returned:  One for the center where the three circles intersect; and one each for the intersections of pairs of circles, where those three intersections do not include the center intersection.

Based on the top answer to this question, I have this code:

from shapely.geometry import Polygon
from shapely.geometry import MultiPolygon
from rtree import index

#read in my polygons separately, not shown, store to array of Polygons

intersections = []
idx = index.Index()
for pos, polygon in enumerate(shapely_arr_polygon):
    idx.insert(pos, polygon.bounds)

for polygon in shapely_arr_polygon:
    merged_polygons = cascaded_union([shapely_arr_polygon[pos] for pos in idx.intersection(polygon.bounds) if shapely_arr_polygon[pos] != polygon])
    intersections.append(polygon.intersection(merged_polygons))

intersection = cascaded_union(intersections)

From this code block, I get "intersections" which is a list of Polygons that are all pair-wise intersections (not what I want), and I get "intersection" which is a MultiPolygon of the union of all of those intersections (not what I want).

I'm thinking there should be a way to use either/both of the above to do what I want, or perhaps a different way to get what I want in that loop, but I can't figure it out ... help?

Full ideal solution: I would have a list of polygon intersections that best represent the overlaps of the original polygons, such that in the triple Venn diagram case, 4 overlaps with that center being separate because it's 3 overlaps instead of 2.

Maybe good enough solution:  Unique overlaps that just don't overlap with each other.  So in the triple Venn diagram case, it could be good enough to just return 3 overlaps, where the center is included in ONE of the pair-wise overlaps.  I think that could avoid a cominatorial situation.


texts: Why does dropping a single point dramatically affect location of polygon centroid?

texts:
I'm a relative novice working with some coordinate data using shapely and have found a weird behavior when trying to identify outlier points. The approach I am using is to iterate over a series of points creating a Polygon using all but the current point in the iteration, and then to calculate the distance from the current point to the centroid of the Polygon created in step 2. Any point significantly further from the centroid of the others than the mean of those distances is considered an outlier.
One case has yielded strange results. The points are as follows, all relatively close together:
from shapely.geometry import Polygon
from shapely import wkt

p = ['POINT Z (28.219403 -9.360669 0)',
     'POINT Z (28.3255635 -9.280321799999999 0)',
     'POINT Z (28.219435 -9.360731700000001 0)',
     'POINT Z (28.2077176 -9.442305599999999 0)',
     'POINT Z (28.3244627 -9.286111200000001 0)',
     'POINT Z (28.2191343 -9.360647500000001 0)']

When you compute the centroid of the polygon with all of them included, the answer is very reasonable.
>>> p = [wkt.loads(x) for x in p]
>>> Polygon(p).centroid.xy
(array('d', [28.250548831170587]), array('d', [-9.362973996290707]))

But when you run the same calculation dropping the third point in the list, the resulting centroid is incredibly far away from all of the other points.
>>> p_prime = ['POINT Z (28.219403 -9.360669 0)',
               'POINT Z (28.3255635 -9.280321799999999 0)',
               # 'POINT Z (28.219435 -9.360731700000001 0)',
               'POINT Z (28.2077176 -9.442305599999999 0)',
               'POINT Z (28.3244627 -9.286111200000001 0)',
               'POINT Z (28.2191343 -9.360647500000001 0)']

>>> p_prime = [wkt.loads(x) for x in p_prime]
>>> Polygon(p_prime).centroid.xy
(array('d', [-6.857533437170832]), array('d', [-191.95569600984862]))

If you drop any of the other points, the result is around (28.3, -9.3). Just this single point seems to throw off the calculation.
Why is this happening and how can I avoid it in the future?


texts: Donating and volunteering geographical knowledge?

texts:
I know a lot about the places I visit, and I want to share it with everyone.

I have been donating and volunteering at OpenStreetMap, and I am also aware of Google Maps.

Where else can I help with collaborative mapping?


texts: I have had no luck snapping points from one layer to another line layer in both QGIS 2.18 and 3.6.3. The snapping plugin is enabled and editing is on for both layers but there are no 'Apply' or 'OK' buttons.


texts: ArcGIS Toolbox Python Script Validation: extracting field names from a multi-value list of feature classes

texts:
I am writing a python script within an ArcGIS toolbox (not a Python Toolbox), that takes as input, a multivalue parameter of feature classes from which I would like to populate a second parameter with a drop down list of fields from the first feature class in the multivalue list. When using a single feature class as input, the field parameter can simply be set using the "Obtained from" functionality in the "Parameters" tab, but not so with a multivalue input. I am attempting to perform the operation using the validation classes instead. Below is my attempt using the "updateParameters" class. The validation doesn't fail, but doesn't produce a list of fields as expected. Any thoughts?

def updateParameters(self):
  """Modify the values and properties of parameters before internal
  validation is performed.  This method is called whenever a parameter
  has been changed."""

  if self.params[0].value:
      self.params[1].filter.list = arcpy.ListFields(self.params[0].valueAsText.split(";")[0])
  return

Edit - I’ve also tried:

[f.name for f in arcpy.ListFields(self.params[0].valueAsText.split(";")[0])]

With no success


texts: Avoid calling tm_shape() multiple times?

texts:
Posting here (rather than Github) as it is not an issue with the tmap package.

Situation:

I am doing simple choropleth maps of varinterest for U.S. counties (N~3000) using tmap (data stored in data object). Each map also has a key column that we will call category. I am rendering a variable of about ~300 counties for a varinterest for a category, but I would also like to fill the counties with missing data as gray.

My current method is as follows:

tm_shape(data) +
     tm_fill() +
tm_shape(data %>% filter(category='cat1')) +
     tm_polygons(varinterest) +

However, this is a bit slow (taking 5-7 minutes to render on a very strong computer). Mapping only the variable of interest (i.e. the second tm_shape(... takes only a few seconds.

I feel like I am missing something obvious that might speed up rendering - any advice?

I will work on a reproducible example and post my times.


texts: Working with year before 1900 and ArcGIS Space-Time Cube?

texts:
I am trying to create a space time cube for 10 years from 1870-1880, but I get this Error:

ValueError: year=1873 is before 1900; the datetime strftime() methods require year >= 1900

I tried to work around with .isoformat   (which allows strftime() to handle years greater than 1900), but I got yet another Error.

import arcpy
import datetime
arcpy.env.overwriteOutput = True
arcpy.env.workspace = r"C:\Users\rin\ArcGIS\Default.gdb"

points_input = r'C:\Users\rin\Analysis\space_time_cube\Export_Output.shp'
outpath = r'C:\Users\rinfall\Analysis\space_time_cube\washington_cube.nc'

with arcpy.da.SearchCursor(points_input, ['YEARCODE']) as cursor:
    for x in cursor:
        arcpy.CreateSpaceTimeCube_stpm(points_input, outpath, x[0].date().isoformat, "FISHNET_GRID")

Is there anther way to work around with the year issue?


texts: Using mapbox to implement lighting effects?

texts:
I use mapbox.gl.js + postgrepsql+mapnik+ mapbox vector tile;
this is my map:

but I want this, The point of overlap is brighter:


texts: Checking if python script is run from ArcGIS (arcmap or on server) or a stand alone python script

texts:
For debugging purposes I have created a hardcoded conditional that checks it the script is run as stand alone or as an ArcGIS tool. Is there anyway to fugure this out at run time? e.g. some environment variable, global variable, etc.

e.g.

dryRun= 1
if dryRun:
    par= 'Hello'
else:
    par= arcpy.GetParameterAsText(0)
# Do something


texts: Batch processing/extracting raw data of one raster using one shapefile (with many polygons)?

texts:
If we have a raster, say of integer elevation data for a country,
and one polygon shapefile, say of 300 river basins in that country, with a unique name for each,  how would we most easily get output like this for them all?

basinID, gridcellelev
a, 320
a, 321
a, 320
b, 17
b, 18
b, 19

The most burdensome way seems to be filtering/converting the single shapefile into 300 shapefiles,
clipping the raster 300 times into 300 uniqueID rasters, reading them back in, generating individual tables for each basin, then combining them all together.

On the other hand the ideal way seems to be skipping the file generation, not saving the xy data, and creating the same table from just using one raster and one shapefile - by iteratively selecting the cells within a basin, stamping them with the basinID, creating a table, losing the coordinates, and keep repeating and appending that table until the 300th basin.

I'm not looking for any statistics, just the raw data listing of the grid cell elevations that would have been part of some standard clip method. I believe the attribute table that comes with the clipped raster from ArcMap is the counts/frequency of the cells. That output table format works for me too.

I don't know to minimally reproduce a raster and a polygon shapefile,
so I'd just be grateful for any tips/libraries/functions/examples. If starting with R, here's a starting point:

library(tidyverse)
library(raster)
library(rgdal)
library(sf)

elev_raster <- raster("spain_elev_meters.tif") #integer raster
basins <- readOGR("spainbasins.shp", "spainbasins") %>% st_as_sf() #unique basin ID column: `basinID`

I'd prefer to do everything in R, but happy to try any batch methods in ArcGIS  as well (I have 10.6 with Spatial Analyst, but not Pro).


texts:
I need to download distribute_setup.py from the following website and the website seems to be broken, has anyone got an idea where else I can download this from?

http://python-distribute.org/distribute_setup.py

I am following this:
https://github.com/AstunTechnology/Loader/wiki/Installation

to install a loader to finally convert an Ordnance Survey (OS) topographic layer (which is in .GZ format) into PostGIS format.


texts: When making a point feature a different point feature appears

texts: I am digitizing a map and I have sample locations and station locations. The program (ArcMap10) worked fine until recently. Now when I try to add a sample it adds a station instead. Any ideas on how to fix this is much appreciated1


texts: Added Column via attribute table in QGIS will not show up in style menu

texts: I have my data within a spatialite DB. I added a line feature to the canvas in QGIS, opened up the attribute table and the field calculator and added a new column while calculating something. When I want to style this line feature layer, I can not choose the new column in the styling menu . Why is that?


texts: I am trying to add a building geometry as GeoJSON to Cesium virtual globe. But the 3d building geometry is seen 2D polygon. Here is a picture from GeoJSON file and Cesium. As can be seen GeoJSON file includes z coordinates.


texts: Open Layers get WMS Layer attributes without single click event

texts:
I am using javascript open layers v6.4.3 to load data from postgreSQL database. Here is my code to create feature layer and the click event to get identify feature:
roads3Layer = new ol.layer.Tile({
    source: new ol.source.TileWMS({
        url: "http://localhost:90/cgi-bin/mapserv.exe?map=d:/code/geoportal/geoportalWMS.map&",
        params: { "LAYERS": "Roads3", "TILED": true }
    })
});

navMap.addLayer(roads3Layer);

navMap.on('singleclick', function (evt)
{
    var viewResolution = navMap.getView().getResolution();
    var url = roads3Layer.getSource().getFeatureInfoUrl
    (
        evt.coordinate, viewResolution, 'EPSG:3857',
        {
        'INFO_FORMAT': 'text/json',
        'QUERY_LAYERS': 'Roads3, Regions'
        }
    );
});

This works fine and identifies the feature without any problem using the click event.
My concern is that I need to get the attributes of roads3Layer without the click event. Something like this line:
var att = roads3Layer.getAttributes();

I need this because in my web application, I want to add the select by attribute functionality, where the user selects feature, then attributes loaded, then user selects one attribute and enter the search value.
I tried to find a url in the capabilities in my mapserver, but could not manage to find one.


texts:
I used to be able to right-click on a layer in the Table of Contents, go to Open Attribute Table, and the attribute table would display.  However, now nothing happens.  The problem occurs regardless of what layer I click on, regardless of what map I'm in.  It does not help to restart the application, or to restart my computer.

Any thoughts?


texts: QGIS not reading elevations from GeoTIFF created by ACADCIVIL?

texts:
I've generated a raster DEM in ACADCIVIL in TIFF format: name.tif & name.tif.aux.xml.
ArcGIS reads the data correctly but QGIS sets the whole data to the same value = minimum height
Looks that QGIS 3.4 doesn't support this format.
Same happens on 2.18

gdalinfo output:

Driver: GTiff/GeoTIFF
Files: Surface1.tif
       Surface1.tif.aux.xml
Size is 405, 528
Coordinate System is:
LOCAL_CS["POSGAR_07_F5",
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (5512262.832999999634922,6109047.309000000357628)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=Autodesk Civil3D 2008
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 5512262.833, 6109047.309)
Lower Left  ( 5512262.833, 6108519.309)
Upper Right ( 5512667.833, 6109047.309)
Lower Right ( 5512667.833, 6108519.309)
Center      ( 5512465.333, 6108783.309)
Band 1 Block=405x100 Type=Float32, ColorInterp=Gray
  Min=38.088 Max=39.473
  Minimum=38.088, Maximum=39.473, Mean=38.565, StdDev=0.301
  NoData Value=3.4028234663852886e+38
  Offset: 38.082,   Scale:0
  Metadata:
    STATISTICS_COVARIANCES=9.08270199017587E-02
    STATISTICS_MAXIMUM=39.473079681396
    STATISTICS_MEAN=38.564526756118
    STATISTICS_MINIMUM=38.088008880615
    STATISTICS_STDDEV=0.30137388124225
    STATISTICS_VALID_PERCENT=52.86

I've changed scale to 1 --> same result
File Properties --> Info
Band

Number Band    No-Data    Min   Max
1    Band 1  3.40282e+38  n/a   n/a


texts: Write accent mark in .shp and .dxf filename

texts:
I'm trying to export with python 2.7 some shapefiles and CAD files that could have specials characters in their filename, like accent marks. In example, I want to have this:

MMT_LT_1215_20101004_05_Crespià_Esponellà.shp

But when I export the files the filename is like this:

MMT_LT_1215_20101004_05_CrespiÓ_EsponellÓ.shp

I've already changed the encoding to UTF-8 with # -*- coding: utf-8 -*-, but it don't result anyway. The fact is that when I create the directory that must contain the files and has the same file name, the accent marks are writed correctly like in the first example. I'm not sure if I have to use the encode() or decode() functions to face that.


texts: How to merge two 'incompatible' polygon layers in ArcGIS?

texts:
I have two layers of polygons with administrative boundaries. Although they come from the same data provider, they seem to have small discrepancies and do not 'fit' each other.

How could I merge these two layers in ArcGIS 9.3 removing overlaps and gaps between them?

Would it be possible to prioritize one layer (grey one in my case) and 'fit' the other one (purple) accordingly?


texts: qgis QgsExpression doesn't work corretly

texts:
I have to do a script that adds a field and changes attributes values utilising an expression so i use QgsExpression.
When i put a expression in QgsExpression where there is only 'field name' it works QgsExpression( """ ((Prof_am*Prof_av)/2)*LargTranch*$length """ ), but when i put in the expression a variable  it doesn't work ' (ratio_roche)*((Prof_am*Prof_av)/2)*LargTranch*$length '

below is my script after adding new fields:

Features = layer.getFeatures()
for feat in Features:
    if feat[field_Prof_am_index]<feat[field_ProfRocher_index] and feat[field_Prof_av_index]<feat[field_ProfRocher_index]:
        #terrain 100% meuble
        exp_Deb_meuble = QgsExpression( """ ((Prof_am*Prof_av)/2)*LargTranch*$length """ )
        exp_Deb_meuble.prepare(layer.pendingFields())
        feat[field_Deb_meuble_index] = exp_Deb_meuble.evaluate(feat)
        feat[field_Deb_roche_index] = 0
        layer.updateFeature(feat)
    else:
        #terrain en partie rocheux
        #ratio de la partie rocheuse
        ratio_roche=0.2
        #evaluer quantite de la partie meuble
        exp_Deb_meuble = QgsExpression( ' (ratio_roche)*((Prof_am*Prof_av)/2)*LargTranch*$length ' )
        exp_Deb_meuble.prepare(layer.pendingFields())
        feat[field_Deb_meuble_index] = exp_Deb_meuble.evaluate(feat)
        #evaluer quantite de la partie rocheuse
        exp_Deb_roche = QgsExpression( ' (1-ratio_roche)*((Prof_am*Prof_av)/2)*LargTranch*$length ' )
        exp_Deb_roche.prepare(layer.pendingFields())
        feat[field_Deb_roche_index] = exp_Deb_roche.evaluate(feat)
        layer.updateFeature(feat)


texts: Analysing point attributes within polygons

texts:
I'm using QGIS 3.4.4.
I have a polygon layer and a points layer. The points have the attributes district and streetname.

How can a check if all the points within a polygon have the same district and streetname?


texts: ArcGIS Javascript API 4.2 basemap toggle thumbnail issue

texts:
I am currently porting an older app I created with ArcGIS JavaScript API 3.14 to API version 4.2.  I had a search widget and a basemap toggle widget in the original application.  I have successfully implemented both widgets in the new app but I can't figure out how to set the default basemaps thumbnailUrl property for my streets basemap.  The current thumbnail looks like the following image after I switch to my imagery basemap.

In the original document I created my own custom basemaps using MapServers by doing the following:

esriBasemaps.Street = {baseMapLayers: [{url:"https://gis.linncounty.org/arcgis/rest/services/PublicBase/baseMapLinnCounty/MapServer"}],
                                                thumbnailUrl: "images/Streets.jpg",
                                                title:"Street"};
esriBasemaps.Imagery = {baseMapLayers: [{
                                            url:"https://gis.linncounty.org/arcgis/rest/services/PublicBase/baseImagery2014/MapServer"},{
                                            url:"https://gis.linncounty.org/arcgis/rest/services/PublicBase/baseImageryLabels/MapServer"}],
                                            thumbnailUrl: "images/Imagery.jpg",
                                            title:"Satellite"};

//Default map extent
//var startExtent = new Extent(5372725,3417182,5502485,3578378,new esri.SpatialReference({wkid:102675}))

//Creates the map object
    map = new Map("map",{
    basemap:"Street",
    //extent:startExtent,
    center: new Point({x:5419598,y:3458787,spatialReference: {wkid: 102675}}),
    zoom: 3
    });

//Basemap toggle widget added to the div #BasemapToggle
    var toggle = new BasemapToggle({
        map: map,
        basemap: "Imagery"},
        "BasemapToggle");
    toggle.startup();

In the new app I am using WebMaps for all my layers and basemaps so I don't know if that is throwing me for a loop.

This is what the new code looks like right now:

var webMapID = "861f2064cd374eafa4bf6379bb451571";
var webmap = new WebMap({portalItem: {id: webMapID}});
var view = new MapView({map: webmap,container: "viewDiv",rotationEnabled:false});
var pItem = new PortalItem({id:"cb956d729be045c1a55dd97b45994cd6"})
var Imagery = new Basemap({ title:"Imagery",
                                    id:"myBasemap",
                                    thumbnailUrl: "images/Imagery.jpg",
                                    portalItem:pItem,
                                    });

toggle = new BasemapToggle({
          view: view, // view that provides access to the map's 'topo' basemap
          nextBasemap: Imagery, // allows for toggling to the 'hybrid' basemap
        });

view.ui.add(toggle, {position:"bottom-left",index: 1});

I can change the icon after the map has loaded by doing the following in Chrome's dev console:

toggle.activeBasemap.thumbnailUrl = "images/Streets.jpg"

If I insert that into my code it doesn't work which leads me to believe that it is loaded asynchronously.

How am I able to set the webMaps basemap icon for the toggle widget? I would think it would be a property I could set when I create the new WebMap instance but I am unsure.


texts: Publishing map from ArcGIS Pro to ArcGIS Server as Map Service?

texts:
As a comment on an Esri blog post entitled Announcing the ArcGIS Pro beta program last May someone named mahenderraj asked:

  I presume that it is not possible to share the map as a mapservice to
  ArcGIS server from ArcGIS Pro. For that we have to use the regular
  ArcMap Desktop product.

I do not seem to have seen an answer to this, and in the ArcGIS Pro 1.0 interface it seems easy to publish maps to ArcGIS Online or Portal for ArcGIS but ArcGIS Server gets no mention.

Is it possible to publish a map from ArcGIS Pro to ArcGIS Server as a Map Service?


texts: Plugging in population data into a map of provinces

texts:
Have shapefiles for provinces in several countries in an Asian region, but currently trying to plug in animal population data for each province with graduated symbols displayed the data. Am trying to work off a Microsoft Excel Chart that I have made with the data (basically with one row for provinces and the other row for populations)- however, importing via a delimited text layer isn't working. I'm pretty new to GIS, and am wondering if there as an elementary error I am making - do I need to enter vector data for each province in my excel file before these will display? Is there an easy way to elicit this from the shapefile(s) that I have that display the provinces in the country?

Thank you so much for your time - I'm trying to learn this as well as I can!


texts: How to open popup on mouseover in Leaflet?

texts:
Possible Duplicate:
  How to show a popup on mouse-over, not on click?

I'm working with Leaflet library. I defined L.layerGroup and add some markers to the group. Can somebody help me with mouseover function on markers. I need to get function which will open popup when user mouseover on marker. How to write code for that function considering this coding I wrote? I'm new with Javascript so please consider that.

Thanks!

var cities = new L.layerGroup();

var city1 = new L.marker([35.896667,14.4625], {icon: myIcon} ).addTo(cities);
var city2 = new L.marker([35.882222,14.522222], {icon: myIcon} ).addTo(cities);
var city3 = new L.marker([35.840854,14.484347], {icon: myIcon} ).addTo(cities);


texts: I am using QGIS to make some boundaries adjustments of the Italian provinces over the years. I have a shapefile with the provincial boundaries in 1961 and another one with the boundaries in 2001. No value is missing in these shapefiles. However, when I use "Union", I get some NULL values. Why does this happen? And how can I handle them?


texts: Perform Tile Index on Newly Created PDFs using PyQGIS

texts:
The python code below successfully performs the tile index algorithm on all PDF files stored in sub-folders named !pdf's within my V:Drive. A list of all paths to each PDF exported from QGIS is created first off, then a processing loop is performed against the variable pdf_path. The output is created as a .gpkg.

import os
import pathlib

my_path = "V:/GIS - Files/1. Client Projects/"
pdf_parent_folder = "!pdf's"
pdf_paths = []

for path, sub, files in os.walk(my_path):
    if pdf_parent_folder in sub:
        for path, sub, files in os.walk(os.path.join(path, pdf_parent_folder)):
            for name in files:
                if os.path.splitext(name)[1] == ".pdf":
                    pdf_paths.append(str(pathlib.PurePath(path, name)))

print(pdf_paths)

out_file = 'C:\\Users\\Username\\Desktop\\Example_Folder\\Output_file.gpkg'#
change to your own file path
for pdf_path in pdf_paths:
    processing.run('gdal:tileindex',
                {'ABSOLUTE_PATH': True,
                'LAYERS': pdf_path,
                'OUTPUT': out_file})

As a one-off this is perfect however, at the end of each working week, I would like to add the tile index of the new PDFs I export to the .gpkg output.

Is there a time constraint I can add to prevent the algorithm inputting PDFs I already have tile indexes for?


texts:
I have a table containing 3D geometry in a PostGIS 2.0 database (the geometry is of type LINESTRING Z).

I need to get the geometry out and into a DXF, DGN or 3D Shapefile. What is the easiest way to do this?

(I was sure I had achieved this via QGIS before, but it seems to remove the Z value and produce only a 2D file when attempting a 'Save As ...' from a PostGIS layer).


texts:
I am trying to convert a GIS LINE shapefile into a .osm.pbf file on Windows.

I've done some research on the topic and some tools that kept coming out are:

Merkaator
shp2osm - which deals mostly with point shapefile
shape2osm.py - which keeps bringing up an error with importing * elementtree.SimpleXMLWriter

But nothing seems to be working.


texts: Snapping points to crossroads using ArcGIS Desktop?

texts:
I have some points (370 points) in a layer . These points must snap to crossroad layer . I could'nt find any tools in arcmap to snap points to crossroad layer.

Is there any way to snap these points to the crossroads ?

I'm using ArcGIS Desktop 10.2.


texts: How to edit font and text styles in Cartodb infowindow?

texts:
I would like to change the font and text styles of the text in my Cartodb infowindow. Specifically, I'd like the name at the top to be larger font and bolded to act as an infowindow 'title' and the rest of the text to be a new font.

I haven't found a way of doing this on the online template, how can this be achieved for a beginner like me?


texts: White page after setting up reverse proxy for geoserver: Why?

texts:
I have installed Geoserver 2.5 in Tomcat6 on a Windows server 2008 R2. I want to have internet access to geoserver so I asked my IT Dept to set up a reverse proxy to internal address using their proxy server in the DMZ. When requesting the geoserver URL I just get a blank page with no source code.

https://egov3.heidekreis.de/geoserver

Can anyone suggest why this is not working?


texts: Exponential notation in attribute table QGIS

texts:
I've imported a spreadsheet from excel, all my spreadsheet data is formated as "number" inside it.
When I import it in QGIS, it changes to exponential notation as shown in the image below. Is there a way to change this?


texts:
I'm running Mac OS X with native Python 2.6.1 installed.

Later today installed GDAL complete pkg, and added the program to the shell path on the terminal with the folowing cmd:

~: echo "export PATH=/Library/Frameworks/PROJ.framework/Programs:$PATH" >>.bash_profile

Resulting on my bash_file this following line:

~: echo "export PATH=/Library/Frameworks/PROJ.framework/Programs:$PATH" >>.bash_profile

Now every time I try on terminal to import the pyproj I receive the error:

ImportError: No module named pyproj

I'm new to Mac's/UNIXs and to OSGgeo in general.


texts: Iterating feature classes in multiple geodatabases using ModelBuilder?

texts:
I have 5 geodatabases. All have a feature dataset that contains a point feature class called "Point" and a line feature class called "Polyline." I want to go through these feature datasets, take only the Polyline feature class and append it to a master feature class.

Is there a way to nest iterators?

I tried created submodels, but that didn't seem to work.


texts:
How can I use WFS-T in openlayers-2?

I'm using geoserver in server side.


texts: ArcGIS 10.2.2 Does "Batch Project" work for different input datum?

texts: I have two vector data sets - one is in NAD27 and the other is in WGS84. Normally I reproject one at a time. But can I just use 'batch project' to reproject both of them to NAD83 at once? I was not sure if the input datum must be the same for all files.


texts: Subdividing polygons using lines in QGIS?

texts:
I have a city that I would like to divide up by its streets with the end result being polygons that are the blocks divided by its streets. In the original data that I have, the city boundary should be the outside boundary. It is a polygon. The city should be divided up by its streets. I have those in two files, I could use either - one is the road right of way which leaves space for the streets, the other is the road centerlines. Those are polylines.

Again, what I would like to do is subdivide the city into blocks. I have tried Union and then polygonize, which yielded an error message. When I checked the geometry the entire boundary came up invalid. When I try to convert lines to polygons, almost the entire city disappears leaving only a couple of blocks.

It's become rather circular at this point. I try one way, run into some problems, and someone tells me to try the other. I try the other way, someone tells me to try the first way. and on and on.

Here is an image of the boundary with the road centerlines.

And here is a conceptual image (not from GIS, only for illustration). In the concept image, the thick blue line would be the polygon corporate boundary and the thin orange lines the road centerlines/rights of way.

What I would like to do is to merge the geometries so that each subdivided segment is its own polygon. In other words, in the end there would be five polygons from the conceptual image.

How can I do this?

Just to be clear I want the result to be polygons not polylines.


texts: Most efficient way to split a polygon with lines, C++ API

texts:
I'm looking for the most efficient way to split a polygon (possibly with holes) with a set of lines (which can intersect) in C++ with GDAL/OGR or any other library which can prove useful.

As shown in the figure below, I would like to split the polygon with the red lines, which would result in multiple polygons (possibly with holes):

Performance is a big issue as I need to make such calculation thousands of times.
For now, I'm doing the following:

Build a MULTILINESTRING with every lines
Buffer it (distance 0.001, precision is not a big issue)
Compute the difference with the POLYGON

Profiling shows that Difference() is a bottleneck.

Some solutions rely on PostGIS (e.g. ST_Split) or other command line/query language tools with no direct binding in the C++ API of GDAL/OGR.


texts: Basic components for building an offline app

texts:
I am very new to GIS app creation and am basically learning from scratch. I am not building the app but researching what would be the best components. These are the essential requirements:

5km (area of interest) ability to zoom out and see whole earth
Refresh and dynamically update overlays
Interactive click to get more info
Display lat/lon on cursor
Switch between road/aerial/topo basemaps
Ability to host all mapping data on an isolated server that isn’t on the internet. The clients are on the same isolated network.
Service that provides standard aerial photographs.
Ability to add our own custom aerial photographs to the server as a new layer
No subscription fees

After some research I think for the basemaps the best options I could find was OpenStreetMaps for the roads layer, USGS Imagery for the aerial layer, and USGS topography map for the topography layer.

OpenLayers as the front end.

GeoWebCache as the tile cacher.

GeoServer as the tile server.

Am I missing anything or is there a better option available? From what I've read all these components can work together and there are methods to get OpenStreetMaps working with GeoServer. I know there will be a lot more to it but I am just trying to get the basics down so I can dive deeper into how they work.


texts: Use the Python port of ogr2ogr in memory

texts:
I found a Python port of ogr2ogr which can be found here
This port allows me to re-project vectors from a Python script instead of having to use the command line:
import pyogr2ogr, tempfile

with fiona.open("B:/Canada-d8-100m/Canada-d8-100m.geojson", "r") as before:
    print(before.crs)
> {'init': 'epsg:3347'}

pyogr2ogr.main(["","-t_srs", "EPSG:4326", "-f", "GeoJSON", "B:/proj_test.geojson", "B:/Canada-d8-100m/Canada-d8-100m.geojson"])

with fiona.open("B:/proj_test.geojson", "r") as after:
    print(after.crs)
> {'init': 'epsg:4326'}

I would like to do the same thing as above, but I want to be able to read/write the output GeoJSON to/from memory instead of to/from the disk.
How can I do this?


texts: Error in passing 'system:index' values from L8_TOA collection to filter L8_SR collection

texts:
I'm new to Google Earth Engine and just starting to learn the tricks. I'm trying to pass 'system:index' values I retrieved from filtering Landsat 8 TOA collection to extract the same images from Landsat 8 Surface Reflectance collection. My code below:
var filtSysInd = filtDateCol.aggregate_array('system:index'); //filtDateCol is the filtered Landsat 8 TOA collection
//print('filtSysInd',filtSysInd.get(1));

/*filtSysInd values to replicate the issue
List (7 elements)
0: LC08_075084_20140223
1: LC08_075084_20150125
2: LC08_075084_20151227
3: LC08_075084_20170303
4: LC08_075084_20190104
5: LC08_075084_20200224
6: LC08_075084_20210125*/

var fImage = L8sr.filterBounds(AOI); //L8sr is the filtered Landsat 8 Surface reflectance collection

var counter = ee.List.sequence(0,filtSysInd.length().subtract(1))
//print('counter',counter);
var filtL8sr = ee.ImageCollection(counter.map(function(i){
    var si = filtSysInd.get(i)

    //return ee.Image(fImage.filterMetadata('system:index', 'equals',si))
    return ee.Image(fImage.filter(ee.Filter.eq('system:index', si)))
}))
print('filtL8sr=',filtL8sr)

The error I get is below:
filtL8sr=
ImageCollection (Error)
ImageCollection.fromImages, argument 'images': Invalid type. Expected type: List<Image<unknown bands>>. Actual type: List<ImageCollection>. Actual value: [<ImageCollection>, <ImageCollection>, <ImageCollection>, <ImageCollection>, <ImageCollection>, <ImageCollection>, <ImageCollection>]

Could someone help me understand what I'm missing?


texts: Adding topography to base maps using Surfer?

texts:
I want to add the topography map with extension (BLN) to my the map which I create it using the Surfer prog, I put only the country board.

How can I add the topography maps as a background using the surfer program?


texts: Merging a vector layer with a Raster layer

texts:
I have a Tiff file that has an aerial view and also a shapefile which has contours of the same area.  How can I put the contours in the tiff file so that you can view the TIFF and see them both in one picture.

I am using QGIS.


texts: Converting ArcPy fieldmap to pandas dataframe

texts:
I want to get a "snapshot" of a fieldMapping / FieldMap.
How do I convert it to a pandas dataframe?
My current approach is to
fieldList = arcpy.ListFields(theFeatureClass)
extracted = []
attributes = []

for field in fieldList:
  attributes = [field.name, field.aliasName, field.type,  field.length, field.defaultValue, field.required]
  extracted.append(attributes)
print (pandas.DataFrame(extracted))

But how do I do this with FieldMappings?


texts: Using Python Script tool in ModelBuilder to check if field exists, if not, create one?

texts:
I'm trying to create a script tool to use in ModelBuilder that will check the incoming dataset for the existence of a particular field. If the field exists it should simply go on to the next calculate field tool which will use that field, however if the field doesn't exist it should create it before running the next calculate field tool so that I don't get an error from the Calculate field tool, because the field it's trying to use doesn't exist.

I added this script to the ModelBuilder, but it keeps giving me the same error as when I try to run the model and the field doesn't exist. So the script tool is apparently not doing its desired job. Can anyone see an issue with my code? Is there a better way of trying to do this

# Import modules
import arcpy
import sys
import traceback

    # Set local variables
    indata = arcpy.GetParameterAsText(0)
    fieldName = "PrintKey_GIS"

    try:
       # check if the PrintKey field exists
       if arcpy.ListFields(indata, fieldName): #if field exists, evaluates to true
           arcpy.AddMessage("PrintKey_GIS field exists")

       else:
           # Add a new field of that name
           arcpy.AddField_management(indata, fieldName, "TEXT", "", "",25, "PrintKey_GIS" , "NULLABLE", "", "",)
           arcpy.AddMessage("PrintKey_GIS field does not exist")

    except Exception as e:
        arcpy.AddMessage(arcpy.GetMessages())


texts: I have a ArcGIS webmap ( http://lvrs.maps.arcgis.com/apps/Viewer/index.html?appid=a2756011923849659bc76ecd5ef7201c) and I am wondering how would I export this and actually be able to run this on google maps as a layer - so i can run it as an app - instead through the website. Also, the end game (this is for a fire department) then I would like to use it for navigation to an address. Any ideas would be helpful.


texts:
I am seeking a new data collection software for a team of biologists who track migrating animals from planes through aerial telemetry.  I currently use Delorme TOPO USA but when collecting point data the time stamp field is maintained as "date modified" causing errors the the analysis of the data.  The typical data collection softwares I have tried in all of the similar questions do not work because they fail to meet one or more of our criteria that must be met. The current preference is for software compatible with Windows 7 based laptops. The criteria are as follows:

Must be able to track current GPS Location with a refresh rate of 3 or less seconds.  (We are in planes and moving quickly).
Must be able to manually place waypoint on map.  We are not collecting current locations.  We are collecting locations of animals based on predicted position of an animal relative to current position.
Must be able to obtain/import offline maps at several zoom levels for an area that covers most of the midwestern US that  are able to refresh at a rate that keeps up with the planes movements.  Working online is not an option.
Must have a topo option for maps that roads/towns/other reference points that can be overlay on map.  USGS does not typically offer enough detail for our purposes.
Must maintain the seconds field in the time stamp for speed calculations.
Must be able to edit/delete points easily with as few operations as possible.  Erroneous points happen and must be able to remove and move on quickly to determine next location within a minute or so.
Must be able to quickly identify a name for asset as there time when there are multiple animals of the same species being tracked and they need to be differentiated.
Editing a point cannot edit time fields, or a "date created" field needs be maintained.
It would be nice if it had some navigation capabilities because we have a ground crew tracking at the same time and its beneficial to have every team on same software.


texts: Extremely slow Nominatim install with Docker

texts:
I am trying to install Nominatim based off of a dockerized version of the tool found here: https://github.com/mediagis/nominatim-docker. I was successful with installing dockerised Nominatim with a smaller country locally. And the process ran pretty fast.

Then, I got an EC2 instance with, 64GB memory, 1TB SSD, and 6 core processors to do an install with the entire planet (48GB). The most of the planet file got processed fast, in about 6hours and when processing the relations gave an error and failed. When I ran the planet file based on the same docker image, the object parsing (specially ways) were extremely slow. About 100 times slower than the initial run. This is the error I got when the planet ran faster:

Processing: Node(5640116k 412.0k/s) Way(625947k 118.89k/s) Relation(628790 712.11/s)2019-12-18 22:48:23.270 UTC [210] LOG:  incomplete message from client
2019-12-18 22:48:23.270 UTC [210] CONTEXT:  COPY place, line 233940695
2019-12-18 22:48:23.270 UTC [210] STATEMENT:  COPY place FROM STDIN
2019-12-18 22:48:23.270 UTC [207] LOG:  incomplete message from client
2019-12-18 22:48:23.270 UTC [207] CONTEXT:  COPY planet_osm_rels, line 598784
2019-12-18 22:48:23.270 UTC [207] STATEMENT:  COPY planet_osm_rels FROM STDIN
Killed
2019-12-18 22:48:23.275 UTC [210] ERROR:  unexpected EOF on client connection with an open transaction
2019-12-18 22:48:23.275 UTC [210] CONTEXT:  COPY place, line 233940695
2019-12-18 22:48:23.275 UTC [210] STATEMENT:  COPY place FROM STDIN
2019-12-18 22:48:23.275 UTC [207] ERROR:  unexpected EOF on client connection with an open transaction
2019-12-18 22:48:23.275 UTC [207] CONTEXT:  COPY planet_osm_rels, line 598784
2019-12-18 22:48:23.275 UTC [207] STATEMENT:  COPY planet_osm_rels FROM STDIN
2019-12-18 22:48:23.275 UTC [207] LOG:  could not send data to client: Broken pipe
2019-12-18 22:48:23.275 UTC [207] STATEMENT:  COPY planet_osm_rels FROM STDIN
2019-12-18 22:48:23.275 UTC [210] LOG:  could not send data to client: Broken pipe
2019-12-18 22:48:23.275 UTC [210] STATEMENT:  COPY place FROM STDIN
2019-12-18 22:48:23.275 UTC [207] FATAL:  terminating connection because protocol synchronization was lost
2019-12-18 22:48:23.325 UTC [210] FATAL:  terminating connection because protocol synchronization was lost
ERROR: No Data
string(7) "No Data"

Finally, I just took a pbf of North America (8GB) and tried the install. Object parsing related to North America went relatively fast, but it is still working on rankings and indexes. The install is going on for over 16 hours.

I didn't do any PostgreSQL optimizations before running the install. When the object parsing is slow, I see the container is only using about 1% of the available memory.

In Nominatim install instructions (here: https://nominatim.org/release-docs/latest/admin/Installation/) I see that with less config (32 GB RAM) a planet install is taking only 2 days. What am I doing wrong and how can I optimize the process?


texts: Geoserver WMS raster too dark when zoomed out

texts:
I'm publishing raster data (3 band geotiff) with Geoserver (for the first time). When i look at the unpublished raster data in QGIS im seeing:

But when i look at the published WMS it appear very dark:

When i zoom in close it lightens a bit.

What do i need to do to make the WMS look more like when i watch the raster in QGIS?

gdal info output:

Driver: GTiff/GeoTIFF
Files: /home/bera/Geodata/Sentinel-2/20150819/2.tif
       /home/bera/Geodata/Sentinel-2/20150819/2.tif.aux.xml
Size is 10980, 10980
Coordinate System is:
PROJCS["WGS 84 / UTM zone 32N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32632"]]
Origin = (600000.000000000000000,6800040.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  600000.000, 6800040.000) ( 10d52' 4.82"E, 61d19'17.77"N)
Lower Left  (  600000.000, 6690240.000) ( 10d48'41.04"E, 60d20'10.72"N)
Upper Right (  709800.000, 6800040.000) ( 12d54'53.65"E, 61d16'40.73"N)
Lower Right (  709800.000, 6690240.000) ( 12d47'47.56"E, 60d17'39.86"N)
Center      (  654900.000, 6745140.000) ( 11d50'51.77"E, 60d48'40.97"N)
Band 1 Block=10980x1 Type=UInt16, ColorInterp=Gray
  Min=125.000 Max=6507.000
  Minimum=125.000, Maximum=6507.000, Mean=2222.441, StdDev=735.814
  NoData Value=0
  Metadata:
    STATISTICS_MAXIMUM=6507
    STATISTICS_MEAN=2222,4410683233
    STATISTICS_MINIMUM=125
    STATISTICS_STDDEV=735,81396086526
Band 2 Block=10980x1 Type=UInt16, ColorInterp=Undefined
  Min=184.000 Max=3718.000
  Minimum=184.000, Maximum=3718.000, Mean=455.960, StdDev=161.904
  NoData Value=0
  Metadata:
    STATISTICS_MAXIMUM=3718
    STATISTICS_MEAN=455,95982458067
    STATISTICS_MINIMUM=184
    STATISTICS_STDDEV=161,90389249049
Band 3 Block=10980x1 Type=UInt16, ColorInterp=Undefined
  Min=397.000 Max=2867.000
  Minimum=397.000, Maximum=2867.000, Mean=705.941, StdDev=123.796
  NoData Value=0
  Metadata:
    STATISTICS_MAXIMUM=2867
    STATISTICS_MEAN=705,94109519
    STATISTICS_MINIMUM=397
    STATISTICS_STDDEV=123,79629499805


texts: pyQGIS - set grid annotation direction 'Vertical Descending'

texts:
I am using pyqgis to create map automatically. However I encounter a problem when I want to set the grid annotation direction to 'Vertical Descending'.
The command I use is the following:

composerMap.setGridAnnotationPosition(QgsComposerMap.OutsideMapFrame, \
                                                  QgsComposerMap.Top)

composerMap.setGridAnnotationDirection(QgsComposerMap.VerticalDescending, \
                                               QgsComposerMap.Top)

However, I get the following attribute error:

  QgsComposerMap has no attributes VerticalDescending

After looking online, I found that this method is deprecated (https://qgis.org/api/deprecated.html) and that I should use grid() or grids() instead. When I look at qgscomposermapgrid.h source code (https://qgis.org/api/qgscomposermapgrid_8h_source.html#l00211) I see that there is an attribute 'VerticalDescending' (line 215) but I cannot find any example on how to use this method.

Any help?


texts: How long should Delete Identical take on 3.25 GB feature class of polygons with spatial index?

texts:
I have a 3.25 GB feature class of polygons with a spatial index. I am using the delete identical tool comparing the shape and one other field which contains a text string of 9 characters.

The tool has been running for 13 days and I suspect it has crashed/got stuck. ArcCatalog says it is still running as does task manager. As there is no progress bar is there a way to estimate how long this task should take?

I am using 64 bit geoprocessing from ArcGIS 10.1 for Desktop with 16GB RAM 3.4 GHz 6 core processor


texts:
I know how to convert a csv file with xyz coords to point25d with ogr2ogr and vrt file
but the coords in the files are the vetices of 3dpolygon (firsth and last ones are the same coords).
I want to convert that files to POLYGON25D with ogr2ogr
How I can do it?

the file has

x,y,z

219980.345,1834052.671,741.279


texts: PostGIS 2.0 and Self Intersecting Features

texts: I have a polygon feature class created by OpenStreetMap that is returning a self intersection error.  Is there a way to filter out these features? Can I just use a SQL statement to find the error polygons?


texts: Viewing feature service view in Collector?

texts:
I have created a hosted feature service view in arcgis online and I have added this to a webmap. Editing is enabled on the hosted feature service view.

Should I be able to edit and view this on collector?