Meteorological Review

1 Introduction

Longitude= -110.86
Latitude = 33.40

source(here::here("scripts","hydro_support.r"))
#For more information
#https://developers.google.com/maps/gmp-get-started\
key_google<-rstudioapi::askForSecret("Google Key")

ggmap::register_google(key = key_google)

#For more information
# https://www.ncdc.noaa.gov//cdo-web/token
key_noaa<-rstudioapi::askForSecret("NOAA Key")

#Get elevation from google maps
Elevation=GM_elev(Longitude,Latitude,key_google)

site.info<-data.frame(Location="Site",
                      Longitude=Longitude,
                      Latitude=Latitude,
                      Elevation=Elevation) 

#Site in looped with variables
site.info_xy<-site.info %>% 
  dplyr::rename(x="Longitude",y="Latitude",z="Elevation") %>% 
  reshape2::melt(id.vars="Location")

pandoc.table(site.info)

Location	Longitude	Latitude	Elevation
Site	-110.9	33.4	1091

#minimum year of information used, to be defined after seen the available info
min.ywi<-10 

2 SOURCES

2.1 NOAA

Local meteorological station records from publicly available National Oceanic and Atmospheric Administration (NOAA) databases: Global Historical Climatology Network Database (GHCND) and Global Surface Summary of the Day (GSOD) were obtained to implement the site data. These stations typically compile precipitation and temperature daily records. A radius of about ~ 100 km around the site was used as an area of interest to compile these meteorological records.

From NOAA, the amount of information in the region of interest is limited. Most of the relevant records present a gap of information after 1990. Meteorological stations with records less than 10years were removed.

On the original analysis, met. stations with information ended before 2000 was removed.

NOAA.file<-here::here("data","NOAA","NOAA Info.rdata")

if(file.exists(NOAA.file)){
  
  load(NOAA.file) 
  
}else{
  
  NOAA.find(longitude = Longitude,latitude = Latitude,radius_deg = 0.7,
            #Get your own key from https://www.ncdc.noaa.gov//cdo-web/token
            #Key from vmunoz@srk.co.uk
            noaa.key = "PacdvUWVYvHOPwXrHRILLOJpvDnORxEB")
  
  selected.NOAA.GHCND<-selected.NOAA.GHCND[selected.NOAA.GHCND$station.distance<50,]
  #work just with the first 10 stations. Ranges can change based on the project
  selected.NOAA.GHCND$maxdate<-if_else(selected.NOAA.GHCND$maxdate>ymd("2020-01-01"),
          as.Date(lubridate::now()),selected.NOAA.GHCND$maxdate)
  
  selected.NOAA.GSOD$END<-if_else(ymd(selected.NOAA.GSOD$END)>ymd("20200101"),
          str_remove_all(as.Date(lubridate::now()),"-") %>% as.integer(),
          selected.NOAA.GSOD$END)
  
  NOAA.capture.GHCND()
  
  NOAA.capture.GSOD()
  
  NOAA.prepare()
  
  
  save(file = NOAA.file,list = c("NOAA.db","NOAA.idx"))  
  
}

NOAA.idx %>% 
  dplyr::select(-WMO_ID,-ICAO) %>% 
  pandoc.table(style="simple",split.table=Inf,
               caption="NOAA STATIONS")

NOAA STATIONS

