Я пытаюсь вычислить площадь (кв. Км или миль) пересечения округов и водоразделов в R, используя растровый пакет и функцию площади.
Мой код выглядит так:
counties <- readOGR('C:\\Shapefiles\\tl_2017_us_county\\tl_2017_us_county.shp')
counties <- spTransform(counties, CRS("+init=epsg:3455"))
huc2_10 <- readOGR('C:\\Shapefiles\\WBD_10_HU2_Shape\\Shape\\WBDHU6.shp')
huc2_10 <- spTransform(huc2_10, CRS("+init=epsg:3455"))
Затем я пересекаю два шейп-файла:
pi <- raster::intersect(huc2_10, counties)
Единицы этой проекции обычно в метрах (я полагаю), поскольку это проекция NAD83 для южной части Южной Дакоты, поэтому функция площади должна вычислять площадь в квадратных метрах. Я пытаюсь вычислить площадь (в квадратных милях) каждого многоугольника, который сформирован в результате этого пересечения, используя функцию площади.
pi$area <- area(pi)/2.589988e6
Однако строка proj4 выглядит следующим образом:
+init=epsg:3455 +proj=lcc +lat_1=44.4 +lat_2=42.83333333333334 +lat_0=42.33333333333334 +lon_0=-100.3333333333333 +x_0=600000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs +ellps=GRS80 +towgs84=0,0,0
В соответствии с этим, единицы измерения - "us-ft". Итак, из-за этого работает ли площадь для каждого полигона в квадратных футах? Кажется, это имеет смысл, но я хотел бы подтвердить это, изменив свой код на:
pi$area <- area(pi)/5280**2
Спасибо.