texts: Converting tile image into raster using rasterio

texts:
I want to convert a tile image (.tif) into a raster data (.tif) and assign lon/lat to the data. (I know how to get the lon/lat of a tile.)
The problem is I don't know how to assign lon/lat to a raster image. The only method I know is copy the format of a raster image to another. But the location of the tile will be totally wrong.
The library I use is rasterio on Python. Here's the code I've tried
import numpy as np
import rasterio
import rasterio.features
import rasterio.warp
from affine import Affine

# copy the format of this raster image
with rasterio.open('SN5_roads_train_AOI_7_Moscow_PS-RGB_chip999.tif') as src:
    transform = src.transform
    #transform = Affine.identity()
    crs = src.meta['crs']
    mask = src.dataset_mask()
    print(rasterio.features.shapes(mask, transform=src.transform))
    # Extract feature shapes and values from the array.
    for geom, val in rasterio.features.shapes(mask, transform=src.transform):
        # Transform shapes from the dataset's own coordinate
        # reference system to CRS84 (EPSG:4326).
        geom = rasterio.warp.transform_geom(
            src.crs, 'EPSG:4326', geom, precision=6)

# read the rgb of the tile image
with rasterio.open('original.tif') as src:
    r, g, b = src.read()
    print("bounds")
    print(src.bounds)
    print("meta")
    print(src.meta)

# generate the raster image
with rasterio.open(
    'generated_raster.tif',
    'w',
    driver='GTiff',
    height=r.shape[0],
    width=r.shape[1],
    count=3,
    dtype=r.dtype,
    crs=crs,
    transform=transform,
    #transform=Affine.identity(),
) as src2:
    src2.write(r, 1)
    src2.write(g, 2)
    src2.write(b, 3)
    print("new.tif meta")
    print(src2.meta)

with rasterio.open('generated_raster.tif') as src:
    transform = src.transform
    crs = src.meta['crs']

    mask = src.dataset_mask()
    print(rasterio.features.shapes(mask, transform=src.transform))

I am new to GIS and raster data.


texts: ArcPy string concatenation with both forward and back slashes

texts:
I am trying to concatenate 2 string variables Var1 and Var2 using python (converting some old ArcGIS VB-based scripts to ArcPy). I need to add a backslash between them.
arcpy.CalculateField_management(TargetDB, "myfield", '!Var1!' + '\\' + '!Var2!', "PYTHON")

I get an error message that says

"SyntaxError: unexpected character after line continuation character"

