Как извлечь определенные значения из DEM (цифровая модель рельефа)?

Question

Я пытаюсь вычислить данные о высоте для пешеходных маршрутов, используя открытые данные (избегая ограничений для лицензирования, таких как Google).

я смог прочитать публичную DEM моей страны (с разрешением до 10 м) с использованием readGDAL (из пакета RGDAL) и proj4string (mygrid) дает мне:

"+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0" 

Начало файл .asc является:

ncols   9000 
nrows   8884 
xllcorner  323256,181155 
yllcorner  4879269,74709 
cellsize  10 
NODATA_value -9999 
978 998 1005 1008 1012 1016 1020 1025 ..... 
..... 
..... 400 Megabytes of elevation values .... 
..... 

Все, что мне нужно, чтобы забрать из этой сетки данных о высоте конкретных узлов маршрута, так чтобы иметь возможность рассчитать прирост высоты, отрицательный наклон, мин/максимальная высота ...

я привожу данные маршруты из OpenStreetMap, используя хороший пакет Осмар, так что таблица данных моего маршрута что-то вроде этого:

RouteId NodeId lat   lon 
1 -13828 -8754 45.36743 7.753664 
2 -13828 -8756 45.36762 7.753878 
3 -13828 -8758 45.36782 7.754344 
4 -13828 -8760 45.36794 7.754541 
.... 

Но я понятия не имею, как преобразовать широты/долготы координаты в системе координат DEM-координат, а затем, как привести соответствующие значения сетки (сделав вид среднего числа ближайших точек?)

Вся документация, которую я нашел в поисковых системах, предназначен для отображения карт сетки, а не для извлечения значений из них.

Любая помощь будет принята с благодарностью!

Cheers, MB

P.S. Второй вопрос должен быть следующим: «Имея несколько сетчатых черепиц, что бы я мог сделать, если маршрут проходит через две или более плитки? Объедините их, укажите оба ...»

Answer 1

Не преобразовывайте DEM. Преобразуйте свои очки. Ваш DEM проецируется по обычной сетке. Если вы перепроектируете его в другое, это может и не быть. Вместо этого сделайте пространственные точки из ваших данных OSM:

require(sp); coordinates(my.OSM.points) <- long + lat 

Тогда преобразование их в исходной системе координат DEM (заметим, возможно, потребуется установить CRS точек первой Таким образом, вы можете сделать это:.

# Set pojection information for points 
proj4string(my.OSM.points) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84") 

# Transform them 
new.points <- spTransform(my.OSM.points , proj4string(mygrid)) 

Затем рассмотрит использование sp::over или raster::extract, чтобы получить значение DEM в местах точек.

Answer 2

большого спасибо Саймона, ваш answser обратился ко мне, чтобы решить мою проблему!

Просто пунктуации: пытаюсь использовать свой код, я получил:

> coordinates(my.OSM.points) <- routesCoord[,lat,lon] 
Error in coordinates(my.OSM.points) <- routesCoord[, lat, lon] : 
object 'my.OSM.points' not found 

Я думал, что ваше заявление работает, потому что помощь «координата» говорит: "установить пространственные координаты в создать в Пространственный объект "

Однако я использовал:

crsString <- "+proj=longlat +datum=WGS84 +ellps=WGS84" 
my.OSM.points <- SpatialPoints(routesCoord[,lat,lon], 
           proj4string=CRS(crsString), bbox = NULL) 
new.points <- spTransform(my.OSM.points , CRS(proj4string(mygrid))) 

Для новичков вроде меня: после этого вы получите хороший кадр данных с данными высоты для всех точек просто с

my.elev <- over(new.points, mygrid) 

P.S. Извините, Саймон, это мой первый вопрос о stackoverflow, поэтому, когда я попытался проголосовать за ваш ответ, Я получил «голосование требует 15 репутации»! Теперь мне 11, как только я набираю еще 4 очка, я вернулся :-)

PS2 Я хотел бы знать, если «over» принимает значение ближайшей точки сетки или вычисляет средневзвешенный ближайших точек ...

PS3 Второй вопрос должен быть: «Имея несколько сеток-плитки, , что я мог сделать, если маршрут через два или более плиток Объединить их, ссылка и ...»

Answer 3

Как сказал Саймон, если вы находитесь в той же проекции вы можете использовать:

library(raster) 
my.pts.elevation <- extract(my.grid,my.point)