25m x 25m 해상도의 래스터를 만들고 싶습니다. 여기서 각 셀에는 셀 중심에서 계산 된 가장 가까운 해안선까지의 거리가 포함됩니다. 이렇게하려면 뉴질랜드 해안선의 모양 파일 만 있으면 됩니다.
나는 R에서 작동하는 Dominic Roye의 튜토리얼 을 따라 노력했다. 약 1km × 1km 해상도로 괜찮지 만 RAM을 높이려고하면 내 PC에서 사용할 수있는 RAM (~ 70GB의 RAM 필요) 또는 내가 액세스 할 수있는 다른 것보다 훨씬 많이 필요합니다. 즉, 이것이 R의 한계라고 생각하며 QGIS 가이 래스터를 생성하는 데 더 계산적으로 효율적인 방법이 있을지 모른다고 생각하지만 새로운 것은 아니며 그것을 수행하는 방법을 알 수는 없습니다.
QGIS를 사용하여 피처까지 거리가있는 래스터 만들기를 시도 했습니까? QGIS에서 생성하지만 다음 오류를 반환합니다.
_core.QgsProcessingException : 입력에 대한 소스 레이어를로드 할 수 없습니다 : C : /…./ Coastline / nz-coastlines-and-islands-polygons-topo-150k.shp를 찾을 수 없습니다
왜 그런지 잘 모르겠습니다.
누구든지 잘못되고 있거나 다른 방법으로 제안 할 수 있습니까?
편집하다:
내가 생산하고자하는 래스터는 약 59684 개의 행과 40827 개의 열을 가지므로 LINZ 의 연간 물 부족 래스터 와 겹칩니다 . 생산되는 래스터가 연간 물 부족 래스터보다 큰 경우 R로 스니핑 할 수 있습니다.
잠재적 인 문제가 될 수 있다고 생각되는 한 가지는 뉴질랜드 해안선의 형태 파일이 섬 사이에 많은 양의 바다를 가지고 있다는 것입니다.이 셀의 해안 거리를 계산하는 데 관심이 없습니다. 토지 일부가 포함 된 셀의 값만 계산하려고합니다. 이 작업을 수행하는 방법을 모르거나 실제로 문제가 있는지 확실하지 않습니다.
답변
PyQGIS 및 GDAL python 라이브러리를 사용하면 그렇게 하기가 어렵지 않습니다. 결과 래스터를 생성하려면 지리적 변형 매개 변수 (왼쪽 상단 x, x 픽셀 해상도, 회전, 왼쪽 상단 y, 회전, ns 픽셀 해상도)와 행 및 열 번호가 필요합니다. 가장 가까운 해안선까지의 거리를 계산하려면 해안선을 나타내는 벡터 레이어가 필요합니다.
PyQGIS를 사용하면 셀 중심으로 각 래스터 포인트가 계산되고 해안선까지의 거리는 QgsGeometry 클래스의 ‘closestSegmentWithContext’메서드를 사용하여 측정됩니다 . GDAL 파이썬 라이브러리는 행 x 열 배열로 이러한 거리 값을 가진 래스터를 생성하는 데 사용됩니다.
다음 코드는 점 (397106.7689872353, 4499634.06675821)에서 시작하여 거리 래스터 (25m x 25m 해상도 및 1000 행 x 1000 열)를 만드는 데 사용되었습니다. 미국의 서쪽 해안선 근처.
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 개만 선택하도록주의하십시오.
이제이 예와 같이 해안까지의 거리를 가진 새로운 점 레이어가 있어야합니다
- 필요한 경우 새 점 레이어를 래스터로 변환 할 수 있습니다 (알고리즘 “래스터 화”).
답변
QGIS 내에서 GRASS 플러그인을 사용해 볼 수 있습니다. 내가 아는 한 R보다 메모리를 더 잘 관리하며 다른 솔루션은 넓은 영역에서 실패 할 것으로 예상됩니다.
GRASS 명령은 r.grow.distance라고하며 처리 도구 모음에서 찾을 수 있습니다. 처음에는 라인을 래스터로 변환해야합니다.
문제 중 하나는 출력 크기 일 수 있으므로 (tif 파일의 경우) BIGTIFF = YES, TILED = YES, COMPRESS = LZW, PREDICTOR = 3과 같은 유용한 생성 옵션을 추가 할 수 있습니다.
답변
나는 다른 방법으로 시도 할 것입니다. NZ의 폴리곤을 사용하는 경우 다각형 모서리를 선으로 변환하십시오. 그런 다음 경계에서 25 미터 거리마다 경계에 버퍼를 만듭니다 (중심은 중지 시점을 결정하는 데 도움이 될 수 있음). 그런 다음 폴리곤으로 버퍼를 잘라낸 다음 해당 폴리곤을 래스터로 변환하십시오. 나는 이것이 효과가 있을지 확신하지 못하지만 분명히 당신은 더 적은 RAM을 필요로 할 것입니다. PostGiS는 성능 문제가있을 때 유용합니다.
그것이 조금이라도 도움이되기를 바랍니다 🙂
답변
나는 원래 내 자신의 질문에 대답하지는 않았지만 (이 사이트를 사용하지 않는) 동료가 나에게 내가 한 일을 수행하기 위해 많은 파이썬 코드를 썼다. 육상 세포에 대해서만 해안까지의 거리를 갖도록 세포를 제한하고 해상 세포를 NA로 남겨 두는 것을 포함한다. 다음 코드는 모든 파이썬 콘솔에서 실행할 수 있어야하며 변경이 필요한 유일한 것은 다음과 같습니다.
1) 스크립트 파일을 원하는 모양 파일과 같은 폴더에 넣으십시오.
2) 파이썬 스크립트에서 shapefile 이름을 shapefile 이름으로 변경하십시오.
3) 원하는 해상도를 설정하고;
4) 다른 래스터와 일치하도록 범위를 변경하십시오.
내가 사용하는 것보다 큰 shapefile에는 많은 양의 RAM이 필요하지만 그렇지 않으면 스크립트가 빠르게 실행됩니다 (50m 해상도 래스터를 생성하는 데 약 3 분, 25m 해상도 래스터를 생성하는 데 10 분).
#------------------------------------------------------------------------------
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)