Raster erstellen, in dem jede Zelle die Entfernung zum Meer aufzeichnet?

10

Ich möchte ein Raster mit einer Auflösung von 25 Metern × 25 Metern erstellen, wobei jede Zelle den Abstand zur nächsten Küste enthält, der aus der Mitte der Zelle berechnet wird. Dazu habe ich nur ein Shapefile der neuseeländischen Küsten .

Ich habe versucht, Dominic Royes Tutorial zu folgen, um es in R zu machen, was funktioniert ... irgendwie. Es ist in Ordnung bis zu einer Auflösung von etwa 1 km × 1 km, aber wenn ich versuche, den RAM-Wert zu erhöhen, wird der auf meinem PC verfügbare Speicherplatz (~ 70 GB RAM erforderlich) oder ein anderer, auf den ich ebenfalls Zugriff habe, deutlich überschritten. Wenn ich das sage, denke ich, dass dies eine Einschränkung von R ist, und ich vermute, dass QGIS eine rechnerisch effizientere Methode zum Erstellen dieses Rasters bietet, aber ich bin neu darin und kann nicht genau herausfinden, wie es geht.

Ich habe versucht, Raster mit QGIS zu erstellen. um es in QGIS zu erstellen, aber es gibt diesen Fehler zurück:

_core.QgsProcessingException: Quellschicht für INPUT konnte nicht geladen werden: C: /..../ Coastline / nz-Küstenlinien-und-Inseln-Polygone-topo-150k.shp nicht gefunden

und ich bin mir nicht sicher warum.

Hat jemand Vorschläge, was schief gehen könnte oder eine alternative Möglichkeit, dies zu tun?

Bearbeiten:

Das Raster, das ich produzieren möchte, hätte ungefähr 59684 Zeilen und 40827 Spalten, so dass es sich mit dem jährlichen Wasserdefizit-Raster von LINZ überschneidet . Wenn das produzierte Raster größer ist als das jährliche Wasserdefizit-Raster, kann ich es in R schneiden ...

Eine Sache, die ich für ein potenzielles Problem halte, ist, dass das Shapefile der neuseeländischen Küste eine große Menge Meer zwischen den Inseln aufweist, und ich bin nicht daran interessiert, die Entfernung zur Küste für diese Zellen zu berechnen. Ich möchte wirklich nur die Werte für Zellen berechnen, die ein Stück Land enthalten. Ich bin mir nicht sicher, wie ich das machen soll oder ob es tatsächlich ein Problem ist.

André.B
quelle
1
Führen Sie dazu ein Skript aus? Oder verwenden Sie die Tools in QGIS? Etwas zu überprüfen, auch wenn es so klingt, wie es sollte - überprüfen Sie, ob die Datei tatsächlich dort existiert, wo Sie es sagen ... Überprüfen Sie auch, ob Sie Lese- und Schreibzugriff auf diesen bestimmten Ordner haben.
Keagan Allan
Ich benutze derzeit die Tools, bin aber sehr daran interessiert, das Skript zu lernen, nur nicht sicher, wo ich anfangen soll. Ich bin sicher, dass die Datei vorhanden ist, da ich die .shp-Datei in QGIS geladen habe und sie als Bild angezeigt wird. Ich sollte auch Lese- / Schreibzugriff haben, da ich ein Administrator auf dem Computer bin und dieser sich nur in meiner Dropbox befindet.
André.B
Versuchen Sie, es aus Dropbox auf ein lokales Laufwerk zu verschieben. Möglicherweise liegt ein Problem mit dem Pfad vor, das dazu führt, dass QGIS ihn ablehnt. Was Sie tun möchten, sollte in QGIS ziemlich einfach sein. Welche Version von QGIS verwenden Sie?
Keagan Allan
1
Ok, versuchen Sie, die Polylinie in ein Raster zu konvertieren. Das Proximity-Tool in QGIS benötigt eine Rastereingabe. Spielen Sie mit den Einstellungen gemäß der Hilfe des Tools: docs.qgis.org/2.8/en/docs/user_manual/processing_algs/gdalogr/… . Beachten Sie, es ist immer noch ein intensiver Prozess, ich teste es jetzt zum Spaß und es läuft seit 30 Minuten und läuft immer noch ...
Keagan Allan
1
Welche Größe des Ausgabe-Rasters in Bezug auf Zeilen und Spalten möchten Sie erstellen? Werden Sie tatsächlich in der Lage sein, mit diesem Raster zu arbeiten, sobald Sie es erstellt haben? Wenn die Dateigröße des Ganzen ein Problem darstellt, können Sie kleinere Kacheln erstellen. Dies können Sie auch aus Gründen der Geschwindigkeit parallel in einem Cluster oder einer Cloud tun.
Spacedman

