#' Plot CCDC (Zhu & Woodcock 2014) time series sample by polygon.
#'
#' Provided a CCDC image and a polygon layer with some category field, this function plots a sample
#' of the pixels with each polygon grouped by category, indicated by the 'catg_col' parameter.
#'
#' @param ccdc_img (SpatRaster,stars) The CCDC image for which to extract coefficients, for one segment only,
#' such as the output from 'gen_latest_ccdc_rast'.
#' @param ccdc_img_date (character) A representation of the CCDC image date that can be interpreted by as.Date()
#' @param zone_poly (SpatVector,sf) Polygons for which to plot CCDC time series sample, by category.
#' @param catg_col (Character) The name of the column that defines the category of each polygon.
#' @param days The number of days for which to plot the CCDC time series, defaults to 730 (2 years)
#' @param band band (character) The name of the band for which to extract CCDC coefficients. One of 'blue', 'green', 'red'
#' 'nir', 'swir1','swir2' or 'therm'.
#' @param N The number of pixels to sample (and plot) for each category.
#' @return (void) Plots the CCDC time series for a sample of pixels stratified by polygon category.
#' @export
<-(ccdc_img,ccdc_img_date,zone_poly,catg_col,days=730,band,N=60)
{
((ccdc_img)=='stars')
{
ccdc_img<-terra::rast(ccdc_img)
}
(((zone_poly)=='sf'))
{
zone_poly<-terra::vect(zone_poly)
}
zonal_vals <- terra::extract(ccdc_img,zone_poly)
zonal_vals$ <- ((zone_poly),catg_col][zonal_vals$ID])
zonal_vals$type_int<-(zonal_vals$)
zonal_vals <- zonal_vals %>%
dplyr::group_by() %>%
dplyr::slice_sample(n=N, = T)
zonal_vals<-(zonal_vals)
zonal_vals<-zonal_vals[complete.cases(zonal_vals),]
i<-1
intp<-(band,"_INTP")
slp<-(band,"_SLP")
<-(band,"_COS")
<-(band,"_SIN")
cos2<-(band,"_COS2")
sin2<-(band,"_SIN2")
cos3<-(band,"_COS3")
sin3<-(band,"_SIN3")
ccdc_img_date<-(ccdc_img_date)
ts_x<-(ccdc_img_date:(ccdc_img_date + days))
((ts_x),(ts_x,
zonal_vals[i,intp],
zonal_vals[i,slp],
zonal_vals[i,],
zonal_vals[i,],
zonal_vals[i,cos2],
zonal_vals[i,sin2],
zonal_vals[i,cos3],
zonal_vals[i,sin3]),
=(0,1.0),
='l',
=zonal_vals$type_int[i],
xlab="Date",
ylab=(band," Surface Reflectance"))
(i (2:(zonal_vals)))
{
((ts_x),(ts_x,
zonal_vals[i,intp],
zonal_vals[i,slp],
zonal_vals[i,],
zonal_vals[i,],
zonal_vals[i,cos2],
zonal_vals[i,sin2],
zonal_vals[i,cos3],
zonal_vals[i,sin3]),
=zonal_vals$type_int[i])
}
(x='topleft',=(zonal_vals$),fill=(zonal_vals$type_int))
}
#' Plot CCDC (Zhu & Woodcock 2014) time series quartiles by polygon.
#'
#' Provided a CCDC image and a polygon layer with some category field, this function plots a sample
#' of the pixels with each polygon grouped by category, indicated by the 'catg_col' parameter.
#'
#' @param ccdc_img (SpatRaster,stars) The CCDC image for which to extract coefficients, for one segment only,
#' such as the output from 'gen_latest_ccdc_rast'.
#' @param ccdc_img_date (character) A representation of the CCDC image date that can be interpreted by as.Date()
#' @param zone_poly (SpatVector) Polygons for which to plot CCDC time series sample, by category.
#' @param catg_col (Character) The name of the column that defines the category of each polygon.
#' @param days The number of days for which to plot the CCDC time series, defaults to 730 (2 years)
#' @param band band (character) The name of the band for which to extract CCDC coefficients. One of 'blue', 'green', 'red'
#' 'nir', 'swir1','swir2' or 'therm'.
#' @param N The number of pixels to sample (and plot) for each category.
#' @return (void) Plots the CCDC time series quartiles (i.e., 25, 50 & 75th percentiles) stratified by polygon category.
#' @export
<-function (ccdc_img,ccdc_img_date,zone_poly, catg_col, days = 730, band)
{
if ((ccdc_img) == "stars") {
ccdc_img <- terra::rast(ccdc_img)
}
ccdc_img_date<-(ccdc_img_date)
zonal_vals <- terra::extract(ccdc_img, zone_poly)
zonal_vals$ <- ((zone_poly), catg_col][zonal_vals$ID])
zonal_vals <- zonal_vals[complete.cases(zonal_vals),]
i <- 1
intp <- (band, "_INTP")
slp <- (band, "_SLP")
<- (band, "_COS")
<- (band, "_SIN")
cos2 <- (band, "_COS2")
sin2 <- (band, "_SIN2")
cos3 <- (band, "_COS3")
sin3 <- (band, "_SIN3")
tStartMin<-(zonal_vals$tStart)[1]
ts_x<-(ccdc_img_date:(ccdc_img_date + days))
<-((ts_x,
zonal_vals[i, intp],
zonal_vals[i, slp],
zonal_vals[i, ],
zonal_vals[i, ],
zonal_vals[i, cos2],
zonal_vals[i, sin2],
zonal_vals[i, cos3],
zonal_vals[i, sin3]),zonal_vals$[i])
for (i (2:(zonal_vals))) {
<-(,((ts_x,
zonal_vals[i, intp],
zonal_vals[i, slp],
zonal_vals[i, ],
zonal_vals[i, ],
zonal_vals[i, cos2],
zonal_vals[i, sin2],
zonal_vals[i, cos3],
zonal_vals[i, sin3]),zonal_vals$[i]))
}
<-()
p25<-(,=(,()]),,probs=(.25))
med<-(,=(,()]),,probs=(.5))
p75<-(,=(,()]),,probs=(.75))
firstplot<-T
ts_x<-(ts_x)
(j (,()]))
{
(firstplot==T)
{
(ts_x,med[med,1]==j,(3:(med)-1)],='l',=j,=(0,1),xlab="Date",ylab=(band," Surface Reflectance"))
firstplot=F
}{
(ts_x,med[med,1]==j,(3:(med)-1)],='l',=j)
}
(ts_x,p25[med,1]==j,(3:(p25)-1)],='l',=j)
(ts_x,p75[med,1]==j,(3:(p75)-1)],='l',=j)
}
(x = "topleft", = (zonal_vals$), fill = ((,()])))
}
#' Uses SAGA-GIS and RSAGA to create N number of random samples with provided polygons, inherits polygon attributes. Points also get attributed
#' with the coefficient from a CCDC image for a given date. Points are returned as a data.frame. Useful function for
#' a large area of interest, i.e., as a alternative to 'plot_ccdc_ts_quartiles'.
#'
#'
#' @param polygons_pth (character) File path to polygon vector file with categories of interest.
#' @param ccdc_img_pth (character) File path to raster file (e.g. GeoTiff) CCDC image for given date (see 'gen_ccdc_at_jdoy')
#' @param pnts_count (integer) Total number of X-Y locations to sample throughout ALL polygons.
#' @param sep_dist (integer) Minimum distance in map units that two sample locations can be.
#' @param bands (character) Vector of band names present in CCDC image (in order), defaults to c('blue','green','red','nir','swir1','swir2','therm')
#' @param env (RSAGA environment) A RSAGA environment to use.
#' @return (data.frame) A data.frame with the sample locations, attributed with the polygon data and CCDC grid coefficient values.
#' @export
<-(polygons_pth,ccdc_img_pth,pnts_count,sep_dist,bands= ('blue','green','red','nir','swir1','swir2','therm') ,env)
{
("Loading Polygons ...","\n")
RSAGA::rsaga.geoprocessor(lib="io_gdal",
module = 3,
param = (FILES=polygons_pth,
SHAPES=((),'shapes')),
env=env)
("Generating Random Points ...","\n")
RSAGA::rsaga.geoprocessor(lib='shapes_points',
module=21,
param = (POINTS=((),"random_pnts"),
EXTENT=3,
POLYGONS=((),'shapes.shp'),
COUNT=pnts_count,
DISTANCE=sep_dist),
env=env)
("Atributing Polygons to Points ...","\n")
RSAGA::rsaga.geoprocessor(lib='shapes_points',
module=10,
param = (INPUT=((),"random_pnts.shp"),
POLYGONS=((),'shapes.shp')),
env=env)
<-(sf::read_sf(((),"random_pnts.shp")))
("Atributing Grids to Points, this may take some time ...","\n")
RSAGA::rsaga.geoprocessor(lib='shapes_grid',module=0,param = (SHAPES=((),"random_pnts.shp"),
GRIDS=ccdc_img_pth,
RESULT=((),"random_pnts.shp"),
RESAMPLING=0),
env=env)
ccdc_img_at_pnts<-sf::read_sf(((),"random_pnts.shp"))
band_names<-()
i<-1
(b bands)
{
(coe ("INTP","SLP","COS","SIN","COS2","SIN2","COS3","SIN3","RMSE"))
{
band_names[i]<-(b,"_",coe)
i<-i+1
}
}
band_names<-(band_names,"tStart","tEnd")
(ccdc_img_at_pnts)[c(:((ccdc_img_at_pnts)-1))]<-band_names
ccdc_img_at_pnts<-(ccdc_img_at_pnts)
("Done ...","\n")
(ccdc_img_at_pnts)
}
#' Function for plotting surface reflectance curves by category (i.e., output of 'sample_ccdc_by_catg').
#'
#'
#' @param ccdc_img_at_pnts (data.frame) A data.frame containing the categories of interest, and associated CCDC coefficients for a given date (i.e., output of 'sample_ccdc_by_catg')
#' @param ccdc_img_date (character) A representation of the CCDC image date that can be interpreted by as.Date()
#' @param catg_field (character) Name of the column indicating the categories of interest.
#' @param band (character) The name of the band for which to forecast CCDC surface reflectance time series. One of 'blue', 'green', 'red'
#' 'nir', 'swir1','swir2' or 'therm'.
#' @param days (integer) Number of days to predict surface reflectance from the CCDC image data ('ccdc_img_date').
#' @param med_only (boolean) When true, only median surface reflectance values are plotted by category
#' @param ylim (numeric) Vector of length 2 defining the Y-Axis limits, defaults to c(0,1)
#' @return (List) A data.frame with all data used for plot generation ('srdata'), a 'ggplot' object ('srplot'), and a data.frame with the legend codes cross-walk ('crosswalk').
#' @export
<-(ccdc_img_at_pnts,ccdc_img_date,catg_field,band,days=730,med_only=F,=(0,1))
{
Rcpp::cppFunction('NumericMatrix predict_sr(int img_date,int days,NumericVector catg, NumericVector coef_intp,NumericVector coef_slp,NumericVector coef_cos,NumericVector coef_sin,NumericVector coef_cos2,NumericVector coef_sin2,NumericVector coef_cos3,NumericVector coef_sin3){
int n = catg.size();
NumericMatrix out(n*days,3);
double pi = 3.1415926535;
int row_idx=0;
int end_date=img_date+days;
for(int i = 0; i < n; i++)
{
double day = img_date;
int current_catg=catg[i];
while(day<(double)end_date)
{
double yhat=coef_intp[i]+coef_slp[i]*day+coef_cos[i]*cos((2.0*pi*day)/365.25)+coef_sin[i]*sin((2.0*pi*day)/365.25)+coef_cos2[i]*cos((4.0*pi*day)/365.25)+coef_sin2[i]*sin((4.0*pi*day)/365.25)+coef_cos3[i]*cos((6.0*pi*day)/365.25)+coef_sin3[i]*sin((6.0*pi*day)/365.25);
out(row_idx,2)=yhat;
out(row_idx,0)=current_catg;
out(row_idx,1)=day;
row_idx++;
day=day+1.0;
}
}
return out;
}')
factor_vector <- ((ccdc_img_at_pnts,catg_field]))
intp <- ccdc_img_at_pnts,(band,"_INTP")]
slp <- ccdc_img_at_pnts,(band,"_SLP")]
<- ccdc_img_at_pnts,(band,"_COS")]
<- ccdc_img_at_pnts,(band,"_SIN")]
cos2 <- ccdc_img_at_pnts,(band,"_COS2")]
sin2 <- ccdc_img_at_pnts,(band,"_SIN2")]
cos3 <- ccdc_img_at_pnts,(band,"_COS3")]
sin3 <- ccdc_img_at_pnts,(band,"_SIN3")]
ccdc_img_date<-(ccdc_img_date)
sr_mat<-predict_sr(ccdc_img_date,days,factor_vector,intp,slp,,,cos2,sin2,cos3,sin3)
sr_mat<-(sr_mat)
str_fact<-(ccdc_img_at_pnts,catg_field])
num_fact<-(((ccdc_img_at_pnts,catg_field])))
lab_tabl<-((num_fact,str_fact))
(lab_tabl)<-("Category","Code")
(sr_mat)<-('catg','jdoy','SR')
sr_mat$catg=(sr_mat$catg)
sr_mat$jdoy<-(sr_mat$jdoy)
median_sr<-sr_mat %>%
dplyr::group_by(jdoy,catg) %>%
dplyr::summarise(sr=(SR,na.rm=T))
p25_sr<-sr_mat %>%
dplyr::group_by(jdoy,catg) %>%
dplyr::summarise(sr=(SR,probs=0.25,na.rm=T))
p75_sr<-sr_mat %>%
dplyr::group_by(jdoy,catg) %>%
dplyr::summarise(sr=(SR,probs=0.75,na.rm=T))
(!med_only)
{
gg<-ggplot2::ggplot()+
ggplot2::geom_line(=p25_sr,ggplot2::aes(x=jdoy,y=sr,color=catg),linetype=2)+
ggplot2::geom_line(=median_sr,ggplot2::aes(x=jdoy,y=sr,color=catg),linetype=1)+
ggplot2::geom_line(=p75_sr,ggplot2::aes(x=jdoy,y=sr,color=catg),linetype=2)+
ggplot2::xlab("Date")+
ggplot2::ylab((band," Surface Reflectance"))+
ggplot2::labs(color = "Category")+
ggplot2::theme_classic()+
ggplot2::()
}{
gg<-ggplot2::ggplot()+
ggplot2::geom_line(=median_sr,ggplot2::aes(x=jdoy,y=sr,color=catg),linetype=1)+
ggplot2::xlab("Date")+
ggplot2::ylab((band," Surface Reflectance"))+
ggplot2::labs(color = "Category")+
ggplot2::theme_classic()+
ggplot2::()
}
((srdata = sr_mat,srplot=gg,crosswalk=lab_tabl))
}
