#####"FRP_statistical" version 1.0###### ###Tested in R 4.1.3 and RStudio 2022.02.01###### ####################CALL LIBRARIES############ Be sure to install and call all the dependent packages before running library(tictoc) tic() #start timer ################CLEAN VARIABLES AND TABS################# rm(list=ls()) #clean global environment graphics.off() #clean previous plots #############SET WORKING DIRECTORY AND READ FIRMS DATASETS################# setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) #set working directory automatically / this can be changed by R-users if needed SNPP<-read.csv(list.files(pattern = "SV(.*)csv$"),sep=",",dec=".",header=TRUE) #read S-NPP data raw_times1=sprintf("%04d", unlist(SNPP["acq_time"]))#extract times and pad 0 if needed raw_dates1=sprintf(unlist(SNPP["acq_date"]))#extract dates datetime_raw1=Map(paste,raw_dates1,raw_times1)#merge dates and times datetime1=strptime(datetime_raw1,format = "%d/%m/%y %H%M")#dates and times as POSIX objects #datetime1=strptime(datetime_raw1,format = "%Y-%m-%d %H%M") #uncomment only if you have problems creating datetime1 FRP1=as.numeric(unlist(SNPP["frp"]))#extract FRP ###IMPORTANT!!! For past eruptions earlier than 2020 there is no NOAA-20 data, we strongly recommend to duplicate part of the SNPP data (3 rows) to avoid errors, another option is to comment all the lines that use the NOOA dataset#### NOAA<-read.csv(list.files(pattern = "J1V(.*)csv$"),sep=",",dec=".",header=TRUE) #read NOAA-20 data raw_times2=sprintf("%04d", unlist(NOAA["acq_time"]))#extract times and pad 0 if needed raw_dates2=sprintf(unlist(NOAA["acq_date"]))#extract dates datetime_raw2=Map(paste,raw_dates2,raw_times2)#merge dates and times datetime2=strptime(datetime_raw2,format = "%d/%m/%y %H%M")#dates and times as POSIX objects #datetime2=strptime(datetime_raw2,format = "%Y-%m-%d %H%M") #uncomment only if you have problems creating datetime1 FRP2=as.numeric(unlist(NOAA["frp"]))#extract FRP par(mfrow=c(2,2)) par(mar=c(3,5,3,3)) ####################PLOT HISTOGRAMS########################## stat_names=c('mean','minimum','5th percentile','1st quartile','median','3rd quartile','95th percentile','maximum') #names of the statistical values stat_vals1=round(c(mean(FRP1),quantile(FRP1,c(0,0.05,0.25,0.5,0.75,0.95,1))),1) #define quantiles of interest and mean stat_legend=Map(paste,stat_names,stat_vals1,sep=': ') #legend formatting for histogram hist1=hist(FRP1,xlim = c(0,stat_vals1[7]),col="orange",breaks = seq(from=min(FRP1),to=max(FRP1)+20,by=10),main=sprintf("Histogram of S-NPP FRP\n%s",paste(min(datetime1),max(datetime1),sep = ' to '))) #create a histogram with breaks of size 20 abline(v=mean(FRP1),lty=2, lwd=1.5) #add vertical black dashed line for the mean abline(v=quantile(FRP1,0.75),lty=3, lwd=1, col="blue") #add vertical blue dotted line for the 3rd quartile legend('topright',legend=stat_legend,cex=1,ncol = 1,bty="n") #legend for the histogram single column, if two columns ncol = 2 stat_vals2=round(c(mean(FRP2),quantile(FRP2,c(0,0.05,0.25,0.5,0.75,0.95,1))),1)#define quantiles of interest and mean stat_legend=Map(paste,stat_names,stat_vals2,sep=': ') #legend formatting for histogram hist2=hist(FRP2,xlim = c(0,stat_vals2[7]),col="purple3", breaks = seq(from=min(FRP2),to=max(FRP2)+20,by=10),main=sprintf("Histogram of NOAA-20 FRP\n%s",paste(min(datetime2),max(datetime2),sep = ' to ')))#create a histogram with breaks of size 20 abline(v=mean(FRP2),lty=2, lwd=1.5) #add vertical black dashed line for the mean abline(v=quantile(FRP2,0.75),lty=3, lwd=1, col="blue") #add vertical blue dotted line for the 3rd quartile legend('topright',legend=stat_legend,cex=1,ncol = 1,bty="n") #legend for the histogram single column, if two columns ncol = 2 ############CUMULATIVE STATISTICAL FRP PARAMETERS EVERY 12 HOURS TO S-NPP######### stat_namesb=c('mean','1st quartile','median','3rd quartile')#names of the statistical values init_time1=strptime(paste(as.Date(min(datetime1)),"00:00:00"),format = "%Y-%m-%d %H:%M:%S")#extract first date of period end_time1=strptime(paste(as.Date(max(datetime1)),"24:00:00"),format = "%Y-%m-%d %H:%M:%S")#extract last date of period times_stats_12a=seq(from=init_time1,to=end_time1,by=43200)#create half-day datetime sequence halfday_cumu_stats_FRPS1=setNames(data.frame(matrix(nrow = length(times_stats_12a)-1,ncol = 5)),c('timesdates12h',stat_namesb))#initialize statistics dataframe every 12 hours for (i in 1:length(times_stats_12a)-1){#calculate and assign cumulative 12 hours statistics statistics_12a=round(c(mean(FRP1[datetime1