Блог

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

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

Чтобы получить шейп-файл береговой линии (обратите внимание, что в 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)
  

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

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

Вот версия, похожая на версию Роберта, оставаясь в 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)
  

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

 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 мили), что совпадает с другим ответом.

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

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

 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 — это долгота / широта, используйте geosphere, а не rgeos, для вычисления длины

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

(т.е. 1162,8 км)

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