ID	NAME	LATITUDE	LONGITUDE	ELEVATION	MINDATE	MAXDATE	DISTANCE
USC00025512	MIAMI	33.4	-110.9	1085	1914-01-01	2017-12-31	1.05
USC00023501	GLOBE #3	33.39	-110.8	1116	2009-01-01	2023-02-07	6.23
USC00023500	GLOBE #2	33.4	-110.8	1091	1975-01-01	2008-12-31	7.2
USC00023498	GLOBE RS	33.38	-110.8	1083	1894-01-01	1975-12-31	7.37
USC00023505	GLOBE	33.38	-110.8	1112	1981-01-01	2006-12-31	8.66
US1AZGL0010	GLOBE 1.8 SSE	33.37	-110.8	1146	2017-01-01	2019-12-31	9.24
USC00020098	ALAMO RS	33.5	-110.8	926.9	1915-01-01	1946-12-31	11.17
USC00026561	PINAL RCH	33.35	-111	1379	1895-01-01	1973-12-31	12.74
US1AZPN0072	MIAMI 8.3 WSW	33.34	-111	1378	2017-01-01	2023-02-07	13.93
USR0000AGLO	GLOBE ARIZONA	33.35	-110.7	1112	1990-01-01	2023-02-07	20.12
USC00028349	SUPERIOR 2 ENE	33.3	-111.1	1266	1974-01-01	1996-12-31	21.98
USC00028348	SUPERIOR	33.3	-111.1	871.7	1920-01-01	2006-12-31	24.69
USC00028351	SUPERIOR SMELTER	33.3	-111.1	851	1948-01-01	1953-12-31	24.94
USC00024345	INTAKE	33.62	-110.9	677	1906-01-01	1952-12-31	25.06
USC00023643	GRAPEVINE	33.63	-111	677	1943-01-01	1950-12-31	31.39
USC00024594	KELVIN	33.1	-111	563.9	1917-01-01	1983-12-31	34.84
USC00027477	SAN CARLOS AP	33.38	-110.5	880.9	1977-01-01	1980-12-31	36.6
USR0000ASCA	SAN CARLOS #1 ARIZONA	33.37	-110.5	846.4	2002-01-01	2023-02-07	37.6
USR0000AROO	ROOSEVELT ARIZONA	33.66	-111.1	664.5	1992-01-01	2009-12-31	38.07
USC00027475	SAN CARLOS	33.35	-110.5	805	1912-01-01	1977-12-31	38.52
USC00027480	SAN CARLOS RESERVOIR	33.18	-110.5	771.8	1900-01-01	2023-02-07	39.43
USC00024590	KEARNY	33.05	-110.9	557.8	1984-01-01	2015-12-31	39.46
USC00027281	ROOSEVELT 1 WNW	33.67	-111.2	672.1	1905-01-01	2023-02-07	40.65
US1AZPN0077	QUEEN VALLEY 0.2 E	33.3	-111.3	627.9	2018-01-01	2023-02-07	41.26
US1AZPN0069	QUEEN VALLEY 1.1 SW	33.29	-111.3	610.2	2016-01-01	2018-12-31	43.16
USC00028940	PARKER CREEK MNTC YRD	33.8	-111	1678	1948-01-01	1951-12-31	45.3
USC00027876	SIERRA ANCHA	33.8	-111	1554	1913-01-01	1979-12-31	45.56
USS0010S01S	Workman Creek	33.81	-110.9	2103	1978-01-01	2023-02-07	45.98
USC00029534	WORKMAN CREEK 1	33.82	-110.9	2126	1948-01-01	1951-12-31	46.68
USR0000AHIL	HILLTOP ARIZONA	33.62	-110.4	1744	2000-01-01	2023-02-07	47.52
USC00024069	HILLTOP	33.62	-110.4	1739	1939-01-01	1948-12-31	49.04
749173_99999	COOLIDGE MUNI	32.94	-111.4	479.8	1944-07-11	2023-02-07	73.89
722786_99999	WILLIAMS GATEWAY	33.3	-111.7	421	2000-01-01	2003-12-31	74.3
722786_23104	WILLIAMS GATEWAY AIRPORT	33.3	-111.7	421.2	1942-03-01	2023-02-07	75.86

#used for MERRA2 - small area coordinates +/- 2.5 degrees
xy <-ggmap::make_bbox(lon = c(0.5*round(Longitude/0.5,1)-2.5,0.5*round(Longitude/0.5,1)+2.5),
                      lat = c(0.5*round(Latitude/0.5,1)-2.5,0.5*round(Latitude/0.5,1)+2.5),
                      f = 0)

