Land Surface Temperature

A short description of the post.

Show code

library(bfast)       # change detection algorithms
library(lubridate)   # handlingdate object
library(osmdata)     # spatial data extraction
library(raster)      # raster data
library(rgee)
library(sf)
library(tidyverse)
library(tmap)

#knitr::write_bib(x = c(.packages(), "distill"), "packages.bib")

• Theoretical and Applied Climatology (3.179)

Abstract

Introduction

The factors affecting LST are complex. latitude, altitude, the NDVI (normalized differential vegetation index), SM (soil moisture), AOD (aerosol optical depth), cloud fraction and atmospheric WV (water vapor), to analyze the driving factors of spatial differences in LST.

climate change India https://doi.org/10.1007/978-981-15-4327-2

• https://public.wmo.int/en/our-mandate/climate/wmo-statement-state-of-global-climate
• https://doi.org/10.1016/j.rse.2012.12.008

Methods

Study Site

• location
• elevation
• climate
• biodiversity

Khol Hi Raitan Wildlife Sanctuary

Show code

## extract boundaries and data from OpenStreetMap
library(osmdata)

## define a bounding box for area of interest
osm_bbox <- c(76.8, 30.6, 77.1, 30.9) # xmin, ymin, xmax, ymax

## extract boundary of the protected area as sf from OpenStreetMap (osm)
protected_areas <- opq(osm_bbox) |>
  add_osm_feature(key = "boundary", value = "protected_area", value_exact = F) |>
  osmdata_sf()

## rename protected areas to short names and save to local file as gpkg
protected_areas$osm_polygons |> select(name, geometry) |>
  mutate(short_name = c("MH", "BS", "SL")) |>
  st_write("data/protected_areas.gpkg")

## extract rivers / streams (waterways) from osm as sf
waterways <- opq(osm_bbox) |> add_osm_feature(key = "waterway") |> 
  osmdata_sf()

## save waterways to local file as gpkg
waterways$osm_lines |> select(name, waterway, geometry) |> 
  st_write("data/waterways.gpkg")

## extract roadways (highway) major roads from osm as sf
highways <- opq(osm_bbox) |> add_osm_feature(key = "highway") |> 
  osmdata_sf()

## save roadways to local file as gpkg
highways$osm_lines |> select(name, highway, geometry) |> 
  st_write("data/highways.gpkg")

## load, merge, and crop pre-downloaded SRTM 30 m elevation data 
elev <- raster("D:/spatial-data/dem/n30_e076_1arc_v3.tif") |>
  merge(raster("D:/spatial-data/dem/n30_e077_1arc_v3.tif")) |> 
  crop(extent(76.85, 77.01, 30.66, 30.75)) # morni hills wls

## prepare and save boundaries for elevation bands
masked_elev <- elev |> mask(mwls)
reclass_matrix <- matrix(data = c(300, 400, 1, 
                                  400, 500, 2, 
                                  500, 600, 3, 
                                  600, 700, 4, 
                                  700, 800, 5), 
                         ncol = 3, byrow = T)
masked_elev |> reclassify(rcl = reclass_matrix) |>
  rasterToPolygons(dissolve = T) |> st_as_sf() |> 
  mutate(elev_zone = c("301-400", "401-500", "501-600", "601-700", "701-800")) |>
  st_write("data/elev_band.gpkg")

## prepare and save hillshade
hillShade(
  # extract slope from elevation
  slope = terrain(elev, opt= "slope"),
  
  # extract aspect from elevation
  aspect = terrain(elev, opt= "aspect"), 
  
  # 45 altitude angle and 315 azimuth angle
  angle = 45, direction = 315) |> 
    writeRaster("data/hill_shade.tif")

Show code

## add hillshade
tm_shape(hill_shade) + 
    tm_raster(palette = "-Greys", n = 100, style = "cont", legend.show = FALSE) +
    
    tm_graticules(lwd = 1.5, col = "grey75", labels.size = 0.85) +  
    
  ## add elevation bands
    tm_shape(elev_band) + tm_fill(col = "elev_zone", alpha = 0.7, palette = "Greens",
                                 title = "Elevation (m)") + 
    
    ## add boundary of protected area
  tm_shape(mwls) + tm_borders(col = "seagreen", lwd = 2.5) +
    
    ## add rivers and streams
  tm_shape(rivers) + tm_lines(col = "#9bbff4", lwd = 6) +
  tm_shape(streams) + tm_lines(col = "#9bbff4", lwd = 3) +
  
  ## add major roads
  tm_shape(tertiary_roads) + tm_lines(col = "white", lwd = 1.5) +
  tm_shape(secondary_roads) + tm_lines(col = "white", lwd = 3) +
  tm_shape(trunk_roads) + tm_lines(col = "white", lwd = 4) +
  
    ## add legend
  tm_add_legend(type = "line", labels = "Roads", col = "white", 
                lwd = 4) +
    tm_add_legend(type = "line", labels = "Rivers", col = "skyblue", lwd = 4) +
    
    
    ## add compass and scale bar
  tm_compass(position = c(0.89, 0.84), size = 2.5, bg.color = "grey90", 
             bg.alpha = 0.7) +
  tm_scale_bar(position = c(0.26, 0.001), breaks = c(0, 5, 10), text.size = 0.75) +
    
    
    ## adjust layout
  tm_layout(frame.double.line = T, compass.type = "arrow", 
            legend.position = c(0.015, 0.015), legend.bg.color = "grey90", 
            legend.bg.alpha = 0.7, 
            main.title = "Khol Hi Raitan (Morni Hills) Wildlife Sanctuary",
            main.title.position = "center",
            legend.text.size = 0.85) 

