In HunterGleason/R.ccdc.tools: Simple Tools for Processing Zhu & Woodcock 2014 CCDC Outputs

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Introduction

This vignette provides a brief tour of the functionality of th R.ccdc.tools package. An important note, this package was designed to work in conjunction with the output raster layers of the Google Earth Engine tools for CCDC API CCDC API. CCDC API leverages Google Earth Engine to apply the Continuous Change Detection and Classification methods of Zhu & Woodcock, 2014 CCDC. A Google Earth Engine script GEE-Script is provided to generate the CCDC outputs, which may then be downloaded locally (from Google Drive) and used with the package provided here, which offers functions for generating derived images and time series plots based on the CCDC output raster of GEE-Script. This package has been tested only with GEE-Script outputs.

Imports

library(R.ccdc.tools)
library(sf)
library(stars)
library(terra)
library(bcdata)
library(ggplot2)
library(RSAGA)

Load example CCDC output from package data, explore data

Below we load in from R.ccdc.tools a CCDC image (exported from GEE-Script) for a small area of interest, we use the R.ccdc.tools::name_ccdc_bands function to rename the bands, as they are lost during export. This CCDC raster contains the CCDC coefficients for each segment, in this case 8 segments were exported using the CCDC API. The Landsat blue, green, red, near-infrared, shortwave-infrared 1, shortwave-infrared 2 and thermal band coefficients are included as blue,green,red,nir,swir1,swir2 and therm, these are the expected band names. In the code below we define a random point within the AOI, and extract the CCDC raster values at that location to illustrate what attributes are present in the output image of GEE-Script. For instance, S1_green_coef_INTP is the green band intercept coefficient for the first segment in the surface reflectance time series at the XY location 1794834, 609232. Notice additional attributes (e.g., S1_tStart) related to change detection and the individual segments are also present in the raster data as additional bands.

ccdc_img = rast(exampleData)
ccdc_img<-name_ccdc_bands(ccdc_img,8)
ccdc_img<-st_as_stars(ccdc_img)


#Arbitrary example coordinates (EPSG:3005) within the image bounds
x_coord<-1794834
y_coord<-609232
epsg<-3005

#Extract the model coefficients at the point from the CCDC image
DF <- data.frame(x = c(x_coord), y = c(y_coord))
DF_sf = st_as_sf(DF, coords = c("x", "y"), 
                 crs = st_crs(epsg), agr = "constant")
extracted <- st_as_sf(st_extract(ccdc_img, DF_sf))

#View the coefficients for each segment, these are used in 'get_ccdc_ts' to generate the time series data
ccdc_coefs <- as.data.frame(extracted)
head(ccdc_coefs)

Plot synthetic surface reflectance at a pixel location

Below we use R.ccdc.tools::get_ccdc_ts to retrieve the CCDC coefficients to generate the synthetic surface reflectance time series at the pixel that intersects the coordinates above for the green and nir bands. Any XY coordinate that intersect the raster data can be specified, for any band exported in the CCDC raster. Notice that the time series has a breakpoint at ~2018 in both the green and nir bands indicating a landcover change at this location occurred then.

#Get the green surface reflectance time series (ts) for all 8-segments at the point.
ts_green<-get_ccdc_ts(ccdc_img,x_coord,y_coord,epsg,'green',8)

#Get the near-infrared (nir) surface reflectance time series (ts) for all 8-segments at the point.
ts_nir<-get_ccdc_ts(ccdc_img,x_coord,y_coord,epsg,'nir',8)

#Plot the time series 
plot(jdoy_to_date(ts_green[,'jdoy_vec']),ts_green[,'ts'],type='l',xlab='Date',ylab='Green Band Surface Reflectance')
plot(jdoy_to_date(ts_nir[,'jdoy_vec']),ts_nir[,'ts'],type='l',xlab='Date',ylab='NIR Band Surface Reflectance')

Generate a synthentic CCDC image for a given date

In the code below we apply the R.ccdc.tools::gen_synthetic_img for generating and plotting synthetic RGB Landsat scenes for several consecutive years on July 15 using the CCDC coefficients found in ccdc_img.Note that these images represent the predicted (i.e., synthetic) surface reflectance values and NOT the measured surface reflectance values, although they should agree closely. Landcover changes with time are apparent in the time series.

dates<-paste0(seq(1986,2020,6),"-07-15")

for(date in dates)
{
  cat(date,'\n')
  synthetic<-gen_ccdc_synthetic_img(ccdc_img,date)
  plotRGB(synthetic,r=3,g=2,b=1,stretch="hist")
  title(date,line=-1)
}

Generate CCDC image with coeffcients that correspond to specfic date only

It can be useful to generate raster with only the CCDC coefficients that correspond to specific point in time. This is accomplished with the R.ccdc.tools::gen_ccdc_at_jdoy function, which returns a raster with all CCDC coefficients but only for the segments that occur during the provided date (within the CCDC time-series, typically 1984~present day). For example, below we generate a CCDC image that contains all the CCDC coefficients but taken only from the segments that occur during July 15, 2021, and we see that the bands no longer have a segment identifier associated with them. We use July 15, 2021 because it represent near present landscape conditions when this document was written, but any date could be used. This output can be used for comparing seasonal surface reflectance curves, or input directly into any classification or regression model. Notice that the start Julian date and end Julian date of the corresponding segment is also included as additional bands. Note for large CCDC rasters this function may be slow and memory intensive, if this is case (not shown here) it is recommended to create tiles (i.e., stars::st_tile) from the CCDC imagery and process the tiles using R.ccdc.tools::gen_ccdc_at_jdoy, a VRT can then be made (i.e., gdalbuildvrt) for working with the output tiles as though they were one image. Just be sure to rename the bands after importing a VRT into R as the original band names will be lost, otherwise the following functions will fail as they look for bands by name and not index.

