Выделение линии из многоугольника и измерение ее длины - PullRequest
1 голос
/ 23 апреля 2019

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

Я следовал этому чрезвычайно полезному уроку по расчету длины береговой линии с использованием линейок разной длины. Тем не менее, он точен только в том случае, если вы хотите измерить всю длину многоугольника, то есть интересующий вас географический объект - это остров.

Чтобы получить шейп-файл береговой линии (обратите внимание, что в вызове ne_countries я специально использую грубую шкалу, чтобы сделать береговую линию более гладкой, и сохраняю только первую возвращенную фигуру - имя "scalerank" не важно):

library(raster)
library(sf)
library(rnaturalearth)
basemap <- rnaturalearth::ne_countries(scale = 110, country = "united states of america", returnclass = "sf")[1]
bbox <- extent(-82, -65, 27, 35) 
cropmap <- st_crop(basemap, bbox)
plot(cropmap)

enter image description here

Возвращает форму, показывающую побережье Южной Атлантики во Флориде. Однако, если я измерим длину этой фигуры, она будет включать все стороны многоугольника, а не только береговую линию. Как изолировать береговую линию (см. Ниже карту того, какая часть многоугольника на самом деле является прибрежной) и просто измерить ее длину в R?

ggplot() +
  geom_sf(data=basemap) +
  geom_sf(data=cropmap, color="blue", fill="blue")

enter image description here

Ответы [ 2 ]

4 голосов
/ 25 апреля 2019

Вот версия, похожая на версию Роберта, оставшуюся в sf.Я преобразовал исходную базовую карту из полигонов в линии, а затем «вырезал» bbox, как вы сделали.Тогда вы можете просто измерить с помощью st_length.

library(sf)
library(rnaturalearth)
basemap <- rnaturalearth::ne_countries(scale = 110, country = "united states of america", returnclass = "sf")[1]

Преобразование из многоугольника в линии, чтобы избежать тех частей многоугольника, которые вам не нужны.

basemap_lines <- basemap %>% st_cast("MULTILINESTRING")
plot(basemap_lines)

Базовая карта какmultilinestring Затем создайте многоугольник печенья, установив crs в широту / долготу

xmin <- -82
xmax <- -65
ymin <- 27
ymax <- 35
bbox <- st_polygon(list(rbind(c(xmin,ymin), c(xmin,ymax), c(xmax,ymax), c(xmax,ymin),c(xmin,ymin)))) %>% st_sfc()
st_crs(bbox) <- 4326

Затем просто обрежьте с помощью st_intersection

crop_lines <- st_intersection(basemap_lines, bbox)
plot(crop_lines)
st_length(crop_lines)

Обрезанная линия

Для ваших bbox размеров я получаю 1162849 метров (~ 723 миль), что соответствует другому ответу.

2 голосов
/ 25 апреля 2019

Это довольно легко, если вы преобразуете полигоны в линии перед кадрированием

Пример данных

library(raster)
library(rnaturalearth)
m <- rnaturalearth::ne_countries(scale = 110, country = "united states of america", returnclass = "sp")
ext <- extent(-82, -65, 27, 35) 

Преобразование и обрезка

m <- as(m, "SpatialLines")   
croplines <- crop(m, ext)

Поскольку CRS являетсядолгота / широта, используйте геосферу, а не rgeos, чтобы вычислить длину

library(geosphere)
lengthLine(croplines)
#[1] 1162849

(т.е. 1162,8 км)

В некоторых случаях вам может потребоваться пересечение с многоугольником вместо экстента (см.raster::drawPoly и raster::intersect) или, возможно, используйте crop, disaggregate и визуально select, чтобы получить нужные вам строки.

...