Skip to content

GitLab

Commit d0a9cc9b authored by Franssen, Wietse's avatar Franssen, Wietse
Browse files

flip y partly solved

parent e0d9994d
Hide whitespace changes
Inline Side-by-side
R/WFRTools.R
View file @ d0a9cc9b
......@@ -82,20 +82,20 @@ r2ggplotDataStructure <-function(data, timestep = 1) {
for(x in 1:length(data$xyCoords$x))
{
#for(x in 1:dim(data$xyCoords$x))
for(y in 1:length(data$xyCoords$y))
{
for(y in 1:length(data$xyCoords$y))
{
lonList[i]=data$xyCoords$x[x];
latList[i]=data$xyCoords$y[y];
i=i+1
}
}
ggplotData<-NULL
ggplotData$lon=lonList
ggplotData$lat=latList
if (length(dim(data$Data))==2) {
ggplotData$data=c(data$Data)
ggplotData$data=c(data$Data)
} else if (length(dim(data$Data))==3) {
ggplotData$data=c(data$Data[timestep,,])
} else if (length(dim(data$Data))==4) {
......@@ -103,7 +103,7 @@ r2ggplotDataStructure <-function(data, timestep = 1) {
}
#ggplotFinal<-data.frame(ggplotData)
ggplotFinal<-na.omit(data.frame(ggplotData))
return(ggplotFinal)
}
......@@ -121,20 +121,20 @@ r2ggplotDataStructure <-function(data, timestep = 1) {
#' @export
ggplotjeold <-function(rdata,lon = NULL, lat = NULL,title=NULL ) {
if (is.null(title)) {
if (!is.na(rdata$Dates$start[1])) {
# title=paste0(rdata$Variable$longName, " [",rdata$Variable$units,"]\n(" ,rdata$Dates$start[1], ")")
title=paste0(rdata$Variable$longName, " [",rdata$Variable$units,"]\n(" ,as.character(format( as.Date(rdata$Dates$start[1]),format="%d %B %Y")), ")")
} else {
title=paste0(rdata$Variable$longName, " [",rdata$Variable$units,"]")
}
}
data<-r2ggplotDataStructure(rdata)
myPalette <- colorRampPalette(rev(brewer.pal(11, "Spectral")), space="Lab")
# myPalette <- colorRampPalette(rev(rainbow(100)), space="Lab")
zp1 <- ggplot(data,
aes(x = lon, y = lat, fill = data)) +
ggtitle(title) + # plot title
......@@ -144,17 +144,17 @@ ggplotjeold <-function(rdata,lon = NULL, lat = NULL,title=NULL ) {
#zp1 <- zp1 + scale_x_discrete(expand = c(0, 0))
#zp1 <- zp1 + scale_x_discrete(limits=c(ggdata$lon[1], ggdata$lon[length(ggdata$lon)]))
#zp1 <- zp1 + scale_y_discrete(limits=c(ggdata$lat[1], ggdata$lat[length(ggdata$lat)]))
#zp1 <- zp1 + coord_cartesian(xlim = c(-20, 40),ylim = c(ggdata$lat[1], ggdata$lat[length(ggdata$lat)]))
#zp1 <- zp1 + coord_cartesian(xlim = c(-20, 40),ylim = c(ggdata$lat[1], ggdata$lat[length(ggdata$lat)]))
#zp1 <- zp1 + ylim(ggdata$lat[1], ggdata$lat[length(ggdata$lat)])
coord_equal() +
theme_bw() +
xlab(expression(Longitude * ' ' * degree * E * ' ')) +
ylab(expression(Latitude * ' ' * degree * N * ' ')) +
theme(legend.title=element_blank())
# xlab(expression(Longitude * ' ' * degree * E * ' ')) +
# ylab(expression(Latitude * ' ' * degree * S * ' ')) +
if (!is.null(lon) && !is.null(lat)) {
zp1 <- zp1 + scale_shape_identity() + geom_point( colour="black", size = 10, aes_string(x = lon, y = lat, shape=3))
zp1 <- zp1 + geom_point( fill="red",colour="black", size = 6, aes_string(x = lon, y = lat, shape=18))
......@@ -164,10 +164,10 @@ ggplotjeold <-function(rdata,lon = NULL, lat = NULL,title=NULL ) {
}
ggplotje <-function(rdata,lon = NULL, lat = NULL,title=NULL, barLimits=NULL, barDiff=FALSE, backgroudColor=NULL ) {
if (is.null(title)) {
if (!is.na(rdata$Dates$start[1])) {
# title=paste0(rdata$Variable$longName, " [",rdata$Variable$units,"]\n(" ,rdata$Dates$start[1], ")")
title=paste0(rdata$Variable$longName, " [",rdata$Variable$units,"]\n(" ,as.character(format( as.Date(rdata$Dates$start[1]),format="%d %B %Y")), ")")
} else {
title=paste0(rdata$Variable$longName, " [",rdata$Variable$units,"]")
......@@ -177,11 +177,11 @@ ggplotje <-function(rdata,lon = NULL, lat = NULL,title=NULL, barLimits=NULL, bar
rdata$Data[rdata$Data <= barLimits[1]]<- barLimits[1]
rdata$Data[rdata$Data > barLimits[2]]<-barLimits[2]
}
data<-r2ggplotDataStructure(rdata)
myPalette <- colorRampPalette(rev(brewer.pal(11, "Spectral")), space="Lab")
# myPalette <- colorRampPalette(rev(rainbow(100)), space="Lab")
zp1 <- ggplot(data,
aes(x = lon, y = lat, fill = data)) +
ggtitle(title) + # plot title
......@@ -189,20 +189,20 @@ ggplotje <-function(rdata,lon = NULL, lat = NULL,title=NULL, barLimits=NULL, bar
#zp1 <- zp1 + scale_x_discrete(expand = c(0, 0))
#zp1 <- zp1 + scale_x_discrete(limits=c(ggdata$lon[1], ggdata$lon[length(ggdata$lon)]))
#zp1 <- zp1 + scale_y_discrete(limits=c(ggdata$lat[1], ggdata$lat[length(ggdata$lat)]))
#zp1 <- zp1 + coord_cartesian(xlim = c(-20, 40),ylim = c(ggdata$lat[1], ggdata$lat[length(ggdata$lat)]))
#zp1 <- zp1 + coord_cartesian(xlim = c(-20, 40),ylim = c(ggdata$lat[1], ggdata$lat[length(ggdata$lat)]))
#zp1 <- zp1 + ylim(ggdata$lat[1], ggdata$lat[length(ggdata$lat)])
coord_equal() +
theme_bw() +
xlab(expression(Longitude * ' ' * degree * E * ' ')) +
ylab(expression(Latitude * ' ' * degree * N * ' ')) +
theme(legend.title=element_blank())
# xlab(expression(Longitude * ' ' * degree * E * ' ')) +
# ylab(expression(Latitude * ' ' * degree * S * ' ')) +
if (!is.null(backgroudColor)) {
zp1 <- zp1 + theme(panel.background=element_rect(fill=backgroudColor))
}
library(scales)
if (is.null(barLimits)) {
if (barDiff==TRUE) {
......@@ -217,7 +217,7 @@ ggplotje <-function(rdata,lon = NULL, lat = NULL,title=NULL, barLimits=NULL, bar
zp1 <- zp1 + scale_fill_gradientn(colours = myPalette(100),limits=barLimits)
}
}
if (!is.null(lon) && !is.null(lat)) {
zp1 <- zp1 + scale_shape_identity() + geom_point( colour="black", size = 10, aes_string(x = lon, y = lat, shape=3))
zp1 <- zp1 + geom_point( fill="red",colour="black", size = 6, aes_string(x = lon, y = lat, shape=18))
......@@ -230,7 +230,7 @@ countryMap <- function(lon,lat){
dum = map('worldHires',xlim = c(lon[1],lon[dim(lon)]), ylim = c(lat[1],lat[dim(lat)]), plot = FALSE)
names(dum) <- c("lon", "lat", "range", "names")
countries = data.frame(dum[c("lon","lat")])
#countries$value = mean(d$value)
#countries$value = mean(d$value)
countries$value = countries$lon
#countries$value = NA
return(countries)
......@@ -244,24 +244,24 @@ ncCheck <-function(ncFile, variable) {
result$dims<-FALSE
result$dimids<-FALSE
result$tArray<-FALSE
lons <-ncFile$dim$lon$vals
lats <-ncFile$dim$lat$vals
times <-ncFile$dim$time$vals
# what is the order of the variable dimensions:
for (i in 1:length(ncFile$var[[variable]]$dimids)) {
result$dims[i]<-ncFile$dim[[i]]$name
result$dimids[i]<-ncFile$dim[[i]]$id
}
if (length(times) > 1) { result$tArray<-TRUE }
if ((lons[2]-lons[1]) < 0) { result$xFlip<-TRUE }
if ((lats[2]-lats[1]) < 0) { result$yFlip<-TRUE }
return(result)
return(result)
}
# @examples
# @examples
# data <- ncLoad( file = paste0(find.package("WFRTools"),"./examples/data/example.nc4"))
# data <- ncLoad( file = "./examples/data/example.nc4")
# data <- ncLoad( file = "~/Desktop/gg/wfd_pr_1974.nc", varName = "pr")
......@@ -279,13 +279,13 @@ ncCheck <-function(ncFile, variable) {
#' @keywords internal
#' @export
ncLoad <-function(file, varName = NULL, lonlatbox = NULL, timesteps = NULL, z = NULL, convertUnit=NULL) {
## first build an empty data structure
data<-rDataStructure()
## Open the netcdf file
ncFile <- nc_open( file )
## if no Variable is given then use the first one in the file
if (is.null(varName)) {
varName <- names(ncFile$var)[1]
......@@ -299,10 +299,10 @@ ncLoad <-function(file, varName = NULL, lonlatbox = NULL, timesteps = NULL, z =
cat(paste0("loading the variable: \"", names(ncFile$var[1]),"\"\n"))
}
}
## Do some sheck on the NetCDF file
ncCheckResult<-ncCheck(ncFile = ncFile, variable = varName)
## SOME CHECKS:
## Check time
if (!is.null(timesteps)) {
......@@ -319,7 +319,7 @@ ncLoad <-function(file, varName = NULL, lonlatbox = NULL, timesteps = NULL, z =
}
}
}
## Check lon and lat
if (!is.null(lonlatbox)) {
minLon <-ncFile$dim$lon$vals[1]
......@@ -332,20 +332,21 @@ ncLoad <-function(file, varName = NULL, lonlatbox = NULL, timesteps = NULL, z =
minLat <-ncFile$dim$lat$vals[1]
maxLat <-ncFile$dim$lat$vals[ncFile$dim$lat$len]
if(lonlatbox[3] < minLat || lonlatbox[4] > maxLat) {
stop(paste0("Latitudes out of range!\n",
"Given: ", lonlatbox[3], " till ", lonlatbox[4], "\n",
"Allowed: ", minLat, " till ", maxLat), call. = FALSE)
##WFF
# stop(paste0("Latitudes out of range!\n",
# "Given: ", lonlatbox[3], " till ", lonlatbox[4], "\n",
# "Allowed: ", minLat, " till ", maxLat), call. = FALSE)
}
}
## Fill lon and lat indexes
if (is.null(lonlatbox)) {
if (is.null(lonlatbox)) {
LonIdx <- c(1: ncFile$dim$lon$len)
LatIdx <- c(1: ncFile$dim$lat$len)
} else {
LonIdx <- which( ncFile$dim$lon$vals >= lonlatbox[1] & ncFile$dim$lon$vals <= lonlatbox[2])
LatIdx <- which( ncFile$dim$lat$vals >= lonlatbox[3] & ncFile$dim$lat$vals <= lonlatbox[4])
}
## Fill time indexes
if (is.null(timesteps)) {
timeIdx<-c(1)
......@@ -354,11 +355,19 @@ ncLoad <-function(file, varName = NULL, lonlatbox = NULL, timesteps = NULL, z =
} else {
timeIdx<-timesteps
}
## Fill lon and lat
data$xyCoords$x <- ncFile$dim$lon$vals[LonIdx]
data$xyCoords$y <- ncFile$dim$lat$vals[LatIdx]
## FLIP DATA
if (ncCheckResult$yFlip) {
data$xyCoords$y <- rev(data$xyCoords$y)
}
if (ncCheckResult$xFlip) {
data$xyCoords$x <- rev(data$xyCoords$x)
}
## Fill time
if(!is.null(ncFile$dim$time)) {
NCtime <- ncvar_get( ncFile, "time" ) [timeIdx]
......@@ -373,27 +382,35 @@ ncLoad <-function(file, varName = NULL, lonlatbox = NULL, timesteps = NULL, z =
data$Dates$end <- format(firstTime + (86400 * (NCtime+1)), format="%Y-%m-%d %T %Z")
}
}
## Fill data
if(is.null(z)) {
if(!is.null(ncFile$dim$time)) {
data$Data <- ncvar_get( ncFile, varName, start=c(LonIdx[1],LatIdx[1],timeIdx[1]),count=c(length(LonIdx),length(LatIdx),length(timeIdx)))
} else {
data$Data <- ncvar_get( ncFile, varName, start=c(LonIdx[1],LatIdx[1]),count=c(length(LonIdx),length(LatIdx)))
}
}
} else {
data$Data <- ncvar_get( ncFile, varName, start=c(LonIdx[1],LatIdx[1],z,timeIdx[1]),count=c(length(LonIdx),length(LatIdx),1,length(timeIdx)))
}
##WFF
if (length(dim(data$Data))==2) {
data$Data<-aperm(data$Data, c(2,1))
if(ncCheckResult$yFlip) data$Data<-apply(data$Data,2,rev)
attr(data$Data,"dimensions") <- rev(ncCheckResult$dims[1:2])
} else if (length(dim(data$Data))==3) {
data$Data<-aperm(data$Data, c(3,2,1))
if(ncCheckResult$yFlip) {
for (iTime in 1:dim(data$Data)[1]) {
data$Data[iTime,,]<-apply(data$Data[iTime,,],2,rev)
}
}
attr(data$Data,"dimensions") <- rev(ncCheckResult$dims)
}
#attr(data$Data,"dimensions") <- rev(ncCheckResult$dims)
data$Variable$varName <- varName
## Fill attributes
attTmp<-ncatt_get( ncFile, varName, "long_name" )
if (attTmp$hasatt == TRUE) {
......@@ -407,21 +424,21 @@ ncLoad <-function(file, varName = NULL, lonlatbox = NULL, timesteps = NULL, z =
} else {
data$Variable$units <- "missing"
}
if(!is.null(convertUnit)) {
data<- convertUnit(data)
}
## Close the file
nc_close(ncFile)
## Print some info
if (.PRINT_INFO == TRUE) {
cat("Variables:\n")
print(names(ncFile$var))
print(ncCheckResult)
}
return(data)
}
......@@ -435,10 +452,10 @@ ncLoad <-function(file, varName = NULL, lonlatbox = NULL, timesteps = NULL, z =
#' @keywords internal
#' @export
makeMask <-function(data) {
## first build an empty data structure
mask<-rDataStructure()
## lon and lats
mask$xyCoords$x<-data$xyCoords$x
mask$xyCoords$y<-data$xyCoords$y
......@@ -446,7 +463,7 @@ makeMask <-function(data) {
if (length(dim(data$Data)) > 2 ) {
mask$Data<-data$Data[,,1]
} else {
mask$Data<-data$Data
mask$Data<-data$Data
}
mask$Data[!is.na(mask$Data)] <- 1
naCells<-sum(is.na( mask$Data ))
......@@ -456,7 +473,7 @@ makeMask <-function(data) {
cat(paste0("valid values: ", naValid, "\n" ) )
cat(paste0("NA values: ", naCells, "\n" ) )
cat(paste0("Total: ", total, "\n" ) )
return(mask)
return(mask)
}
#' ncWrite
......@@ -477,21 +494,21 @@ ncWrite <-function(outFile = "~/out.nc", rData = data) {
dimT <- ncdim_def("time", paste0("days since ",timeString[1]), timeArray, unlim = FALSE)
dimX <- ncdim_def("lon", "degrees_east", round(rData$xyCoords$x,2))
dimY <- ncdim_def("lat", "degrees_north",round(rData$xyCoords$y,2))
FillValue <- NA
FillValue <- 1e20
if (!is.null(rData$Variable$units)) {
datavar <- ncvar_def( rData$Variable$varName, rData$Variable$units, list(dimT,dimY,dimX), FillValue, prec="float")
} else {
datavar <- ncvar_def( rData$Variable$varName, " ", list(dimT,dimY,dimX), FillValue, prec="float")
}
ncid_out <- nc_create(outFile, datavar)
## ADD ATTRIBUTES
ncatt_put( ncid_out, 0, "institution", "Wageningen University and Research centre (WUR)")
ncatt_put( ncid_out, 0, "contact", "Wietse Franssen <wietse.franssen@wur.nl>")
ncatt_put( ncid_out, "lon", "standard_name", "longitude")
ncatt_put( ncid_out, "lon", "long_name", "Longitude")
ncatt_put( ncid_out, "lon", "axis", "X")
......@@ -502,29 +519,29 @@ ncWrite <-function(outFile = "~/out.nc", rData = data) {
ncatt_put( ncid_out, "time", "calendar", "standard")
ncatt_put( ncid_out, rData$Variable$varName, "long_name", rData$Variable$longName)
ncatt_put( ncid_out, rData$Variable$varName, "_FillValue", FillValue)
ncvar_put( ncid_out, datavar, rData$Data )
nc_close(ncid_out)
}
# rm(list=ls())
# library(ncdf4)
# library(fields) # e.g: using the fields library
#
#
# plotje <-function(plottitle) {
# image.plot(NClon,NClat,data, asp = 1, main = plottitle, xlab = '', ylab = '')
# world(add = TRUE)
# }
#
#
# domainName<-c( "GHA", "EU")
# lonmin<- c( 27.75, -24.75)
# lonmax<- c( 49.25, 39.75)
# latmin<- c(-12.25, 33.25)
# latmax<- c( 18.25, 71.75)
#
#
# iDomain<-1
#
#
# ## READ WFD NETCDF
# ncFile <- nc_open( "~/Desktop/gg/wfd_pr_1974.nc" )
# LonIdx <- which( ncFile$dim$lon$vals > lonmin[iDomain] | ncFile$dim$lon$vals < lonmax[iDomain])
......@@ -532,18 +549,18 @@ ncWrite <-function(outFile = "~/out.nc", rData = data) {
# data <- ncvar_get( ncFile, "pr")[ LonIdx, LatIdx, 1]
# landmask<-data
# landmask[!is.na(landmask)] <- 1
# sum( !is.na( landmask ) )
# sum( !is.na( landmask ) )
# nc_close(ncFile)
#
#
# ## READ WFDEI NETCDF
# ncFile <- nc_open( "~/Desktop/gg/wfd_pr_1979.nc" )
# LonIdx <- which( ncFile$dim$lon$vals > lonmin[iDomain] | ncFile$dim$lon$vals < lonmax[iDomain])
# LatIdx <- which( ncFile$dim$lat$vals > latmin[iDomain] & ncFile$dim$lat$vals < latmax[iDomain])
# data <- ncvar_get( ncFile, "pr")[ LonIdx, LatIdx, 1]
# landmask[is.na(data)] <- NA
# sum( !is.na( landmask ) )
# sum( !is.na( landmask ) )
# nc_close(ncFile)
#
#
# ## READ SOIL NETCDF
# ncFile <- nc_open( "~/Desktop/gg/soil_GHA.nc" )
# LonIdx <- which( ncFile$dim$lon$vals > lonmin[iDomain] | ncFile$dim$lon$vals < lonmax[iDomain])
......@@ -551,18 +568,18 @@ ncWrite <-function(outFile = "~/out.nc", rData = data) {
# data <- ncvar_get( ncFile, "stexture")[ LonIdx, LatIdx]
# landmask<-data
# landmask[!is.na(landmask)] <- 1
#
#
# landmask[is.na(data)] <- NA
# sum( !is.na( landmask ) )
# sum( !is.na( landmask ) )
# nc_close(ncFile)
#
#
# ## READ SOIL NETCDF
# nameFileNCin<-"~/Desktop/gg/soil_GHA.nc"
# ncid_in=nc_open(nameFileNCin)
# NCdata <- ncvar_get( ncid_in, "stexture")
# nc_close(ncid_in)
#
#
# ## Plot and add to LandMask
# data<-NCdata[,];plotje(plottitle = "WFDEI")
# landmask[is.na(data)] <- NA
# sum( !is.na( landmask ) )
# sum( !is.na( landmask ) )
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment