texts: Custom Feature Renderer in ArcGIS Server 10.1?

texts:
I've written a custom feature renderer.  It works great in ArcMap (10.1).  Now we want to deploy it to ArcGIS Server (10.1) so we can use it in map services.  But it doesn't seem to work.  Has anyone else successfully deployed a custom feature renderer to ArcGIS Server?

More details:

First of all, this is all running on Windows Server 2008 R2, 64 bit.  And it's ArcGIS Server for Java (and Desktop) 10.1 with SP1 installed.

To build the feature renderer, I followed the steps listed here:

http://resources.arcgis.com/en/help/arcobjects-java/concepts/engine/index.html#/Explore/0001000006ws000000/

And like I said, it works great in ArcMap.

The docs seem to indicate that you can deploy it to ArcGIS Server.  I followed the steps at the following link, namely to build a ".jar" file and copy it to ArcGIS/Server/usr/lib/ext:

http://resources.arcgis.com/en/help/arcobjects-java/concepts/engine/index.html#/Deploying_extensions/000100000671000000/

But when I try to deploy a service that uses the custom symbology, I get an error message (from the "Analyze" step in service deployment) that has these details:

Severity: High
Status: Unresolved
Code: 00006
Description: Layer's symbology is not supported.

I tried copying my jar to a few different places, just to see if it'd pick it up, but to no avail.  (Using ProcessExplorer, I looked at the command line of the various ArcGIS Server processes and saw what entries it was using for the Java classpath and then put my jar in those spots.)

Also, the link above says to deploy the jar at <ArcGIS Server Home>/user/lib/ext.  This path doesn't exist on a stock ArcGIS Server
10.1 install.  But <ArcGIS Server Home>/usr/lib/ext (i.e. without the "e" in "user") does exist.  Just to be sure, I tried both places (by manually creating the directory with the "e").

So I wonder: what additional steps are required to get ArcGIS Server to recognize our custom feature renderer?

Or maybe their documentation is misleading and they don't really support custom feature renderers in ArcGIS Server.


texts: Labeling expressions on AGOL with Arcade for certain features

texts:
I am working on a web map in AGOL for the indoor office map. I have some issues writing custom expressions with Arcade in order to label the types of rooms (office rooms) based on another attribute name.

I have 2 attributes:

attribute "Room_type", which has at least 10 different values, and
attribute "Name", which has the labels I need for different types of rooms.

So what I want to do is to label only those room types that correspond to "Meeting" type, and labels should come from attribute "Name", because it contains the titles of each "Meeting" room type.

I've been trying something like this, but it only returns Null values for me:

a) var room_type = $feature["Room_type"]

var meetlabel = $feature.Name

`If (room_type == 'Meeting') {
return meetlabel;
} else {
return 'None';
}`

b)
iif ($feature["Room_type"] == 'Meeting', $feature.Name, 'None')

c)

`If ($feature["Room_type"] == 'Meeting') {
return $feature.Name;
} else {
return 'None';
}`

d) When( $feature["Room_type"] == 'Meeting', $feature.Name, 'None')

e) Decode( $feature["Room_type"], 'Meeting', $feature.Name, 'None')

I've started thinking to just calculate a new field for each building floor that has the required labels in place with Field Calculator in ArcMap.
But hopefully it's possible to do with Arcade, as it would definitely shorten my workload!


texts: Grouping roads based on existing condition and AADT? (QGIS)

texts:
For QGIS software, I have municipality service agreement roads and an end-plan to group the roads based on close proximity (roads and communities close together or close to each other and on the same side of the main road), with similar need for re-pavement or maintenance - i.e. similar current conditions.... The road shapefile has attributes for roadside development score, AADT (Annual Average Daily Traffic) score, pavement priority points (scores), existing road surface, and existing road condition. I would like to find a way to group the roads based on priority for development/maintenance, such as roads in the worst condition that have the highest traffic and are currently most developed, to roads that need least amounts of development/maintenance. I believe I am able to use the pavement priority score as my main symbology type since this ties in both the roadside development score and AADT, OR it is possible for me to symbolize based on existing condition.

Any ideas of how I could go about this using QGIS?

I know that a buffer would be a good geoprocessing tool to use and will likely be in my end-product... However, I have limited experience in this area, and I am having a hard time finding a good starting point.


texts: Land classification confusion matrix pixel total and toList() pixel total not matching

texts:
Using Landsat imagery and a CART classifier I have successfully carried out land classification and computed a confusion matrix for one full Landsat image, please see: https://code.earthengine.google.com/357ac2012f2275f2aedff34f76a81874. Also here is the confusion matrix:

And according to the tutorials I have read the total of all the boxes is the total number of pixels and therefore 95,591. However in the same piece of code I have also calculated the pixel total of a much smaller classified image/polygon using:
var pixel_List = classified_mosaic.reduceRegion({
    reducer: ee.Reducer.toList(),
    geometry: Tokyo_Pre,
    scale: 30,
});
var class_pixel_array = pixel_List.toArray();
var transposed_array = class_pixel_array.transpose().toList();
print(transposed_array,'prefecture image')

When I print 'transposed_array' the pixel total I get is 2,466,989 which is larger than the confusion matrix's total of 95,591, even though the confusion matrix's image is much larger. Surely the 2,466,989 total is correct, therefore the confusion matrix total is incorrect? If so how do you calculate pixel total from a confusion matrix? Maybe I am missing something here, any help would be hugely appreciated.


texts: Generating Attachment Match Table in Enterprise Geodatabase?

texts:
I've created a file geodatabase, imported a dataset, enabled attachments, generated an attachment match table and then added attachments. Everything worked perfectly and quickly.

When I tried the exact same thing in my enterprise geodatabase, ArcGIS 10.1 just takes forever and I can't seem to create the match table.

Is there anyway to copy the dataset with all the attachments from the file geodatabase to the enterprise geodatabase?


texts: How to run qgis scripts in parallel across remote server(s)

texts:
I have a stand alone python script using qgis that I wrote with the help of the SE question here

The script takes three arguments being; an input file path, an output file path, and another input file path. a la:

myscript.py input1.shp output.shp input2.shp

As per advice from this question the script now happily runs in parallel using gnu parallel and gives the correct output files with the following command:

parallel myscript.py {.}.shp {.}_output.shp /tmp/input.shp  :::  /tmp/tile_{1..443}.shp

I want to do the work in parallel on a remote server (or servers). So I enabled a passwordless login to an available server as per this question, I installed QGIS on the remote server, copied over the script, and made it executable from bash. If I manually copy the required input files over to the server the same parallel command as above runs fine.

The problem I have is when trying to run the remote parallel work from the local machine.

There is an example for running a remote process described here

One could imagine that the equivilent gis type line as per my example might look as follows:

parallel --trc {.}_output.shp -bf /tmp/input.shp -S 192.168.88.222,: \
'myscript.py {.}.shp {.}_output.shp /tmp/input.shp'  :::  /tmp/tile_{1..443}.shp

Unfortunately, there's a fundamental problem with this approach. If a GIS program calls the path to a "shapefile" as an argument; while /path/to/my/file.shp will work most of the time, the "shapefile" actually consists of a collection of files with a common filename prefix. When run, the program actually locates whichever file it needs from the collection from the path to the ".shp" file.

This issue then becomes a problem when the gis script is operating in the remote environment since it doesn't have access to the full collection of files.

The question is then how do we solve this?


texts: Using Google Maps API to selectively highlight railway routes one at a time?

texts:
I am trying to understand the TransitLayer in Google Maps' Javascript API. Utilizing Google's example code, I managed to display a geographically-accurate map of Tokyo Metro, Toei Subway and Yokohama Municipal Subway routes here: https://jsfiddle.net/dsd00wb0/

But what I want to do is highlight specific lines one at a time. For instance, only highlight the Asakusa Line. I can do that in Google Maps, but what I want to know is if I can specifically do it using the Google Maps API.

Moreover, how would you highlight heavy rail/overland lines like the Yamanote Line? The TransitLayer only works on subway lines.

I have been able to extract route relations for specific railway lines in OpenStreetMap (using JOSM), but Google Maps is a lot more confusing.


texts: Fill nodata holes in discrete integer raster with majority value of adjacent pixels QGIS

texts: I have a large land cover raster with classified integer values from 2-22. There are some single or double pixel no-data gaps where I merged 2 adjacent tifs together. If this were a continuous raster I would use the GDAL fill nodata tool to remove these holes but because this is a discrete raster I need to preserve the values I would like a tool that uses the majority value of the pixels adjacent to the nodata pixel instead of the interpolated value.


texts: Why won't MapInfo Pro read my tab files from QGIS 14.4?

texts:
For a whole host of reasons, if I want to load data into SQL Server, I need to convert the shapefiles to tab files, open them in MapInfo Pro and use the Easy Loader tool to create tables on my database server.

When I tried to do that today, my tab files won't open in MapInfo Pro. Version 12.0.2 just crashes, and version 12.5 reports that it cannot open the file.

I have used the Universal Translator on my shapefile and that has created a tab file that MapInfo Pro can read.

I don't know when my work flow stopped working. I haven't had to upload a table for quite a while so I can't remember if it was before or after my latest QGIS version.

Any ideas of what may be happening?


texts: Mosaic overlapping rasters in ArcGIS Desktop?

texts:
I am trying to mosaic 3 data sets, 2 of which are bathymetry and 1 elevation data sets. The bathy data overlap the elevation data where it has NoData and the elevation data overlap the bathy data where there is NoData.

My problem is I keep getting NoData from 1 or the other dataset obscuring the actual data. This normally occurs after I  Build Overviews. Creating footprints also does not work, ArcCatalog just tends to crash.

I have already created mosaics of the individual data sets and then try to mosaic using these. I have simplified my elevation data set and removed al the little islands and small NoData areas within the data set.

See attached picture.


texts: i have successfully installed the ECW support on ubuntu qgis. the main challenge that i have is that whenever i try and load an ecw file, the qgis main window crashes. On the other had i am able to succesfully load all my vector layers with ease. any suggestions?


texts:
I want to count trees in SAGA GIS and I am trying to follow the steps described in Counting trees using segmentation and vectorization in SAGA Gis from dominoc925.blogspot.com (05th, March, 2012).

The problem is while I was able to convert the las file into an .asc DEM, I could not retrieve the .asc DSM from the las file. Hence, I can't accomplish step 1 in the workflow above, because I need to import the DSM in SAGA.

How can I convert a las file into DSM?


texts: Strange behavior with esri.tasks.locator

texts:
I'm facing a very strange issue with respect to the esri.tasks.locator function. If I do a hard-refresh (Ctr+Shift+R), the map starts working without any error. But after that if just try to refresh again (normal F5), the console prints out

  TypeError: esri.tasks.Locator is not a constructor

Then if I again do a hard-refresh, everything works normally. What could be the cause of this behavior? The code I'm using is as follows:

dojo.require("esri.map");
dojo.require("esri.tasks.locator");
.
.
.

var locator = new esri.tasks.Locator("https://geocode.arcgis.com/arcgis/rest/services/World/GeocodeServer"); // this line is throwing error if not hard-refreshed
var params = {address:"200 E Berry St", maxLocations:1};
locator.addressToLocations(params);
locator.on("address-to-locations-complete", function (evt) { ... } );


