texts: Converting from shapefile or raster formats to LAS?

texts:
I am looking for a solution to convert either from ESRI shapefile format to LAS or from a common raster format to LAS. I have access to following software:

ArcGIS 10.2 for Desktop (arcinfo with spatial analyst and 3d analyst extensions),
Global Mapper,
LASTools (unlicensed),
QGIS 2.8,
Whitebox GAT.

I have the functionality in the LASTools shp2las.exe tool, but as I do not have a licence it can only handle 2 million points before it adds XYZ noise.  Due to the large volume of data I am working with, it is impractical to cut it into thousands of smaller pieces in order to process with tiles under 2 million points.

There was the potential of a workaround using the text2las tool, but outputs from this tool are also unacceptable due to a stepping/contouring effect. Despite efforts to ensure no rounding or truncation of input data occurs during transformation.

Is there a way through open source or free software to successfully convert between these formats?

Or am I looking at buying a licence for FME or LASTools or some other software?


texts: Projected coordinate system for South America

texts: Does anybody have any suggestions for a projected coordinate system for the Northern part of South America?  I have a data set that reaches from Venezuela to Southern Peru and I'd like to do some spatial autocorrelation on it.  However, the current WGS84 is obviously a geographic coordinate system and not suitable for this type of analysis.  I could just stick it in UTM 19 or 20 (north or south I guess given it spans the Equator) but I just wanted to know if anyone had any better suggestions.


texts: Moving GPS points with distance + bearing

texts:
I have a list of 73 GPS points that represent animal observations. Each point also lists the bearing (angle in degrees) and distance from the observer to the animal. Instead of displaying each point at the recorded GPS location, I would like to display the location of the animal (the location of the GPS point, offset by the distance and bearing recorded with the point). Does anyone know how to move points in this way?

Other solutions on Stack Exchange suggest using a VB script, which I am not familiar with. But I would be happy to try with a breakdown of the steps.

I am using QGIS 2.18.7.


texts: OpenLayers 3 dynamically add layer from folder contents

texts:
I am new to OpenLayers and have searched far and wide trying to figure out how to programmatically add layers with a default set of parameters from a selected source path, but Google has not produced any meaningful results (at least that I can identify).

I have generated a set of thumbnails to represent each map layer using C# from image files, which I want to use to provide a clickable element that adds the layer to the OL3 map to view/manipulate by the user. The images are dynamically generated, so ideally I want to create a dynamic element (via ASP, JS, or something else) that will feed a full sized image to the map layer. As I see it, it could be done by dynamically changing the URL parameter for a map layer, or adding a whole new layer for each click event on an image. Has anyone tried to do this, or something similar? Perhaps there is a better way to serve images to layers from a folder structure?


texts:
Can anyone explain why for gdal_calc results differ in the two cases below?

case 1: (A>0)-(A<10)+(A>0)-(A<10)
case 2: ((A>0)-(A<10)) + ((A>0)-(A<10))

Mathematically speaking, the two should give the very same result.

With gdal_calc, instead, in case 1 I unexpectedly obtain a raster of zeros and ones.
In case 2 I obtain a raster of zeros and twos, as wished.


texts: Selecting certain fields in layer using ArcPy?

texts:
I have a layer in geodatabase, and I want to use SelectLayerByAttribute to select certain fields in this layer.

Then, use CopyFeature to output those fields to a new shapefile.

for example:

layer = I created in gdb

selection_type = 'NEW_SELECTION'

SelectedFields = """ Field_name = 'XXXX' """

arcpy.SelectLayerByAttribute_management(layer, selection_type, SelectedFields)

But, my where_clause part always invalid.

anyone can help me get the correct SQL expression?


texts:
I have lat/long and when I convert to UTM, my points are in zones 35 and 36.
is it possible to "transform" or "project" (I am not sure which word is used to do this operation) the UTM coordinates from zone 36 to zone 35 ?

Edit: Just to clarify
I have those two coordinates.

point 1: is in 35L 460880.04 E 8327740.29 S
point 2: is in 36L 201956.55 E 8375498.52 S

is there any way to say: point 2 is in 35L xxxxxxx E yyyyyy S ?


texts: Error 99998 when running batch operations?

texts:
I'm running a fairly large script in batch. I've tried it two times now, at it crashes at about 4 hours in and chrashes maybe 20 min into the second batch, giving me the dreaded 99998 error. I've tried to restart both the computer and ArcMap.

ExecuteError: ERROR 999998: Unexpected Error.

I expect it might have something to do with memory, but I had open the task manager while running the script and it only used about ~30 % of the processor and ~30 % memory.

Im not really sure about this, but could it be a problem with arcpy.MakeFeatureLayer_management() or the making of temporary datasets, which I dont delete? An example would be like this:

outCon = Con(viewshed, 1, 0, "value = 0")

I don't delete outCon anywhere, as I thought they might only be temporary.

How would I proceed with that?


texts: No save option for layer in ArcGIS Online

texts:
I am very new to ArcGIS Online. I created a new layer in a map, but there is nor save as button in menu
I tried to create a layer in my map using Map Notes, but after finishing I noticed there is no "save" option either for the layer, nor for the project.


texts: QGIS Function editor calculate area after transform

texts:
In the QGIS function editor, it is possible to transform geometry using transform(geometry,'EPSG:xxx','EPSG:yyy') - however there isn't an area() function to recalculate the area of this.

The reason I am looking to do this is to have area calculated when a map is produced (labelled and as a virtual field in the attribute table); rather than store the area within the database as it is likely to change frequently as users modify the polygons?


texts: Create grid of lat long in polygon at 1 km

texts:
I want all hotel list in that polygon, so I am using the Google Place API to fetch hotel info. But Google Place gives only 60 records limit, so thought to get all lat long with 1 km distance in that polygon, so I can search Google Place with 1 km, and I will get all results.
Solution I tried:

Create polygon in google earth and download KML file
Convert it into GeoJSON file
Now I have polygon lat long co-ordinates
Now I want to draw grid of lat long with distance of 1 KM in same polygon.

I want above yellow lat long which are separated with 1 km in that polygon
Need list of Lat long as below yellow dot.

My data
 {
"type": "FeatureCollection",
"name": "Test",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
"features": [
{ "type": "Feature", "properties": { "Name": "Test", "description": null }, "geometry": { "type": "Polygon", "coordinates":
[ [ [ 76.95921132064926, 28.462199988668662 ], [ 76.963367015913107, 28.454396340060121 ], [ 76.96406171244611, 28.448923074655479 ],
 [ 76.979506480632438, 28.4628132424059 ], [ 76.968779231398202, 28.469906308302541 ], [ 76.961281573952363, 28.472337211248551 ],
 [ 76.95921132064926, 28.462199988668662 ] ] ] } }
]
}


texts: Digitizing lots of points and polygons outside GIS software?

texts:
I have orthophoto of marine bottom from which I want to extract several informations. This information can be for example position of animals (i.e. a point) or the position/contour of a rock (i.e. a polygon). My goal is to produce several vector layers each corresponding to a specific variable (e.g. one layer for seastar position, one layer for rocks, and so on).

I know I can import and georeference the image in QGIS (for example) and there create the layer from scratch adding each single point or polygon. In this way, every time I will add a feature I will be prompted a dialogue box where I should enter some information (e.g. feature ID). The problem in doing this is that having lot of features (e.g. 1000 of points) I will have to do a lot of clicking and moving with the mouse, and that is going to be extremely time consuming.

I was wondering if there is any known workflow/strategy, for example, to just click on the location I want the points to be, then obtain the X-Y position of the point and save their coordinates as a layer (or a csv). Same reasoning with polygons.

TO be noticed is that, while the type of information I want are spatial, I don't necessarily need to have shapefiles georeferenced in any "real word" reference system: I mainly need position of features (points and polygons) relative to each other. For example, I think it could be ok to have all positions in Pixels of the image (subsequently I guess that it could be possible to attach a lat long values to single pixels once the image will be georeferenced eventually, but correct me if I'm wrong).

Therefore I am not constrained to use a GIS software.

It seems that something like this will be possible using ImageJ but I did not figured out yet how to do it and thus I am asking here if anyone knows an already tested and fast strategy/method to achieve my goal.


texts:
I am trying to vectorize a raster with tree coverage information. I classified the data in 3 categories and now I am vectorizing the raster to make calculations and intersecting it with other layers.

The problem is that after vectorizing almost all polygons have invalid geometries because of missing nodes in the intersections. I attached an image where you can see the problem in the arrows. The selected polygon has 2 missing nodes.

I cannot use any vector tool with the layer because of this.

I need to intersect this layer with other complex one, but as a trial I do with a simple temporary layer and I have the same mistake: Input layer A contains invalid geometries (feature 71). Unable to complete intersection algorithm. (I dont know where to check what is feature71).

When I run the Topology checker it gaves me hundreds of mistakes with no valid geometry.

I tried to check some of them and the only thing I find where the missing nodes, that is what I thought that was the problem.


texts: How to find cumulative length at points on a line?

texts:
I am working on a Modelbuilder project to teach myself how to use it. I do not yet have much knowledge in Python. I am currently a little stuck on the project I am working on.

I have a line made from a .csv file with lat/long values. The line retains the attribute data of the original points.

I am now trying to measure the length of the line at each point, starting from the beginning of the line where the value is zero. In other words, I am looking for the cumulative distance from the start of the line to the end, with values at each point along the line.


texts: Openstreetmap - retrieve nearby big cities for a given location

texts:
I was wondering how I can retrieve places near a given location?

For example, if I pass the osm API the location "Paris" is it possible to retrieve all the main biggest cities around Paris within a given distance? (For example: "osm give me please all the biggest cities near Paris, in a range of 50 kilometers").

Is it possible?

Any idea?

thanks.


texts: How to create text labels with background boxes?

texts:
I'm trying to create some custom labels for point features in a QGIS point vector layer, and I want to write a python plugin to do it.  We need our labels to be presented with color tinted rectangles behind them (or else!). The labels are long strings, and are created using labeling engine expressions and conditionals.

I thought the easiest way would be to let the user set up labels for the active layer the way they want, and then let my plugin can grab the label text from the features in the layer and use the text strings to create boxes. My problem is I can't find a way to access those text strings in the QGIS API. How are label strings generated by the labeling engine stored in QGIS and how can I get to them?


texts: Can GDAL Perform Vertical Transformations?

texts:
Can GDAL perform vertical transformations?

If the answer is yes, please indicate which amongst these listed pages can explain how exactly to accomplish vertical transformations?

If there is some other page discussing how exactly to accomplish vertical transformations, please indicate URL.

This question arises after reading Jason's question Vertical Datum Conversion Utilities from a couple of years ago and the partial answer given.


texts: Downloading georeferenced maps from geoportals

texts:
Is it possible to download georeferenced maps from Havering London Borough and Adur & Worthing Councils portals, and if it is, then how?

I tried to use http response (see picture below), however i just got the jpeg image, that is not georeferenced and with the poor resolution.


texts:
The SBAS standard, RTCA-MOPS, contains lots of unclear information. I am trying to implement SBAS corrections by the help of EGNOS SDK which is an open source software by ESA.

When I apply long corrections, results get worse. Fast and ionospheric corrections make it better or worse.

Can you suggest any more detailed documents and source code?


