Urban Heat Island Effects in the FRAME

1) Load packages and set working directory in R.

#set working directory
setwd("/Users/zach/Dropbox (ZachTeam)/Projects/Heat_Island_Effect")

#load packages
library(raster)
library(ggplot2)
library(viridis)
library(rgdal)
library(data.table)

#set global options
knitr::opts_chunk$set(echo = TRUE)

2) Read in Landsat data: bands 4,5, 10, and 11 IMage acquired on 24 October 2019

band4<-raster("/Users/zach/Dropbox (ZachTeam)/Projects/Heat_Island_Effect/Data/LC08_L1TP_014032_20191024_20191030_01_T1/LC08_L1TP_014032_20191024_20191030_01_T1_B4.TIF")

band5<-raster("/Users/zach/Dropbox (ZachTeam)/Projects/Heat_Island_Effect/Data/LC08_L1TP_014032_20191024_20191030_01_T1/LC08_L1TP_014032_20191024_20191030_01_T1_B5.TIF")

band10<-raster("/Users/zach/Dropbox (ZachTeam)/Projects/Heat_Island_Effect/Data/LC08_L1TP_014032_20191024_20191030_01_T1/LC08_L1TP_014032_20191024_20191030_01_T1_B10.TIF")

band11<-raster("/Users/zach/Dropbox (ZachTeam)/Projects/Heat_Island_Effect/Data/LC08_L1TP_014032_20191024_20191030_01_T1/LC08_L1TP_014032_20191024_20191030_01_T1_B11.TIF")

3) Look at example raster

#plot raster
raster::plot(band11, col=viridis(255))

4) Convert Band 10 to Radians: open metadata file from Landsat folder (this is the .MTL file)

1. You are looking for the radiance rescaling multiplier for Band 10; copy that value
   1. NOTE: these are consistent values across Landsat 8 images: Band 10 Radiance Multiplier: 0.0003342 Radiance Add: 0.1 K1: 774.89 K2: 1321.0789
2. Raster Calculator-> 0.0003342 * “LC08_L1TP_014032_20170220_20170301_01_T1_B10.TIF” + 0.1 Grid = RAD10

RAD10<-(band10*0.0003342)+0.1

5) Convert Band 10 to thermal degrees: same metadata file, under TIRS thermal constants

1. Copy K values (as above; these are also constant)
2. Raster Calculator -> (1321.0789/Ln(774.8853/“RAD10” + 1))-273.15 (-273.15 to convert to Celsius) Grid = TB10

TB10<-(1321.0789/log(774.8853/RAD10 + 1))-273.15

6) Convert from At-Satellite Temp to Land Surface Temp

1. Add Band 4 and Band 5 to create NDVI grid
2. Raster Calculator-> Float(Band5-Band4)/Float(Band5+Band4) Grid = NDVI10

NDVI10<-(band5-band4)/(band5+band4)

7) plot NDVI

raster::plot(NDVI10, col=colorRampPalette(c("red", "yellow", "seagreen"))(100))

8) Compute NDVI max and min

NDVImin=min(values(NDVI10),na.rm=TRUE)
NDVImax=max(values(NDVI10),na.rm=TRUE)

9) NDVI grid value range High and Low are the Max and Min

Raster Calculator-> Square((“FebNDVI10”+1)/(1+1)) i.e., ((NDVI10 - NDVImin)/(NDVImax – NDVImin)) squared Grid = Pv

Pv<-((NDVI10-NDVImin)/(NDVImax-NDVImin))^2

10) Compute E10 raster calculation

Raster Calculator -> 0.004 * “FebPv10” + 0.986 Grid = E10

E10<-0.004*Pv+0.986

11) Compute TB10 raster calculation

Raster Calculator-> “FebTB10”/(1+ (10.8 * “FebTB10”/14388) * Ln(“FebE10”)) Grid = LST10 or “FebTB10”/(1+ (0.00115 * “FebTB10”/1.4388) * Ln(“FebE10”))

LST10<-TB10/(1+(10.8*TB10/14388))*log(E10)

12) Extract LST10 Grid by a Mask to trim extreme highs and lows on edges

Grid = LST10clip

#look at range of values
LST10range<-range(values(LST10),na.rm=TRUE) #-0.496812  1.955719


#mask raster with thresholding
LST10clip<-LST10
LST10clip[LST10clip > 1.55]<-NA
LST10clip[LST10clip < -0.5]<-NA

#plotLST10
raster::plot(LST10clip, col=viridis(100,begin=0,end=1))

13) Now repeat previous steps to convert Band 11 to Radians: open metadata file from Landsat folder (this is the .MTL file)

