require(tidyverse);
require(ggplot2); 
require(raster); 
require(sf)

Load in Data

Loading in data from each year. Two from 1995 and two from 2020. A shapefile is also loaded for cropping and masking later in the project. The shapefile is transformed to contain the same coordinate reference systems.

raster_19950727 <- raster('C:/Users/rsaya/Downloads/CapstoneData/LT05_19950727/LT05_L2SP_013032_19950727_20200912_02_T1_ST_B6.TIF')

raster_19950828 <- raster('C:/Users/rsaya/Downloads/CapstoneData/LT05_19950828/LT05_L2SP_013032_19950828_20200912_02_T1_ST_B6.TIF')

raster_20200613 <- raster('C:/Users/rsaya/Downloads/CapstoneData/LC08_20200613/LC08_L2SP_013032_20200613_20200823_02_T1_ST_B10.TIF')

raster_20200820 <- raster('C:/Users/rsaya/Downloads/CapstoneData/LC08_L2SP_013032_20200512_20200820_02_T1_ST_B10.TIF')

nyc_shapefile <- st_read('C:/Users/rsaya/Downloads/CapstoneData/geo_export_706d9c1a-ad08-46de-9dd9-059f649c0e36.shp')
## Reading layer `geo_export_706d9c1a-ad08-46de-9dd9-059f649c0e36' from data source `C:\Users\rsaya\Downloads\CapstoneData\geo_export_706d9c1a-ad08-46de-9dd9-059f649c0e36.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 5 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.25559 ymin: 40.49613 xmax: -73.70001 ymax: 40.91553
## Geodetic CRS:  WGS 84
nyc<- st_transform(nyc_shapefile, crs = st_crs(raster_20200820))

Reprojecting Data

Due to each raster having different dimension to each other, re projection is needed. The raster for August 20, 2020 is used as the master data. The rest of the raster data are reprojected based on the master data dimensions and coordinates reference. The end goal is to be able to stack the rasters.

july95 <- projectRaster(raster_19950727, raster_20200820, method = 'ngb')
august95 <- projectRaster(raster_19950828, raster_20200820, method = 'ngb')
june20 <-  projectRaster(raster_20200613, raster_20200820, method = 'ngb')

Changing Extent

One last step before stacking, the extents need to match. This is done with the simple function, again using August 20, 2020 raster as the master data for reference.

extent(raster_19950727) <- extent(raster_20200820)
extent(raster_19950828) <- extent(raster_20200820)
extent(raster_20200613) <- extent(raster_20200820)

Stack Rasters

Stacking the raster for the respective years, resulting in a RasterStack for 1995 and 2020.

stack1995 <- stack(july95, august95)
stack2020 <- stack(june20, raster_20200820)

Calculations

The averages for both raster stacks are calculated.

average_1995 <- calc(stack1995, fun = mean)
average_2020 <- calc(stack2020, fun = mean)

Rescaling Raster

Since the Landsat data are both Collection 2 and Level 2 products, the only scaling needed is to multiply the scale factor, then adding 149 for the offset. These values are obtained from the USGS Landsat Collection 2 Level 2 Science Product Guide. To convert the values from Kelvin to Celsius, subtract by 273.15.

To compare the land surface temperature changes, the raster from 2020 is subtracted with the 1995.

fun <- function(x) { (x * 0.00341802 + 149) - 273.15 }

lst_1995 <- calc(average_1995, fun)
lst_2020 <- calc(average_2020, fun)

temp_diff <- lst_2020 - lst_1995

Cropping and Masking

To prepare for plotting, the raster data are cropped and masked to match the NYC shapefile extent. This lets the raster keep only the data inside the shapefile boundaries.

They are then converted into a data frame, allowing for visualization, and in this case, plotting the raster. The data frame allows for each value to be tabular, meaning each row and column will represent a variable. The data frame includes x and y columns for the spatial coordinates and a layer column for the temperature variable.

crop_1995 <- crop(lst_1995, nyc)
crop_2020 <- crop(lst_2020, nyc)
crop_diff <- crop(temp_diff, nyc)

mask_1995 <- mask(crop_1995, nyc)
mask_2020 <- mask(crop_2020, nyc)
mask_diff <- mask(crop_diff, nyc)

df_1995 <- as.data.frame(mask_1995, xy=TRUE)%>%drop_na()
df_2020 <- as.data.frame(mask_2020, xy=TRUE)%>%drop_na()
df_diff <- as.data.frame(mask_diff, xy=TRUE)%>%drop_na()

Plotting using ggplot2

ggplot()+
  geom_tile(aes(x=x, y=y, fill=layer), data = df_1995)+
  geom_sf(fill='transparent', data = nyc)+
  scale_fill_viridis_c(option = "C", name='C°', direction = 1)+
  labs(x='Longitude', y='Latitude',
       title='LST in 1995',
       caption='Source: USGS')+
   theme(panel.background = element_blank())

ggplot()+
  geom_tile(aes(x=x, y=y, fill=layer), data = df_2020)+
  geom_sf(fill='transparent', data = nyc)+
  scale_fill_viridis_c(option = "C", name='C°', direction = 1)+
  labs(x='Longitude', y='Latitude',
       title='LST in 2020',
       caption='Source: USGS')+
   theme(panel.background = element_blank())

ggplot()+
  geom_tile(aes(x=x, y=y, fill=layer), data = df_diff)+
  geom_sf(fill='transparent', data = nyc)+
  scale_fill_viridis_c(option = "C", name='Changes in Temperature (C°)', direction = 1)+
  labs(x='Longitude', y='Latitude',
       title='Changes in LST: 1995 vs. 2020',
       caption='Source: USGS')+
   theme(panel.background = element_blank())