Membuat raster di mana setiap sel mencatat jarak ke laut?

10

Saya ingin membuat raster dengan resolusi 25 meter × 25 meter, di mana setiap sel berisi jarak ke garis pantai terdekat, yang dihitung dari pusat sel. Untuk melakukan ini, yang saya miliki hanyalah membentuk garis pantai Selandia Baru .

Saya telah mencoba mengikuti tutorial Dominic Roye untuk melakukannya di R yang berfungsi ... agak. Baik-baik saja hingga resolusi 1 km × 1 km tetapi jika saya mencoba untuk meningkatkan RAM, itu membutuhkan melebihi yang tersedia pada PC saya (~ diperlukan 70 gb RAM) atau yang lain yang saya akses juga. Mengatakan itu, saya pikir ini adalah batasan R dan saya menduga bahwa QGIS mungkin memiliki cara yang lebih efisien secara komputasional untuk menciptakan raster ini, tetapi saya baru dalam hal ini dan saya tidak tahu bagaimana melakukannya.

Saya telah mencoba mengikuti Membuat raster dengan jarak ke fitur menggunakan QGIS? untuk membuatnya di QGIS tetapi mengembalikan kesalahan ini:

_core.QgsProcessingException: Tidak dapat memuat lapisan sumber untuk INPUT: C: / ..../ Garis pantai / nz-garis pantai-dan-pulau-poligon-topo-150k.shp tidak ditemukan

dan saya tidak yakin mengapa.

Adakah yang punya saran tentang apa yang mungkin salah atau cara alternatif untuk melakukan ini?

Edit:

Raster yang saya harapkan akan diproduksi memiliki sekitar 5.9684 baris dan 40.827 kolom sehingga tumpang tindih dengan raster defisit air tahunan dari LINZ. Jika raster yang diproduksi lebih besar dari raster defisit air tahunan saya dapat memotongnya di R meskipun ...

Satu hal yang saya pikir mungkin menjadi masalah potensial adalah bahwa pembentukan garis pantai NZ memiliki sejumlah besar laut di antara pulau-pulau, dan saya tidak tertarik menghitung jarak ke pantai untuk sel-sel ini. Saya benar-benar hanya ingin menghitung nilai untuk sel yang menyertakan beberapa bidang tanah. Saya tidak yakin bagaimana melakukan ini, atau apakah ini benar-benar masalah.

André.B
sumber
1
Apakah Anda menjalankan skrip untuk melakukan ini? Atau apakah Anda menggunakan alat di QGIS? Sesuatu untuk diperiksa, meskipun sepertinya memang seharusnya - periksa file tersebut benar-benar ada di tempat yang Anda katakan ... juga periksa apakah Anda telah membaca dan menulis akses ke folder tertentu.
Keagan Allan
Saat ini menggunakan alat tetapi saya cukup tertarik untuk mempelajari skrip, hanya tidak yakin harus mulai dari mana. Saya yakin file itu ada, karena saya telah memuat file .shp ke dalam QGIS dan muncul sebagai gambar. Saya seharusnya memiliki akses baca / tulis juga karena saya adalah admin di mesin dan itu hanya di dropbox saya.
André.B
Coba pindahkan dari Dropbox ke drive lokal. Mungkin ada masalah dengan jalur yang menyebabkan QGIS menolaknya. Apa yang ingin Anda lakukan harus cukup sederhana di QGIS. Versi QGIS mana yang Anda gunakan?
Keagan Allan
1
Oke, coba ubah polyline menjadi raster. Alat Kedekatan di QGIS membutuhkan input raster. Bermain-main dengan pengaturan sesuai bantuan alat: docs.qgis.org/2.8/en/docs/user_manual/processing_algs/gdalogr/… . Perhatikan, ini masih merupakan proses intensif, saya mengujinya untuk bersenang-senang sekarang dan telah berjalan selama 30 menit dan masih berjalan ...
Keagan Allan
1
Berapa ukuran output raster dalam hal baris dan kolom yang Anda coba buat? Apakah Anda benar-benar akan dapat bekerja dengan raster itu setelah Anda membuatnya? Jika ukuran file dari semuanya adalah masalah, dapatkah Anda membuat ubin yang lebih kecil, yang juga merupakan sesuatu yang dapat Anda lakukan secara paralel pada sebuah cluster atau cloud untuk kecepatan.
Spacedman

Jawaban:

9

Dengan PyQGIS dan GDAL python library tidak terlalu sulit untuk melakukannya. Anda memerlukan parameter transformasi geografis (x kiri atas, resolusi x pixel, rotasi, y kiri atas, rotasi, resolusi ns pixel) dan jumlah baris dan kolom untuk membuat raster yang dihasilkan. Untuk menghitung jarak ke garis pantai terdekat, perlu lapisan vektor untuk mewakili garis pantai.

Dengan PyQGIS , setiap titik raster sebagai pusat sel dihitung dan jaraknya ke garis pantai diukur dengan menggunakan metode 'proximitySegmentWithContext' dari kelas QgsGeometry . Pustaka python GDAL digunakan untuk menghasilkan raster dengan nilai-nilai jarak ini dalam baris x kolom array

