Создание растра, где каждая ячейка записывает расстояние до моря?

Я хочу создать растр с разрешением 25 х 25 м, в котором каждая ячейка содержит расстояние до ближайшей береговой линии, рассчитанное по центру ячейки. Для этого у меня есть только шейп-файл побережья Новой Зеландии .

Я пытался следовать учебному пособию Доминика Роя для того, чтобы сделать это в R, что работает ... вроде. Это нормально до разрешения около 1 км × 1 км, но если я пытаюсь подняться выше, то ОЗУ, которое требуется, значительно превышает доступное на моем ПК (требуется ~ 70 ГБ ОЗУ) или любой другой, к которому у меня есть доступ. Говоря об этом, я думаю, что это ограничение R, и я подозреваю, что QGIS может иметь более вычислительно эффективный способ создания этого растра, но я новичок в этом и не могу понять, как это сделать.

Я пытался следовать Создание растра с расстоянием до объекта, используя QGIS? создать его в QGIS, но он возвращает эту ошибку:

_core.QgsProcessingException: не удалось загрузить исходный слой для INPUT: C: /..../ Coastline / nz-coastlines-and-Islands-polygons-topo-150k.shp не найден

и я не уверен почему.

У кого-нибудь есть какие-либо предположения о том, что может пойти не так или альтернативный способ сделать это?

Редактировать:

Растр, который я надеюсь создать, будет иметь около 59684 строк и 40827 столбцов, чтобы он перекрывался с годовым растром дефицита воды в LINZ. Если создаваемый растр больше, чем годовой растр с дефицитом воды, я могу отсечь его в R, хотя ...

Одна вещь, которая, я думаю, может быть потенциальной проблемой, состоит в том, что шейп-файл береговой линии Новой Зеландии имеет большое количество моря между островами, и я не заинтересован в расчете расстояния до побережья для этих ячеек. Я действительно хочу только рассчитать значения для ячеек, которые включают какой-то кусок земли. Я не уверен, как это сделать, или если это действительно проблема, хотя.

С PyQGIS и библиотекой Python GDAL это не очень сложно сделать. Для создания результирующего растра нужны параметры геопреобразования (верхний левый x, разрешение в пикселях x, вращение, верхний левый y, вращение, разрешение в пикселях ns), а также количество строк и столбцов. Для расчета расстояния до ближайшей береговой линии необходим векторный слой для представления береговой линии.

С помощью PyQGIS вычисляется каждая растровая точка в качестве центра ячейки, а ее расстояние до береговой линии измеряется с помощью метода closestSegmentWithContext из класса QgsGeometry . Библиотека Python GDAL используется для создания растра с этими значениями расстояния в массиве строк х столбцов.

Следующий код был использован для создания растрового расстояния (разрешение 25 м х 25 м и 1000 строк х 1000 столбцов), начиная с точки (397106.7689872353, 4499634.06675821); недалеко от западного побережья США.

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

После выполнения кода выше полученный растр был загружен в QGIS и выглядит так, как показано на следующем рисунке (псевдоцвет с 5 классами и спектральной рампой). Проекция - UTM 10 N (EPSG: 32610)

Может быть решение попробовать:

1. Сгенерировать сетку (Тип «точка», алгоритм «Создать сетку»)
2. Рассчитайте ближайшее расстояние между вашими точками (сеткой) и вашей линией (побережьем) с помощью алгоритма «объединить атрибут по ближайшему». Будьте осторожны, выбирая не более 1 ближайшего соседа.

Теперь у вас должен быть новый точечный слой с расстоянием до побережья, как в этом примере

3. При необходимости вы можете преобразовать ваш новый точечный слой в растр (алгоритм «растеризация»)

В QGIS вы можете попробовать плагин GRASS. Насколько я знаю, он управляет памятью лучше, чем R, и я ожидаю, что другое решение даст сбой на больших площадях.

Команда GRASS называется r.grow.distance, которую вы можете найти на панели инструментов обработки. Обратите внимание, что вам нужно сначала преобразовать вашу линию в растр.

Одной из ваших проблем может быть размер выходных данных, поэтому вы можете добавить некоторые полезные параметры создания, такие как (для файла TIF) BIGTIFF = YES, TILED = YES, COMPRESS = LZW, PREDICTOR = 3

Я бы попробовал другой путь. Если вы используете полигон NZ, тогда преобразуйте ребра полигона в линию. После этого создайте буфер на границе на каждые 25 метров расстояния от границы (возможно, Centorid может помочь определить, когда остановиться). Затем вырежьте буферы с помощью многоугольника и затем конвертируйте эти многоугольники в растр. Я не уверен, что это сработает, но определенно вам понадобится меньше оперативной памяти. И PostGiS великолепен, когда у вас есть проблемы с производительностью.

Надеюсь, это поможет хоть немного :)

Изначально я не собирался отвечать на свой вопрос, но мой коллега (который не использует этот сайт) написал мне кучу кода на Python, чтобы делать то, что я делаю; включая ограничение ячеек расстоянием до побережья только для наземных ячеек и оставление морских ячеек в качестве NA. Следующий код должен быть в состоянии работать с любой консоли Python, единственное, что нужно изменить:

1) Поместите файл скрипта в ту же папку, что и интересующий вас файл формы;

2) измените имя шейп-файла в скрипте Python на любое имя вашего шейп-файла;

3) установить желаемое разрешение и;

4) измените экстент, чтобы соответствовать другим растрам.

Большие шейп-файлы, чем те, которые я использую, потребуют большого объема оперативной памяти, но в противном случае скрипт будет быстро запущен (около трех минут для создания растра с разрешением 50 м и десяти минут для растра с разрешением 25 м).

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

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)