#used for ggsn graphical scale, same ranges but 8% smaller
xy_min <-ggmap::make_bbox(lon = c(xy[1],xy[3]),
                      lat = c(xy[2],xy[4]),f=-0.08) 

#ggmap capture for goggle maps
map <- get_map(xy,scale=2,zoom=7,maptype="hybrid",source="google",
               messaging=FALSE)

NOAA.file.par<-here::here("data","NOAA","NOAA Info par.rdata")

if(file.exists(NOAA.file.par)){

    load(NOAA.file.par) 
  
}else{
  
    par.review<-data.frame(var=c("PRCP"),
                         op=c("sum"),
                         names=c("TOTAL PRECIPITATION"),
                         stringsAsFactors = F)
  
  apply(par.review,1,FUN=function(par){
    
  # par<-par.review[2,]
    
      #variable definition
    par.var<<-par[1]
    par.fun<<-eval(parse(text = par[2]))
    par.title<<-par[3]
    
    #Title (two spaces before \n )
    pander::pandoc.header(par.title,level=2)
    
    # files
    db.name<-paste0("NOAA.",par.var,".db")
    idx.name<-paste0("NOAA.",par.var,".idx")
    syn.name<-paste0("NOAA.",par.var,".syn")
    
    #Prepare of DB
    NOAA.par.db<-NOAA.db %>% 
      dplyr::filter(variable%in%par.var) %>%
      dplyr::select(-NAME,-variable) %>% 
      as.data.frame() %>%  
      split(x=.[c("DATE","value")],f=.$ID)
    
    #Prepare index  
    NOAA.par.idx<-NOAA.idx[NOAA.idx$ID%in%names(NOAA.par.db),]
    
    #Sort db
    NOAA.par.db<-NOAA.par.db[NOAA.par.idx$ID]
    
    #days with information
    st.dwi<-sapply(NOAA.par.db,FUN=function(x){dwi(x)%>%sum})
    
    #1.selection of the stations with information more than 'n' years
    st.sel<-st.dwi>=min.ywi*365.25
    
    #clean the DB based on dwi
    NOAA.par.db<-NOAA.par.db[st.sel]
    NOAA.par.idx<-NOAA.par.idx[st.sel,]
    
    Meteo.analysis(db = NOAA.par.db,
                   FUN.month =par.fun,
                   Min.days.annual = 300,
                   Max.miss.days.month = 5,output = syn.name,EC = F,
                   start.year = 1900) %>% invisible()
    
    # #Mean annual based on a range
    NOAA.par.idx %<>%mutate(ANN=eval(parse(text=syn.name))$"Daily.Comp" %>%
                              # window(start=as.Date("1980-01-01"),
                              #        end=as.Date("2019-12-31")) %>%
                              daily2annual.avg(FUN = par.fun) %>%
                              .$Ann)
    #assign the variables
    assign(x = db.name,value=NOAA.par.db,envir = .GlobalEnv)
    assign(x = idx.name,value=NOAA.par.idx,envir = .GlobalEnv)
    
    #remove the auxiliar variables
    rm(par.var,par.fun,par.title,
       envir = .GlobalEnv)
    
    pander::pandoc.p('')
    
    pander::pandoc.p('')
    
  })
  
  save(file = NOAA.file.par,
       list =c(glue::glue('NOAA.{par.review$var}.db'),
               glue::glue('NOAA.{par.review$var}.syn'),
               glue::glue('NOAA.{par.review$var}.idx')))  
  
  
}

2.2 SRTM

Information obtained indirectly as a CGIAR-SRTM (90 m resolution), this topography covers an extension from 60 deg Lat North / South. For area beyond these limitations, it is recommend to use ArticDEM or another topographical source.

#Download SRTM
SRTM_ras<-SRTM_capture(bbox = xy,file = "SRTM",
                               trim = T,overwrite = T,
                               maxcell =5e5)