Antworten:

9

Mit PyQGIS und GDAL ist die Python-Bibliothek nicht sehr schwierig. Sie benötigen Geotransformationsparameter (x oben links x, x Pixelauflösung, Drehung, oben links y, Drehung, ns Pixelauflösung) und Zeilen- und Spaltennummer, um das resultierende Raster zu erstellen. Zur Berechnung der Entfernung zur nächsten Küstenlinie ist eine Vektorebene zur Darstellung der Küstenlinie erforderlich.

Mit PyQGIS wird jeder Rasterpunkt als Mittelpunkt der Zelle berechnet und sein Abstand zur Küstenlinie mithilfe der Methode 'nextSegmentWithContext' aus der QgsGeometry- Klasse gemessen . Die GDAL- Python-Bibliothek wird zum Erstellen eines Rasters mit diesen Abstandswerten im Array Zeilen x Spalten verwendet.

Der folgende Code wurde zum Erstellen eines Distanzrasters (25 m × 25 m Auflösung und 1000 Zeilen x 1000 Spalten) verwendet, beginnend mit Punkt (397106.7689872353, 4499634.06675821). in der Nähe der Westküste der USA.

from osgeo import gdal, osr
import numpy as np
from math import sqrt

registry = QgsProject.instance()

line = registry.mapLayersByName('shoreline_10N')

crs = line[0].crs()

wkt = crs.toWkt()

feats_line = [ feat for feat in line[0].getFeatures()]

pt = QgsPoint(397106.7689872353, 4499634.06675821)

xSize = 25
ySize = 25

rows = 1000
cols = 1000

raster = [ [] for i in range(cols) ]

x =   xSize/2
y = - ySize/2

for i in range(rows):
    for j in range(cols):
        point = QgsPointXY(pt.x() + x, pt.y() + y)
        tupla = feats_line[0].geometry().closestSegmentWithContext(point)
        raster[i].append(sqrt(tupla[0]))

        x += xSize
    x =  xSize/2
    y -= ySize

data = np.array(raster)

# Create gtif file 
driver = gdal.GetDriverByName("GTiff")

output_file = "/home/zeito/pyqgis_data/distance_raster.tif"

dst_ds = driver.Create(output_file, 
                       cols, 
                       rows, 
                       1, 
                       gdal.GDT_Float32)

#writting output raster
dst_ds.GetRasterBand(1).WriteArray( data )

transform = (pt.x(), xSize, 0, pt.y(), 0, -ySize)

#setting extension of output raster
# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
dst_ds.SetGeoTransform(transform)

# setting spatial reference of output raster 
srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)
dst_ds.SetProjection( srs.ExportToWkt() )

dst_ds = None

Nach dem Ausführen des obigen Codes wurde das resultierende Raster in QGIS geladen und es sieht wie im folgenden Bild aus (Pseudofarbe mit 5 Klassen und Spektralrampe). Die Projektion ist UTM 10 N (EPSG: 32610)

Geben Sie hier die Bildbeschreibung ein

