Skip to content

Instantly share code, notes, and snippets.

@dickoa
Last active February 20, 2017 05:33
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 5 You must be signed in to fork a gist
  • Save dickoa/6078487 to your computer and use it in GitHub Desktop.
Save dickoa/6078487 to your computer and use it in GitHub Desktop.
A R version of the Python script for splines interpolation kindly provided by StormMiner. The code is a straight conversion thus not optimized but it works. Have fun,
require(ncdf4)
loadMesonetData <- function(filename, stationFilename = "station_info.csv") {
data <- read.csv(file.path(data_path, filename))
station_data <- read.csv(file.path(data_path, stationFilename))
list(data = data[,-1], dates = data[,1], station_data = station_data)
}
loadData <- function(filename) {
nc_open(filename)
}
getGrid <- function(data, date, fHour, eMember) {
dateIdx <- which(ncvar_get(data, varid = "intTime") == date)
fIdx <- which(ncvar_get(data, varid = "fhour") == fHour)
eIdx <- which(ncvar_get(data, varid = "ens") == eMember)
ncvar_get(data)[,,fIdx, eIdx, dateIdx]
}
getDailyMeanSumGrid <- function(data, date) {
dateIdx <- which(ncvar_get(data, varid = "intTime") == date)
res <- ncvar_get(data)[, , , , dateIdx]
res <- apply(res, c(1:2, 4), sum)
apply(res, 1:2, mean) * 3600 * 3
}
buildSplines <- function(data, grid, station_data) {
outdata <- vector(mode = "numeric",
length = nrow(station_data))
cat(length(outdata), "\n")
for ( i in seq_along(outdata) ) {
slat <- station_data$nlat[i]
slon <- station_data$elon[i]
nearlat <- which(abs(ncvar_get(data, 'lat') - slat) < 2)
nearlon <- which(abs(ncvar_get(data, 'lon') - slon - 360) < 2)
spline1 <- vector(mode = "numeric",
length = length(nearlon))
for (l in seq_along(nearlat))
spline1[l] <- Spline(grid[nearlon, nearlat[l]],
(slon - floor(slon)))
outdata[i] <- Spline(spline1, (slat - floor(slat)))
}
outdata
}
Spline <- function(y, xi) {
0.5*((2 * y[2]) + (y[3]-y[2]) * xi +
(-y[4] + 4 * y[3] - 5 * y[2] + 2 * y[1]) * xi^2 +
(y[4] - 3 * y[3] + 3 * y[2] - y[1]) * xi^3)
}
data_path <- "./"
data <- loadData(file.path(data_path, "dswrf_sfc_latlon_subset_20080101_20121130.nc")
listsubmit <- loadMesonetData(file.path(data_path, "sampleSubmission.csv"))
dates <- listsubmit$dates
header <- c("Date", names(listsubmit$data))
outdata <- lapply(dates, function(date) {
cat(date, "\n")
grid <- getDailyMeanSumGrid(data, date * 100)
buildSplines(data, grid, listsubmit$station_data)
})
outdata <- do.call(rbind, outdata)
outdata <- cbind(dates, outdata)
colnames(outdata) <- header
write.csv(outdata, "submission.csv", row.names = FALSE)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment