Skip to content

Commit

Permalink
improved plot.ds.pca ++
Browse files Browse the repository at this point in the history
  • Loading branch information
brasmus committed Dec 13, 2024
1 parent e4c4e53 commit 5d4ba7a
Show file tree
Hide file tree
Showing 12 changed files with 311 additions and 83 deletions.
4 changes: 2 additions & 2 deletions DESCRIPTION
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
Package: esd
Version: 1.10.91
Date: 2024-11-04
Version: 1.10.92
Date: 2024-12-11
Title: Climate analysis and empirical-statistical downscaling (ESD) package for monthly and daily data
Author: Rasmus E. Benestad, Abdelkader Mezghani, Kajsa M. Parding, Helene B. Erlandsen, Ketil Tunheim, and Cristian Lussana
Maintainer: Rasmus E. Benestad <[email protected]>
Expand Down
4 changes: 3 additions & 1 deletion R/colorscale.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -81,6 +81,7 @@ col.bar <- function(xleft,ybottom,xright,ytop,breaks,horiz=TRUE,
#ii <- (1:n)%%2 == 1
#text(mids[ii],rep(ybottom+dy,n)[ii],round(breaks,2)[ii],cex=cex.axis, col='grey30')
text(mids[ii]+dm, rep(ybottom,n)[ii]+db, round(breaks,2)[ii],cex=cex.axis, col='grey30')
invisible(list(mids=mids,col=col,breaks=breaks))
}

#' @export
Expand All @@ -101,6 +102,7 @@ colbar <- function(breaks,col,fig=c(0.15,0.2,0.15,0.3),horiz=FALSE,
}
par(fig=par0$fig, xaxt=par0$xaxt, yaxt=par0$yaxt,
mar=par0$mar, las=par0$las, cex.axis=par0$cex.axis)
invisible(list(mids=mids,col=col,breaks=breaks))
}

