Как буферизовать / расширить шейп-файл переписи на 1 милю с помощью rgeos::gBuffer?

Question

Я пытаюсь расширить / увеличить шейп-файл народа чероки, который я загрузил из набора данных переписи здесь: https://www2.census.gov/geo/tiger/TIGER2019/AIANNH/, используя пакет rgeos. Код, который у меня есть для этого сейчас:

library(rgeos)
library(dplyr)

tribe_shp <- readOGR("/location of file/", layer = "tl_2019_us_aiannh")
tribe_shp <- tribe_shp %>% filter(GEOID == "5550R") #filter to cherokee nation

expand_obj <- gBuffer(tribe_shp, byid = F, width = 1000)

>Warning message:
In gBuffer(tribe_shp, byid = F, width = 1000) :
  Spatial object is not projected; GEOS expects planar coordinates

plot(expand_obj)

Результирующий объект теряет фрейм данных исходного SPDF и представляет собой в значительной степени просто круг, который совсем не похож на исходную форму. Что-то мне не хватает?

Answer 1

gBufferиспользует единицы измерения данных. В этом случае данные представлены в градусах широты и долготы, поэтому ширина буфера составляет 1000 градусов. Чтобы буфер в метрах преобразовал в другую систему координат в метрах.

Существует множество систем координат, и вы действительно должны найти ту, которая подходит для вашего местоположения. Я думаю, что в США есть ряд систем, разработанных для каждого штата, так что это, вероятно, будет лучше всего. Но пока я буду использовать EPSG:3857, который Google использует для карт, и на самом деле он не так точен.

Прочтите данные:

tribe_shp <- readOGR("./", layer = "tl_2019_us_aiannh")

Подмножество с использованием выбора - dplyr::filter здесь не работает, но это будет:

tribe_shp = tribe_shp[tribe_shp$GEOID=="5550R",]

Теперь перейдем к другой системе координат:

tribe_shp_trans = spTransform(tribe_shp, "+init=epsg:3857")

и сделать буфер 1км. Если вам нужен буфер в 1 милю, используйте столько метров в миле - тысяча шестьсот с чем-то?

tribe_shp_buf = gBuffer(tribe_shp_trans, width=1000)

Если вы построите их, вы увидите, что буфер больше, чем исходная область:

plot(tribe_shp_trans)
plot(tribe_shp_buf,lty=3,add=TRUE)

Деталь сюжета:

Если вам нужен буфер в формате lat-long, сделайте еще одно преобразование буфера в "+init=epsg:4326".

Вы также можете сделать это с помощью более современных sf пакет, используйте st_read читать данные, st_transform перепроецировать, и st_bufferделать буфер. Это должно быть даже быстрее.