texts:
I have two layers that I'm trying to conditionally join. One layer JUNCTIONS has junction IDs (NAME) and flooding volumes (Flood_cf) for each junction. I'm interested in identifying and symbolizing only the junctions with a flooding volume that exceeds a threshold (400 cubic feet).

The second layer CONDUITS contains pipe features with starting and ending nodes for each pipe. I would like to use the JUNCTIONS layer to identify the pipes that are surcharging (i.e. flooding). These are the pipes with starting nodes (INLETNODE) that match the NAMES from the JUNCTIONS layer which exceed the set threshold. Then, I would like to add a field to the CONDUITS layer in order to tag these pipes for symbology.

I tried doing a subquery (i.e. SELECT FROM (SELECT FROM) ) but I couldn't get the syntax to work. Alternatively, I could create a new feature class for the nodes that exceed the threshold, join the new layer to the CONDUITS layer, and export the results to a new layer. However, this process is tedious and generates a lot of extra and unnecessary files and fields. Is there a faster and a cleaner way to do this? I'll eventually have to create a model for this process because I need to run it on hundreds of files.


texts:
What I typed and what I got:

>>> print 'hello world'
Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
AttributeError: 'geoprocessing sys.stdout object' object has no attribute 'write'

Any suggestions? I have tried restarting to no effect, should I just go straight to re-installing? I am using ArcCatalog 10.4.1


texts: Using Intersect to combine multiple layers: correct and zero results

texts:
I have four rasters with cells of 5x5 arc minutes. I transformed the rasters to vectors/ shapefiles. The cells overlap, so I used intersect to combine the layers. It dawns on me that this was not very efficient, and clip might have done a better/ faster job. Nonetheless, it should have yielded the correct results.

However, for some reason I receive most awkward results: For three of the four layers, seemingly every FID is double. When looking at each of these pairs, one sees that the fourth layer displays a zero one row and the correct value for the other row. I know that the value > zero is correct, because this is the value I get when I check the original shapefile.

What has gone wrong and how can I fix it?


texts: how to extract OSM data to as many layers as possible in DXF( or multiple DXF file) as path for illustrator use?

texts:
So i am just touching QGIS, mainly to explorer way to use map data for my animation. Here's my first try to make it clear to understand what I am aiming for:
https://www.youtube.com/watch?v=jhC-C5f-K68

Right now I can only output 2 layers , "lines any" that include all road, and "multilinestrings any" that include only the main road (seems to be).

I understand the map got a lot of tags and I wish to export them individually, say the highway, small road, bike lane, waterway, building, nature, etc.

googling quite a lot but to no avail.


texts: How to create a line with all attribute information of the original GPX points?

texts:
I am trying to make a shapefile line contain all the data that the GPX file has.

So i load the gpx file into Qgis, (prompts me for tracks, waypoints, trackpoints, routes, etc)

When i choose "Tracks" the attribute table is empty, just one feature, yet when i choose "Track Points" all the data is in the attribute table ( elevation, time, etc )

I would like to be able to make a shapefile line that contains the trackpoint data.  I tried selecting both the track and trackpoint gpx layer together, and then  'save as"  a shapefile.  when i do, it only seems to save one, either just the track, or just the trackpoints as a shapefile.

The conversion to shapefile goes well when saving each gpx layer separately, ... so lets say i have Shapefile.line ( no attributes ) and Shapefile.point ( full attribute table ) .  How can i import the data from .point  into . line  ?

I have tried to merge them, but had no success.

Using Qgis Dufour, on Ubuntu 14.04lts, Gps Garmin Monterra

I have tried the plugin mentioned and have gotten some results, but the resulting line shapefile did not contain any of the information in the attribute table.

Is there not a way to just copy the attribute table from one shapefile, and paste it into another?  I have found a way to copy the table to my clipboard, but have not found a way to paste it to another.

( im working on a road, i want my shapefile to appear as a line but contain attribute table information from trackpoint shapefile, both shapefiles were generated by the same GPX file )

Here is a screenshot of both attribute tables, i would like C Mainline to contain all the data from Mainline points


texts: Creating a folder and setting as workspace using ArcPy?

texts:
I am trying to create a folder and setting that as the workspace for the rest of the script. However I am not sure why it isn't working.

When I hard wire in a location, the env.workspace can be set to that variable. See script below:

import arcpy, datetime

today = datetime.date.today()
date = today.strftime("%Y%m%d_%B_%Y")

folder = r"path\to\folder\\"
arcpy.env.workspace = folder

However when I create a folder using CreateFolder tool, and set that variable as the env.workspace, I keep getting an error. The script doesn't run.

import arcpy, datetime

today = datetime.date.today()
date = today.strftime("%Y%m%d_%B_%Y")

location = r"path\to\location\\"
folder = arcpy.CreateFolder_management(location, date)

arcpy.env.workspace = folder


texts: Zonal statistics using two raster, one with zones and the other with data using Python

texts:
I am trying to calculate the mean of a raster using zone data from another raster. The zone data was created using the Whitebox clump tool where each group of pixels have a unique ID number. There are ~800,000 unique clumps and the rasters are fairly big (27700, 31511) so I am avoiding converting the clumps to vector format because this causes memory issues. I would like the output of this analysis to be another raster where the clump IDs in the original raster are replaced with the mean value of that clump or a table containing the clump IDs and mean value of that clump. This process is the same as ArcMaps Zonal Statistics but I would like to use python and open source packages. I tried this using the code below that I wrote but it is way too slow for the size of data I'm working with.
clumps = raster of clumps
IDs = np.unique(clumps)
values = raster of values
means = clumps
for id in IDs:
    mask = clump
    mask = np.where(mask == id, 1, 0)
    mean_value = mask*values
    mean_value[mean == 0] = np.nan
    mean_value = np.nanmean(mean_value )
    mean_value = mean_value.astype(np.int64)
    means[means == id] = mean_value


texts:
I know that arcpy.Statistics_analysis is able to calculate the sum of Float fields. Now I have a gdb named yyy.gdb that contains a feature class named abc like this:

where:

abc contains 3 fields A, B and field Shape_Area;
field A contains A1, A2, A3...;
field B contains B1, B2, B3....

I want the sum of Shape_Area when A=A1 and B=B1, A=A1 and B=B2... i.e. an output table like this:

import arcpy, os, time

print 'start：' + str(time.ctime())
ws = r"F:\data\yyy.gdb"
arcpy.env.workspace = ws
arcpy.Statistics_analysis('abc', "F:/data/yyy.gdb/abc_1", [["Shape_Area", "SUM"]])
print 'abc is ok'

print 'ok：' + str(time.ctime())


texts: Problem with NotEqualTo filters in SLD and CSS

texts:
I am having a bit of a problem with using a '<>' attribute filter in CSS for GeoServer.  I have parred my CSS down to the absolute minimum for a buildings layer as follows:

[aeroway <> 'terminal']{
  fill: #aa1111;
}

The auto-generated SLD is as follows:

<?xml version="1.0" encoding="UTF-8"?><sld:StyledLayerDescriptor xmlns="http://www.opengis.net/sld" xmlns:sld="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc" xmlns:gml="http://www.opengis.net/gml" version="1.0.0">
  <sld:NamedLayer>
    <sld:Name>Default Styler</sld:Name>
    <sld:UserStyle>
      <sld:Name>Default Styler</sld:Name>
      <sld:FeatureTypeStyle>
        <sld:Rule>
          <ogc:Filter>
            <ogc:PropertyIsNotEqualTo>
              <ogc:PropertyName>aeroway</ogc:PropertyName>
              <ogc:Literal>terminal</ogc:Literal>
            </ogc:PropertyIsNotEqualTo>
          </ogc:Filter>
          <sld:PolygonSymbolizer>
            <sld:Fill>
              <sld:CssParameter name="fill">#aa1111</sld:CssParameter>
            </sld:Fill>
          </sld:PolygonSymbolizer>
        </sld:Rule>
        <sld:VendorOption name="ruleEvaluation">first</sld:VendorOption>
      </sld:FeatureTypeStyle>
    </sld:UserStyle>
  </sld:NamedLayer>
</sld:StyledLayerDescriptor>

According to the documentation both <> for CSS and PropertyIsNotEqualTo for SLD are perfectly valid.  The layer contains data and shows fine if I remove the filter, so the problem does not lie with the data.

Interestingly, PostGIS does not return data whose aeroway data is Null if you send it an SQL as follows:

select * from buildings where aeroway <> 'terminal'

In PostGIS, the above query returns only buildings whose aeroway attribute is both NOT NULL and NOT 'terminal'.  That was a surprise because I was expecting to get absolutely everything that did not have a value of 'terminal' (including the null values).  So, potentially the problem might lie with my expectations of PostGIS's approach to SQL.

Is this normal behaviour and can anybody suggest a work-around please (aeroway = '' does not work)?


texts:
I have a csv with approx 1.3 million places. Some of these have the same coordinates.

My question is:
Is it possible to get visible summary on the point if it has maybe 100 points in the same place. I would like it to be visible as below picture.


texts: What information is stored in a PostGIS geometry and how can I access it?

texts:
Consider the following queries:

SELECT ST_GeomFromEWKT('SRID=32632;POINT(10 20)')::varchar
UNION ALL
SELECT ST_GeomFromEWKT('SRID=4326;POINT(10 20)')::varchar;
UNION ALL
SELECT ST_AsEWKT(ST_GeomFromEWKT('SRID=32632;POINT(10 20)'))
UNION ALL
SELECT ST_AsEWKT(ST_GeomFromEWKT('SRID=4326;POINT(10 20)'))

(Note: the UNION ALL is just so you can copy and paste the whole thing easily)

Results:

0101000020787F000000000000000024400000000000003440
0101000020E610000000000000000024400000000000003440
SRID=32632;POINT(10 20)
SRID=4326;POINT(10 20)

According to the PostGIS documentation, ST_GeomFromEWKT returns a ST_Geometry object. From my first two queries it looks like the resulting objects are identical. They are not, however (see queries 3 and 4).

Now my question is, where is this additional information stored and how can I access the "complete" ST_Geometry object (so all data that are contained in it)?

Version info:

PostgreSQL: PostgreSQL 10.3 (Debian 10.3-1.pgdg90+1) on x86_64-pc-linux-gnu, compiled by gcc (Debian 6.3.0-18+deb9u1) 6.3.0 20170516, 64-bit
PostGIS: POSTGIS="2.4.3 r16312" PGSQL="100" GEOS="3.5.1-CAPI-1.9.1 r4246" SFCGAL="1.3.0" PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 2.2.2, released 2017/09/15" LIBXML="2.9.4" LIBJSON="0.12.1" LIBPROTOBUF="1.2.1" (core procs from "2.4.1 r16012" need upgrade) TOPOLOGY (topology procs from "2.4.1 r16012" need upgrade) RASTER (raster procs from "2.4.1 r16012" need upgrade) (sfcgal procs from "2.4.1 r16012" need upgrade)


texts: Extracting impervious areas (buildings, rails, pavements) from LAS dataset

texts:
How to extract impervious areas (like buildings, rails, pavements etc.) from las  dataset and receive it like a raster? I tried to do it with intensity, but there are many valids with vegetation.