1. You are looking for the radiance rescaling multiplier for Band 11; copy that value
2. NOTE: these are consistent values across Landsat 8 images: Band 11 Radiance Multiplier: 0.0003342 Radiance Add: 0.1 K1: 480.89 K2: 1201.14
3. Raster Calculator-> 0.0003342 * “LC08_L1TP_014032_20170220_20170301_01_T1_B11.TIF” + 0.1 Grid = RAD11

RAD11<-(0.0003342 * band11) + 0.1

# 13.   Convert Band 11 to thermal degrees: same metadata file, under TIRS thermal constants
# a.    Copy K values (as above; these are also constant)
# b.    Raster Calculator (1201.14/Ln(480.89/"RAD11" + 1))-273.15
# (-273.15 to convert to Celsius)
# Grid = TB11
TB11<-(1201.14/log(480.89/RAD11 + 1))-273.15

# 14.   Convert from At-Satellite Temp to Land Surface Temp
# a.    Use NDVI, Pv, and E raster grids calculated above
# Raster Calculator "FebTB11"/(1+ (10.8 * "FebTB11"/14388) * Ln("FebE10"))
# Grid = LST11
# *or*
#   "FebTB11"/(1+ (0.00115 * "FebTB11"/1.4388) * Ln("FebE10"))
LST11<-TB11/(1+ (10.8 * TB11/14388) * log(E10))

#get range of values
LST11range<-range(values(LST11),na.rm=TRUE) #-131.3496   33.7534

# 15.   Extract LST11 Grid by a Mask to trim extreme highs and lows on edges
# Grid = LST11clip
LST11clip<-LST11
#-131.34962   36.84848
LST11clip[LST11clip < -125]<-NA
LST11clip[LST11clip > 30]<-NA

#plotLST10
raster::plot(LST11clip, col=viridis(100))

14) Get mean of raster cells

# 16.   Use Cell Stats to get Mean for LST10 and LST11
LST10mean<-mean(values(LST10clip), na.rm=TRUE) # -0.2090109
print(LST10mean)

## [1] -0.2090109

LST11mean<-mean(values(LST11clip), na.rm=TRUE) # 15.65445
print(LST11mean)

## [1] 15.65445

15) Read in FRAME PA and DE shapefile

FRAMEsites<-readOGR(dsn="/Users/zach/Dropbox (ZachTeam)/Projects/Heat_Island_Effect/Data/GIS_shapefiles/FRAME_SitesPointsRoadsFragments/",layer="FRAME_Sites_DE_PA")

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/zach/Dropbox (ZachTeam)/Projects/Heat_Island_Effect/Data/GIS_shapefiles/FRAME_SitesPointsRoadsFragments", layer: "FRAME_Sites_DE_PA"
## with 38 features
## It has 9 fields

#plot FRAME site polygons
plot(FRAMEsites, col="seagreen")

16) Check and standardize projections

#check projection
crs(FRAMEsites)

## CRS arguments:
##  +proj=utm +zone=18 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
## +no_defs

#set projection of FRAME sites to Landsat raster
crs(FRAMEsites)<- crs(LST11clip)

#look at new projection
crs(FRAMEsites)

## CRS arguments:
##  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0

17) Get Site Names from FRAME shapefile

#get FRAME site names from shapefile
SiteNames<-unique(FRAMEsites@data$Site_Name)
SiteNames.df<-as.data.frame(SiteNames)
SiteNames.df$RowNum<-row.names(SiteNames.df)

unique(SiteNames.df$SiteNames)

##  [1] Andorra             Christina Creek 1   Christina Creek 2  
##  [4] Coverdale           Chamounix           Crows Nest Interior
##  [7] Crows Nest North    Crows Nest West     Chrysler Woods     
## [10] Dorothy Miller      Ecology Woods       Folk               
## [13] Glasgow 1           Glasgow 2           Iron Hill 1        
## [16] Iron Hill 2         Laird               Longwood Square    
## [19] Mount Cuba Interior Mount Cuba Rose     Motor Pool         
## [22] Peacedale East      Peacedale West      Pennypack          
## [25] Phillips            Reservoir           Rittenhouse        
## [28] Stroud Power Line   Stroud South        Sunset Lake 1      
## [31] Sunset Lake 2       Smith Memorial      Walton's Run       
## [34] White Clay 1        White Clay 2        Webb Farm          
## [37] Park Line Drive     Longwood Triangle  
## 38 Levels: Andorra Chamounix Christina Creek 1 ... White Clay 2

18) Crop rasters to extent of FRAME polygons

#crop rasters LST10 and LST11 to extent of FRAME polygons
LST10_crop<- crop(LST10clip, extent(FRAMEsites))
plot(LST10_crop, col=viridis(100))
plot(FRAMEsites, col="red",add=TRUE)

