Grp_P_137
BHY
2023-11-28
Imported median night lights data for an area of the East Coast from all of 2010 using the following GEE code: https://code.earthengine.google.com/3f3fd8fbdada375c2bf2c67deb34f5c9
Imported aver red tailed hawk appearances over a year from Ebird: https://science.ebird.org/en/status-and-trends/species/rethaw/downloads?week=1
Prelims + load data
# set working directory
setwd("~/Desktop/137")
# load packages into library
library(raster)## Loading required package: splibrary(sf)## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUElibrary(viridis)## Loading required package: viridisLitelibrary(RColorBrewer)
# pull .tifs from working directory
nighttime_lights <- raster("median_nighttime_lights.tif")
# print raster and view raster summary to better understand data
print(nighttime_lights) ## class : RasterLayer
## dimensions : 1203, 1587, 1909161 (nrow, ncol, ncell)
## resolution : 0.008983153, 0.008983153 (x, y)
## extent : -81.05499, -66.79872, 35.33074, 46.13747 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : median_nighttime_lights.tif
## names : avg_vissummary(nighttime_lights)## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (5.24% of all cells)## avg_vis
## Min. 0.000000e+00
## 1st Qu. 0.000000e+00
## Median 3.949176e+00
## 3rd Qu. 8.574970e+00
## Max. 1.141667e+03
## NA's 4.256280e+05# pull .tifs from working directory
Hawk <- raster("rethaw_abundance_seasonal_full-year_mean_2022.tif")
# print raster and view raster summary to better understand data
print(Hawk)## class : RasterLayer
## dimensions : 6081, 13314, 80962434 (nrow, ncol, ncell)
## resolution : 2541.849, 2541.868 (x, y)
## extent : -16921198, 16920977, -6996497, 8460601 (xmin, xmax, ymin, ymax)
## crs : +proj=eck4 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
## source : rethaw_abundance_seasonal_full-year_mean_2022.tif
## names : full_yearsummary(Hawk)## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (0.12% of all cells)## full_year
## Min. 0.000000e+00
## 1st Qu. 0.000000e+00
## Median 0.000000e+00
## 3rd Qu. 1.305150e-01
## Max. 1.737937e+00
## NA's 7.487973e+07# Download the shapefile for the state of New York from downloadable-in-R data
ny_shape <- getData("GADM", country = "USA", level = 1)## Warning in getData("GADM", country = "USA", level = 1): getData will be removed in a future version of raster
## . Please use the geodata package insteadny_shape <- subset(ny_shape, NAME_1 == "New York")Match Coordinate Reference Systems
## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs
## WKT2 2019 representation:
## GEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]],
## CS[ellipsoidal,2],
## AXIS["longitude",east,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## AXIS["latitude",north,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]]]## [1] "+proj=longlat +datum=WGS84 +no_defs"## Coordinate Reference System:
## Deprecated Proj.4 representation:
## +proj=eck4 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
## WKT2 2019 representation:
## PROJCRS["unknown",
## BASEGEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]]],
## CONVERSION["unknown",
## METHOD["Eckert IV"],
## PARAMETER["Longitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["False easting",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]Match Resolutions
# get resolutions
res(Hawk) # larger resolutions## [1] 2541.849 2541.868res(nighttime_lights_reprojected) # smaller resolution## [1] 759 998# change smaller res raster to have the larger/more course resolution
NL_LowRes <- resample(nighttime_lights_reprojected, Hawk, method = "bilinear")
# check resolutions to make sure they match
res(Hawk) # larger resolutions## [1] 2541.849 2541.868res(NL_LowRes) # smaller resolution## [1] 2541.849 2541.868Clip and match extents
# Get the extent of the New York vector shape
extent_ny <- extent(ny_shape_reproj)
# Crop and mask rasters to the extent and shape of New York
Hawk_clip <- mask(crop(Hawk, extent_ny), ny_shape_reproj)
NL_clip <- mask(crop(NL_LowRes, extent_ny), ny_shape_reproj)Plot resulting rasters
# Set up a 1x2 layout
par(mfrow = c(1, 2))
data_range_Hawk <- range(values(Hawk_clip), na.rm = TRUE)
# Set up breaks with a slight focus on the lower end of the data
breaks_Hawk <- seq(data_range_Hawk[1], data_range_Hawk[1] + 0.7 * (data_range_Hawk[2] - data_range_Hawk[1]), length.out = 256)
# Plot 2: Nighttime Lights Raster with Vector Overlay using viridis color palette
plot(
Hawk_clip,
main = "Hawk Abundance",
col = viridis(256), # Adjust alpha for transparency
breaks = breaks_Hawk,
asp = 1
)
# Plot 2: Raster plot with vector overlay
# data range
data_range_NL <- range(values(NL_clip), na.rm = TRUE)
# Set up breaks with a slight focus on the lower end of the data
breaks_NL <- seq(data_range_NL[1], data_range_NL[1] + 0.7 * (data_range_NL[2] - data_range_NL[1]), length.out = 256)
# Plot 2: Nighttime Lights Raster with Vector Overlay using viridis color palette
plot(
NL_clip,
main = "Nighttime Lights Raster",
col = magma(256), # Adjust alpha for transparency
breaks = breaks_NL,
asp = 1
)
plot(ny_shape_reproj, add = TRUE)
Merge/join the two data by common long and lat values
# Create a spatial points dataframe by merging with common longitudes and latitudes
# Creating Lat and Long grid to map points onto by using extent
common_points <- expand.grid(
x = seq(from = extent(NL_clip)[1], to = extent(NL_clip)[2], by = res(NL_clip)[1]),
y = seq(from = extent(NL_clip)[3], to = extent(NL_clip)[4], by = res(NL_clip)[2])
)
# Extract values from Hawk_clip and NL_clip
values_hawk <- raster::extract(Hawk_clip, common_points)
values_nl <- raster::extract(NL_clip, common_points)
# Combine values into a dataframe recording Longitude, Latitude, Hawk Abundance, and Nighttime Lights
merged_df_HAWK_NL <- data.frame(
Longitude = common_points$x,
Latitude = common_points$y,
Hawk_Abundance = values_hawk,
Nighttime_Lights = values_nl
)
# view the start of the dataframe as a sanity check
head(merged_df_HAWK_NL)## Longitude Latitude Hawk_Abundance Nighttime_Lights
## 1 -6695340 5064666 NA NA
## 2 -6692798 5064666 NA NA
## 3 -6690256 5064666 NA NA
## 4 -6687715 5064666 NA NA
## 5 -6685173 5064666 NA NA
## 6 -6682631 5064666 NA NAPreliminary clean data of NA values
# remove all rows with NA values because there are many spaces where the rasters both don't cover (from water body removal and state boarder removal, to data flukes/missing points)
merged_df_HAWK_NL_noNA <- na.omit(merged_df_HAWK_NL)Export dataframe as csv
# ile path for the CSV file
csv_file <- "~/Desktop/137/HAWK_NL_merge.csv"
# Export the dataframe, merged_df_HAWK_NL to a CSV file
write.csv(merged_df_HAWK_NL_noNA, file = csv_file, row.names = FALSE)