¿Crear ráster donde cada celda registra la distancia al mar?

Quiero crear un ráster con una resolución de 25 metros × 25 metros, donde cada celda contiene la distancia a la costa más cercana, calculada desde el centro de la celda. Para hacer esto, todo lo que tengo es un archivo de forma de las costas de Nueva Zelanda .

He intentado seguir el tutorial de Dominic Roye para hacerlo en R, que funciona ... más o menos. Está bien hasta una resolución de aproximadamente 1 km × 1 km, pero si trato de subir más la RAM, requiere mucho más que la disponible en mi PC (se requieren ~ 70 gb de RAM) o cualquier otra a la que tenga acceso. Al decir eso, creo que esto es una limitación de R y sospecho que QGIS podría tener una forma más eficiente desde el punto de vista computacional de crear este ráster, pero soy nuevo en él y no puedo entender cómo hacerlo.

He intentado seguir ¿ Crear ráster con distancia a la característica usando QGIS? para crearlo en QGIS pero devuelve este error:

_core.QgsProcessingException: no se pudo cargar la capa de origen para INPUT: C: /..../ Coastline / nz-coastlines-and-islands-polygons-topo-150k.shp no encontrado

Y no estoy seguro de por qué.

¿Alguien tiene alguna sugerencia de lo que podría estar yendo mal o una forma alternativa de hacerlo?

Editar:

El ráster que espero producir tendría aproximadamente 59684 filas y 40827 columnas para que se superponga con el ráster anual de déficit de agua de LINZ. Si el ráster que se produce es más grande que el ráster anual de déficit de agua, puedo cortarlo en R ...

Una cosa que creo que podría ser un problema potencial es que el archivo de forma de la costa de Nueva Zelanda tiene una gran cantidad de mar entre las islas, y no estoy interesado en calcular la distancia a la costa de estas celdas. Realmente solo quiero calcular los valores para las celdas que incluyen algún segmento de tierra. No estoy seguro de cómo hacer esto, o si es realmente un problema.

Con PyQGIS y GDAL, la biblioteca de Python no es muy difícil de hacer. Necesita parámetros de transformación geográfica (arriba a la izquierda x, resolución de píxeles x, rotación, arriba a la izquierda y, rotación, resolución de píxeles ns) y número de filas y columnas para crear el ráster resultante. Para calcular la distancia al litoral más cercano, es necesaria una capa vectorial para representar el litoral.

Con PyQGIS , cada punto ráster como centro de la celda se calcula y su distancia al litoral se mide utilizando el método 'más cercanoSegmentoConContexto' de la clase QgsGeometry . La biblioteca Python GDAL se usa para producir un ráster con estos valores de distancia en filas x columnas.

El siguiente código se usó para crear un ráster de distancia (resolución de 25 m × 25 m y 1000 filas x 1000 columnas) comenzando en el punto (397106.7689872353, 4499634.06675821); cerca de la costa oeste de EE. UU.

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

Después de ejecutar el código anterior, el ráster resultante se cargó en QGIS y se ve como en la siguiente imagen (pseudocolor con 5 clases y rampa espectral). La proyección es UTM 10 N (EPSG: 32610)

Puede ser una solución para probar:

1. Generar una cuadrícula (Escriba "punto", algoritmo "Crear cuadrícula")
2. Calcule la distancia más cercana entre sus puntos (cuadrícula) y su línea (costa) con el algoritmo "unir atributo por más cercano". Tenga cuidado de elegir solo un máximo de 1 vecino más cercano.

Ahora debería tener una nueva capa de puntos con la distancia a la costa como en este ejemplo

3. Si es necesario, puede convertir su nueva capa de puntos en un ráster (algoritmo "rasterizar")

Dentro de QGIS puedes probar el complemento GRASS. Hasta donde sé, administra la memoria mejor que R, y espero que la otra solución falle en áreas grandes.

el comando GRASS se llama r.grow.distance, que puede encontrar en la barra de herramientas de procesamiento. Tenga en cuenta que primero debe convertir su línea en ráster.

Uno de sus problemas podría ser el tamaño de la salida, por lo que puede agregar algunas opciones de creación útiles, como (para un archivo tif) BIGTIFF = YES, TILED = YES, COMPRESS = LZW, PREDICTOR = 3

Lo intentaría al revés. Si está utilizando el polígono de NZ, convierta los bordes del polígono en una línea. Después de eso, cree un búfer en el límite por cada 25 metros de distancia del límite (tal vez centorid podría ayudar a determinar cuándo detenerse). Luego corte los búferes con un polígono y luego convierta esos polígonos en ráster. No estoy seguro de que esto funcione, pero definitivamente necesitarás menos RAM. Y PostGiS es excelente cuando tienes problemas de rendimiento.

Espero que pueda ayudar al menos un poco :)

Originalmente no iba a responder mi propia pregunta, pero un colega mío (que no usa este sitio) me escribió un montón de código de Python para hacer lo que busco; incluyendo limitar las células para que tengan la distancia a la costa solo para las células terrestres y dejar las células basadas en el mar como NA. El siguiente código debería poder ejecutarse desde cualquier consola de Python, siendo las únicas cosas que necesitan modificación:

1) Coloque el archivo de script en la misma carpeta que el archivo de forma de interés;

2) cambie el nombre del shapefile en el script de python a cualquiera que sea el nombre de su shapefile;

3) establecer la resolución deseada, y;

4) cambie la extensión para que coincida con otros rásteres.

Los archivos de forma más grandes que los que estoy usando requerirán grandes cantidades de RAM, pero de lo contrario el script se ejecuta rápidamente (aproximadamente tres minutos para producir un ráster de resolución de 50m y diez minutos para un ráster de resolución de 25m).

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

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)