texts: How to change encoding options for a shp file in ArcGIS

texts: After creating a shp file with attributes containing Turkish characters like ş,ç,ü.., Upon exporting it as another shp file , shp file having problem on displaying those characters, Is there any way in ArcGIS to change encoding option like in QGIS. I can change system encoding from UTF-8 to windows 1254 in QGIS and I have no problem on displays. But how do I do it in ArcGIS???


texts: How to create a csv with lat long coordinates from a WGS84 layer with Grass

texts:
This question seems trivial but I'm having trouble doing this simple task.
I'm working on a WGS84 location in Grass 6.4 and I'm trying to convert a raster layer in a csv with lat long coordinates using r.out.xyz but all I can get is a list of longitudes that range from 180.25 to 539.75 and latitudes from -55.75 to 83.75.
I need them to range from -180 to 180 and from -90 to 90.
r.info output:

  Type of Map:  raster               Number of Categories: 255
  | |   Data Type:    DCELL
  | |   Rows:         294
  | |   Columns:      720
  | |   Total Cells:  211680
  | |        Projection: Latitude-Longitude
  | |            N:     89:06N    S:     58:06S   Res: 0:30:02.44898
  | |            E:       180E    W:       180W   Res:  0:30

g.region -p

  projection: 3 (Latitude-Longitude)
  zone:       0
  datum:      wgs84
  ellipsoid:  wgs84
  north:      89:30N
  south:      58:30S
  west:       180E
  east:       180E
  nsres:      0:30
  ewres:      0:30


texts: Reprojecting Natural Earth "Admin 0 – Countries" shapefile to Web Mercator

texts:
I haven't received a reply posting this same question on the Natural Earth forums, so I figured I'd give it a try here.

I’m trying to reproject the 1:10m "Admin 0 – Countries" shapefile into a Web Mercator projection (EPSG:3857) using QGIS. However, the reprojection fails with the following message:

Export to vector file failed. Error: Failed to transform a point while
drawing a feature of type ‘ne_10m_admin_0_countries’. Writing stopped.
(Exception: forward transform of (3.141593, -1.570796) PROJ.4:
+proj=longlat +datum=WGS84 +no_defs +to +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs
Error: tolerance condition error)

I receive this error for both "ne_10m_admin_0_countries" and "ne_10m_admin_0_countries_lakes".

A post here on gis.stackexchange.com indicates that it might be a problem with points in Antarctica going to infinity during reprojection, but I’m not sure how to fix it if that is indeed the problem.

The 1:50m "Admin 0 – Countries" shapefiles are able to be reprojected successfully using QGIS, which leads me to believe that the root of the problem may lie in the 1:10m shapefiles.

I'm hoping that someone can expand on the answer offered in the linked post and provide answers for the following questions:

Is there an error in the 1:10m shapefile? (Since the 1:50m shapefile reprojects without issue) Or is this sort of reprojection error just to be expected under certain circumstances?
Exactly what steps may I take to fix the shapefile so that it can be successfully reprojected? The linked post says to "delete the objects around the south pole" but I'd prefer not to delete Antarctica. I tried editing the Antarctica feature to chop off the southernmost portion, but it didn't seem resolve the reprojection error.

Please assume I'm a GIS newbie.

Thanks!

RESOLUTION & NEW SHAPEFILES

Second try was the charm for me, and I was able to trim off the bottom of Antarctica using the "Reshape Feature" tool in QGIS. There's a decent description of how to use it here: https://gis.stackexchange.com/a/21651/1770. After removing the southernmost portion of Antarctica, I was able to successfully reproject the map.

To save everyone some time, I've posted the files. These are the "_lakes" files that do NOT include boundary lakes:

Original file in WGS 84 (EPSG:4326), but with the bottom 0.4 degrees of Antarctica trimmed off. You should be able to reproject this file without issue: https://github.com/dmillz/misc/raw/master/shapefiles/ne_10m_admin_0_countries_lakes-trimmed_antarctica.zip
Same file, but reprojected into Web Mercator (EPSG:3857): https://github.com/dmillz/misc/raw/master/shapefiles/ne_10m_admin_0_countries_lakes-EPSG_3857.zip


texts:
I'm trying to Add a Bubble legend with Carto VL to my map but the Bubbles are rather showing as a block and not bubble.
Is the "Bubble-legend" class deprecated or should I style it myself with CSS?


texts:
I am completing a task that requires me to collect info and pictures of gravestones in cemeteries using the Esri Collector app. There are thousands of attachments within the GDB that I would like to be renamed to show the unique section number and also the actual headstone number and the global ID.

