R Преобразование контурных линий в график высот

Я хотел бы иметь возможность создавать график высот из контурных линий в R. Я очень новичок в использовании файлов форм.

На данный момент я скачал данные с здесь который предоставляет файлы .shp для всей Великобритании.

На нем также представлены изолинии, обобщающие топологию Великобритании.

Для графика высот я хотел бы, чтобы data.frame или data.table равномерно расположенных точек (100 м друг от друга) давали выходные данные, дающие значения x, y и z. Где x и y представляют широту и долготу (или восток и север), а z представляют высоту (в метрах над уровнем моря).

Я думаю, что, вероятно, есть некоторые инструменты, которые автоматически выполнят интерполяцию для вас, но я не уверен, как это будет работать с геопространственными данными.

Это мой основной старт...

require(maptools)
xx <- readShapeSpatial("HP40_line.shp")

r sp gis rgdal maptools
person h.l.m    schedule 27.02.2015    source источник
comment
Вы должны напрямую загрузить версию сетки и использовать пакет raster для масштабирования данных и получения SpatialPointsDataFrame.   -  person    schedule 27.02.2015
comment
Не могли бы вы показать мне, как это будет сделано?   -  person h.l.m    schedule 27.02.2015
comment
Я не уверен, что такое версия сетки, но не думаю, что она доступна...   -  person h.l.m    schedule 27.02.2015

Ответы (1)


arrow_upward
3
arrow_downward

Выберите «ASCII Grid and GML (Grid)» в качестве формата загрузки для продукта «OS Terrain 50» и загрузите файл. Это даст вам zip-файл, содержащий множество каталогов zip-файлов, каждый из которых содержит части 50-метровой сетки высот Великобритании (часть, на которую я смотрел, имела 200 x 200 ячеек, что означает 10 км x 10 км). Я зашел в каталог data/su, распаковал туда zip-файл и сделал

library(raster)
r = raster("SU99.asc")
plot(r)

чтобы объединить это в сетку 100 м, я сделал

r100 = aggregate(r) # default is factor 2: 50 -> 100 m

Как упоминалось выше, совет заключается в том, чтобы работать с сетками, поскольку контурные линии получаются из сеток, работа в обратном направлении — это болезненная и большая потеря информации.

Получение значений сетки долготы и широты в виде data.frame можно выполнить двумя способами:

df = as.data.frame(projectRaster(r, crs = CRS("+proj=longlat")), xy = TRUE)

разворачивает сетку в новую сетку по долготе/широте. Поскольку эти сетки не могут совпадать, он минимально перемещает точки (см. ?projectRaster).

Второй вариант — преобразовать сетку в точки и отменить проецирование их на долготу и широту с помощью

df2 = as.data.frame(spTransform(as(r, "SpatialPointsDataFrame"), CRS("+proj=longlat")))

Это не перемещает точки и, как следствие, не создает сетку.

person Edzer Pebesma    schedule 28.02.2015
comment
Я не вижу файл типа .grd... я вижу только .asc .gml и .prj - person h.l.m; 01.03.2015
comment
Также как мне получить data.frame с тремя столбцами широты, долготы и высоты? - person h.l.m; 01.03.2015
comment
вы правы: расширение .asc, а не .grd, я поправил ответ выше; вы можете получить data.frame за as.data.frame(r100, xy = TRUE) - person Edzer Pebesma; 01.03.2015
comment
Спасибо за это, это очень полезно. Итак, я заметил, что при выполнении as.data.frame(r100, xy = TRUE) значения x и y не являются широтой и долготой, как мне преобразовать их в такие... - person h.l.m; 01.03.2015
comment
Хорошо, я добавил это к ответу. - person Edzer Pebesma; 01.03.2015