Как извлечь конкретные значения из матрицы высот (цифровая модель рельефа)?
Я пытаюсь рассчитать данные высоты для пешеходных маршрутов, используя открытые данные (избегая лицензионных ограничений, таких как Google).
Я смог прочитать общедоступную ЦМР своей страны (с разрешением 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, используя симпатичный пакет OSMAR, поэтому таблица данных моего маршрута выглядит примерно так:
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, а затем как вывести соответствующие значения сетки (что-то вроде усреднения ближайших точек?)
Вся документация, которую я нашел, гуглит, что она собирается визуализировать карты сетки, а не извлекать значения из них.
Любая помощь будет оценена!
Ура, мб
PS Второй вопрос должен звучать так: "Имея несколько мозаичных ячеек, что я могу сделать, если маршрут пересекает две или более тайлов? Объединяй их, ссылаясь на оба…"
3 ответа
Не трансформируйте DEM. Преобразуйте свои очки. Ваша ЦМР проецируется на регулярную сетку. Если вы перепроектируете это в другое, это может не быть. Вместо этого сделайте пространственные точки из ваших данных OSM:
require(sp); coordinates(my.OSM.points) <- long + lat
Затем преобразуйте их в систему координат координат матрицы высот (обратите внимание, что вам может понадобиться сначала установить 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 чтобы получить значения матрицы высот в точечных точках.
Как сказал Саймон, если вы находитесь в одной проекции, вы можете использовать:
library(raster)
my.pts.elevation <- extract(my.grid,my.point)
Большое спасибо, Саймон, твой ответчик обратился ко мне, чтобы решить мою проблему!
Просто пунктуация: пытаясь использовать ваш код, я получил:
> 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)
PS Извините, Саймон, это мой первый вопрос по stackru, поэтому, когда я попытался проголосовать за ваш ответ, я получил "для голосования требуется 15 репутаций"! Сейчас мне 11, как только я наберу 4 очка, я вернулся:-)
PS2 Я хотел бы знать, принимает ли значение "over" значение ближайшей точки сетки или вычисляет средневзвешенное значение ближайших точек...
PS3 Второй вопрос должен быть: "Имея несколько плиток сетки, что я мог бы сделать, если маршрут пересекает две или более плиток? Объедините их, ссылаясь на обе…"