I am working with ArcGis Pro 10.1. I generated and reclassified a raster of intensity with range 0-180 => value 1 (values without NODATA), then I created next raster, a nDSM without medium and high vegetation => value 1 (medium, high vegetation - NODATA).

After that, I used a raster calculator and summed this two rasters; value 2 responds like my results. I was doing it with raster cell sizes equal to 0.2m, but the nDSM with this cell size was not satisfactory (the las dataset have 12 pts/m²). For example, buildings were not filled completely, and I created this nDSM with cell size 1m. So, at the moment my results are a little bit generalized, which is not satisfactory too.


texts:
I am trying to open an image file with the extension ".fit". Does anybody know which software to use?

I obtained this image using the WASI software (Water Color Simulator) available on the IOCCG website (http://ioccg.org/resources/software/).

The contents of the image file are described by a additional ".hdr" file. The image file contains several bands, corresponding to the different parameters of the water (constituant concentration, bottom depth, etc) computed by the software.

It would be great if I could open this image with QGIS.

So far I have found this list of image viewers (https://fits.gsfc.nasa.gov/fits_viewer.html), but none of those that I tried have worked yet.


texts:
The L2A_Process script provided by sen2cor expects a full product directory in the .SAFE format. Is there any way to use it for a single tile (a subdirectory in the GRANULES folder)?

If no, is there any conceptual reason preventing it from being implemented this way? As far as I understand, sen2cor in any case processes each tile within the given product on an independent basis.

What I am trying to do is to build and maintain an archive of L2A imagery for a limited area of interest, and I do not want to spend the processor cycles on the tiles that do not overlap with the AOI.


texts:
I'm trying to save filled raster object to geo database by my specified filename.

import arcpy
from arcpy import env
from arcpy.sa import *

env.workspace = r"C:\GIS Data\Cameron_Run\Topo_NED\Projected"

# Set local variables
surfaceRaster = "raw_10m_dem"

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Execute FlowDirection
fillRaster = Fill(surfaceRaster)
arcpy.RasterToGeodatabase_conversion(fillRaster, r"C:\GIS Data\Cameron_Run\Topo_NED
\Projected\fillRasterGDB.gdb")

But this code is saving raster file by its own specified name, for this case Fill_raw_10m1 is the default name of saved raster file.

Is there any way so that I can specify the name of saved raster file?


texts: Using multi-part line features in network analysis?

texts:
Cartographically is there any reason why (pro / con) multi-part line features should not be used when developing road networks for ArcGIS Network Analyst?

For example, say I have a road which is 100 units long and it has two bridges; one located at 20-25 units and one located at 50-55 units. The road is one continuous line feature with no intersecting roads.

So is it better to have one multi-part feature (one database entry) comprised of three features (segment before bridge one, the segment between bridges one and two, and the segment after bridge two), or three feature records in the feature class – one record for each road segment and not have a multi-part feature?


texts: ArcGIS Online - Unable to upload a shapefile with over 250,000 points as hosted feature service

texts:
I have a shapefile with over 250,000 points.

I've tried to upload it as a hosted feature service but it does not populate the attribute table.

It seems to time out despite the option to "keep me always signed in" has been selected.

Anybody know how to upload the file so I can use it in a web map?


texts: R raster: align two rasters with different extent

texts:
I have two rasters, with same dimensions (nrow, ncol, and ncell) and same coordinate system. They differ however by their extent, x min/max points:

Raster 1: has xmin/max -180, 180, i.e. Americas left, Russia right
Raster 2: has 0, 360, i.e. central Africa/Europe left, West Africa/Europe right

I want to align them in R, so that both have -180, 180. I tried several functions such as resample(), projectRaster(), or alignExtent() but all seem to give the right extent, but will basically show only half of the values, i.e. it looks like half of the world disappeared in the new plot.

Which function should I use instead? What is the name of the operation I want to do?

library(raster)

r1 <- raster(nrows=360, ncols=720)
r2 <- raster(nrows=360, ncols=720, xmn=0, xmx=360)

r1[] <- 1:(360*720)
r2[] <- 1:(360*720)

## use projectRaster?
r2_in1 <- projectRaster(r2,r1)

plot(r2_in1)


texts: Why can't osm2po router find vertexIds when I use Graph.SUPPORT_EDGEFLAGS?

texts:
I am using the Osm2po.router in Java, trying to find the clazz of roads from this Get road type (clazz) in osm2po.

When I am using the router without any supports, I am able to find the vertex Ids. That is, the following code returns a path:

Graph graph = new Graph(graphFile);
Properties params = new Properties();
params.setProperty("findShortestPath", "false");
params.setProperty("ignoreRestrictions", "true");
params.setProperty("ignoreOneWays", "true");
params.setProperty("heuristicFactor", "0.0"); // 0.0 Dijkstra, 1.0 good A*

int sourceId = graph.findClosestVertexId(-37.820088f, 144.970264f);
int targetId = graph.findClosestVertexId(-35.281537f, 149.129242f);

Returns a non-negative "sourceId" and "targetId" int.

However, if I change the code to follow the above link, the same target and source nodes return null.

Graph graph = new Graph(graphFile, new Log(), Graph.SUPPORT_EDGEFLAGS, false, null);
Properties params = new Properties();
params.setProperty("findShortestPath", "false");
params.setProperty("ignoreRestrictions", "true");
params.setProperty("ignoreOneWays", "true");
params.setProperty("heuristicFactor", "0.0"); // 0.0 Dijkstra, 1.0 good A*

int sourceId = graph.findClosestVertexId(-37.820088f, 144.970264f);
int targetId = graph.findClosestVertexId(-35.281537f,   149.129242f);

Why? My understanding is that the route should be the same, the only change I've made is about what details about the roads to retain.

I've tried changing the Boolean which is the fourth parameter in the Graph constructor, but that doesn't change the result.


texts: GeoServer - PostGIS database connection issue - default password always appearing

texts:
I am using GeoServer 2.18.2 in a Docker image (Docker version 3.1.0) from my Mac (OS Big Sur 11.2.3). My problem is that when I want to add a PostGIS new store a default password always appears. If I modify it with the proper one, it does not change anything as when I apply or save the connexion the default password comes back.
I have checked with colleagues using the same Docker image also on Mac, and I seem to be the only one with this issue. I disabled all passwords managers but still this default password appears.
Any suggestion?


texts: Connection information provided was for a non-administrative user fails in task scheduler

texts:
When attempting to use the arcpy.DisconnectUser function in a script that runs in task scheduler I receive the error 'Connection information provided was for a non-administrative user'.

I am able to run the script without any hiccups outside of task scheduler using the same connection. Prior to the script throwing the error it had been working properly. The connection uses windows authentication and the SDE DB connection is set up for sys and process admin.

Could this be a SQL permissions related problem?

Below is a sample of the code:

Admin = r'Database Connections\username.sde'
start = datetime.now()
resultsFile.write("Started process at - " + str(start)+"\n")
verList =[ver.name for ver in arcpy.da.ListVersions(Admin) if ver.name.lower() != 'sde.default'and ver.name.lower() != 'gisadmin.qaqc']
print "Created list of versions less default and QaQc at " + str(datetime.now())
resultsFile.write("Created list of versions less default and QaQc at " + str(datetime.now())+"\n")
arcpy.AcceptConnections(Admin, False)
arcpy.DisconnectUser(Admin, "ALL")
print "Disconnected users at " + str(datetime.now())
resultsFile.write("Disconnected users at " + str(datetime.now())+"\n")


texts: Define Workspace for SDE Connection in Python

texts: How do we define SDE Connection for workspace in Python Scripting?


texts: Creating smaller OSM with only basic zoom levels data from larger OSM file?

texts:
From a large .osm/.o5m file (for India), I only want to extract the data that would actually show up when one is looking from zoom levels 1/0 through 8 (which is quite zoomed out.. only major cities, highways, geographical regions and coastlines). See screenshot, of a much smaller area loaded in Maperitive, seen at zoom level 8:

The extracted .osm would then be much smaller in size, and I would be able to load it on Maperitive/other for just basic-zoom-levels rendering. I'm working with India's .osm downloaded from geofabrik.de. The .osm (with author and version info dropped) currently weighs over 3GBs.. impossible to load on my end, and we don't have access to any super-ish computers. Basically I want to get rid of all data that wouldn't show up at zoom level 8 or higher anyway.
I'm working with osmconvert / osmfilter / maperitive, on Windows. Would prefer working with these simple portable apps and not get into postGIS databases etc, in the interest of making this doable for more people.

My understanding is that the data in the OSM file isn't tagged by zoom level; that there are algorithms and map styles that decide what gets displayed and what doesn't. So, it would be great if I could get a list of objects to --keep using osmfilter or similar. I don't want to do any major customization to the regular OSM rendering, and don't need to religiously stick to it either. Example, from the wiki page:

osmfilter norway.osm --keep="highway=primary =secondary waterway=river" >streets.osm

What source should I be looking up to get this --keep list?


texts: Speedup re-projecting rasters in FME Desktop!

texts:
I am trying to re-project some aerial photos (at the resolution 20cm) using Esri re-projector in FME Desktop 2017, and i it takes too much time reprojecting (10-15x slower) comparing to ArcGIS.

Is there any configuration i can do to speed it up?

Mosaicking rasters doesn't take that much, it works normal, but re-projecting takes too much time.


texts:
I can't find any simple way to clone an ol3 style.  There is no clone() method and creating a style from a style (new ol.style.Style(style)) doesn't do it.  JSON.stringify doesn't produce anything usable, either.

Is there any ol3 code I might be overlooking for cloning styles?


texts: Which interpolation method is the most suitable for DEMs?

texts: Which interpolation method is the best suitable for hill/valley areas and normal ground with undulations and rivers? When working with Aster 30mX30m DEMs, which method will give the most realistic elevations? I am trying with different methods and getting different contour distribution and shapes but don't know which will give the most accurate and realistic values.


texts: How to request dynamic Bokeh map and tile updates?

texts:
Google maps in Bokeh has become too difficult with their API keys and registration, so I'm exploring alternatives but really like the Bokeh framework, especially with periodic callbacks.

What's the right way to coerce Bokeh to update the basemap tile better? I assume map tiles are supposed to be updated automatically, and they are, sort of, but it's not working as well as expected.

Here's a simple code clip to convert (lat/lon) pairs into Northing and Easting values to display a region of a basemap.

Zooming and panning works but sometimes leaves large blank tiles that never fill-in.

How can this be addressed?

from bokeh.plotting import figure, show, output_file, curdoc
from bokeh.models.widgets import Slider, Div
from bokeh.layouts import column, layout
from bokeh.tile_providers import get_provider, Vendors
from pyproj import Proj, transform

def LatLon_to_EN(lat, lon):
    try:
      easting, northing = transform( Proj('epsg:4326'), Proj('epsg:3857'), lat, lon) # from WGS-84 to Web Mercator Easting/Northing
      return easting, northing
    except:
      return None, None

tile_provider = get_provider(Vendors.STAMEN_TERRAIN_RETINA)
#tile_provider = get_provider(Vendors.STAMEN_TERRAIN)

# SouthWest and NorthEast corners
SW_corner=(29.191391, -81.049014) # center of campus foot bridge
NE_corner=(29.239598, -80.949180) # out in the ocean

# Convert to mercator Easting and Northing coordinates
SW_corner_EN = LatLon_to_EN( SW_corner[0] , SW_corner[1] )
NE_corner_EN = LatLon_to_EN( NE_corner[0] , NE_corner[1] )

# x_axis_type, y_axis_type may be: "linear", "log", "datetime", "mercator"
plot = figure(x_range=(SW_corner_EN[0], NE_corner_EN[0]), y_range=(SW_corner_EN[1], NE_corner_EN[1]), x_axis_type="mercator", y_axis_type="mercator")
plot.add_tile(tile_provider)

curdoc().add_root(plot)
curdoc().title = "Bokeh example for basemap with non-Google tile providers"


texts: CF compliance - why are certain variables are required if WKT defined?

texts:
I'm trying to make a NetCDF CF-compliant, where the crs variable uses a WKT string. I'm using this compliance checker, testing against CF 1.7. I'm getting errors that certain attributes should be defined, even though the information requested is already in the WKT string.
This seems like a bug, but in example 5.12 in this section of the CF 1.8 conventions, there are attributes defined that are also in the WKT string, e.g. false easting/northing, lon/lat of projection meridian/origin. Is there a reason for this? Or, does this seem like something to bring up with contributors to the CF convention, or developers of the compliance checker?


texts: How to change dialog box with jQuery Mobile on OpenLayers

texts:
I have a problem with dialog box made with jQuery Mobile on OpenLayers. I want to make a popup box instead of a dialog box, but when i click the point the info for the point appears at different window instead of appearing on the map. My code is:

js:

$('#popup').live('pageshow',function(event, ui){
        var li = "";
        for(var attr in selectedFeature.attributes){
            li += "

" + attr + "

"
            + selectedFeature.attributes[attr] + "

";
        }
          //alert("test");

        $("ul#details-list").empty().append(li).listview("refresh");
    });


texts: Matching up fields in different tables using ArcGIS Desktop?

texts:
I work for a University, and we have some utility/site people go out with their GPS and shoot points pretty often when there are new projects, construction, etc. The data dictionary that they have created has some field names that are longer than 10 characters, so when they export to shapefiles, the fields are being cut off to only contain the allowed limit. Once they bring in the shapefile into our existing base map, they copy and then paste the new items into the existing sde feature class, the data that is contained in the truncated fields won't copy. So, then they have to go in and redo some of the fields that didn't copy over.

I know that fields have to match between tables in order to copy/paste data, but is there some kind of script, code, or work around that can match up these fields?

I was hoping that there is some kind of field name match. for example, something that i could use to say that create_dat from the gps export file would equal create_date in the existing feature class.

I have attached an image of both tables and what data will and won't copy over.


texts: Clipping raster layer with *raster* mask layer in QGIS

texts:
Given two input layers

The clip layer - an arbitrary layer (though in my case it's a float 1 band geotiff)
The mask layer - a paletted raster in which every pixel is either NO DATA or a single "color" value, same for every pixel.

I want to create an output raster whose overall extent matches the extent of the mask layer and each pixel is either NO DATA if the corresponding pixel of the mask layer is NO DATA or the value of the corresponding pixel from the clip layer if the corresponding pixel of the mask layer is not NO DATA.
I know I can mask the clip layer with a polygon, but the mask layer is a raster. Polygonizing the raster layer is possible in principle, but would result in a polygon with a lot of very short edges.
Is there a direct way to do this in QGIS?
I suppose there's probably a path that involves creating a new layer and using the raster calculator, but I'm hoping there's something more direct. Further complicating the matter is that the two rasters might not have the same resolution and origin.
I believe this is equivalent to Extract by Mask in ArcMap.


texts: Adding a SQL Server layer always uses one spatial column even if I select a different one

texts:
I have a SQL Server table with two geometry spatial columns (A, B). I have created two rows in geometry_columns table for both of these.

f_table_catalog f_table_schema  f_table_name    f_geometry_column   coord_dimension srid    geometry_type
neil            dbo             Temp_CSVLoad    [A]                 2               27700   POINT
neil            dbo             Temp_CSVLoad    [B]                 2               27700   POINT

If I go into QGIS and add a SQL Server layer, it finds and shows me these two layers as choices correctly, but no matter which one I select it only ever plots A, whether I select to add the row linked to column A or column B. If I look at the properties, it is definitely showing the correct column as part of the properties.

If I null column A and only add a layer for B, I get nothing.
If I null column B and only add a layer for B, I still get A.

Am I doing something obviously wrong?


texts: Comparing a polygon layer within the boundary of a polygon

texts:
I'm pretty new to ArcPy and Python.

In our database, we have a layer containing parcels (polygons) with a neighborhood code and a polygon layer that contains the neighborhood boundary. Within that parcel layer, we have some errors where the neighborhood that the parcel should be in is not. I'm trying to attempt a way to compare each parcel and see if it falls within the neighborhood boundary that it's assigned. An example would be if this parcel is assigned to neighborhood 1, and falls within the boundary of neighborhood 1. If the neighborhood boundary does not match up, I would like that selected so it can be analyzed further. This is what I have so far as my ArcPy code.

I got a runtime error when it could not find the column, but I'm also unsure if a search cursor is what I need. I tried a join using the standard GIS tool, but it didn't return what I needed. I used column 29 because that is the one that contains the neighborhood code, same for column 3.

# Name: ExtractFeaturesByLocationAndAttribute.py
# Description: Extract features to a new feature class based on a spatial
# relationships to another layer, and an attribute query

# Import system modules
import arcpy

# Set the workspace
arcpy.env.workspace = r'C:\ArcMap\Apt_Parcels'

# Select all apartments which overlap the Neighborhood polygon
Apartments_SelectAll = arcpy.SelectLayerByAttribute_management('Parcels_owner', 'NEW_SELECTION',
                                                          "[AcctType] = 'APARTMENT'")

arcpy.CopyFeatures_management(Apartments_SelectAll, 'Apartments_Outside_of_Boundary')

Apartments_neighborhoods = arcpy.SelectLayerByAttribute_management('apartment_neighborhood', 'NEW_SELECTION',
                                                          "[Id] = '0'")

arcpy.CopyFeatures_management(Apartments_neighborhoods, 'Apartment_Neighborhoods')

# Within selected features, further select only those apartments that fall outside
# of the neighborhoods they are supposed to be in

fc = r'N:\My Documents\ArcGIS\Default.gdb\Apartments_Outside_of_Boundary'
fc1 = r'N:\My Documents\ArcGIS\Default.gdb\Apartment_Neighborhoods'

with arcpy.da.SearchCursor(fc,['apt_nbhd']) as cursor1:
    for x in cursor1:
        print x[29]

with arcpy.da.SearchCursor(fc1, ['Nbhd']) as cursor2:
    for y in cursor2:
        print y[3]

arcpy.SelectLayerByAttribute_management(Apartments_Outside_of_Boundary, 'SUBSET_SELECTION',
                                        'cursor1' == 'cursor2')

# Write the selected features to a new featureclass
arcpy.CopyFeatures_management(Apartments_SelectAll, 'Apartments_Outside_of_Boundary_final')


texts: What GIS tools and packages are needed to find out nearby users in android application?

texts:
I'm building a location based android application with following requirements:

Need to find out nearby users(i.e users using the same application)
Need to find out users in the state
Need to find out users in the country

All the processing will be done on the server side. The android application will communicate with the server through web services. What tools and packages do I require to achieve these requirements? I'm new to GIS.
My preferences are java and free software.


texts:
I am using Overpass Turbo, with extended syntax to query places in an area, usually a country, using the following syntax

{{geocodeArea:CH}}->.searchArea;
(
   (node["place"]["name"~"^.+$"]["place"="country"](area.searchArea);)
);
out body;

Now, I would like to use this query against the open API directly.

I did figure out that it is possible but that the API is limited in regard to the extended syntax (the curly braces, it seems). So when I export the query from overpass turbo to text (Export -> Query -> Text), I get the following output:

area(3600051701)->.searchArea;
(
  (node["place"]["name"~"^.+$"]["place"="country"](area.searchArea);)
);
out body;

I can then POST this request to the API and it works.

Now my question:

I see that the text export converted the area code "CH" to the string/id 3600051701.

Is it possible to have a lookup for these IDs? I might need to query many different countries and it would be convenient to have a lookup table somewhere.

I am new to OSM/Overpass.


texts: qgis missing menu bar ver 1.8.0-lisboa on ubuntu linux

texts:
QGIS emits the following warnings when starting, and fails to
display the main menu bar.

$ qgis --nocustomization

Warning: loading of qgis translation failed [/usr/share/qgis/i18n//qgis_en_US]

Warning: loading of qt translation failed [/usr/share/qt4/translations/qt_en_US]

Warning: void DBusMenuExporterPrivate::addAction(QAction*, int): Already trackin
g action "Project Properties..." under id 125

Warning: void DBusMenuExporterPrivate::fillLayoutItem(DBusMenuLayoutItem*, QMenu
*, int, int, const QStringList&): No id for action

Warning: void DBusMenuExporterPrivate::fillLayoutItem(DBusMenuLayoutItem*, QMenu
*, int, int, const QStringList&): No id for action

The "//" in the path "/usr/share/qgis/i18n//qgis_en_US" looks suspect,
and there is no "qgis_en_US" among the 47 files in the i18n directory.

Help and suggestions are welcome.

Thanks.


texts: Filling dictionary with list of row attributes using arcpy.SearchCursor

texts:
I have a dBase database file and I'm viewing and editing the table in ArcGIS Pro. I want to use the arcpy.da.SearchCursor function to fill a dictionary and create a key-value pair by iterating  through multiple fields and row values in my table. I will eventually use this dictionary to compare key-value pairs from another dictionary and see if items in both lists match.

My first step is to understand and create dictionaries with ArcPy and Python. My goal is to create and fill a dictionary where the key is a shapefile and the value is a list of field names thought to be within the shapefile. The code below creates a key and value pair where the key is a shapefile and the value is a list of row attributes from five fields.

import arcpy

arcpy.env.workspace = r"C:\TestFolder"

dbf = "APExistOnly_Edit.dbf"

myfield = "FI_PATH"

FieldNameDict = {}

rows = arcpy.SearchCursor(dbf)
for row in rows:

    FieldNameDict[row.FI_PATH] = [row.AP_ZIP,row.AP_STATE,row.AP_CITY,row.AP_FULLADD,row.AP_BUILDIN]

print(FieldNameDict)

With this code I created a dictionary called FieldNameDict and placed the key as field FI_PATH and assigned the row values from each field to a list. Below is the output from the code. When viewing the output, one can see that the lists contain empty elements with no text. This is because some of the fields contain rows with no value.

Currently, I have messy code and logic and I want to see if someone can help with better code. I want to create the dictionary by using arcpy.SearchCursor and iterate through multiple fields and rows. I want to create lists, associated with a key, that do not include empty row elements '  ' . I would like to exclude rows with no text when iterating through fields with the search cursor object.


texts:
Just trying to run the following code already given in the pyqgis cookbook to create graduated colors symbology (as in example with 2 classes)

I define the source info for the shapefile but when I run execute the code within my test-plugin, it just adds the layer without classifiying. Just adds as usual, all polylines appear same.

Any solution or working code/method?

link to code in cookbook
http://www.qgis.org/pyqgis-cookbook/vector.html#graduated-symbol-renderer

from qgis.core import  (QgsVectorLayer, QgsMapLayerRegistry, QgsGraduatedSymbolRendererV2, QgsSymbolV2,QgsRendererRangeV2)

myVectorLayer = QgsVectorLayer('C:/my_test_shape_file.shp', 'test_shp_file', 'ogr')

myTargetField = 'target_field'
myRangeList = []
myOpacity = 1

# Make our first symbol and range...
myMin = 0.0
myMax = 50.0
myLabel = 'Group 1'
myColour = QtGui.QColor('#ffee00')
mySymbol1 = QgsSymbolV2.defaultSymbol(myVectorLayer.geometryType())
mySymbol1.setColor(myColour)
mySymbol1.setAlpha(myOpacity)
myRange1 = QgsRendererRangeV2(myMin, myMax, mySymbol1, myLabel)
myRangeList.append(myRange1)

#now make another symbol and range...
myMin = 50.1
myMax = 100
myLabel = 'Group 2'
myColour = QtGui.QColor('#00eeff')
mySymbol2 = QgsSymbolV2.defaultSymbol(myVectorLayer.geometryType())
mySymbol2.setColor(myColour)
mySymbol2.setAlpha(myOpacity)
myRange2 = QgsRendererRangeV2( myMin, myMax, mySymbol2, myLabel)
myRangeList.append(myRange2)

myRenderer = QgsGraduatedSymbolRendererV2('', myRangeList)
myRenderer.setMode(QgsGraduatedSymbolRendererV2.EqualInterval)
myRenderer.setClassAttribute(myTargetField)
myVectorLayer.setRendererV2(myRenderer)
QgsMapLayerRegistry.instance().addMapLayer(myVectorLayer)

@Kelly's code works perfectly.

However I noticed that both the code (one I have typed in my first message and yours in the 2nd message) DO NOT work in QGIS v1.7.3 but WORKS in QGIS v1.8.0. I think this was a bug(?) that already resolved in v1.8.0

And one more question;

Do you have any sample code for (natural breaks, Jenks) classification of the "numeric_attribute_field" based on given number of classes (i.e. no of classes will be a variable in the code, say "n", and I will pass that from the plugin GUI i.e. n = spinBox.value() )


texts: enlarge bounding box from a smaller one in GDAL

texts: I have a small raster tiff and I want copy its data in a new raster with a bigger bbox. So is it possible, via a GDAL command, to enlarge a bounding box maintaining same pixel locations from the original raster? In output I expect something like a padded image filled of nodata value around the original bbox.


texts: Arcpy fail to loop between files in different directories

texts:
I want to perform a loop of intersections between different files in different directories. Here is my code:

Intersects = ['Name1','Name2','Name3','Name4','Name5'] #Names of my files
arcpy.env.workspace = "C:\Mydirectory" #Directory of the files above
lands =  "C:\Anotherdirectory\Shapefile_To_Intersect.shp" #Name of the second shapefile to intersect, in another directory
#    Where I want my output to be placed, this folder is already created and empty
out = "C:\Users\Laptop\Documents\ArcGIS\Default.gdb\Intersects"
for files in Intersects:
    Intersects_Names = ['Name1','Name2','Name3','Name4','Name5']
    arcpy.Intersect_analysis([files, lands], os.path.join(out,files))

But this way I get the error:

ERROR 000732: Input Features: Dataset Name1 #;C:\Anotherdirectory\Shapefile_to_Intersect.shp # does not exist or is not supported.

What is the issue in this case? (I checked the file name and directory are correct).


texts: GetFeatureInfo on TMS layer with OpenLayers

texts:
I'm using OpenLayers.Layer.TMS to display external layers that have a different projection than EPSG3857/EPSG900913.
I followed this tutorial.

Now I can't use OpenLayers.Control.WMSGetFeatureInfo or OpenLayers.Control.WMTSGetFeatureInfo controls so how to handle getFeatureInfo request with my TMS layers?


texts: Near real-time sources of lightning strike data

texts:
I am wondering if anyone has any suggestions on sources of near real-time lightning strike data available for download?  For example, Intellicast.com has a map of strikes, but makes no mention of the source of the data they used.

Key attributes would be date/time and lat/long (any additional is just bonus), and strikes that occurred as much as 24 hours previously are logged and available.  It could be in a text/table form (which can be converted to a spatial table), but if it is in a GIS ready format, that would be great too.

I know of a couple of paid sources:

Vaisala
WeatherBank
USPLN, NAPLN, and GLN. (Update: Apparently Intellicast gets its data from WSI Corp, which runs the USPLN).

Other sources I have stumbled across:

blitzortung.org- I haven't figured out how to get the source strike data from the site
According to an 2007 blog entry, GuiWeather.com used to supply weather data to Google Earth, but the site is no longer active.

Is anyone aware of any other paid/free sources of lightning strike data?

UPDATE: Someone posted & deleted a good website for many resources of lightning data.  The NWA Remote Sensing Committee.  It provides a good overview of research and data providers, too.


texts: UpdateCursor is updating empty list correctly but NOT updating data in feature class

texts:
I have looked at other posts and am not able to solve my issue. I'm working with Python 2.7 (must use this ver for work).
I am trying to reformat labels in an attribute column ("ID") and update that column. My current version of the code runs without errors AND when I print the 'labels' list the reformatted edits show. However they do not get populated to my feature class column. I am missing something and/or have a syntax issue.
    mxd = arcpy.mapping.MapDocument(r'D:\JJ_Development\Projects\Basemap\Transportation\Labels.mxd')
fc = arcpy.mapping.ListLayers(mxd, 'RoadCore')[0]

fields = ['ID', 'ADMIN_ORG']
labels = []

def newlabel():
    #newfields = ['ID', 'ADMIN_ORG', 'MAP_ID']
    with arcpy.da.UpdateCursor(fc, fields) as cursor:
        for row in cursor:
            if row[1].startswith ('04'):
                labels.append('{}'.format(row[0].replace("-ADMIN","").replace('-A-ADMIN','').replace('-PARKING','')[2:].lstrip('0').replace('.', '-').lstrip(ascii_letters)))
            cursor.updateRow(row)#Many SO posts show this syntax. Populates the list but does not update the feature class data.

            #cursor.updateRow([row])#per SE post
            #https://gis.stackexchange.com/questions/142965/arcpy-da-updatecursor-error-sequence-size-must-match-size-of-row
            #still errors with TypeError: sequence size mush match size of the row
        #labels.append{row[3]}#Doesn't work
    print(labels)

if __name__ == '__main__':
    newlabel()


texts:
I have done this model calculating least cost path

I'm using "Test input" as set of points with ID as start and also end points. When least cost path is calculated(Cost path tool), it can't use the destination point as start point. Example of correct combination of Start(S)- Destination point(D): S1-D2;S1-D3; S1-D4; S2-D1; S2-D3; S2-D4; S3-D1 etc.
I don't know how to write script for Calculate value (I don't know python).
It should work like this:

IF input_point ID (for destination, in submodel) = input_point ID (in main model)
then use input_point ID+1
else use input_point ID.


texts: Postgis order ST_INTERSECTION(polygon, linestring) distance query results by the intersections order

texts:
Given a query that calculates distances of ST_INTERSECTION of a linestring with many polygons (stored in the polygon_column in the PolygonTable):

SELECT
  SUM(
    ST_LENGTH(
      ST_INTERSECTION(
        po.polygon_column,ST_GeomFromText('LINESTRING(-122 48,-122 38,-122 37,...)', 4326))::geography))
  AS route_length, po.polygon_name
FROM PolygonTable po
GROUP by po.polygon_name;

The results are returned ordered by the PolygonTable polygon_name. How can I order the results in the order the line intersects them?

That is, if I have polygons A,B,C,D, and the line string goes through them in the order of B,A,C,D, I'd like my results to be

  route_length     | polygon_name
------------------+----------
 123.391961613962 |        B
 321.511826869426 |        A
  123.69306414636 |        C
  321.69306414636 |        D


texts: I need to convert a few files which are projected in WGS84 into netCDF. NetCDF requires lat/long. I am a bit confused on how the process is handled by ArcMap especially after looking at the online online Help for Arc10.2. Should I convert my raster to lat long projection before converting the file to netcdf or is the lat long conversion is automatically carried out when using the Raster to netCDF tool?


texts:
i found in this post:
CartoCSS layer visibility
a solution to set a layer not visible for some zoom levels.

But in
https://github.com/gravitystorm/openstreetmap-carto
it seems that a minzoom parameter can be set (according to the project.mml file) :

   {
      "name": "coast-poly",
      "srs-name": "900913",
      "geometry": "polygon",
      "class": "",
      "id": "coast-poly",
      "srs": "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0.0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs +over",
      "Datasource": {
        "type": "shape",
        "file": "data/land-polygons-split-3857/land_polygons.shp"
      },
      "extent": [
        -180,
        -85.05112877980659,
        180,
        85.05112877980659
      ],
      **"properties": {
        "minzoom": 10
      },**
      "advanced": {}
    },

But in Tilemill i can't get it working (lines are drawn in red).

What would be the best/fastest solution to manage minzoom & maxzoom for layers in TileMill?

My TileMill output could be Mbtiles or Mapnik files. Does this output strategy have an impact on the solution?


texts:
I do not seem to be able to paste an image from mac i photo to QGIS print composer.
When i click on add image I get a blank box with no photo.

Also when I click to add a label I do not seem to be able to change the text in the label box.


texts: Open only part of a GeoPackage layer with a certain geometry type in GeoPandas

texts:
I have a GeoPackage (743 MB, zipped)  that has one layer, but two geometry types inside that layer (Point and Polygon). When I open the file in QGIS 3.10, QGIS presents a dialogue to me that lets me choose, which geometry type to add to the project.

How can I achieve the same with GeoPandas? I don't see a way to specify the geometry type:
gpd.read_file('inspire_schlaege_2019.gpkg', layer='inspire_schlaege_2019')


texts:
I am evaluating QGIS with a view to visualize and manage geological exploration data. The information is currently stored in different local coordinate systems that I am looking to overlay with satellite imagery from Google and Bing. I figured out how to load the data I want into QGIS, but I am struggling to find a way to convert them into a common coordinate system?

I presume I should use geographical coordinate system like WGS84. Is anyone can explain what tools are available for conversion from a local projected coordinate system (I believe it is a Mercator based projection) to latitude/longitude using few reference points where I have coordinates in both systems?

I tried to follow suggestion from here:
How to georeference a vector layer with control points?
and downloaded Affine Transformation Plugin. But I cannot figure out how to get the input parameters (scale, translation and rotation) without mastering coordinate projections theory. I have not tried v.transform that was also suggested as I could not see how to make it work for geographical coordinates.


texts:
I am trying to upload a table (shape) in Google Earth Engine. I have never had problems when I tried to upload a table but this time I am having the next problem:

  Error: Encountered 1 bad features during table ingestion.

Does anyone know why is the problem?

I tried to simplify the shape (using less polygons) and does not give me the error but has been carrying the layer for 5 days.


texts: User defined Coordinate Reference System (CRS) in QGIS is not located where expected?

texts:
I have received a shapefile with a user defined projection (i.e., +proj=krovak +lat_0=-32.923433 +lon_0=116.47485 +alpha=111.423433 +k=1.00004 +x_0=44328.529781 +y_0=1345943.211845 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs).

The site is located in Australia but after automatic re-projection into GDA94 (the default CRS) it places the site in northern China.

I have viewed this particular shape file on the client's GIS system (ArcGIS) and it seems to locate fine.

Any suggestions as to what I need to do or am I missing something?

I have attempted to manually change the CRS (did not work) and I also saved the layer to a new SHP file and manually chose the CRS but this also did not work.


texts:
I currently work in the offshore oil and gas industry as an offshore surveyor and I have an undergraduate degree in survey and mapping. I enjoy my job and it is well paid but I want to learn more of GIS that the little I studied at university.

I have been accepted on to a MSc GIS course (free in my country via distance learning). The course seems to be ESRI software only and maybe not cutting edge. I am looking for advice on either taking the course or through self study concentrating on python, databases (PostgreSQL), QGIS and visualization/design.

My end goal would to have the skills and knowledge to change jobs in the future if I decide the sea life is no longer for me.

I would be very interested to hear from everyone who is working within the GIS industry. Thanks for your input.


texts: Finding field with maximum value and field name using ArcPy cursor?

texts:
I was wondering if an equivalent of summary statistics exists for finding such statistics across fields rather than records in a table. For example, of several fields, I wanted to find which one has the highest value for a particular record, and its name. Here, I found a cursor for doing so, but keep running into a problem.

I made few alterations to that cursor, apart from specifying the fields I want to run it on (I have ten such fields, so n = 10) and changing the format and names of some of them where necessary. Below is my version of this cursor:

import arcpy

# input data and fields to find maximum
table = '[the path to my file]'
fields = ["B0R001", "B0R002", "B0R003", "B0R004", "B0R005", "B0R006",   "B0R007", "B0R008", "B0R009", "B0R0010"]
n = len(fields)

# Add fields to input table to store maximum and field name
maxfield = "HIGHEST_M"
maxname = "HGHSTMNAME"
arcpy.AddField_management(table, maxfield, "LONG")
arcpy.AddField_management(table, maxname, "TEXT")

fields2 = fields[:] # shallow copy
fields2.extend([maxfield, maxname])

with arcpy.da.UpdateCursor(table, fields2) as cursor:
for row in cursor:

    # Look at the first n values
    check = row[:n]
    maxval = max(check)

    # Get the index position of the maxval and use it slice into field list
    # which gives us the field name
    # only update the last two rows
    row[-2:] = maxval, fields[check.index(maxval)]

    cursor.updateRow(row)

print "Done!"

The script keeps throwing an error, however. The message is as follows:

  Traceback (most recent call last):   File "[Path to my File]", line
  18, in
      for row in cursor: RuntimeError: A column was specified that does not exist.

Am I right that the "column" refers to the column in which the cursor should start checking maximum values? My 10 columns do not begin with column one. Where in the script should I have specified the column in which the cursor starts looking for maximum values? If I want the cursor to search, for each record, for maximum values only in, say, columns 20 through 29, inclusive, what should I do differently? Sorry about this, I haven't written cursors in a while.


texts: Creating locationtech Geometry using org.geotools.geojson.geom.GeometryJSON

texts:
I am trying to create org.locationtech.jts.geom.Geometry from a json string using org.geotools.geojson.geom.GeometryJSON.
import org.geotools.geojson.geom.GeometryJSON;
import org.locationtech.jts.geom.Geometry;

public class Test {

    public static void main(String[] args) throws Exception {
        GeometryJSON geomJSON = new GeometryJSON();
        Geometry geometry = geomJSON.read("{\"type\":\"Point\",\"coordinates\":{\"longitude\":65.0390625,\"latitude\":62.91523303947614,\"additionalElements\":[]}}");
        System.out.println(geometry);
        geometry = geomJSON.read("{\"type\":\"LineString\",\"coordinates\":[{\"longitude\":80.5078125,\"latitude\":67.06743335108298,\"altitude\":1.0,\"additionalElements\":[]},{\"longitude\":81.2109375,\"latitude\":65.2198939361321,\"altitude\":1.0,\"additionalElements\":[]},{\"longitude\":74.1796875,\"latitude\":63.39152174400882,\"altitude\":1.0,\"additionalElements\":[]}]}");
        System.out.println(geometry);
    }
}

Point is read fine, but LineString is not read and gives error:
Exception in thread "main" java.lang.IllegalArgumentException: Invalid number of points in LineString (found 1 - must be 0 or >= 2)

String is obtained by deserializing a json input as:
"geometry": {
        "type": "LineString",
        "coordinates": [
          [
            80.5078125,
            67.06743335108298,
            1
          ],
          [
            81.2109375,
            65.2198939361321,
            1
          ],
          [
            74.1796875,
            63.39152174400882,
            1
          ]
        ]
      }

How do get it working?


texts:
Where can I find free elevation data for Europe?

The available resolutions and terms of use are important details that I also seek.


texts: What are row.Object and row.Text in Arcpy?

texts:
I'm trying to understand an arcpy script and I'm stuck with the following lines:

cur = arcpy.SearchCursor(Tabelle)
row = cur.next()
while row:
    text1 = row.Objekt
    text2 = row.Text

As far as I understand it I load a table into the variable cur.
the variable row moves to the next line of the table.
But what are text1 and text2? I don't understand the Notation row.Object.

can someone explain me what happens in line 4 and 5?


texts: Is a Helmert transform necessary to convert between WGS84 and NAD83?

texts:
I've read things all over the internet which claim that NAD83 and WGS84 are effectively equivalent (with minor ellipsoid differences). However, if I do geodetic (WGS84) -> ECEF -> geodetic (NAD83) conversions, using only the WGS84 and GRS80 ellipsoid parameters and no Helmert transformation, I can't replicate the results of government sponsored online converters like this one. My geodetic/ecef conversions are done following EPSG note 7.2 and have been rigorously tested for correctness.

Thus my questions are:

Is a Helmert transform required to convert WGS84 to NAD83?
Are there canonical or official coefficient values for doing so? I know there are variants of NAD83 and WGS84 over time, but I'm interested in coefficients that would be used today (I think that means NAD83 (2011) and WGS 84 (G1674)?)


texts: Photogrammetric point cloud classification

texts:
My question is regarding the point cloud classification and editing work for a photogrammetric point cloud generated using drone based imagery.

My objective is the generation of a digital elevation model from the point cloud.

The terrain has high vegetation and relief with few ground points scattered across the terrain. Using auto ground classification tools aren't giving sufficiently accurate or satisfactory results.

I'm looking for specific workflows using open source/proprietary toolkits for classification and editing of the same.

Kindly make some suggestions on the same.

Attached along is the sample snapshot.


texts: Editing feature attribute in layer, which is made from featureCollection, using ArcGIS API for JavaScript?

texts: I know ArgGIS has a JS API to edit attributes of a layer when the layer is published as a feature service but is there a way to edit attribute values if a layer is made from geoJSON as a featureCollection?


texts: convert ArcGrid (.adf) to XYZ using gdal_translate

texts:
I am using:  ">gdal_translate -of XYZ w010101.adf ASCII.txt" in the command window to try to get a file with coordinates and elevations from a raster I downloaded from The National Map.

This is the error message I get: "w010101.adf not recognized as a supported file format" However, on TNM website, they show an example that uses a .adf file as the input file (link).

I have checked to make sure that the gdal_translate and the .adf file are located in the right place.

I have seen other questions related to getting XYZ files from .tif files, but not from .adf, and I'm not sure what I'm missing.


texts: Is it possible to put scale(ex; 1:10000) at print composer on QGIS 1.8.0

texts:
I would like to display a cartographic scale (e.g. 1:10,000) in the print composer.

I can set the scale from map tab.
After setting, I can type by myself as Label.

Is it possible to put scale automatically??


texts:
I have several rectangular polygons that overlap each other and I'd like to label them with the ID's (for atlas page numbers).  The image below shows my initial setup.  I've offset each label from the centroid using the data defined override which I can easily understand how this one works.  But I would like to add more offset to each individual label using a CASE WHEN expression, so here's the issue.

At the bottom of the image are two inputs for "Offset X,Y".  The data defined override for this option says:

Data defined override
Active: yes (ctrl/right-click toggles)
Expected input:
double coord [X,Y] as <= 0.0 =>
Valid input types:
string, int, double
Current definition (expression):

I've used the overide for other individual coordinate values (just X or just Y) and they work fine because they expect just a single number.

But what is the format for the expression to use BOTH X and Y values?

I've tried two comma sperated numbers in square and curly brackets, as well as just numbers and also without the comma separation but I keep getting Expression is invalid from the output preview.  Any answers that I come across only really refer to python which I can't use.


texts: Confusion Matrix in ArcGIS Pro not for accuracy but comparing two raster

texts:
I am wondering how to calculate or process confusion matrix on ArcGIS or other software for two raster data?
I have tried to do it on ArcGIS by create accuracy point (~700 random point) to get the matrix but it is actually not exactly comparing the two raster as it just 700 random point.
I don't know about Python and wondering if there is other way to figure it out?
Both of two raster has been clipped and reclassified into 378 x 378 cell per columns and rows, and also same classes.
I am wondering how can I overlay them and get the confusion matrix?
In the other way, can I just create 378x378 which is 142884 random points and will they sit on every cell of my raster?


texts: I am trying to upload a CSV file with Lat and Long data in it and there are two columns that are not numeric (see picture below). Whenever I try and upload the file, QGIS recognizes it as a string file and loads all the data in one line, which isn't what I want. Is there any way to get this to load as a delimited text layer or should I be loading it in as a different layer? I tried loading it in as a vector layer but then it called the lat long columns "real" not "integer" and I couldn't see them displayed (probably because of this).


texts:
I am a novice in writing scripts.

I have a line geometry. Let's say line.shp.
I want to smooth it with v.generalize, e.g.

v.generalize input=line output=line_chaiken method=chaiken \
threshold=1 error=roads_chaiken_leftover

What do I have to do to get the things done? Do I somehow have to do in three steps, like

Load into GRASS
v.generalize
Write output

How would these steps fit into a script?


texts: Convert EPSG::4326 degrees to km, given coordinates

texts:
I have the following data structure

pointid distance lat long

where pointid is a list of points, distance is distance of point from feature (calculated in QGIS in geographical units under a EPSG::4326 reference system), and lat and long are the coordinates of that point.

How can I convert the distance variable to km?


texts:
I have SharpMap loading a shapefile and displaying the rendered image to a PictureBox.

If I open the shapefile in QGis I can see the SRC:

  +proj=merc +lon_0=100 +lat_ts=-41 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

However, when I tried to convert to Webmercator using SharpMap, the map renders empty.
This is the code I use to re-project (data is the shapeFile object):

    ICoordinateSystem sourceTransform = data.CoordinateSystem;

    ICoordinateSystem webmercator = ProjNet.CoordinateSystems.ProjectedCoordinateSystem.WebMercator;

    CoordinateTransformationFactory ctFac = new CoordinateTransformationFactory();

    ICoordinateTransformation transform = ctFac.CreateFromCoordinateSystems(sourceTransform, webmercator);

    vectorLayer.CoordinateTransformation = transform;

I also tried using different WKT, such as

  PROJCS[""WGS_1984_Web_Mercator"",
                  GEOGCS[""GCS_WGS_1984_Major_Auxiliary_Sphere"",
                      DATUM[""WGS_1984_Major_Auxiliary_Sphere"",
                          SPHEROID[""WGS_1984_Major_Auxiliary_Sphere"", 6378137.0, 0.0]],
                      PRIMEM[""Greenwich"", 0.0],
                      UNIT[""Degree"", 0.0174532925199433]],
                  PROJECTION[""Mercator_1SP""],
                  PARAMETER[""False_Easting"", 0.0],
                  PARAMETER[""False_Northing"", 0.0],
                  PARAMETER[""Central_Meridian"", 0.0],
                  PARAMETER[""latitude_of_origin"", 0.0],
                  UNIT[""Meter"", 1.0]]

... but there's no transformation at all.

How can I create this transformation?


texts: Need new house location added to Google Maps

texts: We have recently built a new house in an existing subdivision.  I am registered with 911 and receiving mail at this address.  Street is not currently shown on Google maps.  How do I add my location to Google Maps so deliveries can find me.  I am also having problems with changing my address (getting invalid address messages) ... apparently since it does not currently show up in Google Maps.


texts: Exporting USDA NASS CDL 2005, 2007 using Google Earth Engine?

texts:
I'm trying to export the USDA NASS CDL masks from 2003 to 2018. However, it seems that years 2005 and 2007 have two versions: 2005a, 2005b, 2007a, 2007b. Does anyone know why there's two, and how to remove or combine them otherwise? I tried checking the FAQ but it doesn't provide much information.

For reference,

https://developers.google.com/earth-engine/datasets/catalog/USDA_NASS_CDL

I tried to remove 2005b and 2007b, but not sure how to filter them out using Google's api.

var dataset = ee.ImageCollection('USDA/NASS/CDL')
              .filter(ee.Filter.date('2002-12-31', '2018-12-31'))
              // This filter doesn't work.
              .filter(ee.Filter.eq("id", "USDA/NASS/CDL/2005b").not());

// Allows me to see number of elements in the query in the console.
print(dataset)

var cropLandcover = dataset.select('cropland');
Map.setCenter(-100.55, 40.71, 4);
Map.addLayer(cropLandcover, {}, 'Crop Landcover');

Edit:

Using the answer below, I found that 2005b adds information about Mississippi and 2007b adds information about California. There isn't any intersection, so it seems like you can just merge 2005a + 2005b and export that as one mask for 2005, and same with 2007.


texts: How to get the previously added layer information for getGetFeatureInfoUrl?

texts:
I am able to fetch the getGetFeatureInfoUrl for the currently selected layer, is it possible to fetch the last added layer FeatureInfo on the map?

My code currently is:

<!DOCTYPE html>

    <title>WMS GetFeatureInfo (Tile Layer)</title>
    <link rel="stylesheet" href="https://openlayers.org/en/v3.19.1/css/ol.css" type="text/css">

      Disaster
      <select id="dict" class="form-control" onchange="loadfloodlayer(this.id);">
        <option>Select disaster</option>
        <option value="guntur">guntur</option>
        <option value="krishna">krishna</option>
        <option value="prakasam">prakasam</option>
      </select>

&nbsp;

&nbsp;


texts: Shortest path calculation from shapefile data using GeoTools?

texts:
I'm using geotools 14.3 with swing to load two shape files. A road network and a node network (road junctions).

Objective

I want to calculate the shortest path from each node to every other node. I've used python with QGIS to successfully do this. But now I want to achieve this programatically using GeoTools. I've been following this tutorial.

I've followed up to graph building, now I want to calculate shortest path. But it's a bit unclear for me on selecting nodes.

//find a source node (usually your user chooses one)
Node start = ..

How could I select a node from my point shapefile?


texts:
In Arcmap (10.1 version)
I create a new field call "gender", and i want the content of "gender" will be changed according to others fields. For example, if a word(boy) is shown in a field call "person" in one row, then the content of "gender" will also change to boy (replace the original content).

Here are my code(VB script):

IF [personA] = "boy" and [personB] = "boy" THEN
    [gender] = "boy,boy"
EndIF

IF  [personA] <> "boy"  THEN
    [gender] = "Non-boy,boy"
EndIF

IF  [personB] <> "boy"  THEN
    [gender] = "boy,Non-boy"
EndIF

however, it cannot run and a message was shown.

A field name was not found or there were unbalanced quotation marks.
ERROR 999999: Error executing function.
Item not found in this collection.
Failed to execute (Calculate Field).


texts: How to process MODIS ndvi data using QGIS/GRASS/R?

texts:
I have MODIS ndvi data (MODIS_Grid_16 days NDVI) from the iberian peninsula (3 tiles), from February of 2000 to August of 2012 (really big data). I want to know how to deal with this kind of files.
My final goal is to extract the ndvi value to points with a associated date, i.e., I have sequencial dated points and I want to know the value of ndvi to that date point.)

How do I (automatically) merge the 3 tiles that compose the iberian peninsula for each 16-days (2000-2012)?
(I think if I use R to identify each filename and use gdal merge, it might work, but maybe there's a simplier way)

How do I create a time series from it and extract the values associated?


texts:
I created a map in Tilemill, exported to mbtiles, and created a cache of tiles using TileStache. Until now, I have been using Apache as a reverse proxy for TileStache to serve the tiles, but I'm wondering if it would be more efficient to simply serve the cached tiles with nginx.

Since tiles are static files, and I don't plan on creating new ones very often, is there any reason not to serve them out of the disk cache with nginx, as opposed to using some other piece of software?


texts: Stand alone Calculate Field script in Python

texts:
I am familiar with using the ArcPy Calculate Field tool in Python but was wondering how this is done in a stand alone Python script with a .csv.
My ArcPy Snippet looks like this. I assume it's a basic procedure but I couldn't find any previous questions. Feel free to direct me to a similar question.

codeblock = """def Reclass(PA,dist):
                    if PA=="M" and dist==1:
                        return "MM"
                    elif PA=="M" and dist==0:
                         return "MW""""

arcpy.CalculateField_management(fc, field="NEW", expression="Reclass(!Field1!,!Field2!)", expression_type="PYTHON_9.3", code_block=codeblock)


texts: QGIS Join table with shp file. Get Null Values on the join fields

texts:
If my problem is that my ID types are not the same, I need help in addressing that.  I started with an excel file, I read that on a mac to I need to save that as a Windows comma Separated file.  I did.

QGIS 2.14
Add Layer - Add Delimited Text Layer -

the csv file has

parcel_num as text
situs_cmptaddr_number as text
acres as double

After join the Shp file has

PARCEL_NUM as string
situs_cmptaddr_number as text
acres as double

I have tried loading the csv file as a vector layer with a like named .csvt file but it did not pick the the details in the csvt file.

I am new to QGIS.

My file names are --- mailingList_w_csv.csv and mailingList_w_csv.csvt.

I am using a mac.  I have since used a PC to create the csv file from an excel file.  the csvt file is being edited on my mac

The csvt file looks like this - “String”,”Text”,”Double”


texts: QGIS 3x calculate distance between two points with coordinates

texts:
I have a problem, very similar to Calculating distance between two points with coordinates using ArcGIS Desktop where I would like to calculate the distance between two polygons (circles) with their coordinates defined by a centroid function in the field calculator.
I found some possible options to calculate the distance between them:

with Python console:  where it looks like I should do the calculation point by point
Distance matrix - not really applicable in y case, as I have the polygons instead of points
NNjoin plugin - it looks like some layers must be merged, I don't want that
GRASS v.distance - I used it, but I am getting an error as you see below.

I also used the following formula:
  distance(($geometry),geometry(get_feature_by_id('Area9_Poles', 1)))

as per the answer below:
https://stackoverflow.com/questions/63231575/calculate-distance-between-two-points-in-qgis-field-calculator
but it gives me the cumulative distance from one point only.

So, the distance is calculated always from P1. If you look at the image above you will get my point.
In the Distance [meters] column I need the calculation:
from P1 to P2 distance (which is correct)
from P2 to P3 distance, which is not correct, because it's calculated from P1 to P3 and gives 47m, whereas it is 27.
What is missing in this formula then?
How can I calculate the distance consequently from the last point instead of the very first one?
My full data attribute table is below:


texts: How can I access custom parameters for WMS layers via ArcObjects, and create a new layer file with the updated parameters?

texts:
I have a set of data in PostGIS related to a large number (7,000) of species, and have a single PostgreSQL view, and corresponding GeoServer layer, that contains data for all species.  Based on input I received for another question, I am using a CQL_FILTER in the custom parameters for an ArcMap layer, to display data for one species at a time, as below:

This has worked great for me.  Now, I am looking for a way to make it easy for my biologists (who will no doubt make mistakes when editing the CQL_FILTER -- it seems simple, but I could see their eyes glazing over when I tried to explain the process) to see data for a given species.  What I would like to do is to programmatically create one new .lyr file for each species, by editing the CQL_FILTER and saving as a new .lyr file, with an appropriate name for the species, so that biologists can just add the appropriate .lyr file for the species.

I have poked around in the ArcObjects API documentation, but cannot find where these custom parameters are described.  As far as I can tell, they only pertain to web mapping service-type layers.  It seems like these parameters should be accessible via the WMSMapLayerClass object, but I'm not finding any description of it in the documentation.


texts:
Is there any way to get a representation of a WMS legend in a form that isn't an image? eg XML

Using an image for the website in our system will mean the text is not sharp, it doesn't resize with the web site, and the format of it is limited to what WMS provides.


texts: Troubleshoot code error for If then else Logical VBScript in ArcMap 10.3.1

texts:
I realize this question has been asked in its basic form before, but I've tried a lot of code samples on answers posed to similar questions on this site to no avail.

I'd like to create a new field in my attribute table (Label) that uses the VBScript comparison operators in Field Calculator to define the value based on text values of a field within the table (STATUS).

If the field "STATUS" is "USFS", I want "Label" to equal "USFS".

If the field "STATUS" is not "USFS", I want the "Label" to equal "Other".

However, I always get a return value of "Other".  Here is the pre-logic script code I'm using:

dim n

if [STATUS] = USFS Then

n = "USFS"

elseif [STATUS] <> USFS Then

n = "Other"

end if

See screenshots.


texts: Counting number of polygons within which points are located using ArcGIS Desktop?

texts:
I have used the network analyst to calculate service areas of varying facilities.  Additionally I have points which represent buildings.

To make it short, I have the following data:

Polygons - which represent service areas
Points - which represent buildings

Is there a possbility to determine which buildings have access to zero, one, two, etc. etc. facilities?

In my understatement that would mean to to count the number of polygons in which a certain point is located.

How can I do that?

I am using ArcMap 10.6 with an advanced desktop license.


texts: Creating a Origin-Destination flow map with Non-Spatial dataset

texts:
I am working on a project where I am trying to learn some geopandas. The dataset that I am using is the LEHD Origin-Destination Employment Statistics that gives you a csv that has a work/home geoid for each block. To simplify this I have done a group by using the first 12 numbers in the geoid to get the block groups and have replicated this process for counties.
    STATEFP     COUNTYFP    TRACTCE     BLKGRPCE    GEOID   NAMELSAD    MTFCC   FUNCSTAT    ALAND   AWATER  INTPTLAT    INTPTLON    geometry
0   29  039     870300  3   290398703003    Block Group 3   G5030   S   15503895    348967  +37.7186844     -093.7752580    POLYGON ((-93.79879 37.69977, -93.79877 37.699...
1   29  039     870300  4   290398703004    Block Group 4   G5030   S   14967570    27779   +37.6998810     -093.8132224    POLYGON ((-93.83355 37.67689, -93.83352 37.677...
2   29  039     870100  5   290398701005    Block Group 5   G5030   S   110381394   634771  +37.8005970     -093.9893682    POLYGON ((-94.07332 37.75846, -94.07332 37.759...
3   29  039     870100  1   290398701001    Block Group 1   G5030   S   210378763   2041441     +37.7770123     -093.7312090    POLYGON ((-93.83839 37.76406, -93.83834 37.764...

I am then trying to take shapefiles and associate 2 centroids to each row (work and home) by using df.centroid, but when I do that I am getting a lot of rows that are returning None as the df['center'] value.
Is there a better way to preform this? Am I not calculating the centroid in the right way?
#!/usr/bin/env python
# coding: utf-8

import urllib
import os
import numpy as np
import pandas as pd
from geopandas import GeoDataFrame, read_file, GeoSeries
from shapely.geometry import Point, LineString, Polygon, MultiPoint
from datetime import datetime, timedelta

import sys
sys.path.append("..")
import movingpandas as mpd

# read the block group file from the census
bg = read_file('tl_2019_29_bg.shp')

bg['center'] = bg.centroid
bg['bg_geoid'] = bg.STATEFP.astype(str)+bg.COUNTYFP.astype(str)+bg.TRACTCE.astype(str)+bg.BLKGRPCE.astype(str)

# use geopandas to read the county file from the census
co = read_file('tl_2019_us_county.shp')
co['center'] = co.centroid
co['co_geoid'] = co['GEOID']

# use geopandas to read the origin-destination file from lehd
od = read_file('mo_od_main_JT00_2017.csv')

# calculate out the block groups and counties for each entry
# blocks id's can be concatenated to blog group / census tracts by shortening the length of the geoid
wh_list = [['w_bg_code', 'w_geocode', 12],
           ['h_bg_code', 'h_geocode', 12],
           ['w_co_code', 'w_geocode', 5],
           ['h_co_code', 'h_geocode', 5]]
for k in wh_list:
    od[k[0]] = od[k[1]].str.slice(0,k[2])

# filter columns and od rows so that it only is looking at MO data geoid = '29'
mo_od = od[['w_geocode', 'h_geocode', 'w_bg_code', 'h_bg_code', 'w_co_code', 'h_co_code', 'S000']][od['w_geocode'].str.startswith('29') & od['h_geocode'].str.startswith('29')]

# group by statement for bg flows (takes a long time to run)
mo_od = mo_od.groupby(['w_geocode', 'h_geocode', 'w_bg_code', 'h_bg_code', 'w_co_code', 'h_co_code']).sum().reset_index()

# filter out non-missouri counties GEOID == '29'
mo_co = co[co['GEOID'].str.startswith('29') & co['GEOID'].str.startswith('29')]
mo_co.loc['center'] = mo_co["geometry"].centroid

join_list = ['w','h']

for item in join_list:
    mo_od.merge(bg[['bg_geoid', 'center']], left_on=f'{item}_bg_code', right_on='bg_geoid')
    mo_od.loc[f'{item}_bg_center'] = bg["geometry"].centroid
    mo_od.merge(co[['co_geoid', 'center']], left_on=f'{item}_co_code', right_on='co_geoid')
    mo_od.loc[f'{item}_co_center'] = co["geometry"].centroid

mo_od[mo_od['w_bg_center'] == None]


texts: QGIS - Create single contour line at specific height from DEM

texts: I wish to create a single contour line of height 4.54 m from a DEM layer. Any idea how I may do this please?


texts: Understanding results of sdo_cs.determine_chain

texts:
I run this query:

select sdo_cs.determine_chain(null, null, 3087, 4326) from dual;

This is the result:

MDSYS.TFM_PLAN(MDSYS.SDO_TFM_CHAIN(3087,-1,4152,-2,4326))

How should I interpret the -1 and -2 in the transformation plan result? According to the docs these are operation IDs but I can't anything in SDO_COORD_OPS that matches.

A few other examples:

3087 -> 4326 via 4152

SDO_CS.DETERMINE_CHAIN(NULL,NULL,3087,4326)
MDSYS.TFM_PLAN(MDSYS.SDO_TFM_CHAIN(3087,-1,4152,-2,4326))

3087 -> 4152

SDO_CS.DETERMINE_CHAIN(NULL,NULL,3087,4152)
MDSYS.TFM_PLAN(MDSYS.SDO_TFM_CHAIN(3087,-1,4152,-2,4152))

4152 -> 4326

SDO_CS.DETERMINE_CHAIN(NULL,NULL,4152,4326)
MDSYS.TFM_PLAN(MDSYS.SDO_TFM_CHAIN(4152,-2,4326))

4326 -> 3087 via 4152

SDO_CS.DETERMINE_CHAIN(NULL,NULL,4326,3087)
MDSYS.TFM_PLAN(MDSYS.SDO_TFM_CHAIN(4326,-2,4152,-3,3087))

4326 -> 4152

SDO_CS.DETERMINE_CHAIN(NULL,NULL,4326,4152)
MDSYS.TFM_PLAN(MDSYS.SDO_TFM_CHAIN(4326,-2,4152))

4152 -> 3087

SDO_CS.DETERMINE_CHAIN(NULL,NULL,4152,3087)
MDSYS.TFM_PLAN(MDSYS.SDO_TFM_CHAIN(4152,-2,4152,-3,3087))


texts: Dissolving creates fragments in the middle

texts:
When dissolving a polygon shapefile I get fragments in the middle (see pictures) in Qgis and R using the functions from rgeos. Has anybody an idea where this comes from and what can be done to deal with this?

original data set showing municipalities

dissolved dataset showing fragments in the middle

I would like to have only the contour of all municipalities without the fragments...


texts: Problem converting georeferenced jpgs to geotiff - image is shifted/warped

texts:
I have hundreds of georeferenced jpgs in ArcGIS 10 that I need to convert to geotiffs.

I'm trying to use GDAL translate to create geotiffs from them.  Everything works great except that the resultant geotiff is slightly shifted/warped (see image 1 and note that the road locations are shifted on the right panel.) Can anyone suggest how to convert georeferenced jpgs to geotiffs without this shift/warp?

My GDAL command to create this geotiff was:

gdal_translate -of GTiff input_file.jpg output file.tif

Gdaltranslate indicates exactly the same coordinate system and bounding coordinates between the jpg and the geotiff (image 2). I'm not sure what is wrong.

Thank you for your help!  (Images are provided as dropbox links because I don't have enough reputation points to post images here.)


texts: Filter combobox to show only Points layers

texts:
I am creating a plugin using QT Designer and QGIS 2.18.  Using Plugin Builder along with examples I have been able to create a combo box that list all current layers in QGIS.  I only want to list point layers.  From hours of research I keep reading similar threads saying to use custom QT widgets like QgsMapLayerComboBox.  But this is where I get confused.

Are these Qgs widgets added into QT Designer in some manner?  Or in my script.py directly? Then I read one has to edit the ‘dialog_base.ui’ every time I make change in QT Designer?  Can someone lead me on how to filter the combobox below to only view point layers?

Here is my code:

def run(self):
    """Run method that performs all the real work"""
    layers = self.iface.legendInterface().layers()
    layer_list = []
    for layer in layers:
        layer_list.append(layer.name())
    self.dlg.comboBox.addItems(layer_list)

    # show the dialog - this is needed
    self.dlg.show()
    # Run the dialog event loop
    result = self.dlg.exec_()
    # See if OK was pressed
    if result:
        # Do something useful here - delete the line containing pass and
        # substitute with your code.
        pass


texts: Updating location field (geometry) in Oracle Spatial using asp.net web application

texts:
I want to update a geometry column in Oracle Spatial database via Asp.net web application, I couldn't find samples for update query or datatypes that can hold geometry data in asp.net.

I use C# for server side , Oracle Spatial 11g ,asp.net web forms and ESRI JsAPI.


texts: Plot spectral signatures to compare with signatures from field spectrometer

texts:
What software I can use to plot the spectral signature of a Landsat 8 or Sentinel 2 pixel to compare it with the spectral signature of a field spectrometer (HR2000+)?

How can I do this?


texts:
My two input feature classes are populated with data that intersect, the output feature class does not exist, and I have try: except: blocks with traceback, but python crashes here and doesn't get to any of that.

    outfc = r'C:\Temp\temp.shp'
    out_buff_fc = r'C:\Temp\out_buf.shp'
    out_line_fc = r'C:\Temp\out_line.shp'

    arcpy.SpatialJoin_analysis(target_features = out_buff_fc,
                               join_features  = out_line_fc,
                               out_feature_class = outfc)

What does it take to get this to run?


texts: Mapbox GL JS with self hosted vector tiles has missing levels?

texts:
Using Mapbox GL JS with self hosted vector tiles in protobuf format where it looks for tiles using {z}/{x}/{y}.pb paths.

I have a scenario where I only have data for zoom levels 9, 11 and 13.

Is there any way to prevent Mapbox GL JS to try loading other zoom levels?

At the moment I get a bunch of 404's on zoom level 12 for example.


texts: How to get random points from large FeatureCollection without exceeding the user memory limit?

texts:
How do you select random points from a large FeatureCollection without exceeding the user memory limit?
var bigData = ee.FeatureCollection('GLIMS/current');
var randomPoints = ee.FeatureCollection.randomPoints(bigData, 100);
print(randomPoints.size()) // Error: User memory limit exceeded

This works on smaller FeatureCollections:
var smallData = bigData.filter(ee.Filter.eq('geog_area','Iceland'))
var smallRandomPoints = ee.FeatureCollection.randomPoints(smallData, 100);
print(smallRandomPoints.size()) // 100


texts: Why is the performance speed between ArcGIS and QGIS so different?

texts:
Ok I am not a programmer but a prolific GIS user. I know that QGIS is written in C++ and ArcGIS in ??? but for most of my tasks lately I always try to use QGIS not only because it is free but due to the fact that its User experience is so good.

All the GIS Gurus out there can you tell me some reasons for the difference in speed between these two systems? Honestly it pains me to use ArcGIS 10 due to its speed and I have a PC with 8 GB RAM.


texts: Rule-based symbology based on current selection in QGIS

texts:
I have two oracle layer (point(p1) and polygon(pl1)). I want to build a rule-based symbology which is based on current selection on two cases:
Case 1- When I select a point from table p1, I want that the polygons from table pl1 to be coloured differently based one the relation that they have with table p1 which is p1.id_b=pl1.id_b (one to many).
Case 2- When I select a polygon from table pl1 other polygons from the same layer that are related to current selection to be coloured differently. Relation between polygons is on the field pl1.id_group, so polygons are grouped by the same pl1.id_group (one to many).
Software used is QGIS.

For illustration I select point one and the related polygons are coloured differently from the others.