Créer un raster où chaque cellule enregistre la distance à la mer?

Je veux créer un raster avec une résolution de 25 mètres × 25 mètres, où chaque cellule contient la distance jusqu'au littoral le plus proche, calculée à partir du centre de la cellule. Pour ce faire, je n'ai qu'un fichier de formes des côtes de la Nouvelle-Zélande .

J'ai essayé de suivre le tutoriel de Dominic Roye pour le faire en R qui fonctionne ... en quelque sorte. C'est bien jusqu'à une résolution d'environ 1 km × 1 km, mais si j'essaie d'aller plus haut, la RAM qu'elle nécessite dépasse largement celle disponible sur mon PC (~ 70 Go de RAM requis) ou toute autre à laquelle j'ai également accès. En disant cela, je pense que c'est une limitation de R et je soupçonne que QGIS pourrait avoir un moyen plus efficace de calculer ce raster, mais je suis nouveau dans ce domaine et je n'arrive pas à comprendre comment le faire.

J'ai essayé de suivre Créer un raster avec une distance à utiliser avec QGIS? pour le créer dans QGIS mais il renvoie cette erreur:

_core.QgsProcessingException: impossible de charger la couche source pour INPUT: C: /..../ Coastline / nz-coastlines-and-islands-polygons-topo-150k.shp not found

et je ne sais pas pourquoi.

Quelqu'un a-t-il des suggestions sur ce qui pourrait mal tourner ou sur une autre façon de procéder?

Éditer:

Le raster que j'espère produire aurait environ 59684 lignes et 40827 colonnes de sorte qu'il chevauche le raster annuel de déficit hydrique de LINZ. Si le raster produit est plus grand que le raster annuel de déficit hydrique, je peux le couper en R ...

Je pense qu'une chose pourrait être un problème potentiel, c'est que le fichier de formes du littoral néo-zélandais a une grande quantité de mer entre les îles, et je ne suis pas intéressé par le calcul de la distance à la côte pour ces cellules. Je veux vraiment seulement calculer les valeurs des cellules qui incluent une tranche de terrain. Je ne sais pas comment faire cela, ni si c'est vraiment un problème.

Avec PyQGIS et GDAL, la bibliothèque python n'est pas très difficile à faire. Vous avez besoin de paramètres de géotransformation (en haut à gauche x, résolution en pixels x, rotation, en haut à gauche y, rotation, résolution en pixels ns) et d'un nombre de lignes et de colonnes pour créer le raster résultant. Pour calculer la distance au littoral le plus proche, il est nécessaire d'avoir une couche vectorielle pour représenter le littoral.

Avec PyQGIS , chaque point raster en tant que centre de la cellule est calculé et sa distance au littoral est mesurée à l'aide de la méthode «le plus procheSegmentWithContext» de la classe QgsGeometry . La bibliothèque python GDAL est utilisée pour produire un raster avec ces valeurs de distance dans un tableau de lignes x colonnes.

Le code suivant a été utilisé pour créer un raster de distance (résolution de 25 m × 25 m et 1000 lignes x 1000 colonnes) commençant au point (397106.7689872353, 4499634.06675821); près de la côte ouest des États-Unis.

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

Après avoir exécuté le code ci-dessus, le raster résultant a été chargé dans QGIS et il ressemble à l'image suivante (pseudocolor avec 5 classes et rampe spectrale). La projection est UTM 10 N (EPSG: 32610)

Peut être une solution à essayer:

1. Générer une grille (Type "point", algorithme "Créer une grille")
2. Calculez la distance la plus proche entre vos points (grille) et votre ligne (côte) avec l'algorithme "join attribute by most". Attention à ne choisir qu'un maximum de 1 voisin le plus proche.

Maintenant, vous devriez avoir une nouvelle couche de points avec la distance à la côte comme dans cet exemple

3. Si nécessaire, vous pouvez convertir votre nouvelle couche de points en un raster (algorithme "rasterize")

Dans QGIS, vous pouvez essayer le plugin GRASS. Pour autant que je sache, il gère mieux la mémoire que R, et je m'attends à ce que l'autre solution échoue sur de grandes surfaces.

la commande GRASS s'appelle r.grow.distance, que vous pouvez trouver dans la barre d'outils de traitement. Notez que vous devez d'abord convertir votre ligne en raster.

L'un de vos problèmes pourrait être la taille de la sortie, vous pouvez donc ajouter des options de création utiles telles que (pour un fichier tif) BIGTIFF = YES, TILED = YES, COMPRESS = LZW, PREDICTOR = 3

J'essaierais autrement. Si vous utilisez un poligon de NZ, convertissez les arêtes du polygone en ligne. Après cela, créez un tampon sur la frontière tous les 25 mètres de distance de la frontière (peut-être que le centoride pourrait aider à déterminer quand s'arrêter). Découpez ensuite les tampons avec un polygone, puis convertissez ces polygones en raster. Je ne suis pas sûr que cela fonctionnerait, mais vous aurez certainement besoin de moins de RAM. Et PostGiS est idéal lorsque vous rencontrez des problèmes de performances.

J'espère que cela pourrait aider au moins un peu :)

Au départ, je n'allais pas répondre à ma propre question, mais un de mes collègues (qui n'utilise pas ce site) m'a écrit un tas de code python pour faire ce que je recherche; notamment en limitant les cellules à la distance de la côte pour les cellules terrestres uniquement et en laissant les cellules basées sur la mer comme des AN. Le code suivant devrait pouvoir s'exécuter à partir de n'importe quelle console python, les seules choses nécessitant une modification étant:

1) Placez le fichier de script dans le même dossier que le fichier de forme qui vous intéresse;

2) changez le nom du fichier de formes dans le script python en quelque soit le nom de votre fichier de formes;

3) définissez la résolution souhaitée, et;

4) modifiez l'étendue pour qu'elle corresponde à d'autres rasters.

Des fichiers de formes plus grands que ceux que j'utilise nécessiteront de grandes quantités de RAM, mais sinon le script est rapide à exécuter (environ trois minutes pour produire un raster de résolution de 50 m et dix minutes pour un raster de résolution de 25 m).

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

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)