xunilk
quelle
Dies ist vielleicht kein Problem, aber eine Sache, über die ich mir ein wenig Sorgen mache, ist, dass das Polygon aus Neuseeland und seinen umliegenden Inseln besteht, was bedeutet, dass es einen großen Teil des umgebenden Meeres umfasst. Ich versuche, mich mit dem Code vertraut zu machen, aber könnte ich mit Ihrem Beispiel den Wert für alle Zellen auf See auf NA setzen? Ich interessiere mich wirklich nur für die Entfernung zum Meer von Punkten an Land.
André.B
Entschuldigung im Voraus, wenn dies eine dumme Frage ist, aber wie wähle ich einen neuen Ausgangspunkt in Neuseeland so, wie Sie die Koordinaten für die Staaten festlegen? Wie behalte ich es auch in EPSG: 2193?
André.B
7

Kann eine Lösung sein, um zu versuchen:

  1. Generieren Sie ein Raster (Typ "Punkt", Algorithmus "Raster erstellen")
  2. Berechnen Sie den nächsten Abstand zwischen Ihren Punkten (Gitter) und Ihrer Linie (Küste) mit dem Algorithmus "Attribut nach nächstem verbinden". Achten Sie darauf, nur maximal 1 nächste Nachbarn auszuwählen.

Jetzt sollten Sie wie in diesem Beispiel eine neue Punktebene mit dem Abstand zur Küste haben Geben Sie hier die Bildbeschreibung ein

  1. Bei Bedarf können Sie Ihre neue Punktebene in ein Raster konvertieren (Algorithmus "rasterize")

Geben Sie hier die Bildbeschreibung ein

Christophe
quelle
2

In QGIS können Sie das GRASS-Plugin ausprobieren. Soweit ich weiß, verwaltet es den Speicher besser als R, und ich erwarte, dass die andere Lösung auf großen Flächen fehlschlägt.

Der Befehl GRASS heißt r.grow.distance und befindet sich in der Verarbeitungssymbolleiste. Beachten Sie, dass Sie Ihre Zeile zunächst in Raster konvertieren müssen.

Geben Sie hier die Bildbeschreibung ein

Eines Ihrer Probleme könnte die Größe der Ausgabe sein, sodass Sie einige nützliche Erstellungsoptionen hinzufügen können, z. B. (für eine TIF-Datei) BIGTIFF = YES, TILED = YES, COMPRESS = LZW, PREDICTOR = 3

Radouxju
quelle
Gibt es eine Möglichkeit, ein Meeresgebiet zu eliminieren, um die Größe / Rechenzeit zu reduzieren?
André.B
Wenn Sie theoretisch den Abstand vom See (alle sehen Pixel mit demselben Wert, dh einem Polygon) anstelle von der Küstenlinie verwenden, sollten Sie Rechenzeit sparen. Die unkomprimierte Größe des Rasters sollte gleich sein, aber die Komprimierung ist effizienter, sodass Sie auch die endgültige Größe reduzieren sollten.
Radouxju
0

Ich würde es anders herum versuchen. Wenn Sie ein Poligon von NZ verwenden, konvertieren Sie Polygonkanten in Linien. Erstellen Sie danach alle 25 Meter Abstand von der Grenze einen Puffer an der Grenze (möglicherweise hilft Centorid bei der Bestimmung, wann angehalten werden soll). Schneiden Sie dann Puffer mit Polygon aus und konvertieren Sie diese Polygone in Raster. Ich bin mir nicht sicher, ob das funktionieren würde, aber auf jeden Fall werden Sie weniger RAM benötigen. Und PostGiS ist großartig, wenn Sie Leistungsprobleme haben.

Hoffe es könnte wenigstens ein bisschen helfen :)

Kumbra
quelle
0