# to a data.frame
SRTM_tbl<-raster::rasterToPoints(SRTM_ras) %>% 
  data.frame()

names(SRTM_tbl)<-c("x","y","value")

SRTM_tbl%<>%
  dplyr::filter(complete.cases(value))

plot(SRTM_ras)

points(Longitude,Latitude)

2.3 ERA5 - Land

As a reference, it was also ERA5-Land. Era-5 is a climatic gridded model from 1950 to the present.

era5_dt<-arrow::read_parquet(here::here("data","rds","bhp_era5land_prcp.par"))

era5_x<-era5_dt$x %>% unique()
era5_y<-era5_dt$y %>% unique()

#MERRA TS
sel_LON<-
  era5_x[(Longitude-era5_x)^2%in%min((Longitude-era5_x)^2)][1]

sel_LAT<-
  era5_y[(Latitude-era5_y)^2%in%min((Latitude-era5_y)^2)][1]

era5_z<-era5_dt %>% 
  dplyr::filter(x==sel_LON, y==sel_LAT) %>% 
  dplyr::select(-x,-y) %>%
  as.data.frame() %>% 
  dplyr::rename(date="variable") %>% 
  mutate(date=as.Date(date)) %>% 
  tk_zoo(silent=T)

plot.zoo(era5_z,yax.flip = T,
         main="Local time series based on era5-land",ylab="mm",xlab=NULL)

2.4 Daymet

library(daymetr)

daymet_info<-daymetr::download_daymet(site = "Daymet",lat = Latitude,lon = Longitude,
                         start = 2000,end=lubridate::year(now())-1,silent = T)


daymet_z<-daymet_info$data %>%
  dplyr::mutate(dates=as.Date(paste0(year,"-01-01"))+lubridate::days(yday)-1) %>%
  dplyr::select(dates,prcp="prcp..mm.day.",tmax="tmax..deg.c.",tmin="tmin..deg.c.") %>%
  dplyr::mutate(tavg=0.5*tmax+0.5*tmin,
                rain=if_else(tavg>0,prcp,0),
                snow=if_else(tavg<=0,prcp,0)) %>%
  tk_zoo(silent=T)

plot(daymet_z$prcp, yax.flip = T,
         main="Local time series based on Daymet",ylab="mm")

3 Total Precipitation

NOAA.PRCP.idx$NAME<-NOAA.PRCP.idx$NAME %>%make.names(unique=T) %>% 
  str_replace_all("\\."," ") %>% 
  str_wrap(string = .,width = 20)


st.names<-glue('{NOAA.PRCP.idx$NAME}\n{round(NOAA.PRCP.idx$ANN,1)} mm/yr')

st.names<-glue('{NOAA.PRCP.idx$NAME}')


dwi.graph(x=NOAA.PRCP.syn$"Daily.Comp",
            station.name = st.names,
            start.year = 1800)

3.1 Station selection

Station selected to be in the band between 1000 and 1200 masl.

Based on this definition these stations were revised:

#Reduced the station in the 1000 to 1200 band
noaa_prcp_sel<-NOAA.PRCP.idx %>% 
  dplyr::filter(ELEVATION>1000,ELEVATION<1200)

noaa_prcp_comb<-merge(
  NOAA.PRCP.syn$Daily.Comp[,noaa_prcp_sel$ID],
  era5_z,
  daymet_z$prcp) 

names(noaa_prcp_comb)<-c(noaa_prcp_sel$NAME,"ERA5","daymet")

pandoc.table(noaa_prcp_sel %>% 
               dplyr::select(-WMO_ID,-ICAO),split.table=Inf,style="simple")

#used for ggsn graphical scale, same ranges but 8% smaller
xy_site <-ggmap::make_bbox(lon = noaa_prcp_sel$LONGITUDE,
                      lat = noaa_prcp_sel$LATITUDE,
                      f=0.1) 

#ggmap capture for goggle maps
map_site <- get_map(xy_site,scale=2,zoom=12,maptype="hybrid",source="google",
               messaging=FALSE)