#' Display a color bar object on an existing plot.
Expand Down Expand Up @@ -160,7 +162,7 @@ colbar.ini <- function(x,FUN=NULL,colbar=NULL,verbose=FALSE) {

## Set breaks and n
if (verbose) {print('Set breaks and n:')}
if (!is.null(colbar$col)) {
if (!is.null(colbar$col) & is.null(colbar$breaks)) {
colbar$n <- length(colbar$col)
if (is.null(colbar$breaks)) {
colbar$breaks <- round(seq(x.rng[1],x.rng[2],length.out=length(colbar$col)+1),nd)
Expand Down
4 changes: 2 additions & 2 deletions R/map.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -795,12 +795,12 @@ map.pca <- function(x,...,it=NULL,is=NULL,ip=1,new=FALSE,projection="lonlat",
if (is.element(FUN,args)) {
z <- map.station(X,new=new,colbar=colbar,
xlim=xlim,ylim=ylim,zlim=zlim,
plot=TRUE,#add=add,fig=fig,
plot=plot,#add=add,fig=fig,
verbose=verbose,...)
} else {
z <- map.station(X,new=new,colbar=colbar,FUN=FUN,
xlim=xlim,ylim=ylim,zlim=zlim,
plot=TRUE,#add=add,fig=fig,
plot=plot,#add=add,fig=fig,
verbose=verbose,...)
}
}
Expand Down
77 changes: 56 additions & 21 deletions R/plot.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -124,8 +124,8 @@
#' @param main main title
#' @param xlab label of x-axis
#' @param ylab label of y-axis
#' @param errorbar if TRUE show errorbar
#' @param legend.show if TRUE show legendp
#' @param errorbar if TRUE show error bar
#' @param legend.show if TRUE show legend
#' @param map.show show map of stations
#' @param map.type 'points' to show stations on map, 'rectangle' to show area
#' @param map.insert if TRUE show map as insert, else show map in new window
Expand Down Expand Up @@ -372,7 +372,7 @@ plot.station <- function(x,...,plot.type="single",new=TRUE,
mar <- c(4.5, 4.5, 4.5, 1.0)
}
}

if (is.null(ylim)) {
ylim <- range(pretty(range(x,na.rm=TRUE)))
ylim <- ylim + diff(ylim)*0.1*c(-0.05,0.05)
Expand Down Expand Up @@ -465,9 +465,9 @@ plot.station <- function(x,...,plot.type="single",new=TRUE,

if(legend.show) {
legend("bottomleft",legend=paste(attr(x,'location'),": ",
round(attr(x,'longitude'),2),"E/",
round(attr(x,'latitude'),2),"N (",
attr(x,'altitude')," masl)",sep=""),
round(attr(x,'longitude'),2),"E/",
round(attr(x,'latitude'),2),"N (",
attr(x,'altitude')," masl)",sep=""),
bty="n",cex=0.65,ncol=2,
text.col="grey40",lty=1,col=col)
}
Expand Down Expand Up @@ -770,7 +770,7 @@ plot.eof.comb <- function(x,...,new=FALSE,xlim=NULL,ylim=NULL,
}

#par(cex.axis=0.75,cex.lab=0.7,cex.main=0.8)

if (length(grep('eof',what))>0) {
#if (!is.null(mfrow)) par(fig=c(0,0.5,0.5,1))
map(x,ip=ip,verbose=verbose,colbar=colbar,...)
Expand All @@ -796,13 +796,13 @@ plot.eof.comb <- function(x,...,new=FALSE,xlim=NULL,ylim=NULL,
apps <- anms[grep('appendix',anms)]
n.app <- length(apps)
if (is.null(ylim)) {
ylim <- range(coredata(x[,n]))
for (i in 1:n.app) {
if(verbose) print(apps[i])
z <- attr(x,apps[i])
zz <- try(coredata(z[,n]))
if (!inherits(zz,'try-error')) ylim <- range(c(ylim,zz),na.rm=TRUE)
}
ylim <- range(coredata(x[,n]))
for (i in 1:n.app) {
if(verbose) print(apps[i])
z <- attr(x,apps[i])
zz <- try(coredata(z[,n]))
if (!inherits(zz,'try-error')) ylim <- range(c(ylim,zz),na.rm=TRUE)
}
}

if(verbose) print(xlim)
Expand Down Expand Up @@ -831,7 +831,7 @@ plot.eof.comb <- function(x,...,new=FALSE,xlim=NULL,ylim=NULL,
plot.zoo(x[,ip],lwd=2,ylab=ylab,main=main,xlim=xlim,ylim=ylim,
cex.main=0.8,bty="n",cex.axis=0.9,cex.lab=1,xaxt="n")
taxis <- range(index(x))

## Plot the common PCs
for (i in 1:n.app) {
z <- attr(x,apps[i])
Expand Down Expand Up @@ -1651,8 +1651,10 @@ plot.ds.pca <- function(x,...,ip=1,
if (is.null(colbar2)) colbar2 <- colbar1
attr(y,'longname') <- attr(y,'longname')[1]

if(is.null(what)) what <- c("predictor.pattern","pca.pattern",
"xval","timeseries")
# if(is.null(what)) what <- c("predictor.pattern","pca.pattern",
# "xval","timeseries")
## REB 2024-12-10
if(is.null(what)) what <- c("eof+pca.patterns","xval","timeseries")
if(is.null(attr(y,'evaluation'))) what <- what[!what %in% c("xval")]

i <- 1
Expand All @@ -1669,7 +1671,7 @@ plot.ds.pca <- function(x,...,ip=1,
figlist <- list(fig1=c(0.05,0.45,0.05,0.975),
fig2=c(0.5,0.975,0.05,0.975))
} else figlist <- NULL

if(new) dev.new()
par(mar=mar, mgp=mgp)
if('pca.pattern' %in% what) {
Expand Down Expand Up @@ -1698,6 +1700,39 @@ plot.ds.pca <- function(x,...,ip=1,
par(cex=1) ## REB 2019-08-06
#title(paste("Predictor pattern #",ip,sep=""), cex=0.8)
}

if('eof+pca.patterns' %in% what) {
## REB 2024-12-10: Default three panels, map with sympols on top
mY <- map.pca(y,ip=ip,plot=FALSE,verbose=verbose,...)
par(new=(add | i>1))
if(!is.null(figlist)) par(fig=figlist[[i]],mar=rep(1,4))
#fig=c(0.55,0.975,0.5,0.975),new=TRUE)
i <- i+1
pp <- attr(y,'predictor.pattern')
attr(pp,"variable") <- paste("EOF.Pattern.",ip,sep="")
attr(pp,"unit") <- "unitless"
colbar2$show <- FALSE
colbar2$breaks <- pretty(c(pp[,,ip],-pp[,,ip]),n=15)
if (varid(y)[1]=='precip') colbar2$pal <- 'precip.ipcc' else
colbar2$pal <- 't2m.ipcc'
colbar2$col <- colscal(length(colbar2$breaks)-1,colbar2$pal)
eofmap <- map(pp,ip=ip,new=FALSE,showaxis=FALSE,
colbar=colbar2,verbose=verbose,
lab=paste("Predictor pattern # ",ip,sep=""))
par(cex=1) ## REB 2019-08-06
## Add rain gauge information
col <- rep('grey',length(mY))
breaks <- attr(eofmap,'colbar')$breaks
## re-scale PCA weights to match the colour scale of the map
mY <- mY*IQR(eofmap)/IQR(c(mY))
#mY <- mY*attr(y,'eigenvalues')[ip]
## Add symbols with colours
breaks <- breaks[-1]
for (j in 1:length(mY)) col[j] <- attr(eofmap,'colbar')$col[breaks >= mY[j]][1]
points(lon(mY),lat(mY),pch=19,cex=1.5,col=col)
points(lon(mY),lat(mY),pch=21,cex=1.5,col='black')
}

if('xval' %in% what) {
if (verbose) print('Evaluation results')
par(new=(add | i>1))
Expand Down Expand Up @@ -2010,7 +2045,7 @@ plot.cca <- function(x,...,ip=1,
# paste(attr(x$Y,'source')[1],attr(x$Y,'variable')[1])),
# col=c("red","blue"),lwd=2,lty=1,
# bty="n",cex=0.5,ncol=2,text.col="grey40")

## KMP 2023-02-22: reset if layout has been used
if(panels>1) par(mfrow=c(1,1))
par(mar=par0$mar, mgp=par0$mgp, xaxt=par0$xaxt , yaxt=par0$yaxt)
Expand Down Expand Up @@ -2098,7 +2133,7 @@ plot.dsensemble.pca <- function(x,...,pts=FALSE,target.show=TRUE,map.show=TRUE,
legend.show=TRUE,verbose=FALSE) {
if (verbose) print("plot.dsensemble.pca")
stopifnot(inherits(x,'dsensemble') & inherits(x,'pca'))

d <- index(x[[3]])
pc <- x[3:length(x)]
pc <- array(unlist(pc), dim = c(dim(pc[[1]]), length(pc)))
Expand Down Expand Up @@ -2277,7 +2312,7 @@ plot.dsensemble.one <- function(x,pts=FALSE,it=0,
if (pts) for (i in 1:d[2]) {
points(year(t),coredata(z[,i]),pch=19,cex=0.3,col=col)#"red")
}

# Produce a transparent envelope
nt <- length(index(z))
t2 <- c(year(t),rev(year(t)))
Expand Down
120 changes: 79 additions & 41 deletions R/station.GHCND.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,42 @@
#' nz <- meta.GHCND(lon=c(5,10),lat=c(60,65),alt=500)
#' esd::map(nz)
#'
#' ## Pull all GHCND-data before they may disappear with new US president
#' ## Make one netCDF file for each country
#'
#' if (length(grep('meta',ls()))==0) meta <- meta.GHCND()
#' meta$country <- gsub(' ','',meta$country)
#' cntrs <- rownames(table(meta$country))
#' save(meta,file='meta.GHCND.rda')
#'
#' done <- list.files(pattern='ghcnd.precip.')
#' done <- sub('.nc','',sub('ghcnd.precip.','',done[grep('.nc',done)]))
#' print(done)
#' cntrs <- cntrs[!is.element(cntrs,done)]
#'
#' for (cntr in cntrs) {
#' print(cntr)
#' m <- subset(meta,cntr=cntr)
#' if (dim(m)[1] > 0) {
#' X <- station.GHCND(m,verbose=TRUE)
#' save(X,file=paste('ghcnd.2024',cntr,'rda',sep='.'))
#' vars <- names(X)
#' for (varid in vars) {
#' if (!is.null(dim(X[[varid]])))
#' if (dim(X[[varid]])[2]>2) {
#' print(paste('ghcnd',varid,cntr,'nc',sep='.'))
#' nv <- apply(X[[varid]],2,'nv')
#' x <- subset(X[[varid]],is=nv > 3625)
#' write2ncdf4(x,file=paste('ghcnd',varid,cntr,'nc',sep='.'))
#' } else {
#' print(paste('One sries - ghcnd',varid,cntr,'rda',sep='.'))
#' x <- X[[varid]]
#' save(x,file=paste('ghcnd',varid,cntr,'rda',sep='.'))
#' }
#' }
#' }
#' }

#' @exportS3Method
#' @export station.GHCND
station.GHCND <- function(x=NULL,cntr=NULL,param=NULL,lon=NULL,lat=NULL,
Expand All @@ -77,51 +113,53 @@ station.GHCND <- function(x=NULL,cntr=NULL,param=NULL,lon=NULL,lat=NULL,
ii <- 1
for (file2get in filenames) {
if (verbose) print(sub(paste0(url,'/'),'',file2get))
ghcnd <- read.table(file2get,sep=sep,header=TRUE)
content <- names(ghcnd)
if (length(grep('PRCP',content))>0) {
precip <- zoo(x=ghcnd$PRCP/10,order.by=as.Date(ghcnd$DATE))
precip <- as.station(precip,stid=ghcnd$STATION[1],loc=ghcnd$NAME[1],
ghcnd <- try(read.table(file2get,sep=sep,header=TRUE))
if (!inherits(ghcnd,'try-error')) {
content <- names(ghcnd)
if (length(grep('PRCP',content))>0) {
precip <- zoo(x=ghcnd$PRCP/10,order.by=as.Date(ghcnd$DATE))
precip <- as.station(precip,stid=ghcnd$STATION[1],loc=ghcnd$NAME[1],
lon=ghcnd$LONGITUDE[1],lat=ghcnd$LATITUDE[1],
alt=ghcnd$ELEVATION[1],cntr=x$country[ii],
param='precip',unit='mm',longname='24-hr precipitation',
src='GHCN',url=url)
if (is.null(Precip)) Precip <- precip else Precip <- combine.stations(Precip,precip)
}
if (length(grep('TMAX',content))>0) {
tmax <- zoo(x=ghcnd$TMAX/10,order.by=as.Date(ghcnd$DATE))
tmax <- as.station(tmax,stid=ghcnd$STATION[1],loc=ghcnd$NAME[1],
lon=ghcnd$LONGITUDE[1],lat=ghcnd$LATITUDE[1],
alt=ghcnd$ELEVATION[1],cntr=x$country[ii],
param='precip',unit='mm',longname='24-hr precipitation',
param='tmax',unit='degC',longname='daily maximum temperature',
src='GHCN',url=url)
if (is.null(Precip)) Precip <- precip else Precip <- combine.stations(Precip,precip)
}
if (length(grep('TMAX',content))>0) {
tmax <- zoo(x=ghcnd$TMAX/10,order.by=as.Date(ghcnd$DATE))
tmax <- as.station(tmax,stid=ghcnd$STATION[1],loc=ghcnd$NAME[1],
lon=ghcnd$LONGITUDE[1],lat=ghcnd$LATITUDE[1],
alt=ghcnd$ELEVATION[1],cntr=x$country[ii],
param='tmax',unit='degC',longname='daily maximum temperature',
src='GHCN',url=url)
if (is.null(Tmax)) Tmax <- tmax else Tmax <- combine.stations(Tmax,tmax)
}
if (length(grep('TMIN',content))>0) {
tmin <- zoo(x=ghcnd$TMIN/10,order.by=as.Date(ghcnd$DATE))
tmin <- as.station(tmax/10,stid=ghcnd$STATION[1],loc=ghcnd$NAME[1],
lon=ghcnd$LONGITUDE[1],lat=ghcnd$LATITUDE[1],
alt=ghcnd$ELEVATION[1],cntr=x$country[ii],
param='tmin',unit='degC',longname='daily minimum temperature',
src='GHCN',url=url)
if (is.null(Tmin)) Tmin <- tmin else Tmin <- combine.stations(Tmin,tmin)
}
if (length(grep('TAVG',content))>0) {
t2m <- zoo(x=ghcnd$TAVG/10,order.by=as.Date(ghcnd$DATE))
t2m <- as.station(t2m,stid=ghcnd$STATION[1],loc=ghcnd$NAME[1],
lon=ghcnd$LONGITUDE[1],lat=ghcnd$LATITUDE[1],
alt=ghcnd$ELEVATION[1],cntr=x$country[ii],
param='t2m',unit='degC',longname='daily average temperature',
src='GHCN',url=url)
if (is.null(T2m)) T2m <- t2m else T2m <- combine.stations(T2m,t2m)
if (is.null(Tmax)) Tmax <- tmax else Tmax <- combine.stations(Tmax,tmax)
}
if (length(grep('TMIN',content))>0) {
tmin <- zoo(x=ghcnd$TMIN/10,order.by=as.Date(ghcnd$DATE))
tmin <- as.station(tmax/10,stid=ghcnd$STATION[1],loc=ghcnd$NAME[1],
lon=ghcnd$LONGITUDE[1],lat=ghcnd$LATITUDE[1],
alt=ghcnd$ELEVATION[1],cntr=x$country[ii],
param='tmin',unit='degC',longname='daily minimum temperature',
src='GHCN',url=url)
if (is.null(Tmin)) Tmin <- tmin else Tmin <- combine.stations(Tmin,tmin)
}
if (length(grep('TAVG',content))>0) {
t2m <- zoo(x=ghcnd$TAVG/10,order.by=as.Date(ghcnd$DATE))
t2m <- as.station(t2m,stid=ghcnd$STATION[1],loc=ghcnd$NAME[1],
lon=ghcnd$LONGITUDE[1],lat=ghcnd$LATITUDE[1],
alt=ghcnd$ELEVATION[1],cntr=x$country[ii],
param='t2m',unit='degC',longname='daily average temperature',
src='GHCN',url=url)
if (is.null(T2m)) T2m <- t2m else T2m <- combine.stations(T2m,t2m)
}
ii <- ii + 1
}
ii <- ii + 1
}
if (is.null(param)) result <- list(precip=Precip,tmax=Tmax,tmin=Tmin,t2m=T2m) else
result <- switch(tolower(param),'precip'=Precip,'tmax'=Tmax,'tmin'=Tmin,'t2m'=T2m)
attr(result,'references') <- c('Menne, M.J., I. Durre, R.S. Vose, B.E. Gleason, and T.G. Houston, 2012: An overview of the Global Historical Climatology Network-Daily Database. Journal of Atmospheric and Oceanic Technology, 29, 897-910, doi.10.1175/JTECH-D-11-00103.1.',
'Durre I., M. J. B.E. Gleason, T. G. Houston, and R. S. Vose, 2010: Comprehensive automated quality assurance of daily surface observations. Journal of Applied Meteorology and Climatology., 49, 1615-1633, doi.10.1175/2010JAMC2375.1.',
'Durre, I., M.J. Menne, and R.S. Vose, 2008: Strategies for evaluating quality assurance procedures. Journal of Applied Meteorology and Climatology, 47, 1785–1791, doi: 10.1175/2007JAMC1706.1')
'Durre I., M. J. B.E. Gleason, T. G. Houston, and R. S. Vose, 2010: Comprehensive automated quality assurance of daily surface observations. Journal of Applied Meteorology and Climatology., 49, 1615-1633, doi.10.1175/2010JAMC2375.1.',
'Durre, I., M.J. Menne, and R.S. Vose, 2008: Strategies for evaluating quality assurance procedures. Journal of Applied Meteorology and Climatology, 47, 1785–1791, doi: 10.1175/2007JAMC1706.1')

invisible(result)
}
Expand Down Expand Up @@ -209,8 +247,8 @@ meta.GHCND <- function(url='https://www1.ncdc.noaa.gov/pub/data/ghcn/daily/ghcnd
if (length(alt)>1) sel <- (meta$altitude >= min(alt)) & (meta$altitude <= max(alt)) else
if (alt <0) sel <- (meta$altitude <= abs(alt)) else
sel <- (meta$altitude >= alt)
if (verbose) print(paste(sum(sel),'selected'))
meta <- meta[sel,]
if (verbose) print(paste(sum(sel),'selected'))
meta <- meta[sel,]
}

if (plot) {
Expand All @@ -222,7 +260,7 @@ meta.GHCND <- function(url='https://www1.ncdc.noaa.gov/pub/data/ghcn/daily/ghcnd
if (verbose) print(dim(meta))
class(meta) <- c("stationmeta","data.frame")
attr(meta,'references') <- c('Menne, M.J., I. Durre, R.S. Vose, B.E. Gleason, and T.G. Houston, 2012: An overview of the Global Historical Climatology Network-Daily Database. Journal of Atmospheric and Oceanic Technology, 29, 897-910, doi.10.1175/JTECH-D-11-00103.1.',
'Durre I., M. J. B.E. Gleason, T. G. Houston, and R. S. Vose, 2010: Comprehensive automated quality assurance of daily surface observations. Journal of Applied Meteorology and Climatology., 49, 1615-1633, doi.10.1175/2010JAMC2375.1.',
'Durre, I., M.J. Menne, and R.S. Vose, 2008: Strategies for evaluating quality assurance procedures. Journal of Applied Meteorology and Climatology, 47, 1785–1791, doi: 10.1175/2007JAMC1706.1')
'Durre I., M. J. B.E. Gleason, T. G. Houston, and R. S. Vose, 2010: Comprehensive automated quality assurance of daily surface observations. Journal of Applied Meteorology and Climatology., 49, 1615-1633, doi.10.1175/2010JAMC2375.1.',
'Durre, I., M.J. Menne, and R.S. Vose, 2008: Strategies for evaluating quality assurance procedures. Journal of Applied Meteorology and Climatology, 47, 1785–1791, doi: 10.1175/2007JAMC1706.1')
invisible(meta)
}
16 changes: 16 additions & 0 deletions R/station.YR.R
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,16 @@
station.YR <- function(is,url="https://api.met.no/weatherapi/locationforecast/2.0/",verbose=TRUE) {
# Load jsonlite package
library(jsonlite)

if (inherits(is,'station')) is <- list(lon=lon(is),lat=lat(is))
## If field, then need to repeat so that there is equal lons and lats
# URL
#url <- "https://api.met.no/weatherapi/locationforecast/2.0/compact?lat=60.10&lon=9.58"

URL <- paste0(url,'compact?lat=',lat,'&lon=',lon)
# Read JSON data from the URL
data <- fromJSON(url)

# View the data
View(data)
}
Loading

0 comments on commit 5d4ba7a

Please sign in to comment.