So for example it would show [(section#) - (Headstone#) (GlobalID)] 1-1  {9E0682F9-36F6-490C-8CDE-57E3E9BCCF33}. I am quite the novice when it comes to Python.

The obviously wrong script that I have been using is:

import arcpy
from arcpy import da
import os

inTable = arcpy.GetParameterAsText(0)
fileLocation = arcpy.GetParameterAsText(1)

with da.SearchCursor(inTable, ['DATA', 'ATT_NAME', 'ATTACHMENTID', 'GLOBALID']) as cursor:
    for item in cursor:
        attachment = item[0]
        filenum = "ATT" + str(item[2]) + "_"
        filename = 'GLOBALID' + filenum + str(item[1])
        open(fileLocation + os.sep + filename, 'wb').write(attachment.tobytes())
        del item
        del filenum
        del filename
        del attachment

Esri support said to put print statements in the code after (item1) and (item2) but I'm not sure what to put in those actual statements. Attached is what was made when I ran the script. I can tell that the script just put globalID in the name and not the actual global ID.

Also attached is the table in AGOL.


texts: Cannot get features of External WFS on GeoServer

texts:
I need to add an external WFS (https://www.wfs.nrw.de/geobasis/wfs_nw_alkis_vereinfacht) to GeoServer.
I am able to create the store, see all the layers available and publish them.
But when trying to get features from one of the layers, I get this error:
    This XML file does not appear to have any style information associated with it. The document tree is shown below.
<ServiceExceptionReport xmlns="http://www.opengis.net/ogc" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.2.0" xsi:schemaLocation="http://www.opengis.net/ogc http://schemas.opengis.net/wfs/1.0.0/OGC-exception.xsd">
<ServiceException> error:Translator error Translator error org.geotools.data.wfs.internal.WFSException: WFS returned an exception. Originating Request: GET_FEATURE[ handle: GeoTools 23.0(c19e36629fec6970bd3c8495b600f64c1c36a50e) WFS 2.0.0 DataStore @953c2fcf91af#8 outputFormat: <application/gml+xml; version=3.2> method: GET onlineResource: <https://www.wfs.nrw.de/geobasis/wfs_nw_alkis_vereinfacht?> final URL: <https://www.wfs.nrw.de/geobasis/wfs_nw_alkis_vereinfacht?PROPERTYNAME=oid%2Caktualit%2Cgeometrie%2Cnutzart%2Cbez%2Cname&REQUEST=GetFeature&RESULTTYPE=RESULTS&OUTPUTFORMAT=application%2Fgml%2Bxml%3B+version%3D3.2&SRSNAME=urn%3Aogc%3Adef%3Acrs%3AEPSG%3A%3A25832&VERSION=2.0.0&TYPENAMES=ave%3ANutzung&SERVICE=WFS> ] [ErrorInPropertyName: [ Request rejected due to errors. InternalExceptionCode: iiOperationParsingFailed. Reason: Could not resolve ambiguous namespace for element namein PropertyName "name" (for FeatureType ave:Nutzung ). ]] WFS returned an exception. Originating Request: GET_FEATURE[ handle: GeoTools 23.0(c19e36629fec6970bd3c8495b600f64c1c36a50e) WFS 2.0.0 DataStore @953c2fcf91af#8 outputFormat: <application/gml+xml; version=3.2> method: GET onlineResource: <https://www.wfs.nrw.de/geobasis/wfs_nw_alkis_vereinfacht?> final URL: <https://www.wfs.nrw.de/geobasis/wfs_nw_alkis_vereinfacht?PROPERTYNAME=oid%2Caktualit%2Cgeometrie%2Cnutzart%2Cbez%2Cname&REQUEST=GetFeature&RESULTTYPE=RESULTS&OUTPUTFORMAT=application%2Fgml%2Bxml%3B+version%3D3.2&SRSNAME=urn%3Aogc%3Adef%3Acrs%3AEPSG%3A%3A25832&VERSION=2.0.0&TYPENAMES=ave%3ANutzung&SERVICE=WFS> ] [ErrorInPropertyName: [ Request rejected due to errors. InternalExceptionCode: iiOperationParsingFailed. Reason: Could not resolve ambiguous namespace for element namein PropertyName "name" (for FeatureType ave:Nutzung ). ]] </ServiceException>
</ServiceExceptionReport>

Edited for more clarity: I can definitely make direct requests to the aforementioned WFS, what does not work is making requests to my GeoServer layer, that is connected to the other WFS.
The problem has to be something in the connection of my GeoServer with the External WFS.
If anybody recognizes the mistake or can help with connections to External WFS, I'd be thankful!


texts:
I spent days looking in all questions here and trying tutorials.
I need to build a graph for a region (It will be my input in simulation)
I downloaded Postgresql and Postgis and osm2pgsql. Now i can load .osm files into a database.
I want to have a complete database with roads ids so i can generate a table of edges and nodes then build my graph

nodes = intersections and edges = roads for a city using python.

Please any detailed help with steps.I saw a lot of keywords (Postgresql/postgis, Networkx, QGIS ..) but i couldn't manage to assemble a good tutorial

I am a very beginner.


texts:
I work on vue postgresql made with table,
so for save I have to use trigger.

My trigger do a lot of thing and I have to send message to the user to make him understand what happend.

But QGIS only show RAISE EXCEPTION and not the RAISE NOTICE, but the
RAISE EXCEPTION stop the trigger and nothing append.

How to send a message from a trigger to the user of QGIS?


texts: Adding labels to a Carto VL visualization?

texts: It was possible to add labels to polygon visualizations in Carto.js 3.x via the CartoCSS; is there a way to add them in Carto VL visualizations? From what I can tell, Carto VL doesn't have the CartoCSS styling options and I don't know what the alternative or workaround would be.


texts: Why does the QGIS Server WMS layer not load in QGIS Desktop?

texts:
I have created a .QGS file with one PostGIS layer which is the alaska.shp file that comes from the samples. It was loaded via the shp2pgsql utility. When looking at my qgs project referencing the postgis layer, I see the alaska shape, everything appears good.

Then I copy my .QGS file to my cgi-bin directory (along with the qgis_mapserve.fcgi file) to set up my WMS.  Then I create a new project in Quantum GIS (1.7.2) and add WMS layer. I connect to my WMS and select the layer I want to import...so far so good. But after I add the WMS layer to my new QGIS project, I see nothing. The layer is checked "on" to render, and "render" is checked in the bottom right of Quantum GIS.

Is there some way to tell if I have something set up wrong?


texts: RESTful Web Services for Geometry Operations

texts: I'm trying to get an idea what web services exist for doing geometry operations. I know that Esri's GeometryService supports several operations, but does anything else exist? I've started a very simple project built on JTS and Google App Engine to support a Code for America project an am trying to decide if I should build it out or not. Any thoughts?


texts:
I have been using base maps in ArcMap 10.3.1 for the past week and all of a sudden today I cannot load anything.

I was using the National Geographic map with no problems (apart from occasional slowing down) but now when I even click on 'Add basemap' the whole application freezes i.e. I need to press Ctrl Alt Delete to exit ArcMAP.

I have gone to Customize->ArcMap Options and ensured that Hardware acceleration was enabled under data view. I've also cleared the display cache. The only licence settings I've altered is changing between basic and advances. My internet connection is fine as well (i.e. I'm on a work network, not at home).

I've read through the previous posts on adding base maps but I haven't seen any issue like this.


texts: Rotate polygon preserving topology in QGIS?

texts:
Is it possible in QGIS? What I need is:

1) Rotate a polygon inside another one (bigger), obviously the polygons are not overlapped.
2) Rotate a given polygon in a topologically correct layer (like a cadastral parcels layer)

What I would like to avoid is make holes, merge polygon pieces, etc. (the "old way").


texts: Error georefencing an aerial image with QGIS (for relatively small area)

texts:
I have an aerial image (.jpg format) of a relatively small site and have been trying to georeference it and generate a corresponding .tif file. I followed several online tutorials, but I keep getting an error "Transform Failed: Failed to compute GCP transform: Transform is not solvable."

Here is a screenshots, of the QGIS Georeferencer:

And for reference, here are the locations & coordinates of the ground control points used:

Multiple transformation types were used (Thin Plate Spline, Projective, Linear), the resampling method is nearest neighbour, and "Set target resolution" is unchecked.

I wonder if the problem lies in the fact that QGIS truncated the ground control points coordinates (as shown in the table in first screenshot) or some sort of issue with the coordinates format I used (Degrees - Minutes - Seconds).

I also tried fixing the ground control points directly against Google Maps Satellite images via the QGIS QuickMapServices (QMS) tool, but got the same error message.

Would anyone have an idea how I could properly georeference this aerial image?

(Note: this is not for commercial use, only a school project).


texts: ST_SnapToGrid and ST_Union resulting in slithers in the data

texts:
I am currently trying to merge polygons (streets and pavements) into one, the resulting merged polygon has slithers (Image 1).

Changing the tolerance (0.1) has an effect on this but results in slithers elsewhere in the data, and when manually merged (in QGIS) I am getting some outlier points as can been seen in image 2.

I am using this workflow for a large dataset that will be updated regularly so manually going in and fixing the problems wont be a way of solving this.

Has anyone had any similar problems?

Example code can be seen below.

Merge using sql (Image 1):

Manualy merged in QGIS (Image 2):

--Drop and create Type 1/2 adopted layer

    DROP TABLE IF EXISTS highway_network.t12_adopted_network;
    CREATE TABLE highway_network.t12_adopted_network AS
    SELECT ST_Union(ST_SnapToGrid(wkb_geometry,0.01)) as wkb_geometry, t12_usrn as usrn, t12_desc as street
    FROM highway_network.highway_adopted_final
    GROUP BY t12_usrn, t12_desc;

-- Add a serial column and make it a primary key
    ALTER TABLE highway_network.t12_adopted_network  ADD COLUMN ogc_fid serial NOT NULL;
    ALTER TABLE highway_network.t12_adopted_network  ADD CONSTRAINT t12_adopted_network_pk PRIMARY KEY (ogc_fid);

--Drop and create Type 1/2 adopted layer components
    DROP TABLE IF EXISTS highway_network.t12_adopted_network_component;
    CREATE TABLE highway_network.t12_adopted_network_component AS
    SELECT ST_Union(ST_SnapToGrid(wkb_geometry,0.01)) as wkb_geometry, t12_usrn as usrn, t12_desc as street, usage_code as component
    FROM highway_network.highway_adopted_final
    GROUP BY t12_usrn, t12_desc, usage_code;

-- Add a serial column and make it a primary key
    ALTER TABLE highway_network.t12_adopted_network_component ADD COLUMN ogc_fid serial NOT NULL;
    ALTER TABLE highway_network.t12_adopted_network_component ADD CONSTRAINT t12_adopted_network_component_pk PRIMARY KEY (ogc_fid);


texts: Gdal+Geos point in polygon test fails under MinGW64

texts:
I have to build a grid of points and remove the points that falls inside a polygon.
Polygons are defined in a ESRI shapefile, points are in a list.

This is how I make the test:

OGRLayer *layer = mShapefile->GetLayer(lr);
layer->ResetReading();
layer->SetSpatialFilter(&point); //getting only the feature intercepting the point
OGRFeature *ftr;
while (( ftr = layer->GetNextFeature()) != 0) {
    if (point.Within(ftr->GetGeometryRef())) {
         n.good = false;
         break;
    }
}

This code works perfectly on Linux (Ubuntu 14.10, GDAL 1.11).

I had to recompile the library under MinGW64, first I forgot to compile with GEOS and the Within() call always returned false (as defined in the API).

If I compile with GEOS, the result is apparently a random response - some point outside the polygon set is filtered out, other inside are accepted.

Is there any known issue, or workaround, or maybe did I forget some other required library for gdal to work properly?

In a previous version I was only testing for results of GetNextFeature(), and it worked perfectly under linux. Under MinGW I found the same strange result.


texts: How can I show a trip on a map, color-coded by how long the trip actually took?

texts:
I have a series of origin and destination addresses in a spreadsheet, with start times and arrival times for each. I'd like to plot them on a map, with a line between each origin and destination, color-coded by how long the trip took.

For example, if the trip from point A to point B took 1 hours and 45 minutes, display a red line between the points, and if the trip from point C to point D took 37 minutes, display the line between them in green.

My preference would be to use a free tool, but if the only way to do this is with a commercial offering, I'd like to know that as well.


texts: QGIS QNEAT3 AttributeError by creating isodistances/isochrones

texts:
Trying to generate isodistances/isochrones in the QNEAT 3 plugin I got the following error:

> Traceback (most recent call last): File
> "C:/Users/User/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\QNEAT3\algs\IsoAreaAsPolygonsFromPoint.py",
> line 249, in processAlgorithm polygon_featurelist =
> net.calcIsoPolygons(max_dist, interval, output_path) File
> "C:/Users/User/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\QNEAT3\Qneat3Framework.py",
> line 469, in calcIsoPolygons band_in = ds_in.GetRasterBand(1)
> AttributeError: 'NoneType' object has no attribute 'GetRasterBand'
> Execution failed after 53.74 seconds

Parameters:

All layers are in a projected CRS. Matplotlib library in the OSGeo4W Shell has been installed.

I am using QGIS 3.6.1

I have tried to play a bit with different parameters but without success (seems it doesn´t have any influence on the error message). Besides I still don´t understand the meaning of the contour interval parameter...?

QNEAT3 infopage

Related thread (but unfortunately didn´t answer my question).

Creating isochrones in QGIS using QNEAT3

Any thoughts why I´m getting this error?

UPDATE:
Changing the parameters so the contour interval >> raster size:

After reprojecting the network layer to a projected CRS via the appropriate geoprocessing tool:


texts: QgsProcessingException in Processing script QGIS3

texts:
I'm trying to use QgsProcessingException within a processing script if the output parameter is empty or corresponds to something.

In the specific case the output defined is QgsProcessingParameterFolderDestination and I would like to display the nice yellow warning in the main dialog if the path chosen is equal to a specific path and not allowing the algorithm to run.

In the processAlgorithm method I tried to use something like:

# set the output folder
Results = self.parameterAsString(parameters, self.OUTPUT, context)
if Results.startswith('/tmp'):
    raise QgsProcessingException(self.tr('Please choose another path'))

it works nice meaning that the algorithm does not start, but the warning message is written in the Log tab and not in the main dialog tab.

Python API are not so useful (https://qgis.org/pyqgis/3.0/core/Processing/QgsProcessingException.html) while C++ API are a little more verbose (https://qgis.org/api/classQgsProcessingException.html)

BTW: using QGIS master (3.1) on a linux debian machine.

To be clear, left the correct handling of the exception for a missing raster parameter, right what happens with the exception with the code of above:


texts: How to copy attribute table values from one feature class to another

texts:
I am looking for a way to copy attribute values from one attribute table to another. I am mainly looking to do this using python, so would prefer a way to do this with a tool or a a stand-alone python script.

A little background about my problem. I have a collection of polylines situated within a collection of polygons (a watercourse network). The crux of what I want to do is identify, spatially, which polygons hold which polylines and record this in the relevant field. The script below shows how far I have got. I have used python to select the relevant polygon and highlight the polylines that are situated inside it (selected features). However, I have so far been unable to copy these selected features into a particular field of another feature class. I have tried using both the Append (got 'warning' 000597) and Join tools but neither have worked.

import arcpy

reachesBoundaries = "C:\Users\maureen\Documents\ArcGIS\EDRN\Catchments\MainRiversCatchmentsPolygonMarch2013.shp"
edrnLink = "C:\Users\maureen\Documents\ArcGIS\EDRN\EDRN_LINK.shp"
reachField = "CATCHNAME"
reachName = "Aire - Lady Beck Catchment"

arcpy.MakeFeatureLayer_management(edrnLink, "edrnLayers9")
arcpy.MakeFeatureLayer_management(reachesBoundaries, "boundaryLayer8", '"' + str(reachField) + '" =' + "'" + str(reachName) + "'")

arcpy.SelectLayerByLocation_management("edrnLayer9", "WITHIN", "boundaryLayers8")

The script works up the point mentioned above, if anyone has any ideas on how to progress it, please let me know.


texts: Calculating albedo values with landsat 7 ImageCollection/FeatureCollection error in Google Earth Engine?

texts:
I am facing a technical difficulty with Google Earth Engine

I am trying to retrieve the albedo value of the state of Idaho by using Landsat7 TOA reflectance imagery over the span of a year. ('2005-01-01','2005-12-31').

The equation to retrieve the albedo is known as:

var albedo = ((0.356*blue)+(0.130*red)+(0.373*nir)+(0.085*swir)+
(0.072*swir2)- 0.018)/ (1.016);

Here is the whole code that I have so far:

var region = table;
var col05 = donnee.filterBounds(region).filterDate('2005-01-01', '2005-12-31');
var calc_albedo = function(donnee){
  var red = donnee.select('B3');
  var blue = donnee.select('B1');
  var nir = donnee.select('B4');
  var swir = donnee.select('B5');
  var swir2 = donnee.select('B7');
  var albedo = ((0.356*blue)+(0.130*red)+(0.373*nir)+(0.085*swir)+(0.072*swir2)- 0.018)/ (1.016);
  return albedo
}
var albedo_2005 = col05.map(calc_albedo);
var clipped05 = albedo_2005.median().clip(region);
strong textMap.addLayer(clipped05,{min: 0.0, max: 1},'albedo_2005');

var region = table is a shapefile of the Idaho state I have upload via the tool google provides.
var donnee is the ImageCollection "USGS Landsat 7 TOA Reflectance (Orthorectified)"

Finally, the error I am receiving is :

  albedo_2005: Layer error: reduce.median, argument 'collection':
  Invalid type. Expected: ImageCollection. Actual: FeatureCollection.


texts: FME - Filter Parameter choices to only projected coordinate systems

texts:
I have an FME workspace in which the user sets several parameters before running.  One of the parameters is a Coordinate System.

I would like the user to choose only Projected Coordinate Systems as the rest of the workbench is dependent on it. I want this to happen before running the workbench.  Kind of like a validation.

I was thinking a Python Startup Script to only make the projected coordinate systems available.  But I'm unsure how to do that.

I'm comfortable with a Python solution if that's the only way.


texts:
Here is my fiddle, which is exact copy of  Openlayers 3 measure example.

What am I doing wrong when my fiddle isn't working? I can't see the blue dot to start measuring my path or area (same problem like in my localhost environment project)...


texts: When I share my CartoDB-map, in other words,  use the link to my map, I can't get the info windows to appear. It says "there has been an error..." Why? How can i fix this?


texts: Turning layers on/off in web map published to ArcGIS Online from ArcMap

texts: I am trying to create a web service from ArcMap, but I am unable to have the layers turn on and off. I haven't been able to see anything using the Service Editor.


texts:
I am using the unsplit line tool and I noticed that the direction of the lines with multiple segments is getting switched the opposite direction while the single sections seem to continue in the proper direction.

I have enclosed two images showing what I mean. The green lines are the split street centerlines which split at each intersection. The pink lines are the unsplit street centerlines.

  (Green lines by themselves)

  (Pink lines over the green lines)

You can see on Anne St, First St and Second St, that the direction was flipped whereas Burton Pl, Bluff St and Wantland St remained in the proper direction.  I need them to remain in the original direction as it relates to our start and end address ranges.


texts: storage in dams on a river system above specific points

texts:
I want to find the storage in dams snapped to a set of polylines defining rivers above specific points on each river.

Using ArcMap 10.3 I have the polylines representing rivers, the points denoting dams, snapped to the river polylines, with a storage variable, and the points on the rivers where I want to calculate the cumulative upstream dam storage.

How to proceed?


texts: Tweaking some aspects of OpenLayers 3 output using QGIS2Web

texts:
I have made a thematic map in QGIS and exported it to OpenLayers 3 format using QGIS2Web plugin containing 4 layers. I need to tweak the following things if possible:

I want to remove check boxes and include radio buttons for 3 layers
I want to fix the world map layer or at least use check box here so that this layer can be fixed.

This would be a great help as I have browsed through the net but cannot find any resource to do this tweaking.


texts: Why can't Android phones render this KML poly?

texts:
Our outage map works well in browsers and most mobile devices. The map can be found here:
Test Outage Map

However, on Android phones the perimeter polygon is rendered as black. Have tried changing the polystyle but that does not work. The outage polygon is displayed properly if the boundry is commented out.

Screenshot of Android render

Any ideas?


texts: In MapBasic, how do I allow a user to choose a folder?

texts: I have used the FileOpenDlg() function to allow the user to choose a file.  I would like to do the same thing to allow them to choose a folder so I can save files to the folder without the user having to interact every time a file saves.  Is there a function for this?


texts:
I was provided the coordinates of a site in Vietnam in VN-2000 grid coordinate, central meridian: 108 15 zone: 3 degrees, ellipsiod: WGS84

Plotting the  WGS84 coordinates (see attached image) in Google Earth resulted to the relatively correct location. However, when I tried adding the data to ArcMap using Add XY Data and the VN-2000 coordinates, the points were plotted way way east somewhere in Malaysia.

What is the proper way of plotting the points in VN-2000 given the information above?


texts: Using QgsLegendFilterButton to filter a Legend in QGIS using pyqgis

texts:
I am trying to get a layout to autogenerate as there is a certain set of information that I have to produce layouts for hundreds of times. I got the map, legend, and rest of the information to generate but the Legend is too large.  It needs to only show the items in the map.  I tried to do that manually but it still shows all the styles for each layer.  Some of the layers used have 10+ styles so this is causing issues.

I was trying to use QgsLegendFilterButton as that is what I used when making a manual layout map but I can't seem to get that working in pyqgis.  I have have been using guides to get the rest setup but couldn't find any on this so I'm totally lost at how to use QgsLegendFilterButton.  Can anyone post a simple example of adding it to the code below? (if possible)

legend = QgsLayoutItemLegend(layout)
legend.setTitle("Legend")
layerTree = QgsLayerTree()
layerTree.addLayer(layer)
layers = QgsProject.instance().mapLayersByName('Cable')
layer = layers[0]
layerTree.addLayer(layer)
layers = QgsProject.instance().mapLayersByName('Attachment')
layer = layers[0]
layerTree.addLayer(layer)
legend.model().setRootGroup(layerTree)
newFont = QFont("Comic Sans MS", 6)
legend.setStyleFont(QgsLegendStyle.Title,newFont)
legend.setStyleFont(QgsLegendStyle.Subgroup, newFont)
legend.setStyleFont(QgsLegendStyle.SymbolLabel, newFont
layout.addLayoutItem(legend)
legend.attemptMove(QgsLayoutPoint(150,75, QgsUnitTypes.LayoutMillimeters))


texts:
I am brand new to the world of GIS and I have never used Python.

I tried running a calculate field command within the geoprocessing toolbox (not python) and I got an error.  I did a bunch of split line functions to split long lengths of sewer mains into shorter stretches. I split them manhole to manhole.  However, in the attribute table (field name - shape_Leng) the length of each split line stays the same as the original length of the longer line. But in a different field (field name - Shape.STLength()) it does update to the actual length of the split line.  I am trying to run a field calculation to make shape_Leng = Shape.STLength() but I am getting an error.

The reason I am trying to make them equal is because I am trying to integrate it with cartegraph and Cartegraph reads the length of main from the shape_Leng field, whereas I want it to read from the Shape.STLength() field which is the updated length after the split.

I am also trying to figure out how to make cartegraph read from a different field, but have not figured that out yet either.

The error I get from the geoprocessing toolbox function is:

  Parameters Input Table   arcgis.SDE.SANITARY_SEWER_MAIN Field
  Name  Shape_Leng Expression   !Shape.STLength()! Expression
  Type  PYTHON_9.3 Code Block    Updated Input
  Table arcgis.SDE.SANITARY_SEWER_MAIN   Messages Start Time: Friday,
  January 12, 2018 8:23:53 AM  ERROR 000539: Error running expression:
  GPVARIANTOBJECT0.STLength()  Traceback (most recent call last):   File
  "", line 1, in  AttributeError: DescribeGeometry:
  Get attribute STLength does not exist Failed to execute
  (CalculateField). Failed at Friday, January 12, 2018 8:24:09 AM
  (Elapsed Time: 15.18 seconds)


texts: Leaflet - Hide marker before zoom shows marker?

texts:
I have a few thousand markers which are grouped in feature groups which i show depending on the zoom level. At the moment i run the the show/hide function for each feature group after the event zoom end.
But when i'm going from a close zoom level to a smaller one, a few hundrets markers will be shown for a few milliseconds until my function hides them. That's because leaflet shows them before i can hide them with the zoomend event.
This behaviour results in a lag for mobile devices.

What can I do to first run the "Hide/show function" before leaflet shows them all?


texts: Using nested geometric functions in ECQL queries in GeoServer

texts:
When doing a WFS Get Feature request, I am using the cql_filter option to do some more advanced spatial filtering/quering. I am looking for the centroids of the polygons in a database that intersect with a polygon that I specify.

However, when I try to do the following query: "INTERSECTS ( centroid( geom ), POLYGON(SomePoly) )

I get the following error:

<ows:ExceptionText xmlns:ows="http://www.opengis.net/ows/1.1">java.lang.IllegalArgumentException: Binary geometry filter, but first expression is not a property name? (it&amp;apos;s a class org.geotools.filter.function.FilterFunction_centroid)

Binary geometry filter, but first expression is not a property name? (it&apos;s a class org.geotools.filter.function.FilterFunction_centroid)

Does anyone know the correct way of formatting such queries?


texts: ArcPy error handling within a C# wrapper

texts:
I'm running arcPy commands from a C# script i.e.

RunPythonScript("arcpy.Delete_management(\"C:\\\\Temp\\\\NW\\\\land.gdb\", \"Table\")");

private static void RunPythonScript(string script)
        {
            try
            {
                Stopwatch stpw = new Stopwatch();
                stpw.Start();

                string pyfile = "C:\\Temp\\cmd.py";

                StreamWriter sw = new StreamWriter(pyfile);
                sw.WriteLine("import arcpy");
                sw.WriteLine("from arcpy import env");
                sw.Write(script);
...

The problem is this makes debugging very difficult without viewing error messages produced in Python. How can I write errors to a log file or otherwise view the errors being produced by my ArcPy commands?


texts: Creating optional (not required) feature set input in ModelBuilder?

texts:
I am creating a model where the user can specify spatial extent by interactively drawing a polygon (a Clip tool where the clip feature is a Feature Set variable). But the model will be mostly used for the entire extent of the map so I would like the user to be able to just ignore this option if they don't want to use it. However, when running the model it is required because I think the Clip tool has to have an input to run.

Is there a way around it?

Picture of my model:

I am not sure if it could be important for this question but once the model is done, it will be published as a geoproccessing service and later embedded in a Web Application created through ArcGIS Web App Builder as a custom widget.

I am using ArcMap 10.4.1.


texts: QGIS Attribute Table Not Displaying Data

texts:
Using QGIS and one minute I could see the attribute table, the next minute I couldn't. Using QGIS on Linux Mint 13. I tried purge/autoclean/reinstall to get a completely new download of all dependencies and no saved config files, but the problem persists. I can't see the attribute tables for either shapefiles or SpatiaLite data (haven't tested PostGIS yet).

Note, the underlying data is there. I can open the DBF in LibreOffice and see the expected data, and if I click on a feature with the Identify Features tool, I do see the attributes associated with that feature. Symbolization of the data (e.g. a simple choropleth map) also works based on the underlying data. The attribute table displays all expected columns, but no rows, as if it was a newly created completely empty layer.

This was a recent install on a new computer, but for the last week I have been able to add data and view attribute tables as expected. Today I set the configuration option to "Open attribute table in a dock window (QGIS restart required)" and then when I tried look at the attribute table, no data was displayed. I restarted QGIS, still the same. I then turned off the docking, no effect, restarted the computer, no effect. Completely reinstalled QGIS as indicated above, no effect.

One obvious gotcha is that "Show selected only" could be checked, but that is not the case here.

Has anyone run into this or know how to fix it?


texts: Rendering large amount of data in ArcMap?

texts:
Is there a way to load 17 big cadastar plans of urban area in ArcMap and make it "invisible" and use data only when printing to PDF?

Total amount of data is too big to work with and I use it only as a base for urban plans. There are 17 cadastar plans covering area of 220 square kilometers. Cadastar plans are in shp format, made of polygons, lines and points. Largest cadastar plan is around 50MB of shape files and total amount of all of thm is 650 MB.

Every cadastar plan is in separate group of layers. I have adjusted scale ranges for every layer in group of cadastar plan. This way when I zoom in certain part of data shows up.

This is still too much data to process for i7 proc, 16Gb RAM, SSD... Should I create geodatabase file for every cadastar plan? Is there another way to reduce this load? I need every detail of every cadastar plan. Should I remove every attribute I don't need? Will that help?

I can't uncheck group of layers in TOC because data driven pages will need data from all layers, since data is all over city area.


texts: Resources for Learning 508 Compliance/Accessibility Requirements for GIS Figures

texts: Does anyone know of any quality resources (books, papers, presentations, websites) for learning about 508 Compliance/Accessibility requirements in relation to cartography (GIS figures and templates)? This is also in relation to colorblindness and visually-impaired end users.


texts: Create a polygon shapefile using point as boundary

texts:
Following have two problems:

I am now have 4 point in a shapefile and I want to set the point's coordinate for boundary of a polygon
The points are xmin, xmax, ymin, ymax that can just be a rectangular polygon.

I want to make it to create a shapefile to make difference with another polygon shapefile

above problems i want to do it with python script
any ideas?
I use QGIS 1.8.0 and python 2.7.3


texts:
In my workflow using rasterio, I'd like to read an overview from a raster and get only a portion of it through a window at the same time.

Is this possible?

I have a pretty complex script so far and I'm trying to check whether the wrong output is due to this not being possible.

Basically, I am asking if it is possible to do something like this:

data = src.read(out_shape=(1,
                           math.ceil(src.height / 64),
                           math.ceil(src.width / 64)),
                window=window)

where data is a decimated read of src to get source raster overview for factor 64, and window has the width and height whose dimensions are relative to the source raster (not the overview).


texts: Understanding intersection of GDAL and QGIS?

texts:
Background

I'm starting a project that will be done entirely through Python, except perhaps at the end for visualisation. So far I've been doing it with ArcGIS and arcpy, but we'd like to move it to an open source platform.

I've been searching and haven't been able to find out how exactly QGIS fits in with GDAL/OGR (and other open source projects like GRASS). It seems that QGIS includes GDAL/OGR but has lots of its own additions. I'm trying to figure out what's the easiest route to perform a lot of simple, high-level operations.

At the moment what I'm doing with arcpy is converting between rasters and points, extracting raster values into existing points, extracting tables and geometry to NumPy, hydrology, finding intersections and nearest distance between features. There is only one operation that I imagine would have to be reworked quite extensively (inserting points every x metres along a set of vectors).

My question

Is the Python interface of QGIS sufficient to handle this, or is there an advantage in using GDAL/OGR and other packages? I'd like it to be as simple and readable as possible, as lots of other people will be using it and modifying it when we're finished.


texts: How to remove the thousands separator from multiple layers in ArcGIS Portal after publishing from ArcGIS Pro?

texts:
I am publishing approximately 40 layers from ArcGIS Pro 2.5 to ArcGIS Portal. They all have the same data source. Each layer has a numeric 'Year' field. In Pro, the numeric field is formatted without the thousands separator. However, when I publish it to Portal, the thousands separator comma formatting appears and the year is "2,019" rather than "2019" in the pop-ups and data table. How do I prevent this from happening upon publishing? I don't want to fix 40 layers. I need to refresh this data regularly.
I am also trying to avoid creating a new text year field and setting it to equal the numeric field because after that I will have to repeat the of hiding the old field and rearranging the order of the new fields step 40 times (once for each layer) before publishing it again to Portal.


texts:
Does anyone know the XYZ accuracy for the Phased Array L-band Synthetic Aperture Radar (PALSAR) DEMs (12.5 x 12.5 px)?

Is there a document for reference?


texts:
I know there is a similar question but I am getting a different error and none of those answers have helped.
I am trying the following query (the GeoJSON itself is from the PostGIS website.

INSERT INTO layer_radar (id, geom)
VALUES
(
    '3bf24920-225b-11e4-8c21-0800200c9a66',
    ST_GeomFromGeoJSON('{"type":"Point","coordinates":[-48.23456,20.12345]}')
)

The error I am receiving is

ERROR:  (null) (at offset 0)

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

ERROR: (null) (at offset 0)
SQL state: XX000

Has anyone seen anything like this before? I can't get the simplest example to work.


texts:
Is it possible to set the color of the strikeout-line to a different color as used by the font-color in the layer labeling settings?

If yes, how?


texts: How to use the field calculator in QGIS to determine the total of differents fields

texts:
I have an attribute table with 7 fields.

I want my attribute " A1 " to be the total of " IP1 to IP5 ". I am not familiar with the request in the field calculator.

How can I calculate the value of one field so it represents the total of other fields ?


texts:
I have a point layer with meteo-stations data. When I tried to obtain filled contours with contour plug-in it seem ok, but when I tried to work with that layer happened a mistake and I don't know how to solve it.

As you can see in the picture below, in the areas pointed with black arrows doesn't have a value. For example, if you click with info button doesn't appear anything. If you point in the areas pointed with blue arrows, both areas appears selected (they have the same value).

I need to intersect this layer with another one, but algorithm doesn't run because there are invalid geometries (I guess will be those which can not be selected).
 How can I give this part of the layer its attributes??


texts:
I want to count the no data values, because when I loaded them with rioxarray they seemed to be oddly many pixel to be non-valid.
When I open the given file in qgis the valid pixel percentage is 89.6 %.
My attempt to measure it with raster was:
import rioxarray as rxa
import numpy as np

raster = rxa.open_rasterio("Z:/path/to/test.tif")
data = raster.data
print(np.sum(data == raster.rio.nodata) /data.size * 100)

output: 10.122128919010978
Background: I wanted to use numpy masked arrays and mask the e.g. the nodata, black pixels, or values above a threshold. But when I use this method, to mask the nodata pixels, almost everypixel is getting masked.


texts: Working with a small spatial extent from a large dataset in ArcMap

texts:
I am currently working with large amounts of forest data where surveyors and ecologists measure and assess the quality of the vegetation in the field and send the result as CAD files on a daily basis to be converted, edited and stored in a large enterprise geodatabase. It's a rather large project with many people involved.

Even though I've created basemap layers from layers that do not need to be snapped against, it's slow to work. I know that in some CAD software (I only know the basics of CAD), you can load just a small extent of data from a database into the view and work only with that extent without having to render the entire database, but is it possible to do the same thing in ArcMap?

I've tested a lot of workarounds, clipped data to a file-based geodatabase, for example, and then copied the objects to the enterprise database after editing, but it gives a rather clumsy workflow.

If I define the spatial extent in the data frame, is it still working on the entire database? I'm not experiencing improved performance when doing so.

The optimum would have been to open an MXD with a WMS background map, draw a rectangle and load data from the database within the extent from the rectangle but I assume that ArcMap isn't built to work that way?

So, what is the easiest way to work with a small spatial extent of a large dataset?


texts:
What's the difference between open source CartoDB version vs the paid one?

this one the open source: https://github.com/CartoDB/cartodb

and this one the "paid" version: https://carto.com/

i don't understand the advantage of the paid version beside i don't need to setup the server and carto myself.

also, if i use the open source version, which i need to setup it to my own server, is there any disadvantage to it (or any security problem, or will i get "limited" version or something) ?


texts: What are the coordinates of the mysterious hole in Siberia?

texts:
I am interested in taking a Google Earth aerial view of the mysterious hole recently discovered in Siberia and perhaps doing some hydrological analysis of the area in a GIS.  However, there appears to be very little information on the actual coordinates of the hole.  What are the coordinates of this interesting geological feature?
More details on the Siberian hole:

Siberian Times
NBC

Source: Associated Press


texts: Adding arbitrary code to ArcGIS ModelBuilder?

texts:
Is there any way to add arbitrary code snippets to an ArcGIS model?

I envision something similar to the "show code block" the expression builder in advanced labelling or field calculator.

For example my project of the moment could benefit greatly from fetching the username of the person running the tool. I know in python this is a simple import os; os.environ.get('USERNAME'). I think I could write that snippet to get_username.py, add it to a toolbox in Arccatalog, then add that tool to the model and use it as per-requisite model parameter for my current tool. That's a lot of work for such a small thing, and next time I might want something else, like inline variable substitution or a compound variable.

I think a code snippet tool or similar could be address all of these needs in a generic, broadly helpful manner.


texts: Is it possible to show migration arrows on Leaflet map?

texts:
I need to develop Leaflet map which will show the whole world. Also, I need to draw some arrows which will show paths of migrations in the world. The important thing to say that those arrows need to be curved. Do you have some advice how can that be shown on the map in Leaflet API?

My ideas:
Is it possible to show some direction arrows with RaphaelJS? I see that Leaflet has Raphael plugin so we can add Raphael object on the map.
Is it possible to draw a polygon which will be arrow shape and overlay it on the map?


texts: One large layer with filtering vs. many small layers for WMS serving partial datasets in GeoServer

texts:
I am working on a small web application which is supposed to display points on a Leaflet map which are served as WMS from GeoServer, backed by a PostGIS datastore.

Users of the application upload the point data via some HTML/JavaScript forms to the webserver (implemented with Python/Flask) which then parses the data (some csv format) and writes it to the PostGIS database from where it can be served by GeoServer as WMS.

The data and geometry format is the same for each user, however, different users should not interfere with each other but only see their own data on the map. I implemented this by assigning each user visiting the page a unique "upload_id" which is stored as a database column along with his data. Hence, the database looks i.e. like:

id (PK), attribute1, attribute2, ... , geometry, upload_id
1, a, b, ..., POINT(...), jhkgf
2, c, d, ..., POINT(...), jhkgf
3, e, f, ..., POINT(...), jhkgd
4, g, h, ..., POINT(...), sggfz
5, i, j, ..., POINT(...), sggfz
6, k, l, ..., POINT(...), sggfz
7, m, n, ..., POINT(...), sggfz

In this example, the first three points have been uploaded by one user and the last four points by a different user, hence the different upload_id.

Now, in order to display the correct data to each user (his upload_id is available in the JavaScript/Leaflet part of the application), I see two possibilities on how to proceed after the data has been uploaded and inserted into the database as described above:

For a user with upload_id xyz, call GeoServer's REST API to create a new layer called xyz based on a SQL view like SELECT * from data_table where upload_id = xyz and then attach the parameter layer=xyz to all of his WMS GetMap requests in Leaflet. This results in one GeoServer layer per data upload and each data set can be identified via the layer's name.
Initialize the application with one GeoServer layer for all incoming data sets based on a SQL view like SELECT * from data_table, similar to the fact that there is one PostGIS table for all data sets. Then, select the data that is to be displayed to each user by sending a CQL filter parameter or similar with each WMS GetMap request that contains the requested upload_id.

In summary, it comes down to either creating lots of small layers, one for each dataset, or rather filtering the data that is to be displayed out of one large layer. Are there any advantages or disadvantages to either of the approaches, performance wise or other? If that makes a difference, traffic on the application will be low (no more than 3-5 users at the same time) and data sets vary in size but none will contain more than approx. 5000 points.


texts:
What is the best workflow to import a roughly 1x1 km area from AHN data (Actueel Hoogtebestand Nederland) into Sketchup Make? My goal is to create contours of existing buildings and structures.

AHN is publicly available height-map of the Netherlands. It available as point cloud or raster data, but raw and terrain only (maaiveld).

Given that a 1x1km area consists of roughly 4,000,000 datapoints, I suspect that the first step is to create mesh with significantly less points (e.g. 100,000 to 400,000) before importing it into Sketchup.

I'm a novice in GIS, and plan to do this workflow only once or twice, so I'm mostly interested in free or donation-ware type of software. Preferably those that can be used on macOS.

I'll list my current solution as a possible answer below, but so far my solutions are not as best as I wanted.


texts:
I am having trouble using the "Select Layer by Location" tool.
I want to select everything from shapefile A (green) that overlays with Shapefile B (red).

Here is the raw data:

And one example of the "Select Layer by location" (all three of them have the same adjustment, just different input):

Running this, it produces the following output.

In a closeup one can see a few selections. But based on the input data (red - as you can see, it basically covers all of the green area) it should be exactly the same. Or am I totally mistaken?

EDIT 7/14/14 4:58pm

This is how it worked out.
Instead of making one feature layer I made three feature layers and three different selection which need to be saved seperately. Another step was copying the features before merging them. Dont be confused, I added another step to the model (calculate field) but this has nothing to do with the trouble shooting of the original problem.

solution http://s14.postimg.org/6hltytbf9d/Untitled.png


texts: Georeferencing Photos With No Spatial Reference?

texts:
I have about 100-200 field ground shot photos I would like to overlay in ArcGIS 10 that have no spatial information which I need to somehow develop using a photogrammetric formula. The camera was attached 2 meters above ground on a RTK GPS pole. One-to-One relationship between GPS points and photos. Accuracy is not a huge concern, but should at least be within 1 meter (which is approx. the field of view of each photo). The goal is to develop a way to auto-georeference each image using a script or model since it is possible to know the footprint size of each image in the coordinate system I'm using.

What I have to work with:

100-200+ Ground-Shot Field Photos:

-camera height 1.228 meters (from tip of gps)

-camera angle ~90 degrees

-camera focal length 0.006 meters

-North orientation

GPS points shapefile (100-200+ records; 1 record=1 point=1 photo):

-XY location for each photo(approximately the centroid of photo)
-NAD83 UTM Zone 14N (meters)

Basic procedures off the top of my head to accomplish without spending a ton of time georeferencing each individual photo:

*1) Determine the dimension (in meters) of one of the photos. I can calculate Photo Scale (PS = f/H) and Representative Fraction (RF = 1/(H/f)), but I'm needing to calculate length DE in order to find the equivalent corresponding length on ground AB (aka. the width (in meters) I will make for rectangular buffers in the next step).

(source: ucsb.edu)

2) Use those dimensions to create rectangular buffers around the GPS points

3) Spatial Join GPS points to the rectangular buffers (to get the picture name corresponding to each point)

4) Split rectangle buffer layer by attribute (using picture directory/name field)

*5) Replace color fill of each rectangle buffer with the corresponding picture

I'm having trouble with the photogrammetry formulas to use in #1. #5 is more of a hack to visualize the pictures to avoid the georeferencing situation. I don't expect anyone to 'solve my problem'. I like to post problems like these to a) collect my thoughts b) to possibly edify others.

Your thoughts are appreciated.

UPDATED DETAILS:

-photos are ground shots

-North-oriented

-while this problem may seem impractical due to image size vs. scale, I have my reasons :D

Update: 6/20/2013
I calculated the geometric dimensions of each photo (1.52m x 1.14m). I then created rectangular buffers with those dimensions at each gps point centroid. I georeferenced one image by hand and it fits perfectly.

Current problem: Using python, I want to scale down each photo to fit the same dimensions of its correpsonding polygon footprint.I looked at using warpfromfile tool in arcgis 10.1, but need to specify to and from ground control points for both the photo and the footprint which seems a bit tedious. Is there a tool that will auto-scale/fit a raster to the dimensions of a polygon?


texts: Updating spreadsheet from polygons in shapefile using GDAL?

texts:
I am a beginner with GIS data, and have been looking for an answer but have not yet found it.

I have a spreadsheet with lat/long points. I also have a shapefile with several numbered polygons over a geographical area. I would like to do the following:

For each row in the spreadsheet (lat/long point) determine which polygon it is within
Update the spreadsheet (or generate a new spreadsheet/shapefile/etc.) with each row having its polygon number in a new column.

I suspect I could possibly use ogr2ogr to do this, or possibly write some Python code to do it, but I really don't even know where to start.


texts:
I want to set null all raster values that fall below the 75th percentile value. E.g. If you classify the raster by quantile method into 4 classes I want to be left with the top 25% percentile. For the raster shown below, 3.154941635 would be used in the set null expression.

I have attempted the following in the Calculate Value tool and a value of 0 is being returned. (source:https://community.esri.com/thread/162022)


texts:
I am using the VG250 dataset from http://www.geodatenzentrum.de. I am trying to create a new layer that contains only the borders of four federal states that are next to each other.

I want the new layer to basically only contain the blue line. (Sorry for the crude drawing)

My goal is to then run a buffer around these borders and then create a selection of municipalities that lie in the buffer area.

One challenge was that the borders of one federal state were stored as one feature.

What I ended up doing is using the 'split feature' tool to split the border lines. Then I deleted the superfluous borders.

After some trial and error, I think another option to get the exact borders would be to create a layer containing the whole border of each federal state and then using the 'intersection' function to extract the borders of interest.

However, this takes quite some time to do manually, so I still seek a simpler suggestion.


texts: Calculating a spherical polygon centroid

texts:
I'd like a general way to calculate centroids for polygons on a sphere.

So far, the best online reference appears to be:

  Tools for Graphics and Shapes by Jeff Jenness.

The method described there suggests decomposing the polygon into multiple spherical triangles, and computing the average of spherical triangle centroids, weighted by the spherical triangle area.

I know that there are several ways to define a spherical polygon centroid, but I'm looking for something analogous to the following definitions for points and polylines:

Points: arithmetic mean of the Cartesian vectors representing the points.
Polylines: weighted mean of Cartesian vectors representing midpoints of each line segment, weighted by the (spherical) length of each segment.

It seems a reasonable continuation to have polygon centroids defined as the weighted mean of the triangular decomposition, weighted by area.

My question is whether the method in the above reference will work regardless of the triangle decomposition used. In particular, it mentions decomposing into triangles relative to an arbitrary point, even external to the polygon, such that some triangles will have negative areas that contribute a negative weight.

Related: How to find the center of geometry of an object?


texts: OpenLayers: Use scales instead of resolutions

texts:
Venturing around in the OpenLayers source code, it seems like everything related to zoom-levels is done using resolutions. Fair enough since this related directly to screen space.

I would like though, to use scales instead, and further more explicitly define the allowed scales.
- Is this possible?

So far I have created a set of scales which I transform into resolutions:

var scales = [500, 1000, 2000, 4000, 10000, 25000, 50000];
var resolutions = [];
for(var i = 0; i < scales.length; i++) {
    resolutions.push(OpenLayers.Util.getResolutionFromScale(scales[i], units));
}

map = new OpenLayers.Map('', {
    ...
    minScale: scales[scales.length - 1],
    maxScale: scales[0],
    resolution: resolutions,
    numZoomLevels: scales.length,
    ...
});

My problem is that the first and the last scale are respected, but in between the current resolution converted into scale is not among the ones I specified.
- What have I missed?


texts: Resizing Map width and height (avoid distortion)

texts:
I am working with OpenLayers library on React

My first question is how to avoid distortion while resizing map width and height,

The second issue is when I am changing with/height each time it creates a new map.
I know that issue is happening because of initializing new Map each time.

How can resize map without initializing new Map?

import React, { useEffect, useState } from 'react';
import * as ol from 'openlayers';
import { layers_collection, vw } from '../map - constants/map-tools';
import { Button, Input, TextField } from '@material-ui/core';
import MapSource from './map-source';

const initMap = () => {
   return (new ol.Map({
      target: 'map',
      layers: layers_collection,
      view: vw,
   }));
};

var map = initMap()

const OpenLayerMap = () => {
   const [currentWidth, setWidth] = useState<string>("400")

   useEffect(() => {
      initMap()
      map.updateSize()
   });

   return (
      <>
         <Button variant="contained" color="primary" onClick={() => {
               setWidth("400")
         }}>
            400
      </Button>
         <Button variant="contained" color="primary" onClick={() => {
               setWidth("800")
         }} >
            800
     </Button>
         <div id="map" style={{width: `${currentWidth}px`}} >
            <MapSource />
         </div>
      </>
   );
};

export default OpenLayerMap;


texts:
I'm trying to produce a map showing water bodies like oceans/rivers/lakes polygons for purposes of plotting vessel voyages using cartopy. Ultimately I'm looking to create a sea mesh grid using the water polygons. However i'm stuck at the first step in generating the polygons for the water bodies. My marine polys look alright but my lakes rivers centerlines look like the polgons are clipped. Does anyone know how I can fix this? I'm pretty sure the rectangularish lakes are incorrectly displayed. My code in python is as follows:

ax = plt.axes(label = 'main_geoaxes',projection=ccrs.PlateCarree())
rivers_lake_centerlines = cf.NaturalEarthFeature(category='physical', name='rivers_lake_centerlines', scale='ne_10m',facecolor='b')
ax.add_feature(rivers_lake_centerlines,zorder=1)


texts: OGR2OGR FileGDB 2 other format with -SQL query change structure

texts:
I'm using GDAL 1.11.1 with ESRI FileGDB Driver, GDB version 10.x.

When I use ogr2ogr, my goal is to sort differently my fields and in some case, to rename them. The problem is when I call FileGDB driver to an other format, my table structure change for every strings.

For example:

Input FileGDB :string field, 3 of length
Output FileGDB:string field, all fields are set to 0 (ogrinfo) or 65536 (ArcGIS 10.1)
Output SHP:string field, every field is set to 80 (by ogrinfo)

here my cmd.
ogr2ogr -f "FileGDB" TEMP.GDB Source.gdb -sql "select * from input_layer"

I have tried to do an ALTER TABLE nation ALTER COLUMN myfield2 TYPE character(15) but it's do not work on FileGDB?

Any ideas?


texts: Looping through all rows in single column and altering string values of each element in shapefile using ArcPy

texts:
I have a shapefile with an attribute table that partially looks like this:

Row | Shape | Date
0   | Poly  | 05-Sep-10 16:12:00
1   | Poly  | 06-Sep-10 17:23:01

Is there a simple method that I can loop through the 'Date' column and add a 20 before the 10 in each cell? So that the output looks like this:

Row | Shape | Date
0   | Poly  | 05-Sep-2010 16:12:00
1   | Poly  | 06-Sep-2010 17:23:01

Date fields are text format.

I have been able to hand edit multiple rows, but let's say I have 10,000 rows, this would take a considerable amount of time.

Attempt using ArcPy in ArcMap:

import arcpy

path='C:/path/to/file.dbf'

field='Date' #name of the second column

twenty = '20'

with arcpy.da.UpdateCursor(path, field) as cursor:

    for row in cursor:
        new_row = row[0][:7] + twenty + row[0][7:]
        row = new_row
        print(row)
        cursor.updateRow(row)

Error I receive with this is:

  TypeError: sequence of size must match the size of the row


texts:
I would like to add a radio button control (like a list) to my leaflet map (which look like basemap dropdown menu).

So users can switch the apperance of some layers by simply selecting corresponding radio button, e.g., display big circles rather than small ones.

How is that possible?

var radius = 600;
var smallRadiusCircle = L.circle([44.710, 6.84], radius, {fill: true, fillColor:'gray', fillOpacity: 0.4, color: 'none', weight: 1});

var radius2 = 1800;
var bigRadiusCircle = L.circle([44.710, 6.84], radius2, {fill: true, fillColor:'gray', fillOpacity: 0.4, color: 'none', weight: 1});


texts: Is it possible to make GeoServer filter features by scale dynamically

texts:
Suppose there are 100,000 points, all of them will be rendered at a small scale with the default style, we can use Rules in SLD style to filter features by properties. But most of the tutorial told me to render, but there are still two problems:

There are no extra properties for the layer

For example, if the feature only have a name property, every feature have the same property, it would be better if GeoServer can filter/simplify the features (maybe discard some of them randomly) according to the current scale.

There are more properties related to the styling.

For example, if two properties(area,typecode) have to be considered during the styling and priority = f(area,typecode), then if two features will overlap at the current scale, feature with the lower priority should be discarded.

Is this possible in GeoServer?


texts:
I am developing a PostGIS query to sample rasters using a point table to detect land cover change. I have read the documentation on ST_value and I am able to get the raster value of a single point using the following query:

SELECT  rid, ST_Value(rast, 1, ST_SetSRID(ST_Point(484375.820,4742079.979), 32617)) as b1
FROM imagery.l8_2015_09_12
WHERE ST_Intersects(rast, ST_SetSRID(ST_Point(484375.820,4742079.979), 32617)::geometry, 1);

However, when sampling using a table of points, I cannot get ST_Value to work. When running the following query in the SQL window of QGIS:

-- sample a raster from a table of points
SELECT ST_Value(r.rast, 3, p.geom) As band3
FROM analysis_results.sample_points AS p, imagery.l8_2015_09_12 AS r
WHERE ST_Intersects(r.rast,p.geom);

I get the following error:

Attempting to get the value of a pixel with a non-point geometry

The table has a multipoint geometry data type, so why is ST_value not playing nicely?


texts:
I want to create a table (table_diff) that is a result of symmetrical difference (ST_SymDifference) between geometries (linestrings) from 2 tables (table1, table2) stored in 2 databases (db1, db2).

So far I have a command that works on 2 tables that are in the same database:

psql -U postgres -d db1 -c "CREATE TABLE table_diff AS SELECT ST_SymDifference(ST_Union(DISTINCT table1.geom),ST_Union(DISTINCT table2.geom)) AS geom FROM public.table1, public.table2

Below is my effort in trying to implement dblink to access table2 from db2 which doesn't work and which I need help with.

psql -U postgres -d db1 -c "CREATE TABLE table_diff AS SELECT ST_SymDifference(ST_Union(DISTINCT table1.geom),ST_Union(DISTINCT dblink('dbname=db2', 'SELECT geom FROM public.table2'))) AS geom FROM public.table1, dblink('dbname=db', 'SELECT geom FROM public.table2')"

The error I am getting:

ERROR:  a column definition list is required for functions returning "record"


texts: +0-180 degrees of mesh layer are displaying as horizontal lines [QGIS]

texts:
I have been trying to display weather data from ECMWF in QGIS on top of a base map, but half of the heatmap rendered from the mesh layer (imported via grib) displays as straight horizontal lines.

I am using QGIS 3.12.1
The base layer are XYZ tiles (Pseudo-Mercator EPSG:3857).
The overlay was imported from ECMWF in .grib format and imported as a mesh layer.
The two layers are clearly aligned, as the windspeeds clearly drop around the shape of South America and match the base layer precisely.
When I also display vector arrows for wind speed, they seem to cover the whole map:

Why is the heatmap imagery not displaying correctly at +0-180degrees, and how can I fix this?


texts: OpenLayers: How to get *real* centroid of a Line Geometry/Feature?

texts:
I have declared a line geometry as:

var lineString = new OpenLayers.Geometry.LineString([source,target]);

then I am trying to get the centroid of this lineString as lineString.getCentroid().

But what's returning is a Start Point of the lineString.

So, is there any way to get the real centroid of a LineString?

P.S: I have also looked through this discussion. But couldn't figure out the solution.


texts: How to show label according to attribute name?

texts:
I try to show label name according to predefined in QGIS 3.0 by using pyqgis. I used that code in QGIS 3.0 but last update broke it down. Here is my code:

    settings = QgsPalLayerSettings()
    settings.fieldName = 'Subbasin'
    labeling = QgsVectorLayerSimpleLabeling(settings)
    wshedLayer.setLabeling(labeling)
    wshedLayer.triggerRepaint()


texts: Editing PostGIS multipoint table with QGIS

texts:
I created an empty multipoint table in PostGIS and I want to fill it with hand-placed points in QGIS. I have connected them and checked that QGIS identifies the new layer as multipoint. I've got a problem:

I am unable to add multipoints neither adding parts to a point nor combining two or more points.


texts:
I'm trying to cut/create polygons in QGIS 3.2.1 using snapping/tracing, but it only works intermittently - it will follow an existing vector line quite happily for a short distance, then drop it, then pick it up again some distance later. So I end up manually plotting all the vertices, which is very time consuming.

I have had no issues using snapping/tracing previously in this version of QGIS, but here I am using a PostGIS layer (in offline mode), so perhaps this is causing the problem?

I've tried tweaking the various snapping options, tracing different vector layers, saving/reloading etc but nothing seems to make any difference.


texts: Creating a GRASS addon that returns something when using read_command()

texts: I am using grass 7.8.4. I'm calling commands externally using the Grass Python Scripting Library in Python. I wanted to know when you create an addon in Python how to end the addon script so that when you call that addon with read_command() it returns something. I would like to have something similar to what happens when calling read_command(g.gisenv) which returns the grass variables settings.


texts: Iteration of ST_GeneratePoints results in identical points

texts:
I am relatively new to doing this type of operation from SQL so perhaps I have made an poor assumption.  The code below results in the exact same points being generated n times.  I am trying to repeat the creation of random points n times.  For the sake of this test the number of points is set to 10 but will vary by sample.

    CREATE OR REPLACE FUNCTION repeating_random_tree_generator(n INTEGER)
    RETURNS VOID
AS $$
BEGIN
    --clean up
    DROP TABLE IF EXISTS random_trees;

    --new table
    CREATE TABLE random_trees (
        run_number int,
        tree bigint
    );
    --ADD POINT GEOMETRY UTM 11 N
    PERFORM AddGeometryColumn ('public','random_trees','geom',26911,'POINT',2);
    FOR i in 1..n LOOP
        --create 10 random points in box
        CREATE TABLE temp_random_trees as (
        SELECT (row_number() over ()) as
            tree,
            geom
         from
            (SELECT
                (ST_Dump
                    (ST_GeneratePoints(
                        ST_GeomFromText('Polygon ((553690 5435710, 553690 5435700, 553700 5435700, 553700 5435710, 553690 5435710))',26911)
                        , 10)
                    )
                ).geom) r
        );
        --add to random tree table
        INSERT INTO random_trees (
            SELECT 1,tree,geom from temp_random_trees
        );
        DROP TABLE IF EXISTS temp_random_trees;

    END LOOP;
RETURN;
END;$$
LANGUAGE plpgsql;


texts: How can I create a GeoServer store from DigitalGlobe Global Basemap tiles?

texts:
I have been given data from a DigitalGlobe product called Global Basemap. I've been tasked with getting this data into a form usable by our map server, which is running GeoServer 2.4.

The data did not include any form of documentation on the Global Basemap format and I have not been able to find any online. (I found plenty of marketing materials, but no technical specs.) The organization that purchased the data for me does not have any documentation.

The existing raster maps we are hosting on our GeoServer instance take the form of ImagePyramid stores, which I created by running gdal_retile.py on GeoTIFF mosaics. This new, Global Basemap imagery I have been given looks like it's already in some sort of image pyramid format, but not the same format required by GeoServer's ImagePyramid plugin. At the root level of this data there are dozens of directories with names like 8_462_181 and 8_463_173. Inside each one of those directories we have directories and files like so:

drwxr-xr-x  10
drwxr-xr-x  11
drwxr-xr-x  12
drwxr-xr-x  13
drwxr-xr-x  14
drwxr-xr-x  15
drwxr-xr-x  16
drwxr-xr-x  17
drwxr-xr-x  9
-rw-r--r--  cacheViewer.html
-rw-r--r--  root.kml
drwxr-xr-x  var

Directories 9 through 17 contain subdirectories. 10 has more subdirectories than 9, 11 has more than 10, etc. and that is what makes me think that this structure represents an image pyramid. Each of those subdirectories contains JPEG files with matching KML, .proj, and world (.jgw) files. For example:

44542.jgw
44542.jpg
44542.kml
44542.prj
44543.jgw
44543.jpg
44543.kml
44543.prj
44544.jgw
44544.jpg
44544.kml
44544.prj

The structure required by GeoServer's ImagePyramid plugin is that the lowest level image tile files go into a directory name 1, the next level up, i.e., fewer tiles, go into a directory named 2, etc.

Question 1: Do I merely need to rename the existing directories 9 through 17 as 8 through 1.

Question 2: There is no over-arching image, world file, or mosaic that covers the entire region. Will I have to use gdal_merge or something similar to create one?


texts: Add secured WMS service to ArcGIS Enterprise web map

texts:
How can I add a secured WMS (username & password) to a Web Map in ArcGIS Enterprise?
I only see these settings:

Where can I specify password and username?


texts: PostGIS: Get closest point between geometry and Multipoint

texts:
I think my question is not that difficult for GIS professionals, but i didn't found an answer to my problem (even not in other questions).

Problem:
I have given two tables (table1,table2) with Points as geometries in my PostgreSQL-database. Now I'm searching for each point in table1 the closest point from all points in table2. after that the closest point (from table2) should be build to a linestring with the current point from table1.

my approach:

CREATE TABLE new_table_shortest_line WITH (OIDS) AS
SELECT a.name,a.type,
       ST_MakeLine(a.the_geom,
              ST_ClosestPoint(a.the_geom,ST_Collect(b.the_geom))
       ) AS closest_line
FROM table1 AS a, table1 AS b
GROUP BY a.name, a.type, a.the_geom, b.the_geom

the result of this SQL-query should be a new table with a linestrings as geometries, which define the closest lines between a point in table1 and the closest point of table2...


texts: How to filter heights using LAScatalog in lidR

texts:
I would like to filter the points that are used in a canopy height model. In particular I would like filter out any points whose height is over 200 feet. I've tried the code below but when I plot the difference between my normalized and filtered tiles, they are exactly the same and there are still z values above 200 feet. I must be using the opt_filter command incorrectly, or have chosen the wrong mechanism.
library(lidR)

indata <- suppressWarnings(readLAScatalog ("c:/test/laz"))

# Normalize point cloud data
opt_output_files(indata) <- "c:/test/norm/norm_{XRIGHT}_{YBOTTOM}"
normalized_las <- normalize_height (indata, tin(), na.rm=TRUE)

# Keep points where height is between 0 and 200 feet
filtered_las <- suppressWarnings(readLAScatalog ("c:/test/norm/"))
opt_filter (filtered_las) <- "Z>=0 and Z<=200"
opt_output_files(filtered_las) <- "c:/test/filt/filt_{XRIGHT}_{YBOTTOM}"
filt_ctg <- catalog_retile(filtered_las)

# Now read in the same filtered and normalized tile
filt1 <- 'c:/test/filt/filt_2885000_2320000.las'
filtlas <- readLAS(filt1)
norm1 <- 'c:/test/norm/norm_2885000_2320000.las'
normlas <- readLAS(norm1)

# Visualize difference between filtered and unfiltered points
par(mfrow=c(2,2))
boxplot(normlas@data$Z, ylim = c(-10,1000), ylab = "Height (ft)", main = "Before Filtering")
boxplot(filtlas@data$Z, ylim = c(-10,1000), ylab = "Height (ft)", main = "After Filtering")
hist(normlas@data$Z, xlim=c(-2,1000), xlab = "Height (ft)", main = "Before Filtering")
hist(filtlas@data$Z, xlim=c(-2,1000), xlab = "Height (ft)", main = "After Filtering")


texts:
I am trying to perform a "t-test" for 4 population rasters for the entire tropical strip for my graduate thesis. Sample years are 2000-2005-2010-2015.
I am using this script:

setwd("C:/Users/MICITT/Desktop/TESIS_MSC_MARCOMM/Poblacion")
require(raster)
require(rgdal)
pob00 <- raster("PRUEBA_rem_pob00_SC.tif")
pob05 <- raster("PRUEBA_rem_pob05_SC.tif")
pob10 <- raster("PRUEBA_rem_pob10_SC.tif")
pob15 <- raster("PRUEBA_rem_pob15_SC.tif")
pob_ttest <- t.test(pob00, pob05, pob10, pob15, alternative = "two.sided", mu = 0, paired = FALSE, var.equal = FALSE,
       conf.level = 0.95)
writeRaster(pob_ttest, filename = "C:/Users/MICITT/Desktop/TESIS_MSC_MARCOMM/Poblacion/ttestprueba", format="GTiff", overwrite=TRUE"

However, when I try to save the raster, it displays the following error:

Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘writeRaster’ for signature ‘"htest", "character"’

I need to get a final raster that gives me the value of "t", "df" and "p-value" per pixel.
Anyone know how to solve this?


texts:
I'm new to pyqgis and I'm tryig to make a combobox with list of fields (from layer also selected in combobox).

    layers = self.iface.legendInterface().layers()
    layer_list = []
    for layer in layers:
        layer_list.append(layer.name())
    self.dlg.comboBox.addItems(layer_list)

    selectedLayerIndex = self.dlg.comboBox.currentIndex()
    selectedLayer = layers[selectedLayerIndex]
    fields = selectedLayer.pendingFields()
    fieldnames = [field.name() for field in fields]
    self.dlg.comboBox_2.addItems(fieldnames)

The problem is, that in second combobox there are fields from first on layerlist layer, not the selected one.


texts: PyQGIS -- write print composer to project file

texts:
Is there a way to write a print composer to a QGS project file via PyQGIS, so that when the project file is opened later, that print composer is available?

I am able to successfully do the following:

load and manipulate a project (load layers, etc)
load a print composer template (QPT file) and manipulate the composer elements
save the modified QGS file
export the print composer map to PDF

However, I can't figure out how to save that print composer to the QGS file.

Essentially that means that if I want to go into the QGS file afterwards and make edits to the map generated from the print composer, I need to load a new print composer from template and redo all the edits from scratch (map title, extents, map notes, etc).

The following code works to load a composer from a QPT template file and export the map as PDF, but does not save the composer with the QGS file:

#initialize QGIS
app = QgsApplication([], True)
app.initQgis()

#set project and template paths
project_path = "/mypath/myproject.qgs"
template_path = "/mypath/mytemplate.qpt"

#get canvas and load project
canvas = QgsMapCanvas()
project = QgsProject.instance()
project.read(QFileInfo(project_path))

#load print composer template
template_file = file(template_path)
template_content = template_file.read()
template_file.close()

document = QDomDocument()
document.setContent(template_content)
composition = QgsComposition(canvas.mapSettings())
composition.loadFromTemplate(document)

#refresh after manipulating composer elements
composition.refreshItems()

#export to PDF and also save a QGS file
composition.exportAsPDF('myfile.pdf')
project.write(QFileInfo('myfile.qgs'))

I believe I must be missing something either in, or prior to, the last line of code:

project.write(QFileInfo('myfile.qgs'))

to add the composer to the QGS file.


texts: Calculating maximum extent of debris flow from point with mean slope value

texts:
At a very large scale I perform a spatial analysis for which I have to apply mean slope criteria. These criteria allow me to calculate the maximum extent of an event (debris flow). I know that below a certain slope threshold, my debris flow stops. In theory, finding this extent seems simple, but in practice I can't find any GIS tools that allow me to get this.

The data I have are:

The starting points (raster or vector)
The average slope value
The flow paths for each starting point
The altitude values (MNT)

I expect a spatial extent (raster or vector) as a result.

Is there a tool to calculate a maximum spatial extent by entering a starting point and an average slope as a parameter?

I work mainly with QGIS and SAGA so I would prefer solutions for these softwares


texts: Understanding why LatLng returned from L.GeometryUtil.closest() are not from array being checked over?

texts:
I am using this function below from this library Leaflet.GeometryUtil to help find the closest distance between two points.
    L.GeometryUtil,closest(map, array, latlng)

It is returning this object below, I can't find this lat and lng in the array I am checking over to see which one is closest to my starting point.
Could someone explain where the lat and lng comes from?
    {lat: -85.0511287798, lng: 32.577928810608014, distance: 3049.2708964603326}

This is what I have done:
var map = L.map("map", { zoomControl: false }).setView(
 [37.8, -96],
 4
);

L.tileLayer('http://{s}.tile.osm.org/{z}/{x}/{y}.png', {
    attribution: '&copy; OpenStreetMap contributors'
}).addTo(map);

const myLocation = {
 lat: 37.0076752995451,
 lon:-120.5956957740688
}

//add starting points
 L.marker([myLocation.lat, myLocation.lon]).addTo(map)

 getAllSpots()

 var spotsArray = []
 var allCoords = []

//call function to get all surf spots
 function getAllSpots() {
    fetch("http://api.spitcast.com/api/spot/all", {
      method: "GET"
    })
      .then(res => res.json())
      .then(json => {
        json.map(item => {
          var markers = {
            lon: item.latitude,
            lng: item.longitude,
            spot_name: item.spot_name,
            county_name: item.county_name
          };
          var coords =
          [ item.longitude,
             item.latitude]
          spotsArray.push(markers);
          allCoords.push(coords);

        });

        findDistances(allCoords)
      })
  }

    function findDistances(spotsArray) {
            console.log("allCoords",allCoords)
            let closest_latlng = L.GeometryUtil.closest(map, allCoords, [myLocation.lat, myLocation.lon])
console.log(closest_latlng)
L.marker([closest_latlng.lng, closest_latlng.lat]).addTo(map)
  }

    var mcg = L.markerClusterGroup();

      var markersArray = [];

 function addMarkers(json) {
json.forEach(function(feature) {

        markers = mcg.addLayer(
          L.marker([feature.lng, feature.lat], {
            icon: L.icon({
    iconUrl: "map-marker.svg",
    iconSize: [64, 74],
    iconAnchor: [32, 74],
    popupAnchor: [0, -70]
  })
          })
        );

        markersArray.push(mcg);
        return markers;
      })
     featureGroup = L.featureGroup(markersArray).addTo(map);
  }


texts: Adding a sequential number using ArcGIS field calculator to string file and keeping the 0 in front of the number?

texts:
I want to know if anyone has a code to count +1 starting at a certain number for Field Calc. ArcGIS. (I.e. 04312, 04313) It also must keep the 0 in front.

rec=0
def autoIncrement():
 global rec
 pStart = 040150117669
 pInterval = 1
 if (rec == 0):
  rec = pStart
 else:
  rec = rec + pInterval
 return rec

    str(autoIncrement( )).zfill(12)

parsing error syntaxerror:invalid token (line 4)


texts: Trouble switching the project CRS in a new project in QGIS

texts:
I'm trying to change the project CRS from WGS 84 to NAD83/ UTM Zone 17N. However QGS will not let me and an error pops up. The text from the error is:

An error has occurred while executing Python code:

_core.QgsCsException: forward transform of (inf, inf) Error: Fallback transform failed
Traceback (most recent call last):
  File "C:/Users/haile/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\go2streetview\go2streetview.py", line 466, in setPosition
    actualSRS = self.transformToCurrentSRS(actualWGS84)
  File "C:/Users/haile/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\go2streetview\go2streetview.py", line 623, in transformToCurrentSRS
    return xform.transform(pPoint) # forward transformation: src -> dest
_core.QgsCsException: forward transform of
(inf, inf)
Error: Fallback transform failed

Python version: 3.7.0 (v3.7.0:1bf9cc5093, Jun 27 2018, 04:59:51) [MSC v.1914 64 bit (AMD64)]
QGIS version: 3.14.16-Pi Pi, df27394552

Python Path:
C:/Users/haile/DOWNLO~1/GEOG36~1/apps/qgis/./python
C:/Users/haile/AppData/Roaming/QGIS/QGIS3\profiles\default/python
C:/Users/haile/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins
C:/Users/haile/DOWNLO~1/GEOG36~1/apps/qgis/./python/plugins
C:\Users\haile\DOWNLO~1\GEOG36~1\apps\Python37
C:\Users\haile\DOWNLO~1\GEOG36~1\apps\Python37\Scripts
C:\Users\haile\Downloads\Geog 366\bin\python37.zip
C:\Users\haile\DOWNLO~1\GEOG36~1\apps\Python37\DLLs
C:\Users\haile\DOWNLO~1\GEOG36~1\apps\Python37\lib
C:\Users\haile\Downloads\Geog 366\bin
C:\Users\haile\DOWNLO~1\GEOG36~1\apps\Python37\lib\site-packages
C:\Users\haile\DOWNLO~1\GEOG36~1\apps\Python37\lib\site-packages\win32
C:\Users\haile\DOWNLO~1\GEOG36~1\apps\Python37\lib\site-packages\win32\lib
C:\Users\haile\DOWNLO~1\GEOG36~1\apps\Python37\lib\site-packages\Pythonwin
C:/Users/haile/AppData/Roaming/QGIS/QGIS3\profiles\default/python
C:/Users/haile/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\go2streetview
C:\Users\haile\AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugins\go2streetview\res
C:\Users\haile\AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugins\HCMGIS/forms

Another error also pops up when I add the layer I want and change that CRS from WGS 84 to NAD83/ UTM Zone 17N:
An alternative, ballpark-only transform was used when transforming coordinates between EPSG:4326 - WGS 84 and EPSG:26917 - NAD83 / UTM zone 17N. The results may not match those obtained by using the preferred operation:
Possibly an incorrect choice of operation was made for transformations between these reference systems. Check the Project Properties and ensure that the selected transform operations are applicable over the whole extent of the current project.
I have tried uninstalling and reinstalling to see if that would help, I have also used both 3.12 and 3.14 to see if the version was the issue. My classmates have been able to switch the same project CRS without any issues on their PCs.