#Get a CCDC image with coefficients for segments occurring during July 15, 2021 (i.e., near present day landscape conditions)
jul_15_2021<-gen_ccdc_at_jdoy(ccdc_img,"2021-07-15")
plot(jul_15_2021)

Assemble a category layer

As a very simple example below we load in basic landcover polygons for use in later functions, these are simple landcover types created from visual inspection of high resolution imagery, wetland polygons are derived using the bcdata package, categories are indicated by the landcover field. Note any categorized polygon layer could be employed, such as BEC Zones.

#Load a example categorical layer (e.g., landcover polygons)
landcover <- landcover
landcover <- vect(landcover)


#Plot category layer
ggplot()+
  geom_sf(data=st_as_sf(landcover),aes(color=landcover))+
  coord_sf(datum=st_crs(epsg))

Comparing seasonal surface reflectance curves by landcover category

A useful aspects of the CCDC is that it can be used to visually compare surface reflectance signatures for different landcover types at a given point in time, using the above jul_15_2021 raster as input, the code bellow illustrate the application of several R.ccdc.tools functions that make creating such plots very simple. The function R.ccdc.tools::plot_ccdc_ts_bypoly selects random sample locations of specified N within each polygon category, and for each sample plots the synthetic surface reflectance from the date of the input image up to any number of days, typically one or two years into the future. The function R.ccdc.tools::plot_ccdc_ts_quartiles is very similar, but uses all available pixels and plots the daily 25, 50 and 75th percentile for each category instead of each individual time series, which is illustrated below for the green and nir Landsat bands.

#Plot CCDC 60 sample time series by polygon type (i.e., landcover) for the green band ('green') over a 730 day period 
plot_ccdc_ts_bypoly(jul_15_2021,"2021-07-15",landcover,'landcover',730,'green',60)

#Same plot as above, only showing quartiles
plot_ccdc_ts_quartiles(jul_15_2021,"2021-07-15",landcover,'landcover',730,'green')

#Plot CCDC 60 sample time series by polygon type (i.e., landcover) for the near-infrared band ('nir') over a 720 day period 
plot_ccdc_ts_bypoly(jul_15_2021,"2021-07-15",landcover,'landcover',730,'nir',60)

#Same plot as above, only showing quartiles
plot_ccdc_ts_quartiles(jul_15_2021,"2021-07-15",landcover,'landcover',730,'nir')

Comparing seasonal surface reflectance curves by landcover category (large areas of interest)

Although we use the same CCDC image here, often the size of the exported GEE-Script raster can be large and many polygons may be used for exploring surface reflectance signatures between categories. The R.ccdc.tools::plot_ccdc_by_catg function achieves essentially the same result as R.ccdc.tools::plot_ccdc_ts_quartiles, but uses a random sample like plot_ccdc_ts_bypoly by first calling R.ccdc.tools::sample_ccdc_by_catg, which uses SAGA-GIS with RSAGA for more efficient computation. It also returns the computed surface reflectance data back as a data.frame. Because SAGA-GIS cannot use R objects directly, they must first be written to disk. The two functions are applied below with 300 random samples locations (total, constrained to polygons) with a minimum distance of 50 m between random locations, the result is a list with the raw-data, plot and legend crosswalk, and are also shown. Plots for the green and nir bands are generated as an example.

#Write jul_15_2021 and landcover to temp directory 
writeRaster(jul_15_2021,file.path(tempdir(),"jul_15_2021.tif"))
writeVector(landcover,file.path(tempdir(),"landcover.shp"))

#Use RSAGA for sampling CCDC coefficients (jul_15_2021) by polygon category, useful for large areas of interest with lots of polygons.  
sampled<-sample_ccdc_by_catg(file.path(tempdir(),"landcover.shp"),file.path(tempdir(),"jul_15_2021.tif"),300,50,rsaga.env())
head(sampled)

#Use he result from above and create a plot from them for green band 
sampled_ts<-plot_ccdc_by_catg(sampled,"2021-07-15","landcover","green")
sampled_ts

#Use he result from above and create a plot from them for near-infrared band 
sampled_ts<-plot_ccdc_by_catg(sampled,"2021-07-15","landcover","nir")
sampled_ts

#Legen crosswalk
sampled_ts['crosswalk']

#Synthetic surface reflectance data from Jul 15, 2021 -- 730 days forward. 
srdata<-as.data.frame(sampled_ts['srdata'])
head(srdata)

Summarise (min, median, max) NDMI and NDWI indicies from CCDC image over specified time period (Currently only NDMI & NDWI).

sumr_2020_ndmi_stats<-ccdc_ndmi_stats(ccdc_img,"2020-06-1","2020-08-30")
plot(sumr_2020_ndmi_stats)

Conclusion

The examples above demonstrate only a few possible applications of the CCDC outputs, but there are many not shown here. The R.ccdc.tools package is designed to simply working with CCDC images in R.

HunterGleason/R.ccdc.tools documentation built on July 1, 2022, 10:37 a.m.