сбор данных о высоте - извлечение геолокации из файла геотифов - PullRequest
Новая
       50

сбор данных о высоте - извлечение геолокации из файла геотифов

1 голос
/ 19 февраля 2020

Я пытаюсь добавить данные высот на график, используя пакет rayshader. Я могу построить область, в которой я хочу найти данные высот, используя пакет leaflet. 

library(whitebox)
library(leaflet)
library(rayshader)
library(rayrender)
library(raster)

# define bounding box with longitude/latitude coordinates
bbox <- list(
  p1 = list(long = -3.6525599, lat = 40.4065001),
  p2 = list(long = -3.7525599, lat = 40.4965001)
)

leaflet() %>%
  addTiles() %>% 
  addRectangles(
    lng1 = bbox$p1$long, lat1 = bbox$p1$lat,
    lng2 = bbox$p2$long, lat2 = bbox$p2$lat,
    fillColor = "transparent"
  ) %>%
  fitBounds(
    lng1 = bbox$p1$long, lat1 = bbox$p1$lat,
    lng2 = bbox$p2$long, lat2 = bbox$p2$lat,
  )

Я слежу за этим блогом: https://wcmbishop.github.io/rayshader-demo/ и в данный момент на этой части:

«Загрузка данных о высоте»

Автор использует USGS данные, полученные по этой ссылке: 

https://elevation.nationalmap.gov/arcgis/rest/services/3DEPElevation/ImageServer/exportImage?bbox=-122.522%2C37.707%2C-122.354%2C37.84&bboxSR=4326&size=600%2C480&imageSR=4326&time=&format=jpgpng&pixelType=F32&noData=&noDataInterpretation=esriNoDataMatchAny&interpolation=+RSP_BilinearInterpolation&compression=&compressionQuality=&bandIds=&mosaicRule=&renderingRule=&f=html

Я заменил лат, длинные координаты в ссылке на mine: 

https://elevation.nationalmap.gov/arcgis/rest/services/3DEPElevation/ImageServer/exportImage?bbox=-3.6525599%2C40.4065001%2C-3.7525599%2C40.4965001&bboxSR=4326&size=600%2C480&imageSR=4326&time=&format=jpgpng&pixelType=F32&noData=&noDataInterpretation=esriNoDataMatchAny&interpolation=+RSP_BilinearInterpolation&compression=&compressionQuality=&bandIds=&mosaicRule=&renderingRule=&f=html

, которая вернула пустое изображение - что неудивительно, поскольку USGS является поставщиком данных США. Поэтому я загружаю следующий файл:

DTM Spain Материк 20м http://data.opendataportal.at/dataset/dtm-spain/resource/38816df0-9d50-476f-832e-0f7f5fa21771

Файл имеет размер 2,4 ГБ и является файлом для всей Испании, но Мне нужны только данные высот для границ моих точек широты и долготы. 

rgdal::GDALinfo("DTM Spain_Mainland (2019) 20m.tif")
rows        48394 
columns     58187 
bands       1 
lower left origin.x        -23380 
lower left origin.y        3901190 
res.x       20 
res.y       20 
ysign       -1 
oblique.x   0 
oblique.y   0 
driver      GTiff 
projection  +proj=utm +zone=30 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
+no_defs 
file        DTM Spain_Mainland (2019) 20m.tif 
apparent band summary:
   GDType hasNoDataValue NoDataValue blockSize1 blockSize2
1 Float32           TRUE      -32767          1      58187
apparent band statistics:
         Bmin       Bmax Bmean Bsd
1 -4294967295 4294967295    NA  NA
Metadata:
AREA_OR_POINT=Point

Мой вопрос: как я могу получить геотиф для длинных и широтных позиций в bbox или отфильтровать текущий DTM Spain_Mainland (2019) 20m.tif Файл до длинных широтных координат?

1 Ответ

1 голос
/ 19 февраля 2020

Вы можете взглянуть на пакет elevatr, который даст вам нужный растр: 

library(elevatr)
library(raster)

bbox2 <- data.frame(x = c(-3.6525599, -3.7525599), y = c(40.4065001, 40.4965001)) 

elev <- get_elev_raster(bbox2, z = 13, clip = "bbox",
                        prj = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")

enter image description here

...