Tworzysz raster, w którym każda komórka rejestruje odległość do morza?

Chcę utworzyć raster o rozdzielczości 25 metrów × 25 metrów, w którym każda komórka zawiera odległość do najbliższej linii brzegowej, obliczoną od środka komórki. Aby to zrobić, wszystko, co mam, to kształtna linia brzegowa Nowej Zelandii .

Próbowałem wykonać samouczek Dominica Roye'a dotyczący robienia tego w języku R, który działa ... w pewnym sensie. Jest w porządku do rozdzielczości około 1 km × 1 km, ale jeśli spróbuję przejść na wyższą pamięć RAM, wymaga to znacznie więcej niż na moim komputerze (wymagane ~ 70 GB pamięci RAM) lub dowolnego innego, do którego mam dostęp. Mówiąc to, myślę, że jest to ograniczenie R i podejrzewam, że QGIS może mieć bardziej wydajny pod względem obliczeniowym sposób tworzenia tego rastra, ale jestem nowy i nie jestem w stanie dowiedzieć się, jak to zrobić.

Próbowałem już tworzyć Tworzenie rastra z odległością do funkcji za pomocą QGIS? aby utworzyć go w QGIS, ale zwraca ten błąd:

_core.QgsProcessingException: Nie można załadować warstwy źródłowej dla INPUT: C: /..../ Coastline / nz-coastlines-and-islands-polygons-topo-150k.shp nie znaleziono

i nie jestem pewien dlaczego.

Czy ktoś ma jakieś sugestie, co może pójść źle lub alternatywny sposób na zrobienie tego?

Edytować:

Raster, który mam nadzieję wyprodukować, miałby około 59684 wierszy i 40827 kolumn, tak aby pokrywał się z rocznym rastrem deficytu wody z LINZ. Jeśli produkowany raster jest większy niż raster rocznego deficytu wody, mogę go wyciąć w R ...

Jedną z rzeczy, które moim zdaniem mogą być potencjalnym problemem, jest to, że plik kształtu linii brzegowej NZ ma dużą ilość morza między wyspami i nie jestem zainteresowany obliczeniem odległości do wybrzeża dla tych komórek. Naprawdę chcę tylko obliczyć wartości dla komórek, które zawierają kawałek ziemi. Nie jestem pewien, jak to zrobić, czy też jest to rzeczywiście problem.

Z PyQGIS i biblioteką Pythona GDAL nie jest to bardzo trudne. Potrzebujesz parametrów transformacji geograficznej (lewy górny x, rozdzielczość x pikseli, obrót, lewy górny y, obrót, rozdzielczość ns pikseli) oraz liczbę wierszy i kolumn, aby utworzyć wynikowy raster. Do obliczenia odległości do najbliższej linii brzegowej potrzebna jest warstwa wektorowa do reprezentowania linii brzegowej.

W PyQGIS obliczany jest każdy punkt rastrowy jako środek komórki, a jego odległość do linii brzegowej jest mierzona przy użyciu metody „closestSegmentWithContext” z klasy QgsGeometry . Biblioteka Pythona GDAL służy do tworzenia rastra z tymi wartościami odległości w tablicy wierszy x kolumn.

Poniższy kod został użyty do utworzenia rastra odległości (25 m × 25 m rozdzielczości i 1000 wierszy x 1000 kolumn), zaczynając od punktu (397106.7689872353, 4499634.06675821); blisko zachodniego wybrzeża 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

Po uruchomieniu powyżej kodu powstały raster został załadowany do QGIS i wygląda jak na poniższym obrazie (pseudokolor z 5 klasami i rampą spektralną). Projekcja to UTM 10 N (EPSG: 32610)

Może być rozwiązaniem, aby spróbować:

1. Wygeneruj siatkę (wpisz „punkt”, algorytm „Utwórz siatkę”)
2. Oblicz najbliższą odległość między punktami (siatką) a linią (wybrzeżem) za pomocą algorytmu „połącz atrybut przez najbliższego”. Uważaj, aby wybrać maksymalnie 1 najbliższego sąsiada.

Teraz powinieneś mieć nową warstwę punktową z odległością do wybrzeża, jak w tym przykładzie

3. W razie potrzeby możesz przekonwertować nową warstwę punktową na raster (algorytm „rasterize”)

W QGIS możesz wypróbować wtyczkę GRASS. O ile wiem, lepiej zarządza pamięcią niż R i spodziewam się, że inne rozwiązanie zawiedzie na dużych obszarach.

polecenie GRASS nazywa się r.grow.distance, które można znaleźć na pasku narzędzi przetwarzania. Zauważ, że najpierw musisz przekonwertować linię na raster.

Jednym z problemów może być rozmiar danych wyjściowych, więc możesz dodać kilka przydatnych opcji tworzenia, takich jak (dla pliku tif) BIGTIFF = TAK, TILED = TAK, KOMPRESS = LZW, PREDICTOR = 3

Spróbowałbym na odwrót. Jeśli używasz wielokąta NZ, przekonwertuj krawędzie wielokąta na linię. Następnie utwórz bufor na granicy dla każdych 25 metrów odległości od granicy (być może centorid może pomóc w określeniu, kiedy zatrzymać). Następnie wytnij bufory za pomocą wielokąta, a następnie przekonwertuj te wielokąty na raster. Nie jestem pewien, czy to zadziała, ale na pewno będziesz potrzebować mniej pamięci RAM. A PostGiS jest świetny, gdy masz problemy z wydajnością.

Mam nadzieję, że to może trochę pomóc :)

Początkowo nie zamierzałem odpowiadać na moje własne pytanie, ale mój kolega (który nie korzysta z tej strony) napisał mi sporo kodu Pythona, aby zrobić to, o co mi chodzi; w tym ograniczenie komórek do odległości od wybrzeża tylko dla komórek lądowych i pozostawienie komórek morskich jako NA. Poniższy kod powinien być uruchomiony z dowolnej konsoli Pythona, a jedyne rzeczy, które wymagają modyfikacji, to:

1) Umieść plik skryptu w tym samym folderze, co plik kształtu;

2) zmień nazwę pliku shapefile w skrypcie python na dowolną nazwę twojego pliku shapefile;

3) ustaw żądaną rozdzielczość, oraz;

4) zmień zakres, aby dopasować do innych rastrów.

Większe pliki shapefile niż te, których używam, będą wymagały dużej ilości pamięci RAM, ale w przeciwnym razie skrypt jest szybki do uruchomienia (około trzech minut, aby uzyskać raster o rozdzielczości 50m i dziesięć minut dla rastra o rozdzielczości 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)