-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathExtract Precipitation NetCDF - Tuolumne.Rmd
111 lines (73 loc) · 3.78 KB
/
Extract Precipitation NetCDF - Tuolumne.Rmd
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
101
102
103
104
105
106
107
108
109
110
111
##This script was written by Vicky Espinoza
##For use by Aditya Sood for CERC-WET
#this script extracts rainfall data from a single grid point (in this case the grid cell closest to 37.9491N, 119.7739W associated with Hetch Hethcy)
#prints csv file with date and water year, month
```{r}
library(raster)
library(ncdf4)
library(maptools)
library(foreign)
library(RNetCDF)
library(rgdal)
library(lubridate)
library(readxl)
#rm(list=ls(all=TRUE)) #start with empty workspace
```
```{r}
#extract_single_grid <- function(path, nc_file, output_path){
precipitation_path <- "C:/Users/gusta/Box/VICE Lab/RESEARCH/PROJECTS/CERC-WET/Task7_San_Joaquin_Model/pywr_models/precipitation"
GCMs <- c("ACCESS1-0", "CCSM4", "CESM1-BGC", "CMCC-CMS",
"GFDL-CM3", "HadGEM2-CC")
rcps <- c("rcp45", "rcp85")
# output_path <- "../../Box/CERC-WET/Task7_San_Joaquin_Model/pywr_models/data/Tuolumne River/hydrology/gcms/MIROC5_rcp85/precipitation"
# Loop through all GCMs for this specific function
for(gcm in GCMs){
for(rcp in rcps){
for(year in 2006:2099){
gcm_path <- paste(precipitation_path,"/",gcm,"/",rcp,"/","rainfall.",
as.character(year),".v0.CA_NV.nc", sep="")
nc_files <- list.files(path = gcm_path)
output_path <- paste0("C:/Users/gusta/Box/VICE Lab/RESEARCH/PROJECTS/CERC-WET/Task7_San_Joaquin_Model/pywr_models/data/Tuolumne River/hydrology/gcms/", gcm, "_", rcp,"/precipitation/")
#if the directory doesn't exist, make it!
if (!dir.exists(output_path)){
dir.create(file.path(output_path), recursive = TRUE)
}
# for(nc_file in nc_files){
# extract_single_grid(gcm_path, nc_file, output_path)
# date.start<- '2006-01-01'
#date.end<-'2099-12-31'
obsdata <- nc_open(paste(gcm_path, sep = "/"))
print(obsdata) # check that dims are lon-lat-time
# location of interest
lon <- 37.9491 # longitude of a
lat <- -119.7739 # latitude of location
# get values at location lonlat
obsoutput <- ncvar_get(obsdata, varid = 'rainfall',
start= c(which.min(abs(obsdata$dim$Lon$vals - lon)), # look for closest long
which.min(abs(obsdata$dim$Lat$vals - lat)), # look for closest lat
1),
count = c(1,1,-1)) #count '-1' means 'all values along that dimension'that dimension'
# create dataframe
#datafinal <- data.frame(dates= obsdatadates, obs = obsoutput)
datafinal <- data.frame(date = seq(from = as.Date(paste(year,"-01-01", sep="")), to = as.Date(paste(year,"-12-31", sep="")), by = 'day'), precip = obsoutput) %>%
rename(`precip (mm)` = precip)
# get dates
#obsdatadates <- as.Date(obsdata$dim$Time$vals, origin = '1950-01-01')
# d <- seq(as.Date(date.start), as.Date(date.end),1)
# datafinal$dates <- d
# datafinal$month <-month(d)
# datafinal$year <-year(d)
# datafinal$day <-day(d)
# datafinal$WY <- datafinal$year + (datafinal$month %in% 10:12) #this method is called vectorization
write.table(datafinal, #vector we want to save
file= paste0(output_path, "/precipitation_Hetch_Hetchy_mm.csv.csv"), #save csv files per basin
append=TRUE, #this gathers all the loop's results, if FALSE we'll have results being overwritten over and over again in the first line
sep=",",
col.names=!file.exists(paste0(output_path,"/precipitation_Hetch_Hetchy_mm.csv.csv")), #if we set as TRUE, we'll have headings repeated per each row, in this way we just have one heading
row.names=FALSE, #no names for rows
quote = FALSE)
#write.csv(datafinal, file = paste0(output_path,"/precipitation_Hetch_Hetchy_mm.csv.csv"))
}
}
}
```