ggmap(map_site)+
  geom_point(data=noaa_prcp_sel,aes(x=LONGITUDE,y=LATITUDE),size=5)+
  geom_label_repel(data=noaa_prcp_sel,aes(x=LONGITUDE,y=LATITUDE,label=NAME),size=2.5)+
  ggsn::scalebar(x.min = -110.9,x.max=-110.73,y.min=33.31,y.max=33.45,dist = 2.5,dist_unit = "km",transform = T,model = "WGS84")

srtm_site_tbl<-raster::crop(SRTM_ras,
                            raster(xmn=-110.92,xmx=-110.73,ymn=33.30,ymx=33.47)) %>% 
  rasterToPoints() %>% data.frame() %>% rename(z=layer)

ggmap(map_site)+
  geom_tile(data=srtm_site_tbl,aes(x=x,y=y,
                                   fill=cut(z,c(800,900,1000,
                                                1100,1200,1500,2000,2500),dig.lab=7)),
            alpha=0.6)+
  geom_point(data=noaa_prcp_sel,aes(x=LONGITUDE,y=LATITUDE),size=5)+
  geom_label_repel(data=noaa_prcp_sel,aes(x=LONGITUDE,y=LATITUDE,label=NAME),size=2.5)+
  scale_fill_brewer(palette = "Spectral")+
  ggsn::scalebar(x.min = -110.9,x.max=-110.73,y.min=33.31,y.max=33.45,dist = 2.5,
                 dist_unit = "km",transform = T,model = "WGS84")+
  labs(fill="Elevation\n[masl]")

3.2 Relationships between stations

noaa_prcp_comb %>% 
  daily2monthly.V2(FUN=sum) %>% 
  fortify.zoo() %>% 
  dplyr::filter(Index>=as.Date("1950-01-01")) %>% 
  reshape2::melt(id.vars=c("Index","ERA5")) %>%
  mutate(yr=year(Index),
         dec=cut(yr,seq(1940,2030,10),dig.lab=6,include.lowest=T)) %>% 
  mutate(value=if_else(value>1000,NaN,value)) %>% 
  ggplot(data=.,aes(x=ERA5,y=value))+
  geom_point(alpha=0.3)+
  stat_poly_eq(use_label(c("adj.R2")),
               formula = x ~ poly(y, 1, raw = TRUE),size=2.5)+
  geom_abline(slope = 1,col="red")+
  facet_grid(variable~dec)+
  theme_light()+
  theme(strip.text.y.right = element_text(angle=0))+
  labs(x="Monthly Precipitation by ERA5-Land [mm/mon]",
       y="Monthly Precipitation by source [mm/mon]")

psych::cor.plot(noaa_prcp_comb)

noaa_prcp_comb %>% 
  daily2monthly.V2(FUN=sum) %>% 
  monthly2annual(FUN=sum,na.rm=F) %>% 
  fortify.zoo() %>% 
  # dplyr::filter(Index>=as.Date("1950-01-01")) %>% 
  reshape2::melt(id.vars=c("Index","ERA5","daymet")) %>% 
  mutate(value=if_else(value>1000,NaN,value)) %>% 
  mutate(yr=year(Index),
         dec=cut(yr,seq(1840,2030,10),dig.lab=6,include.lowest=T)) %>% 
  ggplot(data=.,aes(x=Index,y=value,fill=variable))+
  geom_col(position = "dodge")+
  geom_line(aes(x=Index,y=ERA5),linetype="dotted",size=1)+
  geom_line(aes(x=Index,y=daymet),linetype="dashed",size=1)+
  
  facet_wrap(~dec,scales="free_x")+
  theme_light()+
  labs(y="annual precipitation [mm]",title="annual precipitation per decade",
       subtitle="ERA5 presented as dotted line/ Daymet presented as dashed line",fill="station:")

3.3 Double mass relationships

