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

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.

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)

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

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

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.

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

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 :)

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)