From 23fac335bc71c4f05d5b9e1d7b2c218b1734fc80 Mon Sep 17 00:00:00 2001 From: LaMargna Date: Tue, 11 Jul 2017 11:11:35 +0200 Subject: [PATCH 1/3] Modified generalWrangleRadar.R MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Commented out some code I blieve not yet needed Made projection strings global (for now) Closed open netcdf files when no more needed added a function “makeCellfunction” that as a function of a netcdf file creates a function mapping the Dutch grid (global SOURCE string) to the polar stereographic coordinate in the TARGET string (not yet obtained from the netcdf) --- generalWrangleRadar.R | 107 +++++++++++++++++++++++++++++++++++------- 1 file changed, 90 insertions(+), 17 deletions(-) diff --git a/generalWrangleRadar.R b/generalWrangleRadar.R index d0b35ce..50301b1 100644 --- a/generalWrangleRadar.R +++ b/generalWrangleRadar.R @@ -1,35 +1,54 @@ library(ncdf4) library(data.table) -library(sp) -library(raster) +#library(sp) +#library(raster) library(proj4) +# SOURCE contains the projection string of the external data +# this is currently the Dutch grid. +# SOURCE should eventually be a parameter, if we want international use +# For now, we keep it as a global parameter -#De Bilt 140844, 457994 +SOURCE <- "+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.4171,50.3319,465.5524,-0.398957388243134,0.343987817378283,-1.87740163998045,4.0725 +units=m +no_defs" -pointX<-140844 -pointY<-457994 - -convertToPolar<-function(pointX, pointY){ +# projection string in the polar stereographic projection +# see alsoe http://adaguc.knmi.nl/contents/datasets/productdescriptions/W_ADAGUC_Product_description_RADNL_OPER_R___25PCPRR_L3.html +TARGET <- "+proj=stere +lat_0=90 +lon_0=0.0 +lat_ts=60.0 +a=6378.137 +b=6356.752 +x_0=0 +y_0=0" + -RDH<- "+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.4171,50.3319,465.5524,-0.398957388243134,0.343987817378283,-1.87740163998045,4.0725 +units=m +no_defs" +# for development, we use a typical NetCDF file, containing radar data -pointDF<-data.frame(pointX,pointY) +exampleNetCDFfile <- "./KNMI/PrecipitationRadar/KNMI-Data/RAD_NL25_RAC_MFBS_5min_201501312355.nc" +filepath <- exampleNetCDFfile -coordinates(pointDF)<-~pointX+pointY -crs(pointDF)<-RDH -polar<-"+proj=stere +lat_0=90 +lon_0=0.0 +lat_ts=60.0 +a=6378.137 +b=6356.752 +x_0=0 +y_0=0" +# the location of the KNMI in the Dutch grid +#De Bilt 140844, 457994 -polars<-ptransform(cbind(pointDF$pointX,pointDF$pointY),RDH,polar, silent = F) +pointX<-140844 +pointY<-457994 -polars<-data.frame(polars) -colnames(polars)<-c("X","Y","Z") -#polar is a dataFrame -polars +# the convertToPolar function converts coordinates in the SOURCE projection to +# coordintates in the TARGET projection + +convertToPolar <- function(pointX, pointY) { + + + + # coordinates(pointDF)<-~pointX+pointY + # crs(pointDF)<-SOURCE + + polars <- ptransform(cbind(pointX, pointY, 0), SOURCE, TARGET, silent = FALSE) + + polars <- data.frame(polars) + colnames(polars) <- c("X", "Y", "Z") + + #polars is a dataFrame + polars + } @@ -47,6 +66,7 @@ fileFull<-nc_open(filepath) x<-ncvar_get(fileFull, "x") y<-ncvar_get(fileFull, "y") dataGrid<-ncvar_get(fileFull,variableOfInterest) +nc_close(fileFull) @@ -76,7 +96,60 @@ getMetaDataNetCDF<-function(filepath, variableOfInterest = "image1_image_data"){ scaleFactor <- scaleFactor$value addOffsetFactor <- ncatt_get(fileFull,variableOfInterest,"add_offset") addOffsetFactor <- addOffsetFactor$value + nc_close(fileFull) list(missingValue = missingValue, scaleFactor = scaleFactor, addOffset = addOffsetFactor, variable = variableOfInterest) } + +getMetaDataNetCDF(exampleNetCDFfile) + + +makeCellfunction <- function(filepath) { + + #loading NetCDF file + fileFull<-nc_open(filepath) + + #reading in the variables + x<-ncvar_get(fileFull, "x") + y<-ncvar_get(fileFull, "y") + + #close the NetCDF file + nc_close(fileFull) + + + cellSize <- x[[2]] - x[[1]] #constant size square grid is assumed in this projection + + if (cellSize != abs(y[[1]] - y[[2]])) stop("Grid is not square") # we may drop this test + + # note that the y-values can decrease as we go down from the top line (the first line number is one), + # as we assume y is largest in to of the figure, so y[[1]] - y[[2]] > 0 + # in some cases this may be opposite. + # x values increase if we go to the right + + offsetGrid <- cellSize / 2 #x,y coordinate of the grid are mid points + + xmin <- min(x) + ymax <- max(y) + + function(xvar, yvar) { + + inPolarCoordinates <- convertToPolar(xvar, yvar) + +# inPolarCoordinates$X may be up to "offsetGrid" smaller than xmin to fall in the first column +# inPolarCoordinates$Y may be up to "offsetGrid" larger than ymax to fall in the first row (top row) + + columnX <- ceiling((inPolarCoordinates$X - xmin + offsetGrid) / cellSize) # getting to the NetCDF array cell from the coordinate (X) + + if (y[[1]] < y[[2]]) + rowY<-ceiling((yvar-ymin+offsetGrid)/cellSize) + else + rowY <- ceiling((ymax - inPolarCoordinates$Y + offsetGrid) / cellSize) # and (Y) + list(X = columnX, Y = rowY, xvar = xvar, yvar = yvar) + } +} + +transFormFunction <- makeCellfunction(exampleNetCDFfile) + + +transFormFunction(pointX, pointY) From 118593609c0956025b38264f2af32c7ebf1b4d4b Mon Sep 17 00:00:00 2001 From: LaMargna Date: Tue, 11 Jul 2017 12:51:53 +0200 Subject: [PATCH 2/3] Added ImportPrecipitationRadarData.R This program imports precipitation data from a tree containg NetCDF files. The program assumes all NetCDF files relevant (no other types should be in the path). The program is designed to read from and save to a share pointed to by the environment variable ProjectData --- .../ImportPrecipitationRadarData.R | 237 ++++++++++++++++++ 1 file changed, 237 insertions(+) create mode 100644 KNMI/PrecipitationRadar/ImportPrecipitationRadarData.R diff --git a/KNMI/PrecipitationRadar/ImportPrecipitationRadarData.R b/KNMI/PrecipitationRadar/ImportPrecipitationRadarData.R new file mode 100644 index 0000000..36b175b --- /dev/null +++ b/KNMI/PrecipitationRadar/ImportPrecipitationRadarData.R @@ -0,0 +1,237 @@ +# ImportPrecipitationRadarData.R +# +# The following R script is used to create a key value table tables +# "Time", "X", "Y", "Precipitation" +# from KNMI provided NetCDF files +# http://adaguc.knmi.nl/contents/datasets/productdescriptions/W_ADAGUC_Product_description_RADNL_OPER_R___25PCPRR_L3.html +# +# This program assumes the NetCDF files already to be downloaded (for instance from +# http://opendap.knmi.nl/knmi/thredds/catalog/radarprecipclim/RAD_NL25_RAC_MFBS_5min_NC/catalog.html +# ) and decompressed into the directory pointed to by the environment variable +# "ProjectData" pointing to a share (part of a path at least), and from there in +# directories "KNMI", "PrecipitationRadar", and "NetCDF", for instance ProjectData equals +# "/Volumes/ProjectData", combined it would make /Volumes/ProjectData/KNMI/PrecipitationRadar/NetCDF/ +# If the environment variable is not set, the data is assumed to be decompressed into +# the path ProjectData/KNMI/PrecipitationRadar/ relative to where the program starts. +# The envirionment variable also controls the output location. If not set +# the path is KNMI/PrecipitationRadar/RData/ relative to where the program starts. +# Currently, the zip files offered by KNMI are organized per year and month (and [ignored] hour). +# This structure is assumed by the program: the program first determines a list of all +# directories one level up from "NetCDF", and iterates over each of them. +# Assuming these directories represent data for separate years, the program first checks +# for the presence of a output dataset associated with that year. If it finds one, it +# assumes that year is already processed and skips to the next directory (more about this below). +# If the program assumes the year is not yet processed, it lists all subdirectories one +# level up which it assumes to represent months of the year. In principle, the program +# processes each month in turn, saves the monthly data in a separate file and concatenates +# all per-month data and saves the annual aagregate data as well. If the program finds a +# per-month data file, it loads the data just for concatenation, so in the end, the complete +# aggregate annual data can be saved. This way it is relatively easy to add extra data once +# new data becomes available (could be done more efficiently though). +# If data needs to be replaced, either the output datasets can be deleted, or the function +# check.file.exists below can be changed to return FALSE +# +# Currently the program uses the Parallel package to load the monthly NetCDF files in +# parallel using the forking mechanism, so this might not work on windows machines. + +# It is unclear to me whether the time used in the file name is the beginning of the +# 5-minute interval or the end. May not matter much, but anyway. + +check.file.exists <- function(thePath) file.exists(thePath) + +# uncomment the next line to replace all data +# check.file.exists <- function(thePath) FALSE + + +# +# load required packages +# +suppressPackageStartupMessages(require(dplyr)) +suppressPackageStartupMessages(require(stringr)) +suppressPackageStartupMessages(require(tidyr)) +suppressPackageStartupMessages(require(ncdf4)) +suppressPackageStartupMessages(require(parallel)) + + +# One can use the "ProjectData" environment variable to determine the output location + +if (Sys.getenv("ProjectData") == "") { + outputDirectory <- file.path("KNMI", "PrecipitationRadar", "RData") # could be simplified +} else + outputDirectory <- file.path(Sys.getenv("ProjectData"), + "KNMI", + "PrecipitationRadar", + "RData") + +# +# list the year-directories +# +theDirs <- list.dirs(path = file.path(ifelse(Sys.getenv("ProjectData") == "", + "ProjectData", + Sys.getenv("ProjectData")), + "KNMI", + "PrecipitationRadar", + "NetCDF"), + full.names = TRUE, + recursive = FALSE) + +# debugging code +# theDir <- theDirs[[1]] +# theDir + +# iterate over the years + +for (theDir in theDirs) { + +# Path to potential output file + + pathToYearFile <- file.path(outputDirectory, + paste("PrecipitationRadarData-", basename(theDir), ".Rdata", sep = "")) + +# if file exists, skip to next... + + if (check.file.exists(pathToYearFile)) + next ; + + theMonthDirs <- list.dirs(path = theDir, + full.names = TRUE, + recursive = FALSE) + +# debugging code +# theMonthDir <- theMonthDirs[[1]] +# theMonthDir + + +# create annual RadarData by concatenating all monthly data + + RadarData <- bind_rows( + lapply( + theMonthDirs, + function(theMonthDir) { + +# Path to potential output file + + pathToMonthFile <- file.path(outputDirectory, + paste("PrecipitationRadarData-", basename(theDir), "-", basename(theMonthDir), ".Rdata", sep = "")) + +# if file exists (but not the annual aggregate, which still has to be made), load the data... + + if (file.exists(pathToMonthFile)) + load(pathToMonthFile) + else { + +# if the file does not exist, create it +# first tell the world what happens (and the time we started) + + print(paste(date(), theMonthDir)) + +# make a list of all NetCDF files (assumed *.nc) below the month directory +# this includes per day directories + + theFiles <- list.files(path = theMonthDir, + pattern = glob2rx("*.nc"), + full.names = TRUE, + recursive = TRUE, + no.. = TRUE) + +# debugging code +# theFile <- theFiles[[1]] + + nCores <- detectCores() + my.cluster <- makeForkCluster(nCores) + + RadarDataMonth <- bind_rows( + parLapply( + my.cluster, + theFiles, + function(theFile) { + + timepoint <- as.POSIXct(str_extract_all(theFile, "_[[:digit:]]{8}[[:digit:]]{4}\\.")[[1]], + format = "_%Y%m%d%H%M.", + tz = "UTC") + +# code lended from Andrea's "getMetaDataNetCDF" + + fileH <- nc_open(theFile) + + variableOfInterest <- "image1_image_data" + missingValue <- ncatt_get(fileH, variableOfInterest, "_FillValue") + missingValue <- missingValue$value + scaleFactor <- ncatt_get(fileH, variableOfInterest, "scale_factor") + scaleFactor <- scaleFactor$value + addOffsetFactor <- ncatt_get(fileH, variableOfInterest, "add_offset") + addOffsetFactor <- addOffsetFactor$value + + x <- ncvar_get(fileH, "x") + y <- ncvar_get(fileH, "y") + image_data <- ncvar_get(fileH, variableOfInterest) + # + nc_close(fileH) + + # we obtained a matrix of integers, mostly containing na's + # note that this matrix first dimension is x, and the second is y + # so in terms of i, and j, the y values are columns (j) + + # the next procedure is as follows + # first image_data is converted into a data frame. + # this data.frame has columns V1 to Vn with n the number of columns (this is the Y dimension) + # this data.frame has row.names 1:m with m the number of rows (this is the X dimension) + # we create the data.frame, + # set the X value to the numeric value of the row numbers and + # set the Start value equal to the time point + # + dataFrame <- as.data.frame(image_data) + dataFrame$X <- as.integer(rownames(dataFrame)) + dataFrame$Start <- timepoint + + # next we gather, basically stack, all columns starting with "V", including X and Start + + thisFrameRadarData <- gather(dataFrame, + key = "yLabel", value = "Precipitation", + starts_with("V")) %>% + + # Filter the NA's and the missing values, + + filter(!is.na(Precipitation), + Precipitation > 0, + Precipitation != missingValue) %>% + + # Compute Y as the numerical value of the column name without the "V" + # calculate the actual measurement + + mutate(Y = as.integer(substring(yLabel, 2)), + + # calculate the actual measurement + + Precipitation = scaleFactor * Precipitation + addOffsetFactor) %>% + + # Select the required data + + select(Start, X, Y, Precipitation) + + return(thisFrameRadarData) + })) + + # Refresh the cluster (all leaked data, if any gone) + + stopCluster(my.cluster) + + # save the monthly data + + save(RadarDataMonth, + file = pathToMonthFile) + } + + return(RadarDataMonth) + })) + + # save the annual data + + save(RadarData, + file = pathToYearFile) + +} + +warnings() +sessionInfo() +proc.time() From 74a989b7e6e1a2ac2deefc5be545e8e9addea749 Mon Sep 17 00:00:00 2001 From: LaMargna Date: Tue, 11 Jul 2017 15:29:29 +0200 Subject: [PATCH 3/3] Minor bug fix: one comment deleted --- KNMI/PrecipitationRadar/ImportPrecipitationRadarData.R | 1 - 1 file changed, 1 deletion(-) diff --git a/KNMI/PrecipitationRadar/ImportPrecipitationRadarData.R b/KNMI/PrecipitationRadar/ImportPrecipitationRadarData.R index 36b175b..3f189b3 100644 --- a/KNMI/PrecipitationRadar/ImportPrecipitationRadarData.R +++ b/KNMI/PrecipitationRadar/ImportPrecipitationRadarData.R @@ -197,7 +197,6 @@ for (theDir in theDirs) { Precipitation != missingValue) %>% # Compute Y as the numerical value of the column name without the "V" - # calculate the actual measurement mutate(Y = as.integer(substring(yLabel, 2)),