Whooping Crane Habitat Quality Assessment
Ifeoma Okonye
2024-05-09
PART 1 - BACKGROUND AND GOALS
Introduction
The whooping crane (Grus americana) is a federally endangered bird species found only in North America (Lewis, T. L., Jr., 1997) and it faces critical challenges during its long migratory journey (Armbruster, M.J. 1990). As one of the rarest species, by the late 1940s, only about 15 individuals remained in the wild pushing the species almost to extinction. However, through combined conservation efforts, such as habitat restoration, captive breeding programs, and different protection measures, there has been a significant increase in crane population over time and their numbers have risen to about 600 today (Canadian Wildlife Service and U.S. Fish and Wildlife Service, 2007). Land use change, predation and hunting pressure during its migration are factors which pushed this species to the brink of extinction (Pearse, Aaron T., et al, 2021 & U.S. Fish and Wildlife Service., 2009).
Key to the survival of the Whooping crane is the preservation and restoration of its habitat, particularly wetlands. Wetlands are an important ecosystem for the species’ breeding, foraging and migration activities. These wetlands not only provide suitable grounds for breeding but also provide a variety of food sources such as aquatic invertebrates, fish, and plant matter which are essential parts of its diet (BirdLife International, 2020). Furthermore, wetlands play a crucial role in facilitating the migratory journeys of Whooping Cranes, serving as critical stopover sites where they rest and refuel during their long-distance travels between breeding and wintering grounds (Armbruster, M.J. 1990; Watershed Institute, Inc. 2013; Cowardin, L.M. et al., 1979; Urbanek, R. P. et al., 2018). During migration, the WC crosses 7 states including Kansas and one of its major stopover habitats in Kansas is the Quivira National wildlife Refuge. In line with previous research efforts within the Aransas-Wood Buffalo Population migratory corridor, this study focuses on modeling habitat suitability for whooping cranes in Quivira National Wildlife Refuge (NWR), Kansas, using a multicriteria decision analysis (MCDA) approach. Multicriteria Decision Analysis (MCDA) stands out as a valuable approach, offering a systematic framework for evaluating and ranking various habitat attributes to inform conservation decision-making (Adem E., Blal, & Davide Geneletti, 2018; Belton, L., 2019). By integrating multiple criteria—such as water regime, wetland size, type, density, and proximity to essential resources—MCDA enables conservationists to assess habitat quality comprehensively and identify areas of critical importance for species conservation (Cowardin, L.M. et al., 1979).
Considering the Whooping Crane’s dependence on wetlands for survival, understanding the dynamics of these ecosystems and their conservation status is important. This research aims to understand the intricate relationship between the Whooping Crane and wetlands, using spatial analysis techniques to assess the quality of wetland habitats within Quivira NWR in Kansas, United States. The results of this research will serve as a critical foundation for conservation efforts within Quivira NWR in the Kansas portion of the Aransas-Wood Buffalo Population migratory corridor. This research informs strategies for development avoidance, habitat improvements, and protective measures, supporting the conservation of the endangered whooping crane. Moreover, it contributes valuable insights into habitat suitability modeling for broader wildlife conservation practices.
Research Goals
The primary objective of my research will be to evaluate the quality and suitability of wetland habitats for Whooping Cranes within the Quivira National Wildlife Refuge (NWR) in Kansas, United States. To achieve this overarching goal, my research aims to address the following specific research objectives:
- Habitat Unsuitability Screening: Identify and exclude areas within the refuge that are unsuitable for Whooping Cranes due to proximity to human disturbances such as roads, railroads, power lines, and bridges.
- Habitat Quality Assessment: Evaluate the quality of the remaining wetlands based on key habitat characteristics that are known to influence Whooping Crane habitat selection. These characteristics include water regime, wetland size, wetland type (natural vs. artificial), wetland density, and distance to food.
- Composite Habitat Quality Scores: Calculate composite habitat quality scores for each wetland by summing the scores of the individual habitat characteristics. This will provide a comprehensive measure of habitat quality that takes into account multiple factors.
- Habitat Mapping: Create a spatial map of the composite habitat quality scores to visually represent the distribution of high-quality habitats within the refuge.
- Statistical Analysis: t-test, correlation matrix and linear regression models were used to analyze and validate the model.
- Conservation Recommendations: Based on the findings, provide recommendations for habitat management and conservation actions that could enhance the quality of Whooping Crane habitats within the Quivira National Wildlife Refuge.
Through these goals, my research aims to contribute to the ongoing conservation efforts for the Whooping Crane by providing valuable information on the quality and distribution of their habitats within a key stopover site in their migration route.
Materials and Methods
Study Area:
The study focuses on the Quivira National Wildlife Refuge (NWR) located in central Kansas, United States. Encompassing approximately 22,135 acres, Quivira NWR is renowned for its diverse wetland habitats, including marshes, lakes, and riparian areas, making it an important stopover and wintering site for migratory birds, including the endangered Whooping Crane (Grus americana).
Data Acquisition:
- Spatial Data: Spatial data layers were obtained from various sources, including:
- Quivira NWR boundary shapefile.
- National Wetland Inventory (NWI) wetland classification data.
- Telemetry data on Whooping Crane movements within the study area.
- Road and railroad networks from the United States Census Bureau.
- Electric power transmission lines and non-state bridges data from Kansas Department of Transportation.
- Public Land Survey System (PLSS) quarter sections shapefile.
- Digital Elevation Model (DEM) raster data for to set the extent of each habitat characteristic.
- Remote Sensing Data: Landsat imagery for the year 2019 was used to derive land cover classification data, specifically Cropland Data Layer (CDL) raster dataset.
The data was then preprocessed to ensure compatibility for analysis. This involved selecting relevant variables, setting coordinate systems, creating a 10km buffer around the study area, and preparing study area layers.
PART 2 - ANALYSIS
All analyses were conducted using R statistical software. The sf, sp, terra, raster, tidyverse, mapview, ggspatial, ggthemes, viridis, scales, grid, rgdal, spdep, and rasterVis packages were used for data manipulation, analysis, and visualization. The pacman package was used for package management.
library(pacman)
p_load(sf, sp, terra, raster, tidyverse, mapview, ggspatial, ggthemes, viridis, scales, grid, rgdal, rasterVis, spdep)Data Preparation and Formatting
All required data, as earlier stated, was imported and the data processing was carried out. The data processing and formatting carried out on the data includes:
Coordinate System Standardization: All spatial datasets were standardized to a consistent coordinate reference system (CRS) using the st_transform() function in the sf package for vector data and the projectRaster() function in the raster package for raster data.
Study Area Definition: A 10-kilometer buffer around the Quivira NWR boundary was created using the st_buffer() function to define the study area boundary.
Tidying and Cropping: Spatial datasets were intersected with the study area boundary to retain only relevant features within the study area. Additionally, raster datasets were cropped to the extent of the study area using the crop() function for vector data and the mask() function for raster data.
#######################################################################################
#loading the data ----
#######################################################################################
Quivira_NWR <- st_read("Data/Quivira_NWR.shp")## Reading layer `Quivira_NWR' from data source
## `C:\Users\ifeom\OneDrive - Kansas State University\Courses\3. HSI Research\KDWP_WhoopingCrane_TWI\Scripts\Data\Quivira_NWR.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 11 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -10972210 ymin: 4589451 xmax: -10959910 ymax: 4612157
## Projected CRS: WGS 84 / Pseudo-MercatorNWI_Code_Definitions <- read.csv("Data/NWI_Code_Definitions.csv")
WHCR_Telemetry <- st_read("Data/Telemetry/WHCR_KS_Telemetry_GPS.shp")## Reading layer `WHCR_KS_Telemetry_GPS' from data source
## `C:\Users\ifeom\OneDrive - Kansas State University\Courses\3. HSI Research\KDWP_WhoopingCrane_TWI\Scripts\Data\Telemetry\WHCR_KS_Telemetry_GPS.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 155051 features and 50 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -1617990 ymin: 456357.2 xmax: 212560.3 ymax: 4803673
## Projected CRS: NAD83 / Conus AlbersKS_Roads <- st_read("Data/Tiger_2020_Roads/Tiger_2020_Roads.shp")## Reading layer `Tiger_2020_Roads' from data source
## `C:\Users\ifeom\OneDrive - Kansas State University\Courses\3. HSI Research\KDWP_WhoopingCrane_TWI\Scripts\Data\Tiger_2020_Roads\Tiger_2020_Roads.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 331170 features and 6 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: -102.0517 ymin: 36.99302 xmax: -94.59291 ymax: 40.00315
## Geodetic CRS: WGS 84KS_Railroads <- st_read("Data/Tiger_2020_Railroads/Tiger_2020_Railroads.shp")## Reading layer `Tiger_2020_Railroads' from data source
## `C:\Users\ifeom\OneDrive - Kansas State University\Courses\3. HSI Research\KDWP_WhoopingCrane_TWI\Scripts\Data\Tiger_2020_Railroads\Tiger_2020_Railroads.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 3719 features and 5 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: -102.0487 ymin: 36.99361 xmax: -94.59363 ymax: 40.00219
## Geodetic CRS: WGS 84quivira_PLSS_q <- st_read("Data/Quivera_PLSS_Quarters.shp")## Reading layer `Quivera_PLSS_Quarters' from data source
## `C:\Users\ifeom\OneDrive - Kansas State University\Courses\3. HSI Research\KDWP_WhoopingCrane_TWI\Scripts\Data\Quivera_PLSS_Quarters.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1496 features and 21 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -98.67487 ymin: 37.98419 xmax: -98.34404 ymax: 38.3047
## Geodetic CRS: WGS 84power_lines <- st_read("Data/Electric_Power_Transmission_Lines.shp")## Reading layer `Electric_Power_Transmission_Lines' from data source
## `C:\Users\ifeom\OneDrive - Kansas State University\Courses\3. HSI Research\KDWP_WhoopingCrane_TWI\Scripts\Data\Electric_Power_Transmission_Lines.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 93047 features and 21 fields
## Geometry type: MULTILINESTRING
## Dimension: XY
## Bounding box: xmin: -159.7793 ymin: 13.44758 xmax: 144.8235 ymax: 65.01721
## Geodetic CRS: WGS 84bridges <- st_read("Data/KDOT_NonStateBridges.shp")## Reading layer `KDOT_NonStateBridges' from data source
## `C:\Users\ifeom\OneDrive - Kansas State University\Courses\3. HSI Research\KDWP_WhoopingCrane_TWI\Scripts\Data\KDOT_NonStateBridges.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 19361 features and 135 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -357346.3 ymin: -166685.8 xmax: 300905 ymax: 172732.8
## Projected CRS: Lambert_Conformal_Conic_2SPLULC_2019 <- raster("Data/Cropland_Data/CDL_2019_20.tif")
DEM <- raster("Data/Qui_DEM.tif")
path <- "Data/KS_shapefile_wetlands"
Wetland_shf <- list.files(path, pattern = "\\.shp$", full.names = TRUE)
Wetland_shpf_list <- lapply(Wetland_shf, st_read)## Reading layer `KS_Wetlands_East' from data source
## `C:\Users\ifeom\OneDrive - Kansas State University\Courses\3. HSI Research\KDWP_WhoopingCrane_TWI\Scripts\Data\KS_shapefile_wetlands\KS_Wetlands_East.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 736676 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -215119.6 ymin: 1536759 xmax: 132411.7 ymax: 1903565
## Projected CRS: NAD_1983_Albers
## Reading layer `KS_Wetlands_West' from data source
## `C:\Users\ifeom\OneDrive - Kansas State University\Courses\3. HSI Research\KDWP_WhoopingCrane_TWI\Scripts\Data\KS_shapefile_wetlands\KS_Wetlands_West.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 326781 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -540452.8 ymin: 1539544 xmax: -204909.5 ymax: 1917443
## Projected CRS: NAD_1983_AlbersKS_Wetlands <- do.call(rbind, Wetland_shpf_list)
#removing intermediate datalayers
#rm(path, Wetland_shf, Wetland_shpf_list)
# Check which elements are NULL
which(sapply(Wetland_shpf_list, is.null))## integer(0)#######################################################################################
#Tidying and Cropping the data ----
#######################################################################################
# selecting variables
glimpse(NWI_Code_Definitions)## Rows: 5,665
## Columns: 23
## $ ATTRIBUTE <chr> "E1AB/UBL", "E1AB/UBLx", "E1AB1L", "E1AB1L6", "E1AB1Lh", "E…
## $ SYSTEM_NAM <chr> "Estuarine", "Estuarine", "Estuarine", "Estuarine", "Estuar…
## $ SYSTEM_DEF <chr> "The Estuarine System consists of deepwater tidal habitats …
## $ SUBSYSTEM_ <chr> "Subtidal", "Subtidal", "Subtidal", "Subtidal", "Subtidal",…
## $ SUBSYSTEM1 <chr> "The substrate in these habitats is continuously covered wi…
## $ CLASS_NAME <chr> "Aquatic Bed", "Aquatic Bed", "Aquatic Bed", "Aquatic Bed",…
## $ CLASS_DEFI <chr> "Includes wetlands and deepwater habitats dominated by plan…
## $ SUBCLASS_N <chr> "", "", "Algal", "Algal", "Algal", "Algal", "Rooted Vascula…
## $ SUBCLASS_D <chr> "", "", "Algal beds are widespread and diverse in the Marin…
## $ SPLIT_CLAS <chr> "Unconsolidated Bottom", "Unconsolidated Bottom", "", "", "…
## $ SPLIT_CL_1 <chr> "Includes all wetlands and deepwater habitats with at least…
## $ SPLIT_SUBC <chr> "", "", "", "", "", "", "", "Sand", "", "", "", "", "", "",…
## $ SPLIT_SU_1 <chr> "", "", "", "", "", "", "", "The unconsolidated particles s…
## $ WATER_REGI <chr> "Subtidal", "Subtidal", "Subtidal", "Subtidal", "Subtidal",…
## $ WATER_RE_1 <chr> "Saltwater Tidal", "Saltwater Tidal", "Saltwater Tidal", "S…
## $ WATER_RE_2 <chr> "Tidal salt water continuously covers the substrate.", "Tid…
## $ FIRST_MODI <chr> "", "Excavated", "", "Oligohaline", "Diked/Impounded", "Exc…
## $ FIRST_MO_1 <chr> "", "Special Modifier", "", "Water Chemistry", "Special Mod…
## $ FIRST_MO_2 <chr> "", "", "", "Coastal Halinity", "", "", "", "", "", "", "Co…
## $ FIRST_MO_3 <chr> "", "This Modifier is used to identify wetland basins or ch…
## $ SECOND_MOD <chr> "", "", "", "", "", "", "", "", "", "", "", "", "", "", "",…
## $ SECOND_M_1 <chr> "", "", "", "", "", "", "", "", "", "", "", "", "", "", "",…
## $ SECOND_M_2 <chr> "", "", "", "", "", "", "", "", "", "", "", "", "", "", "",…NWI_Code_Definitions <- NWI_Code_Definitions %>%
select(ATTRIBUTE, SYSTEM_NAM, WATER_REGI, FIRST_MODI)
crs(Quivira_NWR)## [1] "PROJCRS[\"WGS 84 / Pseudo-Mercator\",\n BASEGEOGCRS[\"WGS 84\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4326]],\n CONVERSION[\"unnamed\",\n METHOD[\"Popular Visualisation Pseudo Mercator\",\n ID[\"EPSG\",1024]],\n PARAMETER[\"Latitude of natural origin\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"False easting\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"easting (X)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"northing (Y)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"Web mapping and visualisation.\"],\n AREA[\"World between 85.06°S and 85.06°N.\"],\n BBOX[-85.06,-180,85.06,180]],\n ID[\"EPSG\",3857]]"crs = "epsg:4326"
# Setting coordinate systems
Quivira_NWR <- st_transform(Quivira_NWR, st_crs(KS_Wetlands))
# KS_Wetlands <- st_transform(KS_Wetlands, st_crs(Quivira_NWR))
# KS_Wetlands <- st_transform(Quivira_NWR, st_crs(KS_Wetlands))
WHCR_Telemetry <- st_transform(WHCR_Telemetry, st_crs(KS_Wetlands))
KS_Roads <- st_transform(KS_Roads, st_crs(KS_Wetlands))
KS_Railroads <- st_transform(KS_Railroads, st_crs(KS_Wetlands))
quivira_PLSS_q <- st_transform(quivira_PLSS_q, st_crs(KS_Wetlands))
power_lines <- st_transform(power_lines, st_crs(KS_Wetlands))
bridges <- st_transform(bridges, st_crs(KS_Wetlands))
LULC_2019 <- projectRaster(LULC_2019, crs = crs(KS_Wetlands))
DEM <- projectRaster(DEM, crs = crs(KS_Wetlands))
# making a 10km buffer around the study area
quivira_10km <- st_buffer(Quivira_NWR, dist = 10000)
# preparing study area layers
quivira_wetland <- st_intersection(KS_Wetlands, quivira_10km) %>%
select(-(ORGNAME:CostCenter))
quivira_telemetry <- st_intersection(WHCR_Telemetry, quivira_10km) %>%
select(-c(PROGRAM:REC_TYPE, FREQ_FLAG:QUALITY, TX_FREQ:CostCenter))
quivira_roads <- st_intersection(KS_Roads, quivira_10km) %>%
select(-(ORGNAME:CostCenter))
quivira_railroads <- st_intersection(KS_Railroads, quivira_10km) %>%
select(-(OBJECTID:CostCenter))
quivira_PLSS <- st_intersection(quivira_PLSS_q, quivira_10km) %>%
select(-(OBJECTID.1:CostCenter))
quivira_powerlines <- st_intersection(power_lines, quivira_10km) %>%
select(-c(NAICS_CODE:SUB_2, OBJECTID.1:GlobalID.1))
quivira_bridges <- st_intersection(bridges, quivira_10km) %>%
select(-c(4:8, 10:20, 23:138))
quivira_LULC <- crop(LULC_2019, quivira_10km)
quivira_DEM <- mask(DEM, quivira_10km)
res(quivira_DEM) # Check the resolution of the DEM## [1] 8.03 10.40extent(quivira_DEM) # Check the extent of the DEM## class : Extent
## xmin : -233429.5
## xmax : -202176.7
## ymin : 1663398
## ymax : 1701878# removing unused variables
#rm(bridges, DEM, LULC_2019, KS_Railroads, KS_Roads, power_lines, Quivira_NWR, quivira_PLSS_q)
#Plotting the study area
studyarea <- ggplot() +
geom_sf(data = quivira_wetland, aes(fill = WETLAND_TY), alpha = 0.7, col = "black") +
geom_sf(data = Quivira_NWR, fill = NA, alpha = 0.6, color = "darkred", lwd = 0.8) +
geom_sf(data = quivira_10km, color = "yellow", fill = NA, lwd = 2) +
labs(title = "Quivira National Wildlife Refuge, Kansas, US") +
theme_minimal()
studyarea
#ggsave("StudyArea.png", plot = stdyarea)Habitat Unsuitability Screening
Areas within the refuge that are unsuitable for Whooping Cranes due to proximity to human disturbances such as roads, railroads, power lines, and bridges were identified and excluded from the analysis. This was achieved by creating buffers around areas of human disturbance and assessing the extent to which these overlap with wetland habitats. Finally, due to the size of individual birds (approximately 5 feet tall and 7 feet wide (wing tip to wing tip)) and the fact that the birds prefer habitats of certain size (Watershed Institute, Inc., 2013), all wetlands of 2.5 acres and below were screened out as unsuitable habitats.
# Identifying Wetlands near human disturbances
## Create buffers from human disturbances
### Buffer distance from different road types
calcBufferDist <- function(MTFCC) {
if (MTFCC %in% c("S1500", "S1740", "S1810", "S1830")) {
return(200)
} else {
return(400)
}
}
quivira_roads <- quivira_roads %>%
mutate(BufferDist = sapply(MTFCC, calcBufferDist))
roads_buffered <- st_buffer(quivira_roads, dist = quivira_roads$BufferDist)
### railroad buffer distance
railroad_buffered <- st_buffer(quivira_railroads, dist = 400, dissolve = TRUE)
### powerlines buffer distance
powerline_buffered <- st_buffer(quivira_powerlines, dist = 100, dissolve = TRUE)
### bridges buffer distance
bridges_buffered <- st_buffer(quivira_bridges, dist = 400, dissolve = TRUE)
ggplot() +
geom_sf(data = quivira_wetland, aes(fill = WETLAND_TY), alpha = 0.7, col = "black") +
geom_sf(data = roads_buffered, col = "steelblue", alpha = 0) +
theme_minimal()
## merge unsuitable locations
unsuitable <- rbind(roads_buffered %>% dplyr::select(geometry), railroad_buffered %>% dplyr::select(geometry), powerline_buffered %>% dplyr::select(geometry), bridges_buffered %>% dplyr::select(geometry))
unsuitable_areas <- st_union(unsuitable)
ggplot() +
geom_sf(data = quivira_wetland, aes(fill = WETLAND_TY), alpha = 0.7, col = "black") +
geom_sf(data = unsuitable_areas, col = "gray", alpha = 0, lwd = .5) +
theme_minimal()
## remove wetlands that fall within disturbance buffers
wetlandsErased <- st_difference(quivira_wetland, unsuitable_areas)
## Convert all different parts to one.
wetlandsSinglepart <- st_cast(wetlandsErased, "POLYGON")
ggplot() +
geom_sf(data = wetlandsSinglepart, aes(fill = WETLAND_TY), alpha = 0.7, col = "black") +
geom_sf(data = quivira_10km, col = "yellow", alpha = 0, lwd = 2) +
theme_minimal()
# Calculate wetland areas in acres
wetlandsSinglepart <- wetlandsSinglepart %>%
mutate(ACRES = st_area(wetlandsSinglepart) * 0.000247105, # From square meters to acres
ACRES = str_remove(ACRES, " \\[m\\^2\\]"), # Remove the unit part
ACRES = as.numeric(ACRES)) # Convert the column to numeric
# Select only wetlands with areas greater than the 2.5 acres
screenedWetlands <- wetlandsSinglepart %>%
filter(ACRES > 2.5)
# Plot screened Wetlands
ggplot() +
geom_sf(data = screenedWetlands, aes(fill = WETLAND_TY), alpha = 0.7, col = "black") +
geom_sf(data = quivira_10km, col = "yellow", alpha = 0, lwd = 2) +
theme_minimal()
Habitat Quality Assessment - Model Variables
The quality of the remaining wetlands was evaluated based on key habitat characteristics that are known to influence Whooping Crane habitat selection. These characteristics include water regime, wetland size, wetland type (natural vs. artificial), wetland density, and distance to food. The model variables were created by setting a function for each habitat characteristic and then using these set functions and spatial analysis to score different wetland types. Each characteristic was scored on a scale, and these scores were used to calculate a composite habitat quality score for each wetland.
# Using a MCDSM approach, model variables that will be used to assess the quality of remaining wetlands after screening include:
# 1. Water Regime
# 2. Wetland Size
# 3. Wetland Type (Natural vs. Artificial)
# 4. Wetland density
# 5. Distance to food
#Adding NWI Definitions
screenedWetlands <- left_join(screenedWetlands, NWI_Code_Definitions, by = "ATTRIBUTE")
#######################################################################################
## Water Regime ---------
### define water regime function
calcWaterRegime <- function(WATER_REGI) {
if (WATER_REGI == "Permanently Flooded") {
return(5)
} else if (WATER_REGI == "Intermittently Exposed") {
return(4)
} else if (WATER_REGI == "Semipermanently Flooded") {
return(3)
} else if (WATER_REGI == "Seasonally Flooded") {
return(2)
} else if (WATER_REGI == "Artificially Flooded") {
return(1)
} else if (WATER_REGI == "Temporary Flooded") {
return(1)
} else {
return(9999)
}
}
### calculate water regime score
WaterRegimeScore <- screenedWetlands %>%
mutate(RegScore = sapply(WATER_REGI, calcWaterRegime))
### Convert to raster
WaterRegimeScore_raster <- rasterize(WaterRegimeScore, quivira_DEM, field = "RegScore", fun = max)
#######################################################################################
## Wetland Size -----------
### define Wetland Size function based on Acres
calcWetlandSize <- function(ACRES) {
if (ACRES >= 7) {
return(5)
} else if (ACRES >= 5 && ACRES < 7) {
return(4)
} else if (ACRES >= 3 && ACRES < 5) {
return(3)
} else if (ACRES >= 1 && ACRES < 3) {
return(2)
} else if (ACRES < 1) {
return(1)
} else {
return(9999)
}
}
### calculate Wetland Size score
WetlandSizeScore <- screenedWetlands %>%
mutate(SizeScore = sapply(ACRES, calcWetlandSize))
### Convert to raster
WetlandSizeScore_raster <- rasterize(WetlandSizeScore, quivira_DEM, field = "SizeScore", fun = max)
#######################################################################################
## Wetland Type (Natural vs. Artificial) -----------
### define Wetland Size function based on Acres
calcWetlandType <- function(FIRST_MODI) {
if (FIRST_MODI == "Diked/Impounded") {
return(0)
} else if (FIRST_MODI == "Excavated") {
return(0)
} else {
return(2)
}
}
### calculate Wetland Type score
WetlandTypeScore <- screenedWetlands %>%
mutate(TypeScore = sapply(FIRST_MODI, calcWetlandType))
### Convert to raster
WetlandTypeScore_raster <- rasterize(WetlandTypeScore, quivira_DEM, field = "TypeScore", fun = max)
#######################################################################################
## Wetland Density ---------
### define water density function
calcWaterDensity <- function(count) {
if (count >= 5) {
return(5)
} else if (count < 5) {
return(0)
} else {
return(9999)
}
}
# Summarize within: counting the number of unique wetlands in each quarter section
WaterDensityScore <- screenedWetlands %>%
st_join(quivira_PLSS, join = st_intersects) %>% # Join operation first
group_by(FID_1) %>% # Group by FID_1 which should now be available
summarise(count = n_distinct(ATTRIBUTE),
ATTRIBUTE = first(ATTRIBUTE),
WETLAND_TY = first(WETLAND_TY),
ACRES = sum(ACRES),
SYSTEM_NAM = first(SYSTEM_NAM),
WATER_REGI = first(WATER_REGI),
FIRST_MODI = first(FIRST_MODI),
geometry = first(geometry)) %>%
mutate(DensityScore = sapply(count, calcWaterDensity))
### Convert to raster
WaterDensityScore_raster <- rasterize(WaterDensityScore, quivira_DEM, field = "DensityScore", fun = max)
#######################################################################################
## Distance to Food ---------
### define distance to food function
calcFoodScore <- function(Distance) {
if (Distance > 1500) {
return(1)
} else if (Distance > 1000 && Distance <= 1500) {
return(2)
} else if (Distance > 500 && Distance <= 1000) {
return(3)
} else if (Distance > 0 && Distance <= 500) {
return(4)
} else if (Distance == 0) {
return(5)
} else {
return(9999)
}
}
#Set maximum distance to look for food from each wetland
maxDist <- 1500
dist_to_food <- st_buffer(screenedWetlands, dist = maxDist, dissolve = TRUE)
dist_to_food <- rbind(dist_to_food %>% dplyr::select(geometry))
dist_to_food <- st_union(dist_to_food)
plot(quivira_LULC)
plot(dist_to_food, col = NA, add = T)
# Extract LULC within buffered and screened wetlands
Extracted_LULC2019 <- mask(quivira_LULC, st_as_sf(dist_to_food))
unique(Extracted_LULC2019) #see landuse types in the datalayer## [1] 1 2 4 5 6 24 26 27 28 29 31 36 37 44 53 58 61 111 121
## [20] 122 123 124 131 141 143 152 176 190 195 205 225 236 238good_landuse <- c(1,2,4,5,6,24,26,27,28,29,31,36,37,44,53,58,61,225,236)
cropland <- Extracted_LULC2019
cropland[cropland %in% good_landuse] <- 1
cropland[cropland != 1] <- NA
cropland <- na.omit(cropland)
cropland_poly <- rasterToPolygons(cropland)
cropland_poly <- st_as_sf(cropland_poly)
closest <- st_nearest_feature(screenedWetlands, cropland_poly)
nearest_cropland <- cropland_poly[closest, ]
### calculate food distance score
FoodDistScore <- screenedWetlands %>%
mutate(Distance = st_distance(screenedWetlands, nearest_cropland, by_element = TRUE),
Distance = str_remove(Distance, " \\[m\\]")) %>%
mutate(FoodScore = sapply(Distance, calcFoodScore))
# Convert food score to raster
FoodScore_raster <- rasterize(FoodDistScore, quivira_DEM, field = "FoodScore", fun = max)Habitat Quality Assessment - Composite Habitat Quality Scores
Unweighted Composite habitat quality scores were calculated for each wetland by summing the scores of the individual habitat characteristics, adding no weight to the individual characteristics. This provided a comprehensive measure of habitat quality that takes into account different factors.
CompositeHQS <- WaterRegimeScore_raster + WetlandSizeScore_raster + WetlandTypeScore_raster + WaterDensityScore_raster + FoodScore_raster
glimpse(CompositeHQS)## Formal class 'RasterLayer' [package "raster"] with 13 slots
## ..@ file :Formal class '.RasterFile' [package "raster"] with 13 slots
## ..@ data :Formal class '.SingleLayerData' [package "raster"] with 13 slots
## ..@ legend :Formal class '.RasterLegend' [package "raster"] with 5 slots
## ..@ title : chr(0)
## ..@ extent :Formal class 'Extent' [package "raster"] with 4 slots
## ..@ rotated : logi FALSE
## ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
## ..@ ncols : int 3892
## ..@ nrows : int 3700
## ..@ crs :Formal class 'CRS' [package "sp"] with 1 slot
## ..@ srs : chr "PROJCRS[\"NAD_1983_Albers\",\n BASEGEOGCRS[\"NAD83\",\n DATUM[\"North American Datum 1983\",\n "| __truncated__
## ..@ history : list()
## ..@ z : list()Statistical Analysis
Descriptive statistics, such as min, max, mean, and median, were computed to summarize the distribution of habitat suitability scores and other relevant variables.
To validate the composite map, GPS telemetry points were used to see how many occurrence of whooping cranes will coincide with each habitat score. Also, a t-test is used to test if the distribution of habitat quality scores at crane locations is different from the overall distribution of scores in the raster. Finally, a correlation analysis was done to see if there is a correlation between the time of observation and habitat quality score
habitat_suitability <- CompositeHQS
telemetry_points <- WHCR_Telemetry
### Descriptive statistics
# Compute summary statistics for habitat suitability scores
summary_HSI <- summary(habitat_suitability)
hist_HSI <- hist(habitat_suitability)
########################################################################################################
### Statistical Analysis
########################################################################################################
crane_telemetry_sf <- st_as_sf(telemetry_points, coords = c("longitude", "latitude"), crs = 4326)
# Transform the CRS to match the habitat quality raster
crane_telemetry_sf <- st_transform(crane_telemetry_sf, crs(CompositeHQS))
# Extract the habitat quality scores at the crane locations
crane_telemetry_sf$habitat_quality <- raster::extract(CompositeHQS, st_coordinates(crane_telemetry_sf))
# To see if the distribution of habitat quality scores at crane locations is different from the overall distribution of scores in the raster.
# Perform a t-test
t_test <- t.test(crane_telemetry_sf$habitat_quality, na.omit(values(CompositeHQS)))
# To see if there is a correlation between the time of observation and habitat quality score.
# Convert time to a numeric value
crane_telemetry_sf$time_numeric <- as.numeric(crane_telemetry_sf$TIME_OM)
# Perform a correlation test
cor <- cor.test(crane_telemetry_sf$time_numeric, crane_telemetry_sf$habitat_quality)Validation
To validate the final map, a linear regression model of the habitat score and frequency of crane occurence was carried out. The linear regression model validates the Habitat Suitability Index (HSI) map by quantifying the relationship between the observed telemetry points (i.e., the locations where Whooping Cranes were actually found) and the habitat suitability scores predicted by the HSI map. If the HSI map is a good predictor of Whooping Crane habitat selection, I expect a significant positive relationship between the HSI score and the frequency of telemetry points. That is, areas with higher HSI scores (indicating higher habitat suitability) should have a higher frequency of telemetry points (indicating more frequent crane presence).
By fitting a linear regression model and testing the significance of the HSI coefficient, I will statistically validate whether the HSI map accurately predicts Whooping Crane habitat selection. If the HSI coefficient is significantly different from zero, this suggests that the HSI score is a significant predictor of crane presence, validating the HSI map.
############################################################################################
# Linear Regression Model
############################################################################################
data <- read.csv("CompositeScoreScreened_ExportTableQuivira.csv", header=TRUE, sep=",")
head(data)## OID_ RASTERVALU
## 1 1 12
## 2 2 NA
## 3 3 12
## 4 4 NA
## 5 5 12
## 6 6 12# Create a dataframe of abundance and HSI scores and remove factor levels.
table <- as.data.frame(table(data$RASTERVALU))
colnames(table)[1] <- "HSI"
table[,] <- lapply(table, function(x) {as.numeric(as.character(x))})
# Create a histogram
hist(data$RASTERVALU, main = "HSI Frequency (Screened Habitat Scores)", xlab = "HSI Values", col="lightblue", breaks = 8)
# Create a scatterplot of Freq vs. HSI and add a linear regression line of best fit
plot(table$HSI, table$Freq, main="Abundance vs. HSI", xlab="HSI", ylab="Abundance", pch=19)
abline(lm(table$Freq~table$HSI), col = "red", lty = 4, lwd = 2, add = T)
# Apply the lm() function
model <- lm(table$Freq~table$HSI)PART 3 - RESULTS
Results
A spatial map of the composite habitat quality scores was created to visually represent the distribution of high-quality habitats within the refuge. A second map showing the telemetry point distribution on the composite map was also plotted.
# Set up a 1x2 plot layout
layout(matrix(c(1,2), 1, 2, byrow = TRUE))
# composite habitat quality scores
plot(CompositeHQS, main = "Composite Habitat Assessment Map", bty = "n") # Remove the main title here
plot(quivira_10km, col = NA, border = "yellow", lwd = 2, add = TRUE, bty = "n")
plot(Quivira_NWR, col = NA, alpha = 0.8, add = T, bty = "n")
# composite habitat quality scores showing Telemetry points
plot(CompositeHQS, main = "CM showing Telemetry points", bty = "n") # Remove the main title here
plot(quivira_10km, col = NA, border = "yellow", lwd = 2, add = TRUE, bty = "n")
plot(Quivira_NWR, col = NA, alpha = 0.8, add = T, bty = "n")
plot(quivira_telemetry, pch = "+", col = "blue", add = T)
# Add a main title to the entire plot
mtext("Quivira NWR Composite Habitat Assessment Map", outer = TRUE, line = -1, cex = 1.5)
# Reset plot layout
layout(1)Habitat scores ranged from 6 to 18, with an average score of 12.Based on these scores, we can conclude that habitats within Quivira NWR scoring 12 or higher are considered potentially suitable habitats. Within the study area, the total area of potentially suitable habitats was 22.56 hectares.
Descriptive Statistics
The descriptive statistics of habitat suitability scores within the Quivira National Wildlife Refuge (NWR) revealed a mean suitability value of 13, indicating that all habitat score from 13 and above i.e. 12 - 18 are considered suitable habitat across the study area. Histogram shows that the most common wetland score within the study area is 12.
### Descriptive statistics
# Compute summary statistics for habitat suitability scores
print(summary_HSI)## layer
## Min. 6
## 1st Qu. 10
## Median 13
## 3rd Qu. 13
## Max. 19
## NA's 14256120plot(hist_HSI)
Statistical Analysis
The statistical analysis revealed significant relationships between the location of Whooping Crane telemetry points and the habitat quality scores. This suggests that Whooping Cranes are selecting habitats based on certain characteristics, underscoring the importance of these factors in managing and conserving these habitats.
#Result of t-test
print(t_test)##
## Welch Two Sample t-test
##
## data: crane_telemetry_sf$habitat_quality and na.omit(values(CompositeHQS))
## t = 3.1125, df = 21.006, p-value = 0.005268
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.7681389 3.8610795
## sample estimates:
## mean of x mean of y
## 15.45455 13.13994#correlation result
print(cor)##
## Pearson's product-moment correlation
##
## data: crane_telemetry_sf$time_numeric and crane_telemetry_sf$habitat_quality
## t = -3.4714, df = 20, p-value = 0.00241
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.8222227 -0.2583569
## sample estimates:
## cor
## -0.613175Validation - Linear Regression Model
The validation resulted in a p-value for the HSI coefficient as 0.197.
# Print summary of Linear Regression model results
print(summary(model))##
## Call:
## lm(formula = table$Freq ~ table$HSI)
##
## Residuals:
## 1 2 3 4 5 6 7
## -10.4821 5.6741 3.8304 1.1429 4.2991 -0.3884 -4.0759
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 25.7321 9.6180 2.675 0.0441 *
## table$HSI -1.1562 0.7772 -1.488 0.1970
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.218 on 5 degrees of freedom
## Multiple R-squared: 0.3068, Adjusted R-squared: 0.1682
## F-statistic: 2.213 on 1 and 5 DF, p-value: 0.197PART 4 - DISCUSSION AND CONCLUSION
Discussion
Comprehensive assessment of habitat suitability for Whooping Cranes within the Quivira National Wildlife Refuge (NWR) gives valuable insights into the factors influencing their habitat selection. The findings highlighted in this research shows the importance of various habitat characteristics in shaping the quality of habitats for these endangered birds.
The spatial analysis revealed a range of habitat quality scores across the refuge, with scores ranging from 6 to 18. This indicates significant variability in the suitability of habitats for Whooping Cranes. Habitats scoring 12 or higher are considered potentially suitable, covering an area of approximately 22.56 hectares within the study area.
The statistical analysis further supported the significance of habitat quality scores in influencing the distribution of Whooping Crane telemetry points. The t-test results indicated a significant difference in habitat quality scores between areas where telemetry points were recorded and the overall distribution of scores in the study area. Additionally, the correlation analysis revealed a negative correlation between the time of telemetry observations and habitat quality scores, suggesting that cranes tend to occupy habitats with higher quality scores at certain times.
The validation of the habitat suitability map through a linear regression model provided more information on how the Habitat Suitability Index (HSI) performed against real telemetry points. The validation resulted in a p-value for the HSI coefficient as 0.197, which is greater than 0.05. This suggests that HSI is not a significant predictor of frequency of bird occurrence at the 0.05 level. Although the p-value for the HSI coefficient was not significant at the conventional threshold of 0.05 (p = 0.197), the positive relationship between HSI scores and the frequency of telemetry points suggests that the HSI map can serve as a useful tool for predicting Whooping Crane habitat selection within the refuge.
Conclusion
In conclusion, this study provides a comprehensive assessment of habitat suitability for Whooping Cranes within the Quivira National Wildlife Refuge (NWR) in Kansas, United States. By integrating spatial analysis techniques with ecological knowledge, the research identifies key habitat characteristics that influence the distribution and quality of habitats for these endangered birds.
The findings underscore the importance of conserving and managing wetland habitats within the refuge, as they serve as critical stopover sites for Whooping Cranes during their long migratory journeys. Preservation and restoration efforts should focus on maintaining suitable water regimes, protecting large and natural wetlands, and minimizing human disturbances such as roads, railroads, and power lines in proximity to key crane habitats.
FINAL THOUGHTS
The findings of this study shed light on the intricate relationship between habitat suitability, spatial patterns, and conservation management strategies for Whooping Cranes within the Quivira National Wildlife Refuge (NWR). By integrating spatial analysis techniques and ecological principles, this research has provided valuable insights into the factors influencing habitat quality and distribution patterns of Whooping Cranes within the study area.
One of the key takeaways from this study is the importance of wetland habitats as critical refuges for Whooping Cranes, particularly during migration and wintering periods. The significant spatial autocorrelation and clustering of high habitat suitability areas highlight the presence of distinct spatial patterns in habitat quality, emphasizing the need for targeted conservation efforts to preserve and enhance these vital ecosystems.
Moreover, the identification of hotspots and coldspots of habitat suitability offers actionable information for conservation practitioners and policymakers to prioritize conservation interventions and allocate resources effectively. By focusing on areas of high habitat suitability, conservation initiatives can maximize their impact on improving habitat quality and connectivity for Whooping Cranes, thereby contributing to the long-term viability of this endangered species.
However, it is essential to acknowledge the limitations of this study, including potential data gaps, model uncertainties, and the dynamic nature of ecological systems. Future research endeavors should aim to address these limitations by incorporating additional data sources, refining habitat suitability models, and adopting a dynamic, adaptive management approach to conservation.
References
Armbruster, M.J. 1990. Characterization of Habitat Used by Whopping Cranes During Migration. U.S. Fish and Wildlife Service Biological Report. 90(4). 16 pages.
Adem Esmail, Blal, and Davide Geneletti. “Multi‐criteria Decision Analysis for Nature Conservation: A Review of 20 Years of Applications.” Methods in Ecology and Evolution, edited by Lynn Dicks, vol. 9, no. 1, Jan. 2018, pp. 42–53. https://doi.org/10.1111/2041-210X.12899.
BirdLife International. (2020). Grus americana. The IUCN Red List of Threatened Species 2020: e.T22692167A181363916. https://dx.doi.org/10.2305/IUCN.UK.2020-3.RLTS.T22692167A181363916.en.
Belton, L. (2019). Multicriteria decision analysis in ecology and conservation: Concepts, methods, and applications. Ecological Applications, 29(1), e01883. doi:10.1002/eap.1883.
Canadian Wildlife Service and U.S. Fish and Wildlife Service. (2007). International recovery plan for the Whooping Crane. Ottawa: Recovery of Nationally Endangered Wildlife (RENEW), and U.S. Fish and Wildlife Service, Albuquerque, New Mexico. 162 pp.
Cowardin, L.M., V. Carter, F.C. Golet, and E.T. LaRoe. 1979. Classification of wetlands and deepwater habitats of the United States. U.S. Department of the Interior, Fish and Wildlife Service, Washington, D.C. https://www.fws.gov/wetlands/documents/classification-of-wetlands-and-deepwater-habitats-of-the-united-states.pdf
Lewis, T. L., Jr. (1997). Whooping Crane (Grus americana). In The Birds of North America Online. https://birdsna.org/Species-Account/bna/species/whocrane/introduction
Pearse, Aaron T., et al. “Migrating Whooping Cranes Avoid Wind‐energy Infrastructure When Selecting Stopover Habitat.” Ecological Applications, vol. 31, no. 5, July 2021, p. e02324., https://doi.org/10.1002/eap.2324.
U.S. Fish and Wildlife Service. 2009. Whooping Cranes and Wind Development – An Issue Paper. U.S. Fish and Wildlife Service, Regions 2 and 6. 27 pages.
Urbanek, R. P., Jodice, P. G. R., & Kruse, C. D. (2018). Migratory Connectivity of Whooping Cranes. Waterbirds: The International Journal of Waterbird Biology, 41(3), 325–336. https://doi.org/10.1675/063.041.0310
Watershed Institute, Inc. 2013. Potentially suitable habitat assessment for the whooping crane (Grus americana). June 2013. Watershed Institute, Inc., 1200 SW Executive Drive, Topeka, Kansas 66615.