Figure 1: Location and topography of Khol Hi Raitan Wildlife Sanctuary. The map is prepared in R programming environment using the package tmap (Tennekes 2018)

Show code

## remove all variables from the environment
rm(list = ls())

Satellite data

Land surface temperature (LST) data derived from the Moderate Resolution Imaging Spectro-radiometer (MODIS) instrument mounted on Terra satellite. MODIS has 36 spectral bands. The terra satellite was launched in 1999 and it has transit time 10:30 am. MODIS inversely generates LST by using information from middle infrared and thermal infrared bands.

The generalized split window method and day/night method are the official inversion algorithms of MODIS LST data at present. Through many verification analyses and improvements, the average precision of LST inversion for these two algorithms is approximately 1 K. MODIS LST data have been widely used in many fields, such as climate change, water cycles, environmental assessment and agricultural production, due to the good temporal and spatial resolution and wide coverage. The MOD11C3 LST products are inverted by the day/night method and obtained by projection, splicing, re-sampling and average synthesis. https://doi.org/10.1038/s41598-020-63701-5.

The MOD11A1 Version 6 product provides daily per-pixel Land Surface Temperature and Emissivity (LST&E) with 1 kilometer (km) spatial resolution in a 1,200 by 1,200 km grid. The pixel temperature value is derived from the MOD11_L2 swath product. Above 30 degrees latitude, some pixels may have multiple observations where the criteria for clear-sky are met. When this occurs, the pixel value is a result of the average of all qualifying observations. Provided along with the daytime and nighttime surface temperature bands are associated quality control assessments, observation times, view zenith angles, and clear-sky coverages along with bands 31 and 32 emissivities from land cover types.

Data product: MOD11A1.006 Terra Land Surface Temperature and Emissivity Daily Global 1km (Wan, Hook, and Hulley 2015)

This dataset represents an 8-day composite of 1-2 day observation intervals, with spatial resolution of 1000 m. Data are available from March 5, 2000 to present.

• summary of images used

Show code

## initialise the earth engine 
ee_Initialize()

## read the are of interest and convert to ee object
mwls <- st_read("data/elev_band.gpkg") |> sf_as_ee()

## load the modis lst data
modis_lst <- ee$ImageCollection("MODIS/006/MOD11A2")$
    
    ## filter dates
    filterDate(ee$DateRange("2001-01-01", "2021-01-01"))$
    
    ## select the LST band
    select("LST_Day_1km")

## define a function to convert Kelvin to Celcius and get dates for each image
modC <- modis_lst$map(
  function(image){
    return(
        image$multiply(0.02)$
        subtract(273.15)$
        select("LST_Day_1km")$
        copyProperties(image, image$propertyNames())
        )
  }
)

## Extract pixel values for each elevation band from each image
roi_lst <- ee_extract(modC, mwls)

## Reshape the data to longer format and convert dates strings to date format
roi_long <- roi_lst |> 
    pivot_longer(cols = -c(elev_zone, layer), names_to = "date", 
                 values_to = "lst") |>
    ## select characters 2 to 11 from the date
    mutate(date = str_sub(date, 2, 11))

## save to local file for further analysis    
write.csv(roi_long, "data/mwls_modis_lst.csv")

• terrain extraction: We extracted the terrain factors, including the elevation, slope, and aspect and a shaded relief map, from the DEM data.We extracted the terrain factors, including the elevation, slope, and aspect and a shaded relief map, from the DEM data.
• LST inversion. Te radiative transfer model.
• surface emissivity
• Interference removal.
• Accuracy verification of the LST inversion results

Processing

Statistical analysis

• correlation
• multiple regression
• Sen’s slope
• Theil-Sen Approach (TSA)
• Mann–Kendall
• Time series forecasting
• ARIMA, Support vector machines (SVMs)

Results

Show code

## read the previously extracted LST data
lst <- read.csv("data/mwls_modis_lst.csv") |>
    
    ## convert to date column to Date object
    mutate(date = lubridate::as_date(date, format = "%Y_%m_%d")) |>
    
    ## separate year and months from the date
    mutate(year = lubridate::year(date), 
           month = lubridate::month(date, label = T)) |>
    
    ## select all columsn except layer
    select(-layer)