noaa_prcp_mon<-noaa_prcp_comb %>% 
  daily2monthly.V2(FUN=sum) %>% 
  monthly2annual(FUN=sum,na.rm=F)

lapply(noaa_prcp_mon,FUN=function(x){x %>% na.locf() %>% cumsum()}) %>% 
  do.call(merge,.) %>% 
  fortify.zoo() %>% 
  reshape2::melt(id.vars=c("Index","GLOBE RS")) %>% 
   mutate(yr=year(Index)) %>% 
  ggplot(data=.,aes(x=`GLOBE RS`,y=value,col=yr))+
  geom_point()+
    stat_poly_eq(use_label(c("eq", "adj.R2")),
               formula = y ~ x,size=2.5,rr.digits = 5)+
  geom_smooth(method="lm")+
  geom_abline(slope=1,col="red")+
  facet_wrap(~variable)+
  scale_color_binned(type = "viridis")+
  labs(title="Double mass figure",
       subtitle="comparing values from GLOBE RS",y="compared value [mm]",
       x="GLOBE RS[mm]")+
  theme_light()

lapply(noaa_prcp_mon,FUN=function(x){x %>% na.locf() %>% cumsum()}) %>% 
  do.call(merge,.) %>% 
  fortify.zoo() %>% 
  reshape2::melt(id.vars=c("Index","MIAMI")) %>% 
    mutate(yr=year(Index)) %>% 
  ggplot(data=.,aes(x=`MIAMI`,y=value,col=yr))+
  geom_point()+
      stat_poly_eq(use_label(c("eq", "adj.R2")),
               formula = y ~ x,size=2.5,rr.digits = 5)+
  geom_smooth(method="lm")+
  geom_abline(slope=1,col="red")+
  facet_wrap(~variable)+
  scale_color_binned(type = "viridis")+
    labs(title="Double mass figure",
       subtitle="comparing values from MIAMI",y="compared value [mm]",
       x="MIAMI[mm]")+
  theme_light()

lapply(noaa_prcp_mon,FUN=function(x){x %>% na.locf() %>% cumsum()}) %>% 
  do.call(merge,.) %>% 
  fortify.zoo() %>% 
  reshape2::melt(id.vars=c("Index","ERA5")) %>% 
  mutate(yr=year(Index)) %>% 
  ggplot(data=.,aes(x=`ERA5`,y=value,col=yr))+
  geom_point()+
      stat_poly_eq(use_label(c("eq", "adj.R2")),
               formula = y ~ x,size=2.5,rr.digits = 5)+
  geom_smooth(method="lm")+
  geom_abline(slope=1,col="red")+
  facet_wrap(~variable)+
  scale_color_binned(type = "viridis")+
      labs(title="Double mass figure",
       subtitle="comparing values from ERA5",y="compared value [mm]",
       x="ERA5[mm]")+
  theme_light()

3.4 Monthly Relationship between MIAMI / GLOBE RS / GLOBE #3

names(noaa_prcp_comb)<-make.names(names(noaa_prcp_comb))


noaa_prcp_comb %>% 
  daily2monthly.V2(FUN=sum) %>% 
  ggplot(data=.,aes(y=`GLOBE.RS`,x=MIAMI))+
  geom_point()+
        stat_poly_eq(use_label(c("eq", "adj.R2")),
               formula = y ~ x,size=5,rr.digits = 5)+
  geom_abline(slope=1,col="red")+
  geom_smooth(method="lm")+
  theme_light()+
  labs(x="Globe RS [mm/mon]",y="MIAMI [mm/mon]",
       title="Monthly precipitation")

noaa_prcp_comb %>% 
  daily2monthly.V2(FUN=sum) %>% 
  ggplot(data=.,aes(x=`GLOBE.3`,y=MIAMI))+
  geom_point()+
        stat_poly_eq(use_label(c("eq", "adj.R2")),
               formula = y ~ x,size=5,rr.digits = 5)+
  geom_abline(slope=1,col="red")+
  geom_smooth(method="lm")+
  theme_light()+
  labs(x="Globe #3 [mm/mon]",y="MIAMI [mm/mon]",
       title="Monthly precipitation")

