경도 및 위도별로 고도 추출
뉴질랜드지도에 경도 위도에 따른 큰 좌표 파일이 있습니다. 각 지점에서 대략적인 고도를 찾고 싶습니다. 하나의 영역으로 시작하여 필요한 데이터가있는 래스터를 찾았습니다. 나는 그것을 QGIS에로드했고, 나에게 좋아 보인다. 다음과 같은 정보가 있습니다.
Name NZDEM_SoS_v1-0_27_Dunedin_gf
Path C:\...\elevation\kx-27-dunedin-15m-dem-nzsosdem-v10-GTiff\NZDEM_SoS_v1-0_27_Dunedin_gf.tif
CRS EPSG:2193 - NZGD2000 / New Zealand Transverse Mercator 2000 - Projected
Extent 1372000.0000000000000000,4866000.0000000000000000 : 1492000.0000000000000000,5046000.0000000000000000
Unit meters
Width 8000
Height 12000
Data type Float32 - Thirty two bit floating point
GDAL Driver Description GTiff
GDAL Driver Metadata GeoTIFF
여태까지는 그런대로 잘됐다. 이제 R에서 패키지 래스터를 설치했으며 다음을 수행 할 수 있습니다.
fname = "../elevation/kx-27-dunedin-15m-dem-nzsosdem-v10-GTiff/NZDEM_SoS_v1-0_27_Dunedin_gf.tif"
elev.r <- raster(fname)
이것은 나에게이 경고를 주지만 아마도 그것은 문제가되지 않을 것입니다 :
경고 메시지 : In showSRID (uprojargs, format = "PROJ", multiline = "NO") : CRS 정의에서 GRS80 타원체를 기반으로 알 수없는 데이터가 삭제되었지만 + towgs84 = 값은 유지됨
그럼 할 수 있어요
extract(elev.r,1000,1000)
그리고 그것은 아마도 고도 인 1125.455 값을 반환합니다.
내 경도, 위도를 추출 함수가 이해할 수있는 x, y로 변환하려면 어떻게해야합니까?
여기에서 래스터를 다운로드했습니다. https://koordinates.com/my/downloads/2000967/download/?dl
long = 170.605375
lat = -45.859668
xy <- cbind(lat,long)
colnames(xy) <- c('x', 'y')
xy <- as.data.frame(xy)
coordinates(xy) <- ~ x + y # telling R these are spatial points
crs(xy) <- crs(elev.r) # set the same crs as in your_raster
crs(xy)
extract(elev.r, xy)
xy에있는 것처럼 보이는 모든 NA를 반환합니다.
답변
내 경도, 위도를 추출 함수가 이해할 수있는 x, y로 변환하려면 어떻게해야합니까?
포인트는 다음 중 하나 여야합니다 (x, y 포인트는라는 객체에로드되고 xy고도 래스터는 라는 객체에로드된다고 가정합니다 your_raster).
- 2 열 행렬 또는 data.frame으로 표시되는 포인트. 즉, 다음 레이아웃으로 형식을 지정해야합니다.
library(raster)
# read your raster here
your_raster <- raster("path/to/the/elevation/raster")
# creating a data.frame with the x,y data
xy <- data.frame(x = seq(1400000, 1450000, by = 10000),
y = seq(4900000, 5000000, by = 100000))
print (xy)
# x y
# 1 1400000 4900000
# 2 1410000 5000000
# 3 1420000 4900000
# 4 1430000 5000000
# 5 1440000 4900000
# 6 1450000 5000000
extract(your_raster, xy)
- SpatialPoints 또는 SpatialPointsDataframe
library(sp)
library(raster)
# read your raster here
your_raster <- raster("path/to/the/elevation/raster")
# creating a data.frame with the x,y data
xy <- data.frame(x = seq(1400000, 1450000, by = 10000),
y = seq(4900000, 5000000, by = 100000))
print (xy)
# x y
# 1 1400000 4900000
# 2 1410000 5000000
# 3 1420000 4900000
# 4 1430000 5000000
# 5 1440000 4900000
# 6 1450000 5000000
coordinates(xy) <- ~ x + y # telling R these are spatial points
crs (xy) <- crs(your_raster) # set the same crs as in your_raster
extract(your_raster, xy)
게시 한 점 longlat/WGS84은 래스터가 다른 투영에 있는 반면에 표시됩니다. 추출하기 전에 점을 래스터의 crs로 변환 할 수 있습니다.
xy = data.frame(x=170.605375, y=-45.859668)
coordinates(xy) <- ~ x + y
crs(xy) <- CRS("+proj=longlat +datum=WGS84")
xy <- spTransform(xy, crs(your_raster))
extract(your_raster, xy)