-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathpreprocess_ECMWF.R
105 lines (89 loc) · 4.89 KB
/
preprocess_ECMWF.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
101
102
103
104
105
# Author: Kate Doubleday
# Last updated: 2/20/2020
# ----------------------------------------------------------------
# This script pre-processes ECMWF NETCDF files to extract their 51
# ensemble forecasts of GHI from the SSRD (surface solar radiation
# downwards) ECMWF parameter. Estimates for the 7 SURFRAD sites are
# interpolated from the nearest 4 grid points, given a 0.2x0.2 degree
# lat/lon grid. GHI forecast time resolution is hourly average.
# -----------------------------------------------------------------
# Load dependencies
library(here)
library(ncdf4)
library(pracma)
library(abind)
# -----------------------------------------------------------------
# Define constants
site_names <- c("Bondville",
"Boulder",
"Desert_Rock",
"Fort_Peck",
"Goodwin_Creek",
"Penn_State",
"Sioux_Falls")
site_coord <- list(Bondville=list(N=40.05192, W=88.37309),
Boulder=list(N=40.12498, W=105.23680),
Desert_Rock=list(N=36.62373, W=116.01947),
Fort_Peck=list(N=48.30783, W=105.10170),
Goodwin_Creek=list(N=34.2547, W=89.8729),
Penn_State=list(N=40.72012, W=77.93085),
Sioux_Falls=list(N=43.73403, W=96.62328))
input_directory <- here("ECMWF_files", "Raw")
output_directory <- here("ECMWF_files", "SURFRAD_sites")
dir.create(output_directory, showWarnings = FALSE)
month_days <- c(31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31)
day_indx <- 1:365
# Parameters for new NetCDF file
dims <- mapply(ncdim_def, name=c('Day', 'IssueTime', "LeadTime", 'Member'), units=c('', 'Hour', "Hours", ""),
vals=list(day_indx, seq(0, 18, by=6), 1:24, 1:51),
longname=c("Day number starting Jan 1, 2018", "Forecast time of issue", "Forecast lead time", "Ensemble member"),
SIMPLIFY = FALSE)
pvar <- ncvar_def("irradiance", "W/m^2", dims, missval=NA, compression = 9)
# -----------------------------------------------------------------
#' A helper function to import each daily netCDF, spatially interpolate, and transform to irradiance
#' @param day Day of year index
#' @param site Site name
#' @return an array of forecasts
get_daily_irr <- function(day, site) {
month <- min(which(day <=cumsum(month_days)))
day_of_month <- day - ifelse(month>1, sum(month_days[1:(month-1)]), 0)
xp <- site_coord[[site]]$N
yp <- site_coord[[site]]$W
# Extract closest spatial data from NetCDF file
nc <- nc_open(file.path(input_directory, paste("north_america_", ifelse(month<10, paste(0, month, sep=''), month), ifelse(day_of_month<10, paste(0, day_of_month, sep=''), day_of_month), ".nc", sep='')))
lon_indx <- which(round(nc$dim$longitude$vals, 1)==-ceiling(yp*5)/5)
lat_indx <- which(round(nc$dim$latitude$vals, 1)==ceiling(xp*5)/5)
ssrd_grid <- ncvar_get(nc, varid="ssrd", start=c(lon_indx,lat_indx,1,1,1), count=c(2,2,-1,-1,-1))
nc_close(nc)
nc <- nc_open(file.path(input_directory, paste("north_america_ctrl_", ifelse(month<10, paste(0, month, sep=''), month), ifelse(day_of_month<10, paste(0, day_of_month, sep=''), day_of_month), ".nc", sep='')))
ssrd_ctrl_grid <- ncvar_get(nc, varid="ssrd", start=c(lon_indx,lat_indx,1,1), count=c(2,2,-1,-1))
nc_close(nc)
ssrd_ctrl_grid <- array(ssrd_ctrl_grid, dim=append(dim(ssrd_ctrl_grid), 1, after = 2))
ssrd_grid <- abind(ssrd_ctrl_grid, ssrd_grid, along=3)
# Spatially interpolate to SURFRAD coordinates.
x <- nc$dim$latitude$vals[c(lat_indx+1, lat_indx)] # lat coordinates must be nondecreasing
y <- nc$dim$longitude$vals[c(lon_indx, lon_indx+1)]
# interp2 requires Z to length(y) x length(x) rather than length(x) x length(y).
ssrd <- apply(ssrd_grid, MARGIN = c(3, 4, 5), FUN=function(z) (interp2(x, y, matrix(z[c(3,4,1,2)], ncol=2), xp, -yp, method="linear")))
# Calculate irradiance from ssrd
irr <- apply(ssrd, MARGIN=c(1,3), FUN=diff)/3600
# Clean: truncate small negative values to 0
irr[irr<0] <- 0
return(irr)
}
# -----------------------------------------------------------------
#' A helper function to gather the daily files into a single irradiance array for a given site and export to a new NetCDF
#' @param site Site name
export_site_netcdf <- function(site) {
# Get irradiance array for this site
irr <- sapply(day_indx, FUN=get_daily_irr, site=site, simplify="array")
# Reformat to dimensions: [day]x[issue time]x[forecast horizon]x[members]
irr <- aperm(irr, c(4,3,1,2))
# Export new NetCDF file
nc_new <- nc_create(file.path(output_directory, paste(site, ".nc", sep='')), pvar)
ncvar_put(nc_new, pvar, irr, count=pvar[['varsize']])
nc_close(nc_new)
}
# -----------------------------------------------------------------
# Cycle through the sites, exporting their site-specific irradiance NetCDF's
for (site in site_names) export_site_netcdf(site)