## Comparison on actual report records and Globe RS

load(here::here("data","rds","Backup Daymet Local review_2021-05-14.rdata"))

merge(noaa_globe3=zoo(globe3_raw.dt$value,globe3_raw.dt$date),
      patch1_globe3=globe3.adj,
      patch2_globe3=globe3b.adj) %>% 
  # window(start="1980-01-01") %>% 
  plot.zoo(yax.flip = T,main="Globe 3 comparison - Total Prcp [mm/day]")

prcp_report_z<-merge(noaa_globe3=zoo(globe3_raw.dt$value,globe3_raw.dt$date),
      patch1_globe3=globe3.adj,
      patch2_globe3=globe3b.adj)


time(prcp_report_z)<-as.Date(time(prcp_report_z))

merge(prcp_report_z,
      GLOBE.RS=noaa_prcp_comb$GLOBE.RS) %>% 
  window(start=as.Date("1950-01-01"),end=as.Date("1969-12-31")) %>% 
  daily2annual.avg(FUN=sum) %>% 
  pandoc.table(round=0,missing="-",split.table=Inf)

Station	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec	Ann
noaa_globe3	-	-	-	-	-	-	-	-	-	-	-	-	-
patch1_globe3	-	-	-	-	-	-	-	-	-	-	-	-	-
patch2_globe3	45	29	42	15	5	8	55	89	27	25	28	59	427
GLOBE.RS	39	26	36	13	6	10	62	87	25	28	24	51	407

merge(prcp_report_z,
      GLOBE.RS=noaa_prcp_comb$GLOBE.RS) %>% 
  daily2monthly.V2(FUN=sum,na.rm = T) %>% 
  monthly2annual(FUN = sum,na.rm = F) %>% 
  window(end=as.Date("2019-12-31")) %>% 
  fortify.zoo() %>% 
  reshape2::melt(id.vars="Index") %>% 
  mutate(yr=lubridate::year(Index),
         bi_dec=cut(yr,seq(1800,2020,20),dig.lab=8)) %>% 
  ggplot(data=.,aes(x=yr,y=value,col=variable))+
  geom_point()+
  geom_line()+
  facet_wrap(~bi_dec,scales="free_x")+
  theme_light()+
  labs(x="yrs",y="mean annual prcp [mm]")+
  theme(legend.position="bottom")

merge("SRK.2021"=prcp_report_z$patch2_globe3,
      GLOBE.RS=noaa_prcp_comb$GLOBE.RS) %>% 
  daily2monthly.V2(FUN=sum,na.rm = T) %>% 
  monthly2annual(FUN = sum,na.rm = F) %>% 
  window(start=as.Date("1900-01-01"),
         end=as.Date("1979-12-31")) %>% 
  fortify.zoo() %>% 
  # reshape2::melt(id.vars="Index") %>% 
  mutate(yr=lubridate::year(Index),
         bi_dec=cut(yr,seq(1900,1980,20),dig.lab=8),include.lowest=T) %>% 
  dplyr::filter(complete.cases(bi_dec)) %>% 
  ggplot(data=.,aes(x=SRK.2021,y=GLOBE.RS))+
  geom_point()+
          stat_poly_eq(use_label(c("eq", "adj.R2")),
               formula = y ~ x,size=2.5,rr.digits = 5)+
  geom_smooth(method="lm")+
  geom_abline(slope=1,col="red")+
  facet_wrap(~bi_dec)+
  theme_light()+
  # labs(x="yrs",y="mean annual prcp [mm]")+
  theme(legend.position="bottom")+
  labs(title="Annual precipitation",subtitle="information divided by bi-decade",
       x="Annual precipitation - SRK 2021 [mm/yr]",
       y="Annual precipitation - Globe RS [mm/yr]")