when the code block runs. I believe the cause of the error is that Var2 contains a forward slash character. It isn't file or folder path related, it's an ID field that acts as a primary key in a database (which is frustrating, but something I can't do anything about).
Given that I can't really modify the values of Var2 because of the database situation, does anyone have any recommendations of a way around this in Python 3 to work in ArcGIS Pro 2.6? The old VB code worked, so I am assuming this is a python specific issue.


texts: ImportError: cannot import name 'CRS' from 'pyproj'?

texts:
I have installed pyproj version 1.9.6. using the command prompt window. However, when trying to import CRS (from pyproj import CRS), I receive the following ImportError:

ImportError: cannot import name 'CRS' from 'pyproj' (C:\Users\snn490\AppData\Local\Continuum\miniconda3\envs\CoInf\lib\site-packages\pyproj__init__.py)

How can I solve this problem?


texts: Convert SmallWorld GIS topology to pgRouting topology

texts:
Is there any way to convert SmallWorld GIS topology to the format that pgRouting use?
My costumer have use SmallWorld GIS and make his topology on it. but when he export his information as shape file, the topology is gone!

How can the topology be preserved?


texts: Location is not accurate when replicating an existing feature and zooming in

texts:
I am trying to resolve an issue with the accuracy of a selected point using queryRenderedFeatures. The accuracy is fine when zoomed out, but when zooming in it appears off.

I am using the coordinates from the results of the queryRenderedFeatures to set the coordinates of a separate selection layer. This is the "blip" underneath the "properties" layer.

Here is the basic structure of my code:

const features = map.queryRenderedFeatures(e.point);

const selected = features[0].toJSON() // "properties" layer is at index 0

map.addSource('propertySelect', { // source for the selected property
                "type": "geojson",
                selected
            });

map.addLayer({"id": "propertySelect", // layer for the selected property. The "blip"
              "source": "propertySelect",
              "type": "circle"
              "circle-radius-transition": {duration: 0},
              "circle-opacity-transition": {duration: 0},
              "circle-radius": 0,
            }, 'properties');

TLDR

Setting a features coordinates based on an existing feature does not seem to match the same coordinates when zoomed in. See image below.


texts: Correct way to perform point-in-polygon queries on streaming point data?

texts:
Short version: What is the "correct" Esri-based solution to the presumably-common scenario of receiving dynamic asset location data (eg. aircraft), checking to see whether those points are within any of a set of polygons (eg. flight zones), and updating the point record with the identifier of the polygon which contains it (eg. "Aircraft A is within Zone 1")?

Long version: We provide a service which tracks assets and reports whether or not they are inside an area of interest (AOI). We have written our own point-in-polygon code to perform this check, but it's not scaling. We're assuming that a dedicated GIS will be better at it than we can make it.

Currently our system is something like this:

Asset locations are stored as records in MS SQL Server 2008 with X & Y fields for lon & lat
AOIs are stored as polygons with geom in a simple geometry-type field, and minX,minY,maxX,maxY coordinates as attributes
Every few minutes, we check each new record against each of the existing AOI polygons by first checking if the lat/lon is within the bounding coordinates, and if so, we apply our point-in-poly code to determine whether it's really inside the polygon
If it is, we update the point record with the identifier of the polygon it falls inside.

At the moment we get ~10 point updates per second and check them against ~<100 polygons. In the near future, we need to be able to scale up to checking hundreds of point updates per second against up to 10,000 polygons. Our tests indicate that our current system's not going to be fast enough, and that our point-in-polygon code is the bottleneck.

It's been proposed we use Esri's ArcGIS Server software to perform this check. We have access to the software to do a proof-of-concept run. The thing I'm sticking on is exactly how to integrate it into our workflow. Would we hit the GeometryService using a HTTP request? If so, would we first have to translate our point and polygon geometry to JSON and interpret the response? Would we have to do that for every point, and every polygon? If so, it seems this overhead would likely outweigh any gain in processing speed from using Esri's geometry service (but that's an assumption).

Alternatively, would we set up the server to have full awareness of our existing database and use some kind of map/geoprocessing service perform the query? On the surface, that seems like it would be roughly equivalent to building a tesla coil to kill a fly.

Generally, I've got that "surely it can't be that hard" feeling. I'm assuming this is a bread-and-butter use case for GIS, but I'm having trouble finding implementational details. I suspect there's some obvious thing that I'm missing, but it'd also be useful to hear "you've got the right idea, it's not as overkill as it seems". So, any thoughts would be appreciated.

It's worth making explicit that although we have access to the Esri suite to test, and are somewhat interested in the other functionality it provides, if there's an easier/more efficient non-Esri solution to our immediate need it will win out, so I'd appreciate that advice too.


texts: Reading shapefile faster using Python basemap (and pickle)

texts:
I am trying to plot a simple 'merc' map with boundary from shapefile. The total size of the shapefile ne_10m_admin_0_countries_lakes.shp is just 8mb. The simple mslp surface plot from GFS data took about more than 28 sec which I think is too much. After investigation I found that it is the reading of shapefile consuming around 10 sec. A simple code to demonstrate the issue is shown below:-

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

south = 0
north = 5
west = 70
east = 85

m = Basemap(projection='merc', llcrnrlat=south, urcrnrlat=north,
        llcrnrlon=west, urcrnrlon=east, resolution='c')
m.readshapefile('data/gis-data/world_countries/'
    'ne_10m_admin_0_countries_lakes', 'ne_10m_admin_0_countries_lakes', linewidth=0.7)

plt.savefig('map.png')

To show the problem I am facing, I have run the above code with and without commenting the path to shapefile. Here is time the script took after reading from shapefile:-

$ time python test.py
real    0m18.234s
user    0m17.832s
sys 0m1.020s

Here is the result without reading from shapefile:-

$ time python test.py
real    0m2.506s
user    0m2.360s
sys 0m0.324s

I hear that pickle is one of the most suitable solution. But I could not find any of the code with respect to reading shapefile.


texts:
I have an Esri WMS server serving raster data.
I can query this using OpenLayers TileWMS, all appears beautifully on my basemap.

I need to be able to query the WMS based on a polygon and display my rasters returned but for the life of me I can't work out how.

I can query on my extent using bbox nice and easily (out of the box).
I can query the WMS using my polygon (hand drawn polygon) and get a list of object ID's back.
I can use those object ID's and make a WMS request to get back JSON that contains href and the extent of the raster tile.

I am worried I am being stupid and perhaps overthinking this, but how can I get OpenLayers (5 currently) to display said image on my basemap from this GeoJSON?

To confirm, I want to display all rasters that appear geospatially inside a hand drawn polygon within my extent.


texts: How to find the smallest-area ellipse circumscribing a convex polygon?

texts:
We are currently dealing with the problem of:
How to find an ellipse circumscribing a convex polygon?
Recent googling and lots of reading could not help me so much :(


texts: Addresses, Geolocation and driving times

texts:
I have a list of addresses (Germany and Austria) that I

want to display on a map (easy)
from a random point on the map obtain a list of all the addresses that are within a driving time of x mins (hard).

for 2) can this be done with ArcGIS (I have a student account), and are there data sets that estimate driving times between different points (Germany and Austria) that are readily available and can be used for this?

Any other ideas on how to approach this problem? Is there another tool that comes to mind to solve this problem?


texts: Creating Multipoint Lines from CSV File Containing Addresses using ArcGIS Pro

texts:
I am new to network & route analysis and I am trying to use ArcGIS Pro to analyze my companies operations. Each operation has a unique ID,  start and end addresses and can have up to 25 intermediate addresses. I would like each operation to be visualized as a line composed of multiple points (the addresses / stops) and eventually I will need to preform network and route analysis on the operations. The data is structured similar to the bellow table.
+------+-----------------+-----------------+-----+---------------+
| ID   | start_address   | stop1_address   | ... | end_address   |
+------+-----------------+-----------------+-----+---------------+
| OP1  | start_address_1 | stop1_address_1 | ... | end_address_1 |
| OP2  | start_address_2 | stop1_address_2 | ... | end_address_2 |
| OP3  | start_address_3 | stop1_address_3 | ... | end_address_3 |
+------+-----------------+-----------------+-----+---------------+
The problem I have been having is that I seem to only be able to add at most two different addresses / coordinates using the XY to line tool.
Could someone point me to the correct tools or methodology to use?
I am also comfortable with Python if a custom solution is required.


texts:
I was trying to merge stream lines from many HUC polygons into one. Currently these stream lines are stored in one geodatabase as separate feature classes with each feature class corresponding to one HUC polygon. It's polyline feature and all lines have exactly the same attributes. All feature classes have the same schema.
I gave several tries by using arcpy.Merge_management() and arcpy_Append_management(), all ended up only merging several HUCs, but not all. (see the attached picture: try to merge stream lines from all green polygons, but only got stream lines from several upper polygons merged). I also tried to repair geometry first before merge, define extent before merge.., none of these seem to fix the issue.
Still new to python, anything I missed? I use ArcMap 10.7.1. Here is the code I tried:
fcList = arcpy.ListFeatureClasses("D:/test.gdb")
arcpy.Merge_management(fcList, "D:/output.gdb/merged")

emptyFC="merged"
schemaType="NO_TEST"
fieldMapping=""
subtype=""
template="example"
spatial_reference=arcpy.Describe(template).spatialReference

has_m="DISABLED"
has_z="DISABLED"

try:
    arcpy.CreateFeatureclass_management(outWorkspace, emptyFC, "POLYLINE", template, has_m, has_z, spatial_reference)

    listFC=arcpy.ListFeatureClasses("", "POLYLINE")
    arcpy.Append_management(listFC, outWorkspace+os.sep+emptyFC, schemaType, fieldMapping, subtype)

except Exception as err:
        print(err.args[0])


texts: Tried to Mosaic Sentinel 2 images but my final image appeared with some overlapping Null data

texts:
I'm working on an area that required me to use three sentinel images to fully cover it, I used the "mosaic to new raster" tool in ArcGIS to merge the three images and was surprised to find that the output image actually has the Null pixels from raster (3) instead of the fully covered imagery in raster (2) (pic related)
How can I fix this so it can look basically like "pic related 2"?


texts: When exporting MXD to PDF, map layout moves to left?

texts: When I export an MXD to PDF, the layout has moved to the left. Exporting to JPEG doesn't cause this problem. If, for example, the map in layout view is of quad 2 and quad 3, the pdf will be of quad 1 and quad 2.


texts:
Are there any free public services where I can send one point in a coordsystem and it returns the point reprojected in a specified coordinate system?

I will not be doing many requests, for now it is just for a demo.


texts: Vector Shapefile with no known coordinate system. How can I find out and transform it?

texts:
I was looking for a shapefile for watersheds in Ethiopia and I found the following data set
https://datacatalog.worldbank.org/dataset/watersheds-ethiopia/resource/0e6b847c-b9e6-464e-9141-90172ed49842
I can open it up in a new project shapefile, but QGIS assumes that it is WGS 84, and I can view it.

However, I can't work out what the real-CRS is, and I want to overlay it with other data sets whose CRS I do know. Is there a way to work out what the CRS is, so I can transform it to a CRS that is aligned with my other shapefiles?


texts: Buffering polygon until it hits neighbor using PyQGIS?

texts:
I have a shapefile made up of a number of polygons.
Neighboring polygons often share an attribute.

What I need to do is buffer these polygons based on the common attribute. However, I would like to expand said buffers until they meet, at which point they can stop.

I have a hastily put together mock-up:

All I have code-wise is the following:

from qgis.core import *
from qgis.analysis import *

layer=qgis.utils.iface.activeLayer()

QgsGeometryAnalyzer().dissolve(layer, "D:/test_dissolved.shp", onlySelectedFeatures=False, uniqueIdField=3, p=None)

Truthfully I am not very well versed in the QGIS API documentation and despite my best efforts I cannot find a suitable method for this.
QgsGeometryAnalyzer().dissolve allows me to dissolve the buffer but that isn't what I am after.

Should I perhaps be looking in another direction here?


texts:
Just trying to update an old project but for some reason it was set up in  EPSG:4326, WGS 84. I have no idea what base maps were originally used to create it as they are no longer on the system, and even with 'on the fly' enabled whenever I load a map in it says the extents are too large, if I reproject them they end up off the coast of Angola.
If I try loading Bing first then add the layers, it simply wont zoom to them and if I re-project the layer at that point it ends up off the coast of Angola and not in Somerset (where it belongs).
Saved each layer with EPSG 3857, WGS 84 Pseudo and then did it, again Angola.
Can anyone guess why? Or a suitable CRS that will mean it goes back to Somerset. Or do I accept and start anew?


texts:
I'm having issues exporting a georeferenced map from QGIS using version 3.10.4-A coruna. I've tried different (newer) versions of QGIS and also different computers and get the same error. Hoping some of you could help me figure this out.
What I'm trying to do:
I'm trying to export a limited view (approximately 50x50km) of a larger data set (vector data of the whole of Sweden) to a georeferenced PDF in order to view it in map software like "Avenza maps" etc.
My issues are the following:

The pdf won't display in pdf readers. Adobe reader gives me an error message reading: "An error exists on this page. Acrobat may not display the page correctly." image of error: https://i.imgur.com/5YgyIeb.png

Slim pdf reader gives me an error message reading: "#0: Code: PDL::ERR_PDF_SYNTAX Function: Investintech::PDL::DAL::Rendering::Renderer::ProcessOperator_sc_Common Line: 7900". image of error: https://i.imgur.com/poTTky6.png

Opening the pdf in "Avenza maps" displays the pdf correctly but it is not georeferenced. image of error: https://i.imgur.com/kKWV4e8.jpg

Dataset:
Here's a link to the dataset I'm using if you want to try for yourself. Its vector data covering Sweden.
https://mega.nz/file/WohwlKZA#Twhb7EALJfPQGnybRKt-urGpwSNxilr-MsWJDG2b8Wk
Here's what I know/ think so far:

QGIS does not report any errors when exporting the pdf (successfully exported layout).

My export settings are set to default, however I've tried different settings. These are my default settings: https://i.imgur.com/ALZq9qG.png

PDF export of much smaller vector datasets functions correctly in version 3-10.4 .

Opening the PDF in Acrobat reader sometimes loads some features (though far from all of them). You can see an example of this here: https://i.imgur.com/9Olx14f.png

I do not believe my computer is the cause since I've tried multiple and get the same error.

Will update if I think of more things.
If any of you know what could be the issue I'd be very grateful for a response!
Update 1* Added link to pdf: https://mega.nz/file/CtpGSDYY#HGR47VxTgpOQt2KRA6XRDlBYSRv9Sm3hwnkJ4s1k9WU
All the best!


texts: Getting different pixel depths from output raster when done in ArcMap versus Arcpy?

texts:
When I use 'Plus' in Spatial Analyst to add together two 8-bit unsigned rasters (pixel values 0-100) the output raster is 8-bit unsigned (pixel values 0-200) but when I try and do the same thing running a python script the output is a 32 signed raster (pixel values 0-200) . How can I adjust my code to have the output be an 8-bit raster?

# Import needed modules
import arcpy
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial")
raster1 = arcpy.Raster("C:\Users\XXX\Desktop\Raster_test\Confidence_2008_56m.tif")
raster2 = arcpy.Raster("C:\Users\XXX\Desktop\Raster_test\Producers_2008_56m.tif")
# Execute Plus
outPlus = Plus(raster1, raster2)
# Save the output
outPlus.save("C:\Users\XXX\Desktop\Raster_test\script_plus2.tif")


texts: Is there way to color on the outside of polygons like this? I have tried to use the inverted polygon style in QGIS but it didn't work.


texts: Clipped Imagery Larger File Size Than Input

texts:
I have jp2 imagery over a larger area. For our consultant to use this imagery, they require the file size to be smaller. To reduce the file size, I am using ArcMap to clip the imagery to only the area of interest. But when I check the size of the output file (also jp2) they are bigger than the input!

Why is this and how might I get around it?


texts: Creating Offline Maps using QGIS and OpenStreetMap?

texts:
I am working with QGIS 2.10.1, having only little knowledge with the software so far.

I am mapping certain findspots in northern Bavaria and use the OCMlandscape map as basemap. Everything works fine as long as I am in the internet, but I need the whole thing offline.

The Offline plugin Database->offline editing->convert into an offline project doesn't work. The program tells me "unable to open database file".

Afterwards I tried to download the part via Vektor->Open Street Map -> Download Data, but the download gets always interrupted.

I then downloaded the Northern Bavaria map from geofabric (file: .osm). With Vektor -> Open Street Map -> topology imported from XML and afterwards ->topology export to spatial line I get a huge set of data and a blue points or lines or polygons that has the shape of northern bavaria but doesn't look like a topographic card :)

How do I know which tags I have to export?

Or do you have any idea how to get this ground map offline?


texts: Error when importing data source from web using GDAL/OGR in GRASS 7.4.1

texts:
I am running the latest bundle of OSGeo4W with GRASS 7.4.1

I am following the tutorial in the following URL:
Analysing the ECA&D climatic data, where I am trying to import data from web using this code:

INFILE="/vsizip/vsicurl/https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_admin_0_countries.zip"

However, I get the error:

INFILE is not recognizable as an internal or external command, a program or an executable file

I am not sure this error is because of using wrong syntax, or because GDAL/OGR isn´t installed/activated properly

If you need more information please comment bellow and I will gladly provide it


texts: Removing overlapping/duplicate polygons in same layer using QGIS?

texts: I am working in QGIS on a vector layer where areas are classified using an attribute "typo". The problem I'm encountering is that many polygons are duplicated or overlapping, with the result that some areas are classified as two or more "typo". This is an error. In order to make some statistical analysis I need to clean this layer removing the overlapping/duplicates so that every inch of territory is classified as one and only one "typo"; which one is indifferent.


texts: I am searching for a single file of all roughly 11 million US census blocks. Where can I find this?


texts: Seeking algorithm to place maximum number of points within constrained area at minimum spacing?

texts:
I have a polygon layer that describes a constraint; I wish to add points within this area.  I want to add as many points as possible, but they must have a minimum spacing between them.  Is it possible to do this with GIS?

To clarify, it would be best if an ordered grid could be generated, as this would guarantee the most points.  However the constraint would rarely allow this, and it may be preferable to remove points to allow an offset to better fit within the constraint.


texts: Plotting Sentinel 2a images on Earth Engine map shows missing chunks of images interrupting continuous sweeps

texts:
I am by no means knowledgable in this field so I apologise for the lack of a better description in the title, but (shown below) it seems that a lot has been masked out. This was taken from a short 5 day period and only start to fill in after several months of composites. The code is just importing and representing:

var image = ee.ImageCollection('COPERNICUS/S2_SR')
  .filterDate('2018-07-01','2018-07-5')
  .select(['B4','B3','B2'])
  .median()
  .divide(10000);

Map.addLayer(image, {min:0, max:0.3});

What is causing this masking? Shouldn't the image be continuous along the stripes?

Even the scene classification is masked so I can't see whether the pixels were defective, and when I add the SCL band in the .select(), I get this error most of the time:

Layer 1: Tile error: Expected a homogeneous image collection, but an image with an incompatible band was encountered. Mismatched type for band 'SCL':
Expected: Type<Short<0, 255>>.
  Actual: Type<Integer<0, 65535>>.
          Image ID: 20180609T103019_20180609T103214_T31TGK
This band might require an explicit cast.

I didn't know if this was correlated so I thought I should add it in; apologies if it is unrelated.


texts: Cannot start config manager on intramaps roam

texts:
I get the following error when trying to start Intramaps roam on a windows 7 Motion tablet.

Errors Occurred: See the logfile 'C:\Intramaps Roam\ Roam Config Manager.exe.log' for details

The logfile is shown below.

Traceback (most recent call last):
  File “__main__py”, line 23, in <module>
  File “configmanager\ui\configmanagerdialog.pyc”, line 12, in <module>
ImportError: No module named events
Traceback (most recent call last):
  File “__main__py”, line 23, in <module>
  File “configmanager\ui\configmanagerdialog.pyc”, line 12, in <module>
ImportError: No module named events
Traceback (most recent call last):
  File “__main__py”, line 23, in <module>
  File “configmanager\ui\configmanagerdialog.pyc”, line 12, in <module>
ImportError: No module named events`

It has worked in the past on the same tablet, until I upgraded the software version. I have since downgraded but the error still appears.

I have attempted installing ROAM and QGIS on another tablet but I still get the same error.

I have also attempted re-installs and reboots to no avail. I have also edited the roam.config file to ensure the path is correct for QGIS to 'C:/Program Files/QGIS Essen/bin/qgis.bat'.


texts: Hosted ArcGIS Online feature does not return any data?

texts:
I would like to share a table as a feature service hosted in ArcGIS Online.  I am able to publish the table as a feature service but when using the REST API interface with for ArcGIS Server I am unable to successfully return data for a simple 1=1 for my where clause. I am receive an error which states Unable to perform query. Please check your parameters.
No return fields were requested.

What could be the cause of this?  This data is shared to the public, so a token should not be necessary.


texts: Leaflet popups only displays columns where cell information NOT NULL?

texts:
I've got the following map on JSFiddle.

The popup information comes from an Arcgis Server with 23 columns. When you click around the markers you can see that most of the popup information are not filled.

How could I add code to display only info where the info is NOT NULL?


texts: Difference between debian and ubuntugis repositories?

texts:
I just spent some time struggling to upgrade to QGIS 2.18.

I originally had these lines in my /etc/apt/sources.list:

deb http://qgis.org/ubuntugis xenial main
deb-src http://qgis.org/ubuntugis xenial main
deb http://ppa.launchpad.net/ubuntugis/ubuntugis-unstable/ubuntu xenial main

but when I tried to install I was getting error messages like:

  Some packages could not be installed. This may mean that you have
  requested an impossible situation or if you are using the unstable
  distribution that some required packages have not yet been created
  or been moved out of Incoming.
  The following information may help to resolve the situation:

  The following packages have unmet dependencies:
   python-qgis : Depends: libqgis-analysis2.18.14 but it is not going to be installed
                 Depends: libqgis-server2.18.14 but it is not going to be installed

...etc.

When I changed the sources.list entries to the following the problem was fixed:

deb http://qgis.org/debian xenial main
deb-src http://qgis.org/debian xenial main

My question is, why?


texts:
I tried to use tool of QGis, Line to Polygon but it does'nt work.

I tested with close lines in the same direction and I tested with intersection of 4 lines but the result is the same, that's generate an empty Shapefile.

I tested with Saga GIS tool in toolbox but the result is the same.

I've the version 2.14.4 of QGis and this is the link for my file : https://we.tl/l8zPIY8V4S


texts: What exactly is the min/max of an extent? (QGIS heatmap plugin)

texts:
I'm having a hard time making the heatmap plugin do precisely what I want (namely, to make a heatmap showing not only points but an attribute of those points), and now I'm thinking that I may be misunderstanding what "min/max of full extent" means.

To my mind, the minimum and maximum values QGIS pulls out of the attribute field and displays in the band rendering section should be the lowest and highest numbers that appear in that column. Is that actually correct?

(In that case, my max should be 6500 or so. Instead, it's 10.35.)


texts: Automatically snap lines at angles when offsetting

texts:
I used the geometry generator to offset lines.

  offset_curve($geometry,"offset")

Is there a way to get the lines automatically snapped at angles when they are connected?
I can use the rendering levels for this case, but that won't be so easy with a lot of lines.


texts: Seeking pretty icons for GIS development?

texts:
I am developing GIS software and I want to build a pretty GUI for it.

Where can I found some pretty icons for my GUI?

I mean for example, icons for adding a vector or raster.
I mean icons like QGIS icons.

I searched the web and found Maki icons but it was not helpful.

Similar questions:

Seeking free map icons sets?
Seeking Good Symbology for Interactive Web Mapping?


texts:
I have a series of point data from a .csv file that I would like to display as a heatmap raster, the .csv data have various attributes with them and I would like to adjust the color ramp based on which of those attributes I select.
My first thought was to apply a color ramp to a vector grid based on how many points are in the grid cells.

However, I'm a beginner at both GIS and this is my first time using QGIS. I'm happy to thrash around on my own. I've tried to create a vector grid, but pretty sure I've done it wrong because although a layer was added to the TOC I saw nothing on the canvas.


texts: Generate a JSON instead of an image for the legend with OpenLayers

texts:
I'm not clear on how I can pass the params to getLegendUrl. My aim is to obtain the JSON and for this scope I use the code below:
  var graphicUrl = wmsSource.getLegendUrl(
    resolution,
    'params': 'application/json',
  );
  console.log(graphicUrl);

In the console I see only the error below concerned to 'params': 'application/json',:

Uncaught SyntaxError: missing ) after argument list


texts: How to sum pixel values of a floating point raster in ArcMap?

texts:
I have a floating point raster with values from 60 up to 95 (with decimals), it has 7077 pixels of 250x250m size.

Is there a way to sum all pixel values of this raster? For example, if the raster had 3 pixels with values of 60.20, 71.43 and 86.59, the sum would be 218.22.

I tried with zonal statistics and cell statistics, but it does not work. I am using ArcMap 10.1.


texts: What is the sIgnificance of equator in GPS offset calculation?

texts:
I am going through a SAE document where it says that the GPS offsets (NodeLL) should be calculated in such a way that the range should be within one of these:

ranges of +- 0.0008191 degrees
ranges of +- 90.571389 meters at the EQUATOR

The way it is calculated in the code written was just a mere subtraction of latitudes and longitudes. If this is correct what is the significance of equator here? Or is there anything wrong we are doing?


texts: See the source fields of a shapefile in Geopandas?

texts:
I want to access the source fields(like in QGIS)

shp(right click properties->Source fields)

but in Geopandas.

The logic is to access the source fields (maybe export it in txt) of a valid shapefile and programmatically compare it with the corresponding(source fields) of a new similar shapefile, to test for differences that are not appropriate. For example (column names, data types, length of characters.)


texts: How to convert coordinates to shapefile with R package's dataset?

texts:
I need construct spatial weight matrix for data which syear=1995 (house$s1995==1)in housedataset in spdeppackage,So I tried writeOGR function as the link suggeust to get shapefile enter link description here.

The shapfile was built successfully，but when I use readShapePoly to read,it gets error.

There is my code:

library(spdep)
library(rgdal)
data(house)
house=house[which(house$s1995==1),]
class(house)
writeOGR(obj=house, dsn="test", layer="house", driver="ESRI Shapefile")


texts: R - calculate costDistance between two nearest points

texts:
I have two sf points datasets: cities and nearest neighbours, and a conductance transition layer generated from a DEM in the extent of the point datasets using the gdistance package. I want to calculate costDistance between each pair of points, that is each city and its neighbour, but avoid generating a distance matrix for all cities and all nearest neighbours (there are 100,000 points in each dataset, so it gets big). I tried running a loop but it was very slow. Any suggestions on how can I get this done faster?
for (i in 1:nrow(cities)){
  city <- as(cities, "Spatial")[i,]
  nearest <- as(nearest, "Spatial")[i]
  cities$Costdist_1ktown[i] <- costDistance(conductance,
             fromCoords= insc,
             toCoords=nearest)
  }


texts: How to handle temp files created during Extract by Mask (Spatial Analyst)?

texts:
I have three scripts running simultaneously that processes raster data in several directories. The problem I am encountering is that ArcGIS is creating a temporary grid file called 'Extract_recl1' when the Extract by Mask (Spatial Analyst) tool is called. The temporary file is created in a directory where it is looking for raster data. Seeing that I have several scripts running simultaneously, there is occasionally conflict where one script trys to run Extract By Mask and encounters 'Extract_recl1', which was created by Extract by Mask running in a different script. This is the error I receive:

ExecuteError: ERROR 010429: Error in GRID IO: CellLyrCreateInternal: Directory f:\mydirectory\reclass_new_classes\extract_recl1  already exists.
ERROR 010302: Unable to create the output raster: F:\mydirectory\reclass_new_classes\Extract_recl1
ERROR 010067: Error in executing grid expression.
Failed to execute (ExtractByMask).

I have tried to handle this situation by saving the raster object with a unique name to a temp directory and deleting the file upon completion:

# Extract by mask
reclass_clp = ExtractByMask(reclass, mask_raster_temp)
reclass_clp_temp = random_name_generator(temp_ws, "reclass_clp_temp", ".tif") # This is a function that generates a unique output name
reclass_clp.save(reclass_clp_temp)

# Do some processing...

arcpy.Delete_management(reclass_clp_temp)

I have also set the overwriteOutput environment to True.

arcpy.env.overwriteOutput = True

However, neither of these seem to have any effect. What is the best way to handle temporary files created by Extract By Mask (Spatial Analyst)?

In other words, is there a way to force Extract by Mask to generate a unique name for the temp grid 'Extract_recl1'?


texts: Coordinate System choice for Nearest Neighbor

texts:
Here's an interesting scenario I came upon yesterday and I wanted to throw it out as a question to try help me visualize what is happening.

I have a single point (a person, say) who is looking for the nearest transit station. The person's location is in Lat/Long and the transit station is a dataset in UTM.

To do a nearest neighbor I need to reproject either the person into UTM (to match the transit dataset) or the transit dataset into Lat/Long (to match the person).

If I reproject the person to UTM, the nearest neighbor tells me the nearest transit station is X

If I reproject the transit stations to Lat/Long, the nearest neighbor tells me the nearest transit station is Y

So... does it make sense that the nearest neighbor can be different depending on which coordinate system you are working in? My gut feeling is yes, because Lat/Long is a geographic coordinate system and UTM isn't.

Am I right? I think so but am just having problems visualizing why! In particular, which one of the results is "correct"?


texts: Calculate distance between two points in projection EPSG:900913

texts:
I have two points in 900913 projection and I want to calculate the distance between the points.

I don't want to transform points to 4326 and then calculate.

How do I do this?


texts: TileMill: How to access geojson properties with CartoCSS

texts:
When trying to display a geojson layer in TileMill, is it possible to access the 'properties' of each feature?

My geojson looks like this:

        "type": "Feature",
        "geometry": {
            "type": "Point",
            "coordinates": [2.22614, 48.8765]
            },
        "properties": {
            "type": "stop",
            "name": "SOME_NAME",
            "color": "#ff0000"
            }

and I would like to display the name in properties/name, using the color in properties/color.

Is it possible? Properties are clearly accessible since one can write this and it works:

#id[type='stop'] {
  /* various css rules */
  }

But is it possible to use properties in the css rules themselves? This for instance does not work:

#id[type='stop'] {
  marker-line-color: [color];
  }

(or any variation with or without quotes inside or outside brackets for 'color').

If it's not possible, in what other format should I transform the geojson in order to access properties of each element from the css?


texts: Combining raster file with country boundaries into a .csv dataframe in R

texts:
I have a raster file from the FAO's Global Agro-ecological Zone database ( http://gaez.fao.org/Main.html#) which gives the yields of specific crops in 5' grid cells. I would like to assign a country name to each of these grid cell (obviously drop the cells that have no yield data or lie in the oceans). I was thinking of using GADM but I don't know how to reconcile these two sources.

Ultimately, I want to create a data frame with the following columns:

Grid cell ID | Longitude | Latitude | Yield (from raster) | Country

I'm quite clueless as to where to even begin. For some broader context, I'm trying to follow this paper: https://economics.mit.edu/files/7536 by summing up the crop production in each country from the grid cells.


texts:
In the past, by following this tutorial by Tom MacWright, I was able to in a few steps display any jpeg images as an interactive map with Mapbox.

The steps are summarized below:

Run a python script via the command line to convert a jpeg image to a GeoTIFF file.
Using what is now called Mapbox Studios Classic, I make a new layer with the aforementioned GeoTIFF and then upload it to Mapbox Editor.
From Mapbox Editor I can then add features like markers and polygons and other interactive elements to the custom image map.

Unfortunately, the second step now fails ever since the Mapbox Studios desktop app (now referred to as Mapbox Studios Classic) became superseded by Mapbox Studios the web app. I can no longer upload any of my converted GeoTIFFS, not even the sames that I uploaded without a problem in the past. Instead I get this error:

  Error: Error creating Mapnik Datasource: Invalid raster: Invalid pixelsize in geotransform array

I haven't the slightest clue what that means. Googling around hasn't turn up a single hint (or at least a hint I could recognize).


texts:
I am interested in working with the freely available ATDI MAPS. The maps are provided for their proprietary software ICS manager or ICS telecom. I am involved in a project to develop a software in c++ or Python to read these GIS maps. There are available libraries in python for GIS maps, but if I am not wrong they work with other file types.

The downloaded maps are zip files containing RGE, RSO and RBL files.

Spain_25m.zip for example contains spain_25m.RGE, spain_25m.RSO, spain_25m.RBL

As I am completely new in this field I dont know how to convert these files to standard GIS file formats.


texts: How to convert a raster output in a stack in R into ArcGIS

texts:
I am modelling the response of a bird (skylark) to various climatic and topographical variables in the program R. Using the predict function I predicted the abundance per pixel (see below for code).
I would now like to export my three rasters into ArcGIS 10.2?
How do I do that?

Can i add the ref coords. using ArcGis or is it easy to do in R?
Is there a faster way predict the density? (I am happy to map on density without the errors initially?

library(raster)
 # Load the raster layers of the predictors for the entire landscape as a stack.
nsw.stack <- stack(list.files(pattern="tif$", full.names=FALSE))
names(nsw.stack) <- c("riv","lake","ndvi","temp","prec","dem")
summary(nsw.stack)

class       : RasterStack
dimensions  : 2129, 1905, 4055745, 6  (nrow, ncol, ncell, nlayers)
resolution  : 0.008926991, 0.008926991  (x, y)
extent      : 55.99709, 73.003, 43.99541, 63.00097  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
names       : riv,     lake,    ndvi,    temp,    prec,      dem

Then I used the predict function to map the data using my best ranked model (t5)

skylark.psi <- predict(t5, type="state", newdata=nsw.stack)

plot(skylark.psi, ylab="psi - GLM")

> skylark.psi
class       : RasterStack
dimensions  : 2129, 1905, 4055745, 4  (nrow, ncol, ncell, nlayers)
resolution  : 0.008926991, 0.008926991  (x, y)
extent      : 55.99709, 73.003, 43.99541, 63.00097  (xmin, xmax, ymin, ymax)
coord. ref. : NA
names       :    Predicted,           SE,        lower,        upper
min values  : 2.987641e-05, 1.275741e-05, 1.293780e-05, 6.899167e-05
max values  :    779.48502,     85.93247,    675.88093,    898.97033


texts:
I have a database of northings and eastings of boundary markers of land parcels that are in a form that I cannot pin on any format. I would like to perform some conversions on them so that they can be used on Google Earth, Maps, Waze etc. Example land parcel boundary markers:

23463.817 E, -54819.146 N
23449.900 E, -54808.635 N
23445.504 E, -54841.976 N
23435.597 E, -54834.499 N
23434.750 E, -54828.132 N

One of the marker is pointing "close" to 3.184487, 101.718079 on Google Earth.

I'm having trouble determining what format are these because it seems the values are way out of range for many known formats. I have no access to the people providing this information at the moment. Can anyone help me figure this out?

The area of concern is in Kuala Lumpur, West Malaysia.


texts:
Rasterio supports processing using window in order to avoid loading the entire image (as can be seen here). It also supports resampling images when using read (as can be seen here). However, I am facing problems when using both approaches together.
My goal is to proccess two images of different spatial (10m and 20m), downsampling 10m to 20m and processing both at 20m using window to avoid load the entire image.
The proccess flux consists in a loop that opens the first image block, downsample the second image (downsample), open the same image block, perform an operation and write the results.
I wrote the following script with small image cases. However I am doing somethink wrong on the resampling/load window for the second image, due to it not present the same shape of the first one. Is it possible to Resample and load a specific block using rasterio considering that the resample block order is important in this case?
import numpy
import rasterio
from rasterio.enums import Resampling

resample_factor = 1/2 # Downsample to 1/2 of the resolution

path_to_img_10m = 'path/to/img10m.tif'
# http://brazildatacube.dpi.inpe.br/Repository/Archive/S2_MSI_L2_SR_LASRC/2019-07/S2B_MSIL1C_20190730T142759_N0208_R053_T20LLN_20190730T161125/S2B_MSIL1C_20190730T142759_N0208_R053_T20LLN_20190730T161125_sr_band8.tif
path_to_img_20m = 'path/to/img20m.tif'
# http://brazildatacube.dpi.inpe.br/Repository/Archive/S2_MSI_L2_SR_LASRC/2019-07/S2B_MSIL1C_20190730T142759_N0208_R053_T20LLN_20190730T161125/S2B_MSIL1C_20190730T142759_N0208_R053_T20LLN_20190730T161125_sr_band8a.tif
path_to_output = 'path/to/out.tif'

# Load 20m profile and block sizes
with rasterio.open(path_to_img_20m) as dataset:
    profile = dataset.profile
    blocks = list(dataset.block_windows())
    height, width = dataset.shape
    result = numpy.full((height, width), dtype=profile['dtype'], fill_value=profile['nodata'])

# Loop on blocks
for _, window in blocks:
    print('window')
    row_offset = window.row_off + window.height
    col_offset = window.col_off + window.width

    # Open image block
    with rasterio.open(path_to_img_20m) as dataset:
         values1 = dataset.read(
            1,
            masked=True,
            window=window
        )
    print(values1)
    print(values1.shape) # (256, 256)

    # Open image block that needs resample
    with rasterio.open(path_to_img_10m) as dataset:
        values2 = dataset.read(
            out_shape=(
                dataset.count,
                int(window.height * resample_factor),
                int(window.width * resample_factor)
            ),
            resampling=Resampling.average,
            masked=True,
            window=window
        )
    print(values2)
    print(values2.shape) # (1, 128, 128)

    # Fictional operation using blocks
    result[window.row_off: row_offset, window.col_off: col_offset] = values1 + values2

# Write result on disc
with rasterio.open(path_to_output, 'w', **profile) as dataset:
    dataset.write_band(1, result)


texts: Why do layers with different CRS overlap fine but cannot be merged by location?

texts:
Perhaps someone could clear this up for me? I have a shapefile with the CRS Selected CRS (EPSG:29902, TM65 / Irish Grid)

I also have a number of XY points which take the CRS, Selected CRS (EPSG:4326, WGS84).

Although different CRS, both of the points map on to the shapefile perfectly.

However, when I wish to join these attributes by location I am unable to as they have different CRS, thus blank points.

When I try to make sure both CRS match, the shapefile and coordinates no longer match up.

I do this through layer properties > data > general function as detailed in many examples online.

Would anyone be able to tell me why the CRS align correctly when different, yet don't when the same?

Also how could I correct this to merge by location?

Edit - Note this is when both CRS are different

While this is when the CRS are the same e.g. (EPSG:29902)

Is it the case that my XY values are long lat, while the shapefile is metres?

Edit
Perhaps access to the data may help?

The shapefile was downloaded from the CSO

http://census.cso.ie/censusasp/saps/boundaries/ED_SA%20Disclaimer1.htm

the file being - 2013 Constituency boundaries:

These are some of the coordinates I'm working with.

 Long   Lat
-7.24   52.64


texts: How do I convert a GeoTiff's neatline metadata from vector format to pixels using GDAL?

texts:
I have a GeoTIFF that produces the following when run through gdalinfo:

Driver: GTiff/GeoTIFF
Files: PAA_1.tif
Size is 14400, 6600
Coordinate System is:
PROJCS["PAA_1_Asia_South_Lambert_Conformal_Conic",
    GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",-13.3333333],
    PARAMETER["standard_parallel_2",46.66666666],
    PARAMETER["latitude_of_origin",-15],
    PARAMETER["central_meridian",125],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
    AUTHORITY["EPSG","9001"]]]
GeoTransform =
  -7964169.743710292, 1177.758187333239, 322.1610457363416
  6367638.541192908, 322.2290109665665, -1177.712898399229
Metadata:
  AREA_OR_POINT=Area
  CREATION_DATE=D:20141215110023-05'00'
  CREATOR=ESRI ArcMap 9.3.1.4050
  NEATLINE=POLYGON ((-5337534.907769799232483 -155620.899077162524918,-5337535.365463367663324 10233860.720439316704869,10253410.032338818535209 10233861.175854712724686,10253410.596936786547303 -155621.599452315276721,-5337534.907769799232483 -155620.899077162524918))

My goal is to use the given neatline in the metadata to crop the image.  I do not want to use gdalwarp because I do not want to rotate or warp the image, just crop it.  Down the road, there may be other imaging tools (other than gdal) that crop the image, so I would prefer converting the given neatline coordinates into pixels relative to the top left corner and then pass those pixels off to something like imagemagick or PIL for image cropping.

Based on information from other answers, I put this together:

import gdal
from geometry import MapToPixel

ds = gdal.Open('PAA_1.tif')
transform = ds.GetGeoTransform()

# parses the neatline string into a list of (x,y) tuples
neatline = ds.GetMetadata()['NEATLINE'][10:-2]
neatline = [(pt.split(' ')[0], pt.split(' ')[1]) for pt in neatline.split(',')]
neatline = [(float(pt[0]), float(pt[1])) for pt in neatline]

# MapToPixel: https://code.google.com/p/metageta/source/browse/trunk/metageta/geometry.py#497
pixels = [MapToPixel(*pair, gt=transform) for pair in neatline]
print pixels

#RESULTS:
# [(665.3036244770069, 5720.95251596183), (2910.358041930786, -2486.5316409600555), (15226.421624842182, 883.213682625752), (12981.36704208365, 9090.698775704966), (665.3036244770069, 5720.95251596183)]

All this does is parse the NEATLINE metadata into x,y coordinate pairs and use the image's GeoTransform data, along with MapToPixel to convert the coordinates into pixel values.  However, none of these values map to anything that looks like a neatline area of the image, and some of them are beyond the width and height bounds of the image.

If I understand correctly, it should be possible to apply some simple linear function to a set of coordinates to convert them into pixel values using the image's GeoTransform metadata.  MapToPixel should be providing such a transformation.  Please correct me if I'm wrong.  I can also settle for hearing about GDAL utilities that could crop the image using NEATLINE without warping or rotating it (somehow preventing gdalwarp from doing this or using a different utility altogether).


texts:
I have a layer in Geoserver called bln_blocks. The layer is in EPSG:3857.

I am doing the follwoing WFS request with a CQL Filter:

  var geoserverURL = geoserverURL + '/ows?';
  var service = 'service=WFS&';
  var version = 'version=1.1.0&';
  var request = 'request=GetFeature&';
  var typeName = 'typeName=' + layerName;
  var special = "&outputFormat=text/javascript&format_options=callback:getJson";
  var filter = '&CQL_FILTER=WITHIN%28geom%2C+POLYGON%28%281489248.8973085624165833+6895049.90960700158029795%2C1489248.8973085624165833+6895577.80283359717577696%2C1489945.19086309685371816+6895577.80283359717577696%2C1489945.19086309685371816+6895049.90960700158029795%2C1489248.8973085624165833+6895049.90960700158029795%29%29%29';
  var urlWFS = geoserverURL + service + version + request + typeName + filter +special;

The polygon in my CQL_Filter is also in EPSG:3857. Is it possible that I specify my polygon in EPSG:4326?


texts: Tableau: Using UK Shapefiles in PostGIS Database

texts:
I'm trying to use Tableau to display data. Some of this data is suitable for mapping, and therefore I'd like to create maps in Tableau to display this data.

Unfortunately, neither the basic view nor the excellent work done for UK Mapping in tableau that is available freely available quite meets my needs. The right combination of Districts, Counties, Unitaries and metropolitans doesnt - as far as I can tell - exist in tableau format. However I believe I've got all the data and bits I need to DIY... but am struggling putting it all together.

I have Tableau Professional 9.3; Postgres 9.4 (PostGIS 2.1) and QGIS 2.12.2.

I have all of my data held in a PostGres database that I've connected to Tableau (and previously to QGIS).

I have all of my shapefiles downloaded from the ONS and imported into a second PostGres database. I'm particularly interested in this one: (https://geoportal.statistics.gov.uk/Docs/Boundaries/Local_authority_district_(GB)_2014_Boundaries_(Generalised_Clipped).zip). These are EPSG:27700 (1936/British National Grid)

I can query either database from tableau or QGIS. I have tried the method here (https://community.tableau.com/docs/DOC-5831) but end up with a blank attribute table in QGIS after performing the first step to convert to WPGS 84 (i.e. it loses all data when it converts). This therefore doesn't work for me; and I'm wary of the method as it strikes me there should be a more elegant solution at the database level.

Tableau wants:

a path;
a poloygon id;
a name;
a latitude;
a longitude;

Therefore, in order to get this information in a suitable form, in the PostGIS database I believe I need to:

Extract the nodes from the 'geom' type in the table created by the imported shapefile (do I use (ST_DumpPoints(GEOM)).Geom)?)
Extract the path from the 'geom' type in the table created by the imported shapefile (do I use (ST_DumpPoints(GEOM)).Path)?)
Extract the Easting and Northing data from the Nodes
Convert the Easting and Northing data from EPSG:27700 to Lat\Long (i.e. WPS:84) (I know I need to use something like the following set long=st_x(st_transform(geom,4326)), lat=st_y(st_transform(geom,4326))
Put these into a new table\view and link to this in Tableau

Here's where I've got to:

--Drop Table IF Exists ctyua_dec_2014_ew_bgc_tblu;

With CTE_PointsAndNodes as
(
SELECT gid
    , ctyua14cd
    , ctyua14nm
    , ctyua14nmw
--  , geom
    , (st_dumppoints(geom)).geom as points
    , (st_dumppoints(geom)).path as path
  FROM (Select * from ctyua_dec_2014_ew_bgc Where ctyua14nmw is Null) as tbltest
  )

Select gid
    , ctyua14cd
    , ctyua14nm
    , ctyua14nmw
--  , geom
    , points
    , path
    , st_x(st_transform(points,4326))
    , st_y(st_transform(points,4326))

INTO ctyua_dec_2014_ew_bgc_tblu

From CTE_PointsAndNodes

I do this; link to it in Tableau... but I get a blank sheet!

So - am I right... and what am I doing wrong...


texts: How to use an existing data dictionary with a new shapefile in pathfinder office?

texts:
when importing a new shapefile in pathfinder office to later be exported to the trimble geoxt, and selecting the option to import with an existing data dictionary does there need to be some common unique ID? What connects the the data dictionary to the new shapefile? Do I need to enter all of the attributes from the existing data dictionary into my shapefile? This would be a tedious task when my data dictionary has more than  40 attributes.

I have already tried to import --> properties --> and chose the 3rd option to import with an existing data dictionary but when I export and then check the file in Terrasync on the trimble, no data is associated.


texts: Shapely WKT dumps and loads does not always preserve rounding_precision

texts:
Shapely's WKT loads and dumps can be useful for rounding overprecise geometries, as pointed out here (Rounding all coordinates in shapely?) and here (Geopandas [Shapely] spatial difference: TopologyException: no outgoing dirEdge found).
I needed that rounding, because I also ran into problems with GeoPandas overlays:
"TopologyException: found non-noded intersection between LINESTRING ...")
The Shapely dumps -> loads trick is very useful, but I'd like to ask about an irregularity in preserving the rounding.  I'm using GeoPandas 0.8 and Shapely 1.7.0.
In the code below, I am tacking on an outer dumps to see easily see the output of loads.  (If this is not the best way to do so, I apologize.  Seems like other options (e.g., list(zip(*geom.exterior.coords.xy)) always give back 1 decimal)
I can dump with expected behavior:
from shapely.wkt import dumps, loads

>>>type(geom)
<class 'shapely.geometry.polygon.Polygon'>

>>>dumps(geom)
'POLYGON ((222595.7000000000116415 -759823.8000000000465661, 222595.8999999999941792...
>>>dumps(geom, rounding_precision=2)
'POLYGON ((222595.70 -759823.80, 222595.90 -760114.10, ...

But sometimes the original precision is still there, as I can see by loading it back:
>>>dumps(loads(dumps(geom, rounding_precision=2)))
'POLYGON ((222595.7000000000116415 -759823.8000000000465661, 222595.8999999999941792...`

This is confusing, as the output of interior dumps is definitely a string:
>>>type(dumps(geom))
<class 'str'>

So it seems like the decimals should be gone (they can't be hiding in the string... can they?).  How does 'loads' get them back?
I am out of my depth, but in my limited testing, rounding_precision=0 seems to be the most robust.  I.e., as opposed to the example with 2 above, the 0 rounding sticks:
>>> dumps(loads(dumps(geom, rounding_precision=0)))
'POLYGON ((222596.0000000000000000 -759824.0000000000000000, 222596.0000000000000000`

Any ideas about what I'm doing wrong, and how I can round and keep the rounding?


texts: Export a QGIS print composer file (.qpt) as PDF from the command line

texts:
I have a workflow where my R script produces csv's and shapefiles. Those csv's and shapefiles then serve as layers in a QGIS project. Every time the R script is run, it overwrites the previous data, and the QGIS project is "dynamically" updated. The QGIS project is all set up such that it never needs to be updated or modified, since it dynamically loads the layers output by the R script.

My current workflow:

Run R script
Open QGIS project file
Open Print composer manager
Select and Open existing print composer layout
Export (print) as PDF
Quit qGIS
View PDF

Desired workflow:

Run R script
View PDF

I can get part of the way there in R by opening the QGIS file in QGIS, by calling

# if on Windows:
shell.exec("my_qGIS_file.qgs")
# or if on MAC OS X:
system(paste("open", "my_qGIS_file.qgs"))

but then I can't figure out how to programmatically (from R) produce the pdf and quit qGIS.
I've read and reread the PyQgis cookbook section on the print composer, as well as the other examples on the internet (including this one),  but those are both beyond my python comprehension, and seem to be way more complicated than I need. It seems like all the examples are programmatically redefining all the print composer every time the script is run. In this workflow, I already have a print composer layout stored in the project that I had created through the GUI, and have no desire or ability to recreate it with code.

If this can't be done directly in R, I think it's possible to invoke a python script from R via a system call, something like:
source("my_python_script_that_exports_my_PDF.py"), where the python script is self-standing.

Also, I have some QGIS plugins installed that bypass having to open the print composer and just export existing composer layouts cut down the workflow by a few clicks. One of the plugins is "Maps Printer", another is "multiPrint"
Perhaps a python or R script can invoke commands from these plugins?

Most critically I would like any proposed solution to run on windows, but bonus points if the suggested solution:

works on both windows and mac
does NOT actually load the full QGIS GUI (i.e., the person running the script does not see the QGIS loading splash screen or anything QGIS related; they just see the PDF magically show up in a specific directory, perhaps with a file name that is "file_name_" + timestamp.PDF).

I am running QGIS 2.8


texts: How to find ZCTA for the 2010 Census block code?

texts:
I'm trying to find ZIP Code Tabulation Areas for the list of 2010 Census block codes.

How to get ZIP Code Tabulation Areas for Census block code?


texts: Openlayers 2.12 + Firefox 15 = OSM Layer not showing

texts:
I know it's due to the new Access-Control-Allow-Origin in OL 2.12, but it is still not functional to me. I've followed instructions on here on stackoverflow.

Firefox version

15.0.1

Apache2 config

In /etc/apache2/sites-enabled/000-default I've set:

<Directory /home/fradeve/public_html/ark-oia/*>
    Header set Access-Control-Allow-Origin *
</Directory>

OL OSM layer config

var osm_layer = new OpenLayers.Layer.OSM('OSM Layer', null, {
    transitionEffect: 'resize',
    tileOptions: {crossOriginKeyword: null},
    resolutions: [...],
    serverResolutions: [...]
})

But OSM layer is fine on Safari, Chromium 20.0.1132.47, but not in Firefox. What am I missing?


texts: Failure to re-project SpatialPolygonsDataFrame to Equidistant Cylindrical in R (units problem?)

texts:
I am trying to re-project a shapefile to match the projection and units o Black Marble nightlights data (link below) which comes in a raster with Equidistant Cylindrical projection, units in meters. However, when I do this re-projection, the extent of the shapefile is way, way too big... larger than the entire world. I think this is because the shapefile is originally in decimal degrees, and the nightlights raster is originally in meters. However, I'm not sure how to alter my re-projection process, to fix the problem.

Specifics: my nightlights data has this projection:

+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0

My district shapefile has this projection:

+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

When I re-project my shapefile (using spTransform), the extent of the shapefile becomes:

extent : 8911228, 9818899, 2933088, 3389212  (xmin, xmax, ymin, ymax)

This is FAR, far larger than the extent of my nightlights data, which coves the entire world, with an extent of -180, 180, -90, 90. So of course, one can't map them together or extract the nightlights data to the shapefile.

Can anyone advise on how to better do this re-projection? (I also tried reprojecting the raster itself using projectRaster -- and precisely the opposite happened, the resulting raster extent was so small as to be a tiny speck within the shapefile.)


texts: Any .NET library to calculate polygon area?

texts: Is there a .NET library to calculate area enclosed by polygon made up of geo coordinates (lat longs)?


texts: Making the categories layer style which shows features overlapping another layer in a different color

texts:
I have tried with this new QGIS function from expression engine but it is not working.
Here is the query I have tried:
overlay_touches('Master_3307To3339',filter:= within(geometry(@parent),@map_extent))

I don't want it to check all parent layer features every time I create, so make it check only map extend visible parent layer features. Is it possible?


texts: Add and Calculate Fields of Iterated Annotation Classes

texts:
I'm have difficulty with my script Arcmap v10.5.1.  I get a syntax error on line 8.  My question is:  Can I calculate the Added field from an existing field within an iterated value?  If so, what's my syntax error preventing me from doing so.  If not, how should I proceed?  my workspace is a feature dataset

import arcpy, os, sys

arcpy.env.workspace = arcpy.GetParameterAsText(0)

for fc in arcpy.ListFeatureClasses():
    arcpy.AddMessage("Adding fields and updating " + fc + "...")
    arcpy.AddField_management (fc, "KM_FontSize", "DOUBLE")
    arcpy.CalculateField_management(fc, "KM_FontSize", "!FontSize!")
    arcpy.AddField_management (fc, "KM_VerAlign", "SHORT")
    arcpy.CalculateField_management(fc, "KM_VerAlign", "!VerAlign!")
    arcpy.AddField_management (fc, "KM_HorAlign", "SHORT")
    arcpy.CalculateField_management(fc, "M_HorAlign", "!HorAlign!")
    arcpy.AddField_management (fc, "KM_XOffset", "DOUBLE")
    arcpy.CalculateField_management(fc, "KM_XOffset", "!XOffset!")
    arcpy.AddFieldmanagement (fc, "KM_YOffset", "DOUBLE")
    arcpy.CalculateField_management(fc, "KM_YOffset", "!YOffset!")
    arcpy.AddMessage("Finished updating " + fc + "...")

And here's the error message:

Traceback (most recent call last):
  File "H:\Python Scripts\KM_Utility\Create Anno Fields.py", line 10, in <module>
    arcpy.CalculateField_management(fc, "KM_FontSize", "!FontSize!")
  File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\management.py", line 3661, in CalculateField
    raise e
ExecuteError: ERROR 999999: Error executing function.
Syntax error
Failed to execute (CalculateField).


texts: Adjust labels in CartoDB choropleth legend

texts:
CartoDB allows the user to change the legend values for the left label and the right label. However, I would like to add values to the three middle color values. Is there a way to do use within the the CSS or SQL options in CartoDB?

Alternatively, CartoDB allows the user to select 3, 5 or 7 buckets to define the color intensity for a group of values. Is it possible for the color ramp to be on a continuous scale? I.e. the color of each polygon isn't allocated to one of the 3, 5, or 7 colors but rather is determined based on where the value falls within the entire value range?


texts: How to cache online tiles on a development machine?

texts:
I am developing a slippy map web application using openlayers. Right now my tiles come straight away from OSM and google, etc. For the moment I always work on the same area for development purposes and I would like to have those tiles cached on my dev machine so I can work on the web app offline yet still see some map and to reduce network consumption.

What I would like is to set up some sort of proxy which would just cache tiles I watch and if they are not on the local machine download them or show a dummy tile if I am offline.

I thought tilecache was the tool for this, but I get quite confused trying to set it up. It seems more appropriate for mapserver generating tiles from map data. I know firefox does some caching, but I would like some more control and be able to clear that cache to make sure I see the last update of my work.

What is the most straightforward way to achieve this?


texts:
I am in need of a geocoding api to use in my application. I tried Google maps api but it restricts 10 requests per second. The thing is place names will be there in my csv file and it should geocode for those place names in csv file. Since there are around 900 names in my csv file, only 10 were being geocoded using Google maps api.

I would like to know if using geonames or other opensource geocoding apis would be an option, with respect to number of requests per second and maximum number of requests per day.


texts: Select rasters from a folder that intersect the tiles of a polygon shapefile

texts:
Ideally I want to create a model in ArcGIS ModelBuilder that will iterate a folder of rasters and output the ones that intersect polygon tiles in another shapefile.

I have the polygon shapefiles saved to a file location and can read them into the model, I have an iterator that is reading through the folder of rasters.
What tool (if existing) will allow me to output the rasters that intersect the polygons in the shapefile?

Tile Index

Selection from Tile Index: Desired Rasters to pull from folder

Attribute table of selection from Tile Index.  Name column highlighted matches naming convention of Raster DEMs

Raster DEM folder, expanded in Arc Catalog

EDIT:
    I was able to run a batch process since the names of the polygon tiles matched the correlated names of the raster DEMs:
    Create a file list to be used for batch processes
Create a file list (export attribute table as dbf in ArcMap)
Open up file list in excel (delete all columns except common name field)

Column B - type in copy “ ” type in location of files you want to copy

Ex: copy Y:\Data\RasterDEMs\

Column C - copy/paste special Column A  AS VALUES
        Delete column A contents

Column D - .*

Column E - type in location for files to be saved

Column F - concatenate columns B - E
Ex: =CONCATENATE(B1,C1,D1," ",E1)

Sheet 2, Column A - copy/paste special Column F  AS VALUES

Save Sheet 2 as tab delimited .txt
close excel and open Sheet2.txt
Save As with extension .bat (batch file)
double click that file and it should open the CMP box and start running

[ “ “ = Space]

For those who wanted to see the data I've added the images.

I'm running ArcGIS 10.3.1
extensions:
Spatial Analyst
3D Analyst


texts: Change Greenland's continent using rworldmap

texts:
I am using the answer provided here and exactextractr::exact_extract to find the average of my raster file by continent.  I particularly like the rworldmap option because I am looking to analyze by North America/Eurasia (rather than North America/Europe/Asia),  and rworldmap combines Europe and Asia into Eurasia. However, I'd like Greenland to be part of North America and not Eurasia.  How can I change its continent?
library(exactextractr)
library(raster)
library(sf)
library(rworldmap)
library(rgeos)
library(maptools)
library(cleangeo)

sPDF <- getMap()
#### I assume the changing continent part would go here
sPDF <- clgeo_Clean(sPDF)  ## Needed to fix up some non-closed polygons
cont <-
  sapply(levels(sPDF$continent),
         FUN = function(i) {
           ## Merge polygons within a continent
           poly <- gUnionCascaded(subset(sPDF, continent==i))
           ## Give each polygon a unique ID
           poly <- spChFIDs(poly, i)
           ## Make SPDF from SpatialPolygons object
           SpatialPolygonsDataFrame(poly,
                                    data.frame(continent=i, row.names=i))
         },
         USE.NAMES=TRUE)
cont <- Reduce(spRbind, cont)
noant <- cont[cont$continent!="Antarctica",] # antarctica breaks spTransform
myraster<-raster('filepath')
trans <-spTransform(noant,crs(myraster))
nam_eur<- st_as_sf(trans)
raster_mean<-exact_extract(myraster,nam_eur,'mean')

plot(cont, col=heat.colors(nrow(cont)))

shows Greenland as part of Eurasia:


texts: Is it possible to distribute custom python expressions along with .qgs file?

texts:
I have a QGIS project I'm working on with a few other users. The print composer and map tips use custom python expressions in the expression builder.

I know the expressions are store in .qgis2/python/expressions but I want to find a reasonable way to distribute and update those expressions so that when someone alters on it syncs to all the others. It seems odd to me that these can't be embedded somehow in the .qgis file (I understand there's security concerns, but a model similar to macros in microsoft products could be used...)


texts: Requesting tiles without using XML file (GDAL C++)

texts:
I have been able to request a specific tile with my XML through GDAL lib. However, is there is another way of requesting the tiles from a web server
without using XML?

I want to do it directly through my C++ code. I will appreciate any directions or examples.


texts: How to apply cdo mask command to loop over multiple .nc files in Python?

texts:
I want to use the cdo command for masking the .nc file (the command below):
cdo div input.nc maskfile.nc outputfile.nc
The command above should be applied over multiple .nc files (looping) in python, can anyone assist in how to apply it?
Also the same case with weighted mean command from cdo,
cdo fldmean input.nc fldmean.nc


texts:
I'm trying to make a test application using ArcGIS Web AppBuilder, I built the app downloaded it,and unzipped it, tried to launch it using IIS not working and the same thing by nodejs's http-server, it loads for moments then I get this error message "GET /headless-content.js.map" Error (404): "Not found", I tried to locate this file but never found it. Also in the browser I noticed that it stop working after a message 'redirecting http to https' so I edited all urls in the config.json to http but still having same issue.
Any tips how to find out what causes the problem?


texts: Reclassify and download raster in Google Earth Engine

texts:
Playing with Google Earth Engine i'm trying to download Tree Cover Loss data for a country using the Hansen dataset ("UMD/hansen/global_forest_change_2019_v1_7").
I would like to download tree cover loss data from 2014 to 2019 for Cameroon using this script:
// select the Global forest Watch dataset
var gfc = ee.Image('UMD/hansen/global_forest_change_2018_v1_6')

// Create a variable and extract just the "lossyear" band from the GFW dataset
var tclYr = gfc.select('lossyear')

// Create a variable and select only the years from 2014 to 2019
var tclAll = tclYr.gte(14).and(tclYr.lte(19))

// Clip the tclALL on ivory coast (geometry created manually as rectangle)
var tcl_cameroon = tclAll.clip(geometry)

// Add layer to the map, without visualization parameters
Map.addLayer(tcl_cameroon)
// Map.addLayer(ESA_forest)

// Export image to Google Drive
Export.image.toDrive({
  image: tcl_cameroon,                 //  <----- put the image you want to download
  description: "tcl_cameroon",    //  <----- download with this filename
  scale: 30,                        //  <----- set the pixel size in meters
  //region: geometry
});

I realized that it takes too much and maybe could be faster if:

i reclassify years 14 to 19 as one class
i reclassify NODATA for the rest.
Do you think the process will be faster?  and filesize lighter?
How could i do it?


texts: Downloading NASA Earth Exchange Global Daily Downscaled Projections for single lat lon

texts:
I have a lat lon - 23.14, 83.14 and I need to download the NASA Earth Exchange Global Daily Downscaled Projections (NEX-GDDP) data from the below website only for this lat lon. I am looking to download the Tmax.

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

I have not worked with Earth Engine before.

How do I go about extract and saving the daily time series of Tmax for this particular lat lon?


texts: Renaming columns using geopandas - issues with data table in Arc

texts:
What I'm trying to achieve:

Read in an existing ESRI polygon shapefile
Drop an attribute
Dissolve by another attribute
Rename columns
Reorder columns
Write out as new shapefile

What works: Drop column and dissolve all work - I've run these steps and written to file. This opens fine in ArcMap, I can identify things, open the attribute table etc. No issues.

Problem: When I rename (not got to the reorder step without error yet) the columns, everything runs fine and I write to file. When then opening the shapefile in ArcMap (10.3.1) I can query layers using the Identify cursor which shows my renamed attributes as they should be BUT when trying to open the attribute table, I get the following message in a window:

Loading table data

Could not load data from the data source. If you can correct the problem, press the refresh button to reload data. Possible problems can include bad network connection, invalid file, etc.

A column was specified that doesn't exists.

A column was specified that doesn't exists.

So, it looks like something is missing or my changes are not being mapped across the shapefile components.

Given a mock version of the data:

# Make a mock geopandas dataframe
import geopandas as gpd
from shapely.geometry import Point, Polygon
from fiona.crs import from_epsg

data = gpd.GeoDataFrame()

# Coordinates of some locations in British National Grid (epsg: 27700)
coordinates1 = [(19.950899, 60.169158), (20.953492, 60.169158), (20.950958, 61.169990), (19.953510, 61.170104)]
coordinates2 = [(21.950899, 60.169158), (22.953492, 60.169158), (22.950958, 61.169990), (21.953510, 61.170104)]
coordinates3 = [(19.99, 60.3), (20.23492, 60.3), (20.23492, 60.8), (19.99, 60.8)]

# Create a Shapely polygon from the coordinate-tuple list
poly1 = Polygon(coordinates1)
poly2 = Polygon(coordinates2)
poly3 = Polygon(coordinates3) # (sits inside polygon 1)

# Add polygon to dataframe
data.loc[0, 'geometry'] = poly1
data.loc[1, 'geometry'] = poly2
data.loc[2, 'geometry'] = poly3
data.loc[0, 'value'] = 1
data.loc[1, 'value'] = 2
data.loc[2, 'value'] = 3
data.loc[0, 'CLASS'] = 'A'
data.loc[1, 'CLASS'] = 'B'
data.loc[2, 'CLASS'] = 'A'
data.loc[0, 'Hazard_Pot'] = 'High'
data.loc[1, 'Hazard_Pot'] = 'Low'
data.loc[2, 'Hazard_Pot'] = 'High'
data["contaminat"]="something"

# Define CRS
data.crs = from_epsg(27700)

## Plot it
#import matplotlib.pyplot as plt
#ax1=data.plot()
#data.loc[[2],'geometry'].plot(ax=ax1, color='red') # plot on top of others
#plt.show()

We'll now dissolve the polygons based on Hazard_Pot and then rename and reorder things. The problem is with the renaming... What am I missing?

# Drop a column (if writing out file here, does as expected)
data=data.drop(['CLASS'], axis=1)

# Do the dissolve (if writing out file here, does as expected)
data=data.dissolve(by='Hazard_Pot', as_index=False)

# Rename <<< won't work in ArcGIS
data=data.rename(index=str, columns={"Hazard_Pot":"Potential",
                        "contaminat":"Contam."}) # <<< EDIT: notice the punctuation (noticed following the accepted answer)

# Reorder
data=data[['value',
           'Contam.',
           'Potential',
           'geometry']]

# write to file (creating a new CRS based on the old as per https://gis.stackexchange.com/questions/204201/geopandas-to-file-saves-geodataframe-without-coordinate-system?noredirect=1&lq=1)
ofile='./rename_shape.shp'
data.to_file(ofile, driver='ESRI Shapefile')


texts: Trying to get the max value of two bands and put in to array or variable in GEE

texts:
I'm trying to do an array with the maximum value of bands obtained by Reducer.max() with the objetive to get the max value of two bands with Math.max() but in the moment that I tried to make this, obtain a NaN in the console.

//MAXIMO B1 FECHA 1

var IMGLandsat5= ee.ImageCollection ('LANDSAT/LT05/C01/T1_TOA')
.filterDate ('2000-10-01', '2001-03-31')
     .filterBounds (geometry)
  .filterMetadata ('CLOUD_COVER', 'Less_Than', 20);
var Landsat5Filtro = ee.Image(IMGLandsat5.median());
var Landsat5Clip = Landsat5Filtro.clip (geometry);

var b1_1max = Landsat5Clip.select('B1');
var b11max = b1_1max.reduceRegion({
  reducer: ee.Reducer.max(),
  scale: 30,
  maxPixels: 1e10
});

var b11Max = b11max.get('B1');
print (b11Max);

//MAXIMO B1 FECHA 2

var IMGLandsat52= ee.ImageCollection ('LANDSAT/LT05/C01/T1_TOA')
.filterDate ('2010-10-01', '2011-03-31')
     .filterBounds (geometry)
  .filterMetadata ('CLOUD_COVER', 'Less_Than', 20);
var Landsat5Filtro2 = ee.Image(IMGLandsat52.median());
var Landsat5Clip2 = Landsat5Filtro2.clip (geometry);

var b1_2max = Landsat5Clip2.select('B1');
var b12max = b1_2max.reduceRegion({
  reducer: ee.Reducer.max(),
  scale: 30,
  maxPixels: 1e10
});

var b12Max = b12max.get('B1');
print(b12Max);

// ARRAY PARA OBTENER MAXIMO DE LAS DOS FECHAS
var num = [b11Max,b12Max];
var max_total = Math.max.apply(null, num);
print(max_total);


texts: Geodataframe not plotting properly in geopandas

texts:
I am playing around with the gdf.plot method which is supposed to facilitate plotting geodataframes using geopandas. The reason is that I would like to avoid jumping from QGIS to python while developing scripts and just focus on the programming.

Unfortunately for a simple polygon (granted it has a hole inside) this method seems to have problems. On QGIS I get the following plot:

On matplotlib I get the following plot:

I am quite sure that the QGIS version is correct.  The geojson file for which this all applies contains the following:

{
"type": "FeatureCollection",
"features": [
{ "type": "Feature", "properties": { }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 1.0, 2.0 ], [ 1.0, 1.5 ], [ 0.5, 1.5 ], [ 0.0, 1.5 ], [ 0.0, 0.0 ], [ 1.0, 0.0 ], [ 1.0, 0.5 ], [ 1.5, 0.5 ], [ 2.0, 0.5 ], [ 2.0, 2.0 ], [ 1.0, 2.0 ] ], [ [ 1.0, 1.5 ], [ 1.0, 1.0 ], [ 1.5, 1.0 ], [ 1.5, 1.5 ], [ 1.0, 1.5 ] ] ] } }
]
}

Why does geopandas not recognize that there must be a hole in there?


texts: Is it possible to define a No SQL model for spatial data?

texts:
Possible Duplicate:
  Is there any way I can use a Key-Value store for geospatial data?

I read somewhere about No SQL being a solution for large databases management, the kind that google uses for their operations, are there any examples of a spatial model data with NO SQL?


texts: Automating WFS to Shapefile Export using PyQGIS?

texts:
I'm trying to write a simple script to automate the download of several hundred WFS services and create shapefiles. All of these services have the same URL, e.g.: https://fbinter.stadt-berlin.de/fb/wfs/data/senstadt/s_Induktionsschleife while only the last word is changing from file to file. The shapefiles are supposed to be named with a certain structure based on that last word of the URL, e.g. Senat_Induktionsschleife_25833.shp.

Unfortunately, it only creates empty shapefiles, the WFS services fails to load. If I print and copy the created URI and add the WFS manually, it works, so the URI is created correctly.

id = "Induktionsschleife"
uri = r"https://fbinter.stadt-berlin.de/fb/wfs/data/senstadt/s_" +id + "?srsname=EPSG:25833&typename=fis:s_" +id + "+&version=2.0.0&request=GetFeature"
outputPath = '/Users/XXXX/Desktop/WFS-Download/' + 'Senat_' + id + '_25833.shp'

lyr = QgsVectorLayer(uri, "Temp", "WFS")

pfad = outputPath
QgsVectorFileWriter.writeAsVectorFormat(lyr, pfad, 'utf-8', lyr.crs(), 'ESRI Shapefile')

name = pfad.split('/') [-1].replace('.shp', '')
layer = QgsVectorLayer(pfad, name, "ogr")
QgsProject.instance().addMapLayer(layer)

I'm very new to Python and am not able to find out what's the problem with my code.

Do you have any suggestions?

I'm using QGIS 3.6 on MacOS 10.14.


texts: Replacing polygon within composite relationship without deleting tabular record?

texts:
I have a geodatabase that contains a table, a point feature class, and a polygon feature class.

Each record in the table has both a related point and a related polygon, so all one-to-one relationships. The geodatabase is currently set up with composite one-to-one relationships. My issue with this is that we sometimes need to replace the polygons with better data. If a user deletes the old polygon before adding the new one, the corresponding table record is also automatically deleted.

It is important that there is referential integrity in this dataset, but I'm wondering if it might make more sense to maintain the relationship classes as simple, and then explicitly test for referential integrity during the QAQC process. We will eventually be allowing partners to edit replicas of this database, and my concern is that this will cause some problems with users.

Is there an efficient way to replace a polygon within a composite relationshp without deleting the tabular record?


texts: Invert the colour gradient inside a buffer

texts:
I'm struggling to invert the direction of the gradient inside a buffer. I've been playing around with the rotation inside the symbology menu but it's leading me nowhere.

Basically, I need the gradient to go from shore to the ocean. The shore is materialised by the brownish line. What I've done so far is to compute the symbology menu like this:

But the gradient always start from the centre to the outskirts of the buffer. Is there a way to have it homogeneously go from the "top" of the buffer (i.e. from the line), to the bottom of the buffer? Best option would be to have the gradient go from dark close to the shore to clearer close to the bottom of the buffer.


texts: What is the best way to break up a shapefile into discretely sized tiles?

texts:
TL;DR version:  I need to break up road vectors in QGIS to 4,033km2 tiles to import into an UnrealEngine landscape.
I am evaluating a plugin the imports GIS data into UnrealEngine for a simulation. One of the best things about this plugin is that it can import shapefiles to create landscape splines (that's roads and rivers and such for all the non-game programmer types, since this S.E. is about GIS, not UE). It's really great, but for landscapes greater than 16km2 or so in UE, you need to import the data into landscape tiles.
I'm using QGIS to grab the overall data area I want, but I need it broken up into tiles. The topography DEM information was easy enough to split into tiles, but I am having a hard time breaking up the shapefile data. Does anyone here know of any tools for taking a larger-area shapefile and breaking it up into chunks of a specific size?


texts: GeoEvent streaming service labeling in ArcGIS Online or Portal

texts: Is it possible to label streaming services within AGOL or Portal?  There is not a create label option available after adding a streaming service to AGOL or portal.


texts: DEM to RGB-GeoTIFF (land-/waterarea colored separate)

texts:
I want to create a colored basemap out of a ETOPO dataset.

Is there a way to color the land (masked by a landpolygon) and water area separate and combine them later. The result should be one GeoTIFF (RGB) with land and water areas in different color gradients.
To use one color gradient with the break at 0 meters isn't enough.

The software doesn't really matter. I prefer QGIS, but ArcMap is also ok.
I know how to do that with GMT (Generic Mapping Tools).

Is there maybe a possibility to create GeoTIFFs out of the PS-Files with same resolution as the origin ETOPO dataset?


texts: How to set all raster pixels within a shapefile to null?

texts:
I'm working with ArcGIS Pro 2.4.0. I have a raster dataset consisting of several plots of crops. Within each plot, there is a small ground target that has differing reflective values than the surrounding crop. Since I am trying to extract the mean pixel value of the plots, I don't want the higher reflective values of the targets to skew the mean.

Here is a plot I am working with:

"A" denotes the main plot boundary, whereas "B" denotes the boundary around the reflective target.

How can I set all the pixel values within B to null? I've tried a SetNull equation, but typically the SetNull procedure involves setting all pixel values above/below a set value to null. If the solution lies within building a conditional statement, I'm not quite sure how to do that.


texts: Imported WMS and WMTS did not overlap in ArcMap

texts:
I've created WMS and WMTS of the same city then imported them to ArcMap via service capabilities copy link location but they did not overlap.

P.S. the one at the top of the layers is the tif itself and perfectly overlaps with WMS.

What's wrong ?


texts: Shifting multipart labels for an entire layer results into disappearing labels

texts: When I attempt to move a multipart label, all of the labels asides the one I am moving disappear. How can I manually move a label and preserve all of the corresponding multipart labels using QGIS 3.60? Bear in mind i am doing this for an entire layer.


texts: Can a table (.dbf) be encoded in two different systems?

texts:
I am having a *.dbf-file.

When I load it into QGIS (delimited text-file) I am having trouble getting all the characters right.

Either I load it in encoding=windows-1252 - then the headers will be shown correctly, but the features wont or I load it in utf-8 and then the features will be shown correctly, but the headers wont.

I am talking about "æøå" -characters.


texts: Is there a way to label a layer using its file name or layer name?

texts:
Is there a way to label a layer using its file name or layer name (if I'd change it in Layers window TOC)?
I tried some functions but I'm not good on that.

I have lots of shapes (lines) and it would be easy to identify them this way. Of course, this name could be a table attribute but I'm trying a alternative way (maybe more dynamic).


texts: How to compare two sets of points in QGIS?

texts:
I have two point layers with XYZ data; one set are points measured in the field and the second set is a Lidar data set.

I want to check the accuracy of the Lidar in comparison to the measured data.
Is there a way to do so without creating grids? That is, comparing the vectors and not the rasters?

Note: I found similar questions but either there was no answer or it referred to ArcGIS for Desktop.


texts: Why "new" dataset table is deleted when trying to update existing dataset table in Esri Insights workbook?

texts:
I have an established Insights workbook, with designated cards.

The data being used comes from "Table A" when I add "Table B" (updated data) and attempt to 'Update' Table A - which should be a very simple process according to all of Esri's documentation - Table B is deleted from the table of content and Table remains the same.

Both tables are formatted identically and contain the exact same fields.

Why is a "new" dataset table deleted when trying to update an existing dataset table in an Esri Insights workbook?


texts: Images exported from ArcGIS Online as File Geodatabase inaccessible to ArcGIS for Desktop?

texts:
I am a new user of ArcGIS Online.

I have a system that includes features and images of the objects and I use Esri ArcGIS for Desktop to perform data collecting.

Before exporting the whole database as FGDB I use attachment management to download the images as FGDB zip file.

How can I reach the images in the exported FGDB file using ArcGIS for Desktop?


texts: pyQGIS: Which attribute change signal should I listen to?

texts:
I want to receive a signal when changes have been made to a layer, regardless of the type of editing (add features, remove features, change attributes). When I receive this signal I export the layer's data to an excel sheet. I've looked at these signals:

layer.dataProvider().dataChanged: Well, nothing seems to happen here when I start and end an editing session in the layer.

layer.attributeValueChanged: seems to be emitted before changes are committed to the data provider. When I export to excel, the columns for edited features are shifted to the right: If the un-edited data is in columns A..D, the edited data is in columns E..H instead. Is this the "editing buffer representation"?

layer.beforeCommitChanges: same outcome as above: edited features don't occupy the same columns as the un-edited ones.

What I want is probably something like layer.afterCommitChanges, but I can't find anything like that. There are step-by-step signals such as committedAttributesAdded, committedFeaturesAdded and so on but I don't think  that catching all of those is the way to go.

I'm using XlsxWriter to write the excel file. Maybe the problem is related to this choice?


texts: How can I optimize performance for QGIS 2.14 vector processing?

texts:
I am attempting to do a Count points in polygon routine for a clean looking 'density' grid.  I have around 800,000 points over a large polygon grid.

My workstation has 24 cores and 128gb RAM, and yet the process will often take twelve hours or longer to complete.  A look at the processor affinity shows 4 CPUs out of the bunch being hit at any one time for a moment, punctuated by extended periods where nothing appears to happen, and the output filesize does not increase.

My data are in the same SpatiaLite db, with indexes built.

The output does eventually complete, and gives me what I expect.  However, are there any tweaks or options I could explore in QGIS to improve this speed or distribute the processing across more CPUs?


texts: Disconnect processing.run() from source data in QGIS standalone OR How do I remove a file system object from OGR?

texts:
I am using processing.run("qgis:rastersampling") and need a way to disconnect the connection to the corresponding input vector file after each sampling operation. This is what I am thinking of:
#settings
QGIS_install="C:/QGIS"
import sys
sys.path.append(QGIS_install+'/apps/qgis-ltr/python/plugins/processing')
import qgis.core as gc
import processing
from processing.core.Processing import Processing
gc.QgsApplication.setPrefixPath(QGIS_install, True)
qgs = gc.QgsApplication([], False)
qgs.initQgis()
Processing.initialize()

#loop
for i in range(10000):
    InputVector="C:/vector"+str(i)+".gpkg"
    InputRaster="C:/raster.tif"
    OutputFilePath="C:/output"+str(i)+".gpkg"
    processing.run("qgis:rastersampling", { 'COLUMN_PREFIX' : 'xy', 'INPUT' :  InputVector , 'OUTPUT' : OutputFilePath, 'RASTERCOPY' : InputRaster })

    #***at this point in the loop, I am looking for something that closes/deletes the connection of the input vector file to QGIS***

#some hours of other processing

#exit
qgs.exitQGIS()

The problem I am facing is that after several iterations of the loop my RAM is insufficient and Spyder crashes.
Somehow the loop stuffes my memory full of data, and using .exitGIS() would make no sense at this part of the script as it a loop, so this is only possible at the very end of the code.
I do not need the vector input anymore for the further calculations.
The .gpkg-wal and .gpkg-shm files of the opened input vectors accumulate in my folder, so I assume the processing.run() leaves connections opened to all vector files once they have been used as Input for the sampling. gc.QgsProject.instance().removeAllMapLayers() does not help in this case.
=> I am using Python 3.8 in Spyder 5, working with QGIS 3.16
Edit:
It seems as if it has to do something with the "ogr file system", as I get an error if I try to use the same input vector twice for the rastersampling. The error is:
OGR error: A file system object called 'C:/vector1.gpkg' already exists.

So the question could be "How do I remove a file system object from OGR?"


texts: Append expression to 'Select by expression' window using PyQGIS

texts:
I have a PyQGIS script that creates expressions and I would like to append them to the  'Select by expression' window. The script is executed from the Python console. At the moment the script copies expressions to the clipboard and then I paste the expressions manually. However, I wonder if it's possible to append text directly from the script.

The script:

# Script that is supposed to help with selection of features
# with certain attributes
# 1. Make your selection
# 2. Run script
# 3. Point at the feature of interest, observe that
#    marker has changed.
# Script will then look for an attribute that is ~NULL.
# This attribute value is then copied to clipboard as
# an expression ready to use for further processing

from qgis.gui import QgsMapToolIdentifyFeature
from PyQt5.QtGui import QColor
from qgis.PyQt.QtWidgets import QApplication

attribs = ['VAGBETECKN', 'HUVLEDNAMN', 'NYLEDNAMN','VAGNR', 'GATUNAMN', ]
attribs.append('RLID')  # default, always exists

def onFeatureIdentified(f):
    #print('Feature RLID: {}'.format(f['RLID']))
    for attrib in attribs:
        val = f[attrib]
        if val != NULL:
            if not isinstance(val, str):
                val = str(val)
            #print('Value found. {}:{}'.format(attrib, val))
            break

    #print('Using: {} = {} for RLID {}'.format(attrib, val, f['RLID']))
    s = "(\"" + attrib + "\"='" + val + "')"
    clipboard = QApplication.clipboard()
    clipboard.setText(s)
    print('To clipboard: {}'.format(s))

layer = iface.activeLayer()
mc=iface.mapCanvas()
mapTool = QgsMapToolIdentifyFeature(mc)
mapTool.setLayer(layer)
mc.setMapTool(mapTool)
mapTool.featureIdentified.connect(onFeatureIdentified)


texts: Insert empty VIEWS into QGIS from PostgreSQL using PyQGIS

texts:
The Aim:
I am trying to insert a number of empty VIEWS into QGIS using PyQGIS. These  VIEWS have been created and stored in a PostgreSQL database.
I want to also ensure i can successfully re-open the project at a later date whilst some of the layers are still empty or have become empty. These layers can at times be empty depending on what has been happening recently in the database.
The Problem:
QGIS does not like it when these layers are empty on inserting them, or when attempting to re-open the project. It either fails to insert the empty layer, or it informs the user that they are 'bad layers' and need to be deleted in order to continue. The layers can be inserted successfully when they contain data.
Python code used:
Note: "PlatformID" is the primary key.
uri = QgsDataSourceUri()
uri.setConnection(host, port, db_name, user, password)
uri.setDataSource("Data", "FieldOfViewIndex", "geom", "", "PlatformID")
vlayer = QgsVectorLayer(uri.uri(False), "FieldOfViewIndex", "postgres")
QgsProject.instance().addMapLayer(vlayer)

Common error messages:

WARNING    Geometry type and srid for empty column geom of "Data"."FieldOfViewIndex" undefined.
WARNING    invalid PostgreSQL layer

In the PostgreSQL geometry_columns VIEW, their SRID is 0 and the geometry type is GEOMETRY.


texts: Generating a Map (like Google Maps) using HGT files

texts:
I'm doing a programming project. It consists of generating a Map (like Google Maps) from height files .HGT

Currently, I can read each file .hgt (each .hgt file have 1201 * 1201 elevation data). For each point I assign a color depending on the elevation. And from this I can generate an image.

My problem Is Not reading the hgt file or generate the image.

To be exact, my problem is to generate the structure of the map and the square projection (by levels). The image below makes it clear what I want to do.

Each tile of each level will have a size of 256px * 256px.

How can I use the .hgt files to create a square projection? (As I said: I'm not an geographic information systems expert)

I have the idea of generate a big image with all .hgt files. After that, I will cut it into into equal pieces of increments of 2^(2x). But the main problem is the square projection...

Additional information: I am not allowed to use third-party software. It is a complete programming project. I have a cluster of 30 machines to do all the work, so you don't have to worry about the calculation power.


texts: Question about coordinates from Google Earth Pro to ArcMap

texts: I have a project where I must georeference an archaeological site. I have only a picture from Google Earth on which there are drawn some tombs. I am trying to goereference the image but the picture is coming from an old version of Google Earth and thus the landscape (when I view it in a basemap) has changed a lot. Because of this I cannot find some certain points to georeference the picture. I managed to find old satellite photos from Google Earth Pro and find some points (some trees for example) and extract the coordinates. But, as a noob I am confused about the projections on ArcMap. I can't figure out what coordinate system Google Earth Pro uses in order to use the same with my layers on ArcMap and, in result, to locate the exact points and georeference my picture. Maybe this method is totally wrong but I don't know what else to do. I have no coordinates regarding the site, only this old satellite picture.


texts: Can QGIS2 / GDAL convert to .ECW without corruption?

texts:
I have qgis 2.0.1 and it can read *.ecw files that I download.

However, when I convert say a *.jpg or *.tiff to *.ecw the output is corrupt.

Has anyone seen this problem?