lst |> head() |> gt::gt()

X	elev_zone	date	lst	year	month
1	301-400	2001-01-01	15.21948	2001	Jan
2	301-400	2001-01-09	12.74996	2001	Jan
3	301-400	2001-01-17	19.48637	2001	Jan
4	301-400	2001-01-25	22.93128	2001	Jan
5	301-400	2001-02-02	22.66557	2001	Feb
6	301-400	2001-02-10	24.59091	2001	Feb

Show code

## Overall simple regression
ggplot(data = lst, aes(x = date, y = lst)) +
    geom_point(size = 3.5, alpha = 0.5, color = "grey90") +
    geom_line(color = "grey85") +
    geom_smooth(alpha = 0.2) +
    ggpubr::stat_cor(show.legend = F, label.y = 52, p.accuracy = 0.001, 
                     r.accuracy = 0.01, color = "black") +
    ggpubr::stat_regline_equation(label.y = 48) +
    ylim(0, 55) +
    theme_light()    

Figure 2: Overall simple linear regression for LST

Show code

## simple linear regression
ggplot(data = lst, aes(x = date, y = lst)) +
    geom_point(size = 3.5, alpha = 0.5, color = "grey90") +
    geom_line(color = "grey85") +
    geom_smooth(alpha = 0.2) +
    facet_wrap(.~month) +
    ggpubr::stat_cor(show.legend = F, label.y = 52, p.accuracy = 0.001, 
                     r.accuracy = 0.01, color = "black") +
    ggpubr::stat_regline_equation(label.y = 5) +
    ylim(0, 55) +
    theme_light()

Figure 3: Monthly trend for LST

Show code

lst |> drop_na() |> group_by(year) |> 
  summarise(avg_lst = mean(lst), max_lst = max(lst), min_lst = min(lst)) |>
  pivot_longer(cols = -year, names_to = "vars", values_to = "lst") |>
  ggplot(aes(x = year, y = lst, color = vars)) +
  geom_point(size = 3.5) +
  geom_smooth(aes(fill = vars), method = "lm", alpha = 0.1) +
  ggpubr::stat_cor(show.legend = F, p.accuracy = 0.001, label.y = c(35, 55, 8),
                     r.accuracy = 0.01) +
  ggpubr::stat_regline_equation(show.legend = F, label.y = c(32, 52, 5)) +
  ylim(5, 55) +
  theme_bw()

Figure 4: Annual trends of LST

Show code

par(mfrow=c(3,1))
boxplot(lst ~ year, data = lst, cex = 1.2)
boxplot(lst ~ month, data = lst)
boxplot(lst ~ elev_zone, data = lst)

Figure 5: Variation in LST for year, month and elevation zones

Time series analysis

Show code

## calculate monthly average LST for each year
monthly_avg_lst <- lst |> drop_na() |> 
    group_by(year, month) |> 
    summarise(avg_lst = mean(lst))

## create a monthly ndvi time series object, frequency = 12 for 12 months
lst_ts <- ts(monthly_avg_lst$avg_lst, start = c(2001, 1), end = c(2020, 12), 
             frequency = 12)

par(mfrow = c(2, 1))
# classical seasonal decomposition by moving averages
lst_ts |> decompose() |> plot()

Figure 6: Time series decomposition for LST

Show code

# seasonal decomposition of time series by loess
lst_ts |> stl(s.window = "periodic") |> plot()

Figure 7: Time series decomposition for LST

Show code

library(bfast)

## BFAST algorithm for break
lst_ts |> bfast01() |> plot()

Figure 8: Change point detection using bfast

Show code

## classify the break
lst_ts |> bfast01() |> bfast01classify() |> gt::gt(caption = "Classification of breakpoint")

Table 1: Classification of breakpoint

flag_type	flag_significance	p_segment1	p_segment2	pct_segment1	pct_segment2	flag_pct_stable
5	1	0.01475315	0.1348009	0.4952384	0.2477474	NA

Show code

## detecting all breaks
lst_ts |> bfast(formula = response ~ trend + harmon, order = 3, breaks = "BIC") |> 
    plot()

Figure 9: Iterative breakpoint detection

#> NULL

Show code

lst_ts |> bfastmonitor(start = 2015, formula = response ~ trend + harmon) |> 
  plot(ylim = c(-7, 45))

Figure 10: Monitoring changepoint and LST time series

Show code

lst_ts |> bfastlite() |> plot()

Discussions

Conclusion

Authors’ contribution

Abhishek Kumar Conceptualisation, Methodology, Software, Validation, Formal Analysis, Investigation, Resources, Data Curation, Writing - Original Draft, Writing - Review & Editing, Visualization, Funding acquisition Meenu Patil Formal Analysis, Writing - Review & Editing, Funding acquisition Anand Narain Singh Supervision, Project administration