dwi_tbl<-dwi(prcp_report_z$patch2_globe3) %>% 
  fortify.zoo() %>% 
  dplyr::rename("days_with_information"=".",
                "year"="Index") %>% 
  mutate(year=as.numeric(as.character(year)),
         day_in_year=as.numeric(ymd(paste0(year,"-12-31"))-ymd(paste0(year,"-01-01"))+1),
         missing_days=day_in_year-`days_with_information`,
         completeness=`days_with_information`/day_in_year*100) %>% 
  dplyr::filter(year>=1910)

pandoc.table(dwi_tbl,round=1)

year	days_with_information	day_in_year	missing_days	completeness
1910	0	365	365	0
1911	0	365	365	0
1912	0	366	366	0
1913	0	365	365	0
1914	332	365	33	91
1915	363	365	2	99.5
1916	366	366	0	100
1917	365	365	0	100
1918	365	365	0	100
1919	365	365	0	100
1920	366	366	0	100
1921	365	365	0	100
1922	365	365	0	100
1923	365	365	0	100
1924	361	366	5	98.6
1925	365	365	0	100
1926	365	365	0	100
1927	365	365	0	100
1928	366	366	0	100
1929	365	365	0	100
1930	365	365	0	100
1931	362	365	3	99.2
1932	366	366	0	100
1933	365	365	0	100
1934	365	365	0	100
1935	365	365	0	100
1936	355	366	11	97
1937	365	365	0	100
1938	348	365	17	95.3
1939	365	365	0	100
1940	366	366	0	100
1941	365	365	0	100
1942	365	365	0	100
1943	365	365	0	100
1944	362	366	4	98.9
1945	365	365	0	100
1946	362	365	3	99.2
1947	365	365	0	100
1948	366	366	0	100
1949	365	365	0	100
1950	365	365	0	100
1951	362	365	3	99.2
1952	362	366	4	98.9
1953	365	365	0	100
1954	365	365	0	100
1955	365	365	0	100
1956	366	366	0	100
1957	365	365	0	100
1958	359	365	6	98.4
1959	365	365	0	100
1960	362	366	4	98.9
1961	365	365	0	100
1962	358	365	7	98.1
1963	365	365	0	100
1964	366	366	0	100
1965	365	365	0	100
1966	363	365	2	99.5
1967	365	365	0	100
1968	366	366	0	100
1969	365	365	0	100
1970	365	365	0	100
1971	365	365	0	100
1972	366	366	0	100
1973	365	365	0	100
1974	365	365	0	100
1975	365	365	0	100
1976	365	366	1	99.7
1977	363	365	2	99.5
1978	365	365	0	100
1979	365	365	0	100
1980	366	366	0	100
1981	365	365	0	100
1982	365	365	0	100
1983	365	365	0	100
1984	365	366	1	99.7
1985	365	365	0	100
1986	365	365	0	100
1987	365	365	0	100
1988	366	366	0	100
1989	365	365	0	100
1990	365	365	0	100
1991	365	365	0	100
1992	366	366	0	100
1993	365	365	0	100
1994	365	365	0	100
1995	365	365	0	100
1996	366	366	0	100
1997	365	365	0	100
1998	365	365	0	100
1999	365	365	0	100
2000	366	366	0	100
2001	365	365	0	100
2002	365	365	0	100
2003	365	365	0	100
2004	366	366	0	100
2005	365	365	0	100
2006	365	365	0	100
2007	365	365	0	100
2008	365	366	1	99.7
2009	365	365	0	100
2010	365	365	0	100
2011	365	365	0	100
2012	366	366	0	100
2013	365	365	0	100
2014	365	365	0	100
2015	365	365	0	100
2016	366	366	0	100
2017	365	365	0	100
2018	365	365	0	100
2019	365	365	0	100
2020	137	366	229	37.4

4 Output

file.image<-here::here("data","rds",paste0("Backup Meteorological review_",as.Date(now()),".rdata"))

save.image(file = file.image)