-
Notifications
You must be signed in to change notification settings - Fork 0
/
Zhang_extraction.R
100 lines (82 loc) · 4.06 KB
/
Zhang_extraction.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
library(tidyr)
library(dplyr)
library(ncdf4)
library(stars)
library(raster)
library(here)
library(ggplot2)
library(ggthemes)
library(stringr)
library(zoo)
library(SPEI)
library(viridis)
rm(list=ls())
OutDir <- "C:/Users/arunyon/3D Objects/Local-files/RCF_Testing/HAVO_data"
tif <- read_stars("D:/HI_Data/Zhang/Processed/HI_Downscaling_Share/DynDS_RainfallTIF_AllFolders/DynDS_6MonthlyTIF_State_250m/rcp85/State_RF_rcp85_2080_01_250m.tif",along=band)
# boundary <- st_read('./data/HALE/HALE_boundary.shp')
boundary <- st_read('C:/Users/arunyon/OneDrive - DOI/Documents/GIS/HAVO_Kilauea_Summit_Wet_Dry_Zones/HAVO_Kilauea_Summit_Wet_Dry_Zones.shp')
boundary <- st_transform(boundary, st_crs(tif))
rm(tif)
# ggplot() + # Resolution is course
# # geom_raster(data = topo_df ,aes(x = x, y = y,alpha=HYP_HR_SR_W_1), show.legend=FALSE) +
# geom_stars(data = Zhang[1,], alpha = 0.8) +
# # facet_wrap("time") +
# # scale_fill_viridis() +
# #coord_equal() +
# geom_sf(data = boundary, aes(), fill = NA) +
# theme_map() +
# theme(legend.position = "bottom",
# legend.key.width = unit(6, "cm"),
# legend.key.height = unit(.3, "cm"),
# legend.justification = "center",
# plot.title=element_text(size=12,face="bold",hjust=0.5)) +
# # plot.background = element_rect(colour = col, fill=NA, size=5)) +
# # labs(fill = "Water Balance") +
# scale_colour_manual(values = c(rgb(207, 31, 46, maxColorValue = 255)), "#ffda85")
## Monthly Rainfall
Monthly_dir_RF<-"D:/HI_Data/Zhang/Processed/HI_Downscaling_Share/DynDS_RainfallTIF_AllFolders/DynDS_6MonthlyTIF_State_250m/"
scens<-c("present","rcp45","rcp85")
df <- data.frame()
#Create list of file in directory
for(i in 1:length(scens)){
file.list = list.files(path = paste0(Monthly_dir_RF,scens[i],"/"), pattern = '.tif', full.names = TRUE)
monthly_RF<-read_stars(file.list,along="band")
monthly_RF<-setNames(monthly_RF,"Rainfall_mm")
date<-as.Date(paste0(gsub("_","-",str_extract(file.list,"\\d{4}[_]\\d{2}")),"-01"),format="%Y-%m-%d")
mystar_RF = st_set_dimensions(monthly_RF, 3, values =date, names = "date")
RF_crop <-st_crop(mystar_RF,boundary)
saveRDS(RF_crop, file = paste0(here::here('data/Output/Data-files//'),'RF.monthly-',scens[i]))
# RF_avg<-st_apply(RF_crop, c("x", "y"), mean) #mean of all time periods
# RF_max <- st_apply(RF_crop, c("x", "y"), max)
# RF_range <- st_apply(RF_crop, c("x", "y"), range) #min and max bands for each pixel
pcp.time <- st_apply((RF_crop %>% dplyr::select(Rainfall_mm)), c("date"),mean,na.rm=TRUE, rename=FALSE)
df1 <- data.frame(pcp.time)
df1$scen=scens[i]
df<-rbind(df,df1)
rm(df1)
}
write.csv(df, file=paste0(here::here('data/Output/Data-files//'),'RF.monthly.csv'))
## Monthly Temperature
Monthly_dir_Tmean<-"D:/HI_Data/Zhang/Processed/HI_Downscaling_Share/DynDS_TemperatureTIF_AllFolders/DynDS_6MonthlyTIF_State_250m/"
df <- data.frame()
#Create list of file in directory
for (i in 1:length(scens)){
file.list = list.files(path = paste0(Monthly_dir_Tmean,scens[i],"/"), pattern = '.tif', full.names = TRUE)
monthly_Tmean<-read_stars(file.list,along="band")
monthly_Tmean<-setNames(monthly_Tmean,"TmeanK")
date<-as.Date(paste0(gsub("_","-",str_extract(file.list,"\\d{4}[_]\\d{2}")),"-01"),format="%Y-%m-%d")
mystar_Tmean = st_set_dimensions(monthly_Tmean, 3, values =date, names = "date")
Tmean_crop <-st_crop(mystar_Tmean,boundary)
Tmean_crop %>% mutate(TmeanF = (9/5*(TmeanK - 273.15)+32)) ->Tmean_crop
Tmean_crop["TmeanF"]->Tmean_crop
saveRDS(Tmean_crop, file = paste0(here::here('data/Output/Data-files//'),'Tmean.monthly-',scens[i]))
# Tmean_avg<-st_apply(Tmean_crop, c("x", "y"), mean) #mean of all time periods
tmean.time <- st_apply((Tmean_crop %>% dplyr::select(TmeanF)), c("date"),mean,na.rm=TRUE, rename=FALSE)
df1 <- data.frame(tmean.time)
df1$scen=scens[i]
df<-rbind(df,df1)
rm(df1)
}
write.csv(df, file=paste0(here::here('data/Output/Data-files//'),'Tmean.monthly.HALE.csv'))
rm(monthly_RF,mystar_RF,monthly_Tmean,mystar_Tmean,pcp.time,RF_avg,RF_crop,RF_max,RF_range,Tmean_crop,tmean.time)
gc()