1

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"))
1
  • 1
    do you want to st_drop_geometry() |> group_by(DESIGN_ENG) |> summarise(total = sum(GIS_AREA))?
    – Elio Diaz
    Commented Jan 11, 2022 at 17:54

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.