Landscape Metrics and Percent Land Use of Polygons

To keep things simple, I’m going to use a shapefile with only one property in it. I used to work at Rockefeller State Park Preserve (RSPP) and thought it would be interesting to see what the different land cover classifications come up. The property includes fields, forests, and various wetlands.

rspp <- readOGR("RSPP.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/stevendifalco/Documents/GitHub/website_copy/content/post/2021-12-14-landscapemetrics/RSPP.shp", layer: "RSPP"
## with 1 features
## It has 22 fields

map <- ggplot(rspp, aes(x = long, y = lat, group = group)) +
  geom_polygon(color = "black", size = 0.1, fill = "lightgrey") +
  ggtitle("RSPP Map")+
  coord_equal() +
  theme_minimal()

print(map)

NLCD 2016 raster will be used for the land classifications.

nlcd <- raster("NLCD_Clipped.tif")

## get information about the raster itself
#check_landscape(nlcd) 

## to speed up the processes, you can clip the raster to the extent of the shapefile. I already did this step previously but the code would be.
#nlcd_cropped <- crop(nlcd, *area wanted to crop to*)

To mitigate problems with the projections not lining up, use spTransform to match datasets. Afterwards, the next steps will calculate the percentage of each land cover in the raster. I pivot the table and rename the columns so it’s easier to understand the output of this calculation.

#match the spatial projectsion between shapefile and NLCD dataset 
rspp <- spTransform(rspp, proj4string(nlcd))

#calculates % land cover for each type in polygon extent
rspp_nlcd <- sample_lsm(landscape = nlcd, y= rspp, what= c("lsm_c_pland"))

head(rspp_nlcd)

## # A tibble: 6 × 8
##   layer level class    id metric   value plot_id percentage_inside
##   <int> <chr> <int> <int> <chr>    <dbl>   <int>             <dbl>
## 1     1 class    11    NA pland   3.32         1              99.9
## 2     1 class    21    NA pland   1.14         1              99.9
## 3     1 class    22    NA pland   0.409        1              99.9
## 4     1 class    23    NA pland   0.0455       1              99.9
## 5     1 class    41    NA pland  77.0          1              99.9
## 6     1 class    42    NA pland   0.0910       1              99.9

LandscapeMetric analyzes each polygon seperately for each class in the raster. Since the shapefile was a multipolygon, the tool automatically assigned a random number to each of my polygons (plot_id above). Seeing this, I’ll go back and add my own ID so I know which is which polygon from the analysis. I have to convert the shapefile to be an sf object to use st_cast() to extract the multipolygon to singular polygons. Then assign a unique # to each polygon.

# convert shapefile to sf object
rspp_sf <- as(rspp, "sf")

# cast the multi to singular polygons
rspp_cast <- st_cast(rspp_sf, "POLYGON")

#add unique ID to each polygon
rspp_cast$ID <- 1:nrow(rspp_cast)

#lets view the differnce between this and the first map
map <- ggplot() + 
  geom_sf(data = rspp_cast, size = 1, color = rspp_cast$ID, fill = rspp_cast$ID) + 
  ggtitle("RSPP Map") + 
  coord_sf() +
  theme_minimal()

Unfortunately, this shapefile had multiple small slivers in it. There are multiple roads surrounding this property and it is made up of multiple tax parcels, which caused this issue. Not really a bad problem, we’ll work with it as it is for now.

print(map)

To calculate the percentage of each land cover, per unit of the park, we’ll use the code below. Pivoting the table is necessary to reconfigure the data so that each row is only one of the park units. The map is colored by the percentage of ‘Deciduous’ forests of each area.

#calculates % land cover for each type in polygon extent
rspp_nlcd <- sample_lsm(landscape = nlcd, y= rspp_cast, plot_id = rspp_cast$ID, what= c("lsm_c_pland"))

#pivot to wide format consolidating by unique ID
nlcd_wide <- pivot_wider(data= rspp_nlcd, names_from = class, values_from = value)

#rename land cover types
nlcd_wide<- nlcd_wide %>%
  rename(
    Open_Water = '11',
    Dev_OS = '21', 
    Barren = '31', Deciduous = '41', Evergreen = '42', MixedForest = '43',
    Shrub_Scrub = '52', Grassland = '71', Pasture_Hay = '81',
    WoodyWet = '90', EmergentWet = '95'
  )

#join to unique IDS
nlcd_rspp <- inner_join(rspp_cast, nlcd_wide, by= c("ID"="plot_id"))

nlcd_rspp <- nlcd_rspp[nlcd_wide$percentage_inside > 0]

deciduous<- ggplot(nlcd_rspp) +
  geom_sf(aes(fill = Deciduous)) +
  scale_fill_viridis_c(option = "B")

print(deciduous)