LST11_crop<- crop(LST11clip, extent(FRAMEsites))
plot(LST11_crop, col=viridis(100))
plot(FRAMEsites, col="red",add=TRUE)

19) Extract values from rasters within FRAME polygons

#extract values from within each FRAME polygon
FRAME_LST10 <- extract(LST10clip, FRAMEsites)
FRAME_LST11 <- extract(LST11clip, FRAMEsites)

20) Create data.frame of Sites with extracted cell values

#create data.frame of extracted cell values
lst10.df = rbindlist(lapply(FRAME_LST10, function(x) data.frame(t(x))),
               fill = TRUE)
lst10.df$RowNum<-row.names(lst10.df)

#now merge with SiteNames.df
lst10.merge<-merge(lst10.df,SiteNames.df, by="RowNum",all.x=TRUE)
#unique(lst10.merge$SiteNames)

#add metric column
lst10.merge$Metric<-"LST10"

#now for LST11
lst11.df = rbindlist(lapply(FRAME_LST11, function(x) data.frame(t(x))),
                     fill = TRUE)
lst11.df$RowNum<-row.names(lst11.df)

#now merge with SiteNames.df
lst11.merge<-merge(lst11.df,SiteNames.df, by="RowNum",all.x=TRUE)
#unique(lst11.merge$SiteNames)

#add metric column
lst11.merge$Metric<-"LST11"

#now combine lst10.merge and lst11.merge
lstComb<-plyr::rbind.fill(lst10.merge, lst11.merge)

#now melt to get ready to plot
lst.melt<-melt(lstComb, id.vars=c("SiteNames","Metric","RowNum"))

## Warning in melt(lstComb, id.vars = c("SiteNames", "Metric", "RowNum")): The melt
## generic in data.table has been passed a data.frame and will attempt to redirect
## to the relevant reshape2 method; please note that reshape2 is deprecated, and
## this redirection is now deprecated as well. To continue using melt methods from
## reshape2 while both libraries are attached, e.g. melt.list, you can prepend the
## namespace like reshape2::melt(lstComb). In the next version, this warning will
## become an error.

#make Site and Metric factors
lst.melt$SiteNames<-as.factor(lst.melt$SiteNames)
lst.melt$Metric<-as.factor(lst.melt$Metric)

head(lst.melt,10)

##          SiteNames Metric RowNum variable      value
## 1          Andorra  LST10      1       X1 -0.1918848
## 2   Dorothy Miller  LST10     10       X1 -0.2157781
## 3    Ecology Woods  LST10     11       X1 -0.2130466
## 4             Folk  LST10     12       X1 -0.2055094
## 5        Glasgow 1  LST10     13       X1 -0.2037947
## 6        Glasgow 2  LST10     14       X1 -0.2084698
## 7      Iron Hill 1  LST10     15       X1 -0.1954404
## 8      Iron Hill 2  LST10     16       X1 -0.2003476
## 9            Laird  LST10     17       X1 -0.2040477
## 10 Longwood Square  LST10     18       X1 -0.1877210

21) Now plot LST10 and LST11 as boxplot among sites

#now plot in ggplot
heatIslandPlot<-ggplot(data=lst.melt, aes(x=SiteNames, y=value))+
  geom_boxplot(aes(fill=Metric),alpha=0.7)+
  scale_fill_manual(values=c("royalblue2","tomato"))+
  theme(axis.text.x=element_text(angle=90),
        )
#heatIslandPlot

heatIslandFacet<-heatIslandPlot+facet_wrap(.~Metric, scales="free_y",ncol=1)
print(heatIslandFacet)

## Warning: Removed 9630 rows containing non-finite values (stat_boxplot).

22) Now look at each FRAME site up close with LST10 and LST11 clipped rasters

#look at FRAME sites up close
siteList<-unique(lst.melt$SiteNames)

#LST10
for(i in 1:length(siteList)){
  
  newSite<-subset(FRAMEsites, Site_Name==siteList[i])
  
  siteCrop<- crop(LST10clip, extent(newSite))
  siteMask<- mask(siteCrop, newSite)
  plot(siteMask,col=viridis(100),main=paste("LST10 - ", as.character(siteList[i]),sep=""),
       xlab="utm X",ylab="utm Y")
  plot(newSite, border="red",add=TRUE)

}

#LST11
for(i in 1:length(siteList)){
  
  newSite<-subset(FRAMEsites, Site_Name==siteList[i])
  
  siteCrop<- crop(LST11clip, extent(newSite))
  siteMask<- mask(siteCrop, newSite)
  plot(siteMask,col=viridis(n=100,option="magma"),main=paste("LST11 - ", as.character(siteList[i]),sep=""),
       xlab="utm X",ylab="utm Y")
  plot(newSite, border="red",add=TRUE)
  
}