I'm interested in calculating the area of overlap of different protected area types in a polygon I make of an animal's range using a Brownian Bridge. The shape files come from here https://www.protectedplanet.net (If you want to replicate my example you'd need the South Africa and Mozambique data).
I can calculate the area of the Brownian bridge, the total area of overlap, the areas of the different parks but not the area of overlap with different protected area type. That is I'd like to group by parks$DESIG_ENG
and get those areas of overlap with the brownian bridge.
The output would be an area for National Park and UNESCO-MAB Biosphere Reserve in this case.
This answer looks very close to what I want but I can't quite crack it.
Full code and data below:
#' load packages
library(tidyverse)
library(move)
library(sf)
library(amt)
#' create a move object
loc <-
move(
x = track$x_,
y = track$y_,
time = track$t_,
proj = CRS( "+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"),
data = track,
animal = track$id
)
#' Now create a dBBMM object
dbbmm <-
brownian.bridge.dyn(
object = loc,
location.error = 20,
window.size = 31,
margin = 11,
raster = 10000
)
#' load in the cleaned version which is in Albers Equal Area
#' this is the taken from the protectedplanet and projected to albers equal area
merged_Africa_tranform2021 <- read_sf("shapefile//merged_Africa_protected_clean.shp")
st_crs(merged_Africa_tranform)
#' area of bb
st_area(hr_isopleths(dbbmm$layer, level = 0.95)) / 1e6 #' 11342.2
#' what's the area of overlap with the protected areas?
intersection_1_95 <-
st_intersection(hr_isopleths(dbbmm$layer, level = 0.95),
merged_Africa_tranform_2021$geometry)
sum(st_area(intersection_1_95)) / 1e6
#' intersection of the home range with the protected area
parks <- st_intersection(hr_isopleths(dbbmm$layer, level = 0.95),merged_Africa_tranform_2021)
parks$area / 1e6 # this gives the total area of overlap
parks$DESIG_ENG # this gives the groups I'd like to measure
parks$GIS_AREA # this gives the total areas of the protected areas
Here's a sample of my data:
track <- structure(list(x_ = c(712247.593901045, 712249.63599448, 712255.771966198,
712102.501818369, 711909.388675105, 712218.823688024, 712221.896519337,
712221.899749721, 712215.750856712, 712218.830148765, 712220.862551095,
712218.813996915, 711592.258625023, 711214.431947071, 705232.573529737,
704484.660505673, 704574.525571441, 704586.833330272, 704577.582526676,
703839.386231732, 703841.666220748, 705268.479685813, 703601.261667574,
703466.806636679, 703228.446726644, 688054.891827718, 688049.142018659,
688048.879997468, 688047.436483215, 688049.738642599, 688288.470079595,
688298.728382841, 688283.331520927, 688287.233689746, 688286.411520356,
693251.060727369, 691895.94267286, 695934.025526456, 704723.981466042,
704653.085222894, 709478.368513605, 724297.15313828, 737135.325658038,
743123.408321197, 743107.492642517, 743122.890739821, 743120.844240733,
743110.6083413, 743101.349741506, 743280.634206655, 741686.447716499,
740751.050795124, 723565.872650844, 726424.620739435, 742113.896661416,
763891.653299725, 771306.099600916, 784410.706252694, 784420.984524396,
784651.673216051, 784652.197601293, 784650.759834527, 786713.465194893,
776131.160177772, 791978.497428501, 771910.815530644, 777319.113293987,
823508.490325819, 805173.793218746, 772941.448517303, 765688.241223667,
765329.9021302, 765336.071269917, 765329.306439016, 765326.653694197,
765331.229368852, 765963.052544134, 764395.637977457, 764053.609010882,
771278.483883322, 760497.163530225, 726089.048424295, 715548.19436321,
686220.606351535, 682643.839100749, 682584.383499343, 682640.734139815,
682639.730992916, 682643.834403012, 682647.950340522, 682640.765457917,
682639.718465694, 675928.456251057, 666858.343134657, 657227.261135399,
622646.808919555, 609555.364548704, 619466.706883248, 619493.70597879,
619488.561824205), y_ = c(-2897619.44515289, -2897621.63440962,
-2897624.92738795, -2897498.93857344, -2897683.94074959, -2897649.92496668,
-2897649.93406015, -2897648.84246531, -2897649.91587326, -2897647.74177681,
-2897653.20581315, -2897653.19975086, -2896163.41378075, -2895770.36567857,
-2888073.58655877, -2888910.76395045, -2888920.36769288, -2888914.94243702,
-2888925.83792632, -2888108.83180636, -2883175.87994117, -2874770.83015247,
-2858136.04006769, -2858222.11584159, -2858106.49242823, -2827409.03751134,
-2827628.37797414, -2827630.43368107, -2827627.64191529, -2827628.74527147,
-2827418.47954615, -2827414.12170738, -2827423.94881294, -2827420.88895306,
-2827421.76403575, -2832216.16545687, -2822458.43781159, -2822065.22206298,
-2804439.46921934, -2804855.57191824, -2772317.88886974, -2747820.05265006,
-2736444.12926355, -2739205.61405997, -2739214.94124829, -2739207.26710426,
-2739206.15769513, -2739201.71373307, -2739212.7161063, -2739053.04418682,
-2723051.92607288, -2723903.6891992, -2699199.77940061, -2706451.23477482,
-2682789.33824453, -2661478.09448947, -2655059.92404033, -2657267.70538243,
-2657259.97781768, -2657576.71095908, -2657573.10939881, -2657575.4825813,
-2657295.17114333, -2635355.85939375, -2639186.06073835, -2638488.25825561,
-2668239.53429157, -2666989.13461943, -2689928.96222006, -2677759.03523352,
-2669617.1115024, -2667523.0504592, -2667517.52994131, -2667524.36967127,
-2667512.84632221, -2667520.20547296, -2666056.86632241, -2645303.64873537,
-2650660.65859159, -2663045.04796073, -2691436.8095449, -2725953.49148845,
-2714426.12496965, -2714588.2568858, -2711445.14971251, -2711445.53371262,
-2711455.63848273, -2711447.90058841, -2711446.80721079, -2711441.29383709,
-2711444.58849718, -2711452.32058259, -2715498.19840585, -2715635.30879854,
-2702724.33528911, -2696489.11677258, -2663022.66842587, -2632540.6036741,
-2632608.40346369, -2632615.0522259), t_ = structure(c(1480132500,
1480135560, 1480138680, 1480141680, 1480144020, 1480146300, 1480148760,
1480151160, 1480153620, 1480156080, 1480158840, 1480161780, 1480164840,
1480167900, 1480169880, 1480173480, 1480218360, 1480220940, 1480223640,
1480226160, 1480229760, 1480233360, 1480236960, 1480240560, 1480244160,
1480247760, 1480251360, 1480254960, 1480258560, 1480262160, 1480303740,
1480305660, 1480307640, 1480311240, 1480314840, 1480318440, 1480322040,
1480325640, 1480329240, 1480332840, 1480336440, 1480340040, 1480343640,
1480347240, 1480390500, 1480393200, 1480395960, 1480398780, 1480401300,
1480404900, 1480408500, 1480412100, 1480415700, 1480419300, 1480422900,
1480426500, 1480430100, 1480433700, 1480475340, 1480478040, 1480481640,
1480485240, 1480488840, 1480492440, 1480496040, 1480499640, 1480503240,
1480506840, 1480510440, 1480514040, 1480517640, 1480521240, 1480563480,
1480567080, 1480570680, 1480574280, 1480577880, 1480581480, 1480585080,
1480588680, 1480592280, 1480595880, 1480599480, 1480603080, 1480606680,
1480647600, 1480650780, 1480654380, 1480658100, 1480660020, 1480662000,
1480663920, 1480667520, 1480671120, 1480674720, 1480678320, 1480681920,
1480685520, 1480689120, 1480692720), tzone = "UTC", class = c("POSIXct",
"POSIXt")), id = c("ID2", "ID2", "ID2", "ID2", "ID2", "ID2",
"ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2",
"ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2",
"ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2",
"ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2",
"ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2",
"ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2",
"ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2",
"ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2",
"ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2",
"ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2", "ID2",
"ID2", "ID2", "ID2", "ID2")), row.names = c(NA, -100L), class = ("data.frame"))
st_drop_geometry() |> group_by(DESIGN_ENG) |> summarise(total = sum(GIS_AREA))
?