Kode berikut digunakan untuk membuat raster jarak (resolusi 25 m × 25 m dan 1000 baris x 1000 kolom) mulai dari titik (397106.7689872353, 4499634.06675821); dekat dengan garis pantai barat 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

Setelah menjalankan kode di atas, raster yang dihasilkan dimuat dalam QGIS dan terlihat seperti pada gambar berikut (pseudocolor dengan 5 kelas dan jalan Spectral). Proyeksi adalah UTM 10 N (EPSG: 32610)

masukkan deskripsi gambar di sini

xunilk
sumber
Ini mungkin bukan masalah tetapi satu hal yang saya agak khawatirkan adalah poligon adalah dari Selandia Baru dan pulau-pulau sekitarnya yang berarti itu mencakup sejumlah besar laut di sekitarnya. Saya mencoba menyiasati kode, tetapi dengan contoh Anda, apakah saya dapat menetapkan nilai untuk semua sel di laut menjadi NA? Saya benar-benar hanya tertarik pada jarak ke laut dari titik di darat.
André.B
Permintaan maaf sebelumnya jika ini adalah pertanyaan bodoh tapi bagaimana cara memilih titik awal baru di Selandia Baru dengan cara Anda mengatur koordinat untuk negara bagian? Juga bagaimana cara menyimpannya di EPSG: 2193?
André.B
7

Mungkin menjadi solusi untuk dicoba:

  1. Hasilkan kisi (Ketik "titik", algoritma "Buat kisi")
  2. Hitung jarak terdekat antara titik Anda (kotak) dan garis Anda (pantai) dengan algoritma "gabungkan atribut dengan terdekat". Berhati-hatilah untuk memilih hanya maksimal 1 tetangga terdekat.

Sekarang Anda harus memiliki layer titik baru dengan jarak ke pantai seperti dalam contoh ini masukkan deskripsi gambar di sini

  1. Jika perlu, Anda dapat mengubah layer titik baru Anda menjadi raster (algoritma "rasterize")

masukkan deskripsi gambar di sini

Christophe
sumber
2

Di dalam QGIS Anda bisa mencoba plugin GRASS. Sejauh yang saya tahu itu mengelola memori lebih baik daripada R, dan saya berharap solusi lain gagal pada area yang luas.

perintah GRASS disebut r.grow.distance, yang dapat Anda temukan di toolbar pemrosesan. Perhatikan bahwa Anda perlu mengonversi baris Anda menjadi raster pada awalnya.

masukkan deskripsi gambar di sini

Salah satu masalah Anda mungkin ukuran output, sehingga Anda dapat menambahkan beberapa opsi pembuatan berguna seperti (untuk file tif) BIGTIFF = YA, TILED = YA, KOMPRES = LZW, PREDICTOR = 3

radouxju
sumber
Apakah ada cara agar saya bisa menghilangkan area laut untuk mengurangi ukuran / waktu komputasi?
André.B
secara teori, jika Anda menggunakan jarak dari tahta (semua melihat piksel dengan nilai yang sama, itu adalah poligon) alih-alih dari garis pantai Anda harus menghemat waktu perhitungan. Ukuran raster yang tidak terkompresi harus sama, tetapi kompresi akan lebih efisien, jadi Anda juga harus mengurangi ukuran akhir.
radouxju
0

Saya akan mencoba sebaliknya. Jika Anda menggunakan poligon NZ maka konversikan tepi poligon menjadi garis. Setelah itu buat buffer pada batas untuk setiap 25 meter jarak dari batas (mungkin centorid dapat membantu menentukan kapan harus berhenti). Kemudian potong buffer dengan poligon dan kemudian ubah poligon itu menjadi raster. Saya tidak yakin ini akan berhasil tetapi pasti Anda akan membutuhkan lebih sedikit RAM. Dan PostGiS hebat ketika Anda memiliki masalah kinerja.

Semoga ini bisa membantu setidaknya sedikit :)

Kumbra
sumber
0

Saya awalnya tidak akan menjawab pertanyaan saya sendiri tetapi seorang kolega saya (yang tidak menggunakan situs ini), menulis kepada saya banyak kode python untuk melakukan apa yang saya inginkan; termasuk membatasi sel untuk memiliki jarak ke pantai hanya sel terestrial dan meninggalkan sel berbasis laut sebagai NAs. Kode berikut harus dapat dijalankan dari konsol python, dengan satu-satunya hal yang perlu diubah adalah:

1) Letakkan file skrip di folder yang sama dengan file bentuk yang diinginkan;

2) ubah nama shapefile dalam skrip python menjadi apa pun nama shapefile Anda;

3) mengatur resolusi yang diinginkan, dan;

4) ubah sejauh sesuai dengan raster lainnya.

Shapefile yang lebih besar daripada yang saya gunakan akan membutuhkan RAM dalam jumlah besar tetapi jika tidak skrip cepat dijalankan (sekitar tiga menit untuk menghasilkan raster resolusi 50m dan sepuluh menit untuk raster resolusi 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)
André.B
sumber