Ich wollte ursprünglich nicht meine eigene Frage beantworten, aber ein Kollege von mir (der diese Site nicht nutzt) hat mir eine Menge Python-Code geschrieben, um das zu tun, wonach ich suche. einschließlich der Begrenzung der Zellen auf den Abstand zur Küste nur für terrestrische Zellen und des Zurücklassens der Zellen auf Meeresbasis als NAs. Der folgende Code sollte von jeder Python-Konsole aus ausgeführt werden können. Die einzigen Dinge, die geändert werden müssen, sind:

1) Legen Sie die Skriptdatei in denselben Ordner wie die gewünschte Formdatei.

2) Ändern Sie den Namen des Shapefiles im Python-Skript in den Namen Ihres Shapefiles.

3) stellen Sie die gewünschte Auflösung ein und;

4) Ändern Sie das Ausmaß, um es an andere Raster anzupassen.

Größere Shapefiles als das, was ich verwende, erfordern viel RAM, aber ansonsten lässt sich das Skript schnell ausführen (etwa drei Minuten, um ein Raster mit einer Auflösung von 50 m und zehn Minuten für ein Raster mit einer Auflösung von 25 m zu erstellen).

#------------------------------------------------------------------------------

from osgeo import gdal, ogr
import numpy as np
from scipy import ndimage
import matplotlib.pyplot as plt
import time

startTime = time.perf_counter()

#------------------------------------------------------------------------------

# Define spatial footprint for new raster
cellSize = 50 # ANDRE CHANGE THIS!!
noData = -9999
xMin, xMax, yMin, yMax = [1089000, 2092000, 4747000, 6224000]
nCol = int((xMax - xMin) / cellSize)
nRow = int((yMax - yMin) / cellSize)
gdal.AllRegister()
rasterDriver = gdal.GetDriverByName('GTiff')
NZTM = 'PROJCS["NZGD2000 / New Zealand Transverse Mercator 2000",GEOGCS["NZGD2000",DATUM["New_Zealand_Geodetic_Datum_2000",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6167"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4167"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",173],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",1600000],PARAMETER["false_northing",10000000],AUTHORITY["EPSG","2193"],AXIS["Easting",EAST],AXIS["Northing",NORTH]]'

#------------------------------------------------------------------------------ 

inFile = "new_zealand.shp" # CHANGE THIS!!

# Import vector file and extract information
vectorData = ogr.Open(inFile)
vectorLayer = vectorData.GetLayer()
vectorSRS = vectorLayer.GetSpatialRef()
x_min, x_max, y_min, y_max = vectorLayer.GetExtent()

# Create raster file and write information
rasterFile = 'nz.tif'
rasterData = rasterDriver.Create(rasterFile, nCol, nRow, 1, gdal.GDT_Int32, options=['COMPRESS=LZW'])
rasterData.SetGeoTransform((xMin, cellSize, 0, yMax, 0, -cellSize))
rasterData.SetProjection(vectorSRS.ExportToWkt())
band = rasterData.GetRasterBand(1)
band.WriteArray(np.zeros((nRow, nCol)))
band.SetNoDataValue(noData)
gdal.RasterizeLayer(rasterData, [1], vectorLayer, burn_values=[1])
array = band.ReadAsArray()
del(rasterData)

#------------------------------------------------------------------------------

distance = ndimage.distance_transform_edt(array)
distance = distance * cellSize
np.place(distance, array==0, noData)

# Create raster file and write information
rasterFile = 'nz-coast-distance.tif'
rasterData = rasterDriver.Create(rasterFile, nCol, nRow, 1, gdal.GDT_Float32, options=['COMPRESS=LZW'])
rasterData.SetGeoTransform((xMin, cellSize, 0, yMax, 0, -cellSize))
rasterData.SetProjection(vectorSRS.ExportToWkt())
band = rasterData.GetRasterBand(1)
band.WriteArray(distance)
band.SetNoDataValue(noData)
del(rasterData)

#------------------------------------------------------------------------------

endTime = time.perf_counter()

processTime = endTime - startTime

print(processTime)
André.B
quelle