Commit dee3e76f authored by Burkhardt Rockel's avatar Burkhardt Rockel

Version 2.1 2018/01/15

Removed: *plot_colourbar.R*, *plotmap.old.R*

Replaced: *ncdf_times.R* by *nccf_times.R*

In **colfun.R**

* Fixed errors in *get.ethz.palette* and  *get.ncar.palette*

In **plot_colourbar.R**

* new option *sea.lab* , *border*, *triangle_ends*

In **plotmap.R**

* new options *sea.lab* and *lakes*, changed option name *scale_factor* to *scale.factor*

In **polymask.R**

* added support for multiple polygons

In **windvec.R**

* added new option *legend.scale*
parent 008f7eea
Package: ncdf4Utils
Type: Package
Title: provides plot functions to use with NetCDF data
Version: 2.0
Date: 2017-06-29
Author: Jonas Bhend, with contributions
from Burkhardt Rockel <burkhardt.rockel@hzg.de>
Version: 2.1
Date: 2018-01-16
Title: Provides plot functions to use with NetCDF data
Authors@R: c(
person("Burkhardt", "Rockel", "email=burkhardt.rockel@hzg.de", role = c("aut", "cre")),
person("Ronny", "Petrik", "email=ronny.petrik@hzg.de", role = "aut")
Author: Burkhardt Rockel [aut,cre],
Jonas Bhend [aut],
Ronny Petrik [aut]
Maintainer: Burkhardt Rockel <burkhardt.rockel@hzg.de>
Depends: R (>= 1.8.0), ncdf4, maps, mapdata, sp, akima
Suggests: udunits
Suggests: PCICT
Description: This package contains numerous functions to plot
spatial data on both rotated and regular lonlat grids, as
well as an extension to the functions in ncdf allowing for
NetCDF times to be read into R Date objects.
spatial data on rotated and other grids following the CF conventions.
License: GPL
LazyLoad: yes
LazyData: true
Packaged: Thu Jun 29 12:45:52 CEST 2017; rockel
Packaged: 2018-01-16 16:09:23 UTC; rockel
......@@ -161,12 +161,13 @@ function(palette.name, url="https://wiki.c2sm.ethz.ch/pub/Data/VisNCLCosmoLibrar
{
file<-paste(url,palette.name,".ct",sep="")
urlcolors<-readLines(file)
n.levels<-as.numeric(urlcolors[3])
#levels<-as.numeric(strsplit(x = urlcolors[4], split = " ", fixed = T)[[1]][])
levels<-as.numeric(scan(text=urlcolors[4], what="",quiet=TRUE))
urlcolors<-urlcolors[!startsWith(urlcolors,';')]
n.lines<-length(urlcolors)
n.levels<-as.numeric(urlcolors[1])
levels<-as.numeric(scan(text=urlcolors[2], what="",quiet=TRUE))
colors<-c(1:(n.levels+1))
for (i in 1:(n.levels+1)) {
colors.rgb <- as.numeric(scan(text=urlcolors[i+4], what="",quiet=TRUE))
colors.rgb <- as.numeric(scan(text=urlcolors[i+2], what="",quiet=TRUE))
colors[i] <- rgb(colors.rgb[1],colors.rgb[2],colors.rgb[3], maxColorValue=255)
}
out<-list(colors=colors, levels=levels)
......@@ -176,13 +177,13 @@ get.ncar.palette <-
function(palette.name, url="http://www.ncl.ucar.edu/Document/Graphics/ColorTables/Files/")
{
file<-paste(url,palette.name,".rgb",sep="")
file<-paste(url,palette.name,".rgb",sep="")
urlcolors<-readLines(file)
urlcolors<-urlcolors[urlcolors!="" & !startsWith(urlcolors,'#')]
n.colors<-as.numeric(strsplit(x = urlcolors[1], split = "=", fixed = T)[[1]][2])
colors<-c(1:n.colors)
for (i in 1:(n.colors)) {
colors.rgb <- as.numeric(scan(text=urlcolors[i+2], what="",quiet=TRUE))
# print(c(i,colors.rgb))
colors.rgb <- as.numeric(scan(text=urlcolors[i+1], what="",quiet=TRUE))
colors[i] <- rgb(colors.rgb[1],colors.rgb[2],colors.rgb[3], maxColorValue=255)
}
out<-list(colors=colors)
......
geo2multirot <- function(loncent, latcent, lon, lat, false_easting=0, false_northing=0, earth_radius=NULL){
if (is.null(earth_radius)){
rearth <- 2.E7/pi
} else {
rearth <- earth_radius
}
zpir18 <- pi/180
zrpi18 <- 180/pi
latpol <- asin( cos(latcent*zpir18) ) * zrpi18
if (latcent > 0) {
lonpol <- loncent -180
longn <- 0
} else {
lonpol <- loncent
longn <- 180
}
tmp <- geo2rot(lonpol, latpol, lon, lat, longn)
lonh <- zpir18*tmp$x
lath <- zpir18*tmp$y
zarg1 <- -cos(lonh)
zarg2 = tan(lath)
zarg2[zarg2 == 0] <- 1.0e-20
ypol <- abs(atan(zarg1/zarg2))
yr = pi/2. - ypol
ind <- lath < 0 & !is.na(lath)
yr[ind] <- -yr[ind]
zarg1 = - sin(lonh) *cos(lath)
zarg2[ind] = -sin(ypol[ind])*cos(lath[ind])*cos(lonh[ind]) + cos(ypol[ind])*sin(lath[ind])
ind <- lath >= 0 & !is.na(lath)
zarg2[ind] = -sin(ypol[ind])*cos(lath[ind])*cos(lonh[ind]) - cos(ypol[ind])*sin(lath[ind])
zarg2[zarg2 == 0] <- 1.0e-20
xr = abs(atan(zarg1/zarg2))
ind <- lonh < 0 & !is.na(lonh)
xr[ind] <- -xr[ind]
#... radians on earths surface into Cartesian distances
xc = xr * rearth + false_easting
yc = yr * rearth + false_northing
out <- list(x=xc, y=yc)
out
}
##########################################################
multirot2geo <- function(loncent, latcent, xc, yc, false_easting=0, false_northing=0, earth_radius=NULL){
if (is.null(earth_radius)){
rearth <- 2.E7/pi
} else {
rearth <- earth_radius
}
zpir18 <- pi/180
zrpi18 <- 180/pi
xr <- (xc - false_easting) / rearth
yr <- (yc - false_northing) / rearth
lath <- asin ( cos( xr ) * sin( yr ) ) *zrpi18
lonh <- atan2( tan( xr ) , cos( yr ) ) *zrpi18
latpol <- asin( cos(latcent*zpir18) ) *zrpi18
if (latcent > 0) {
lonpol <- loncent - 180
longn <- 0
} else {
lonpol <- loncent
longn <- 180
}
out <- rot2geo(lonpol, latpol, lonh, lath, longn)
}
`nccf_times` <-
function(nc) {
# this function converts netcdf times in a netCDF file
# with CF conventions to the
# R date-time format
timei <- which(names(nc$dim) %in% c("time", "TIME", "tim", "TIM"))
units <- nc$dim[[timei]]$units
vals <- nc$dim[[timei]]$vals
tmp <- ncatt_get(nc, names(nc$dim)[timei], "calendar")
if (tmp$hasatt) {
calendar <- trimws(tmp$value)
} else {
calendar <- "standard"
print(paste("Warning: Calendar is missing in", nc$filename))
print(paste(' "standard" calendar is used'))
}
# if (calendar == "proleptic_gregorian") calendar <- "standard"
#... valid calendar names
cals <- c("standard", "gregorian", "proleptic_gregorian", "360_day", "360", "365_day", "365", "noleap")
itype.cals <- c(0, 0, 0, 1, 1, 2, 2, 2)
if (identical(calendar %in% cals, FALSE)){
print(paste("Warning: ",calendar," is not a valid calendar name" ))
print(paste("valid names are: ", cals))
print(paste(' "standard" calendar is used'))
calendar <- "standard"
}
units.dt <- c("seconds since", "minutes since", "hours since","days since","months since","years since")
mul.dt <- c(1, 60, 3600, 86400, 1, 1)
#... find the index of matching units.td in units
index.dt <- match(1,sapply (units.dt,function(x) charmatch(x, units, nomatch=0)),nomatch=0)
if (index.dt == 0) {
print(paste("ERROR: ",units,"is no valid unit"))
print(paste("accepted units are:"))
print(units.dt)
stop
# quit()
}
ref.time <- trimws(gsub(units.dt[index.dt],"",units))
ref.time.pcict <- as.PCICt(ref.time, calendar)
if (index.dt == 5) {
ymd<-format(as.Date(ref.time), "%Y-%m-01")
ref.time.pcict.month<-as.PCICt(ymd,calendar)
y<-sapply (vals,function(x) seq(ref.time.pcict.month, by = "month", length = x+1))
z<-lapply(seq(1,length(y),1), function(x) tail(y[[x]],-(length(y[[x]])-1)))
zz<-do.call(c,z)
zzz<-as.PCICt(zz, calendar) - ref.time.pcict.month
vals<-do.call(c,lapply(seq(1,length(zzz),1),function(x) zzz[[x]]))
}
if (index.dt == 6) {
ymd<-format(as.Date(ref.time), "%Y-01-01")
ref.time.pcict.year<-as.PCICt(ymd,calendar)
y<-sapply (vals,function(x) seq(ref.time.pcict.year, by = "year", length = x+1))
z<-lapply(seq(1,length(y),1), function(x) tail(y[[x]],-(length(y[[x]])-1)))
zz<-do.call(c,z)
zzz<-as.PCICt(zz, calendar) - ref.time.pcict.year
vals<-do.call(c,lapply(seq(1,length(zzz),1),function(x) zzz[[x]]))
}
mul<-mul.dt[index.dt]
times <- vals * mul
time <- ref.time.pcict + times
time
}
`ncdf_times` <-
function(nc, as.Rdate=TRUE, force=TRUE, tz="UTC") {
# this function converts netcdf times to the
# R date-time format or to the udunits dates
# you can choose to switch to uduints format
# for gregorian calendars by setting as.Rdate
# to FALSE
# non-gregorian calendar dates are output using
# udunits date format
# you can force to get back an R date format, even
# if the calendar used is not gregorian using
# force=T (may return udunits dates if conversion
# is not successful)
# WARNING: time zones are not fully supported yet
# check whether udunits is available
.udunitsInclude <- FALSE
if (any(.packages() == "udunits") & class(try(utInit(), silent=T)) != "try-error"){
.udunitsInclude <- TRUE
}
timei <- which(names(nc$dim) %in% c("time", "TIME", "tim", "TIM"))
units <- nc$dim[[timei]]$units
vals <- nc$dim[[timei]]$vals
tmp <- ncatt_get(nc, names(nc$dim)[timei], "calendar")
if (tmp$hasatt) {
calendar <- tmp$value
} else {
calendar <- "standard"
print(paste("Warning: Calendar is missing in", nc$filename))
}
if (calendar == "proleptic_gregorian") calendar <- "standard"
if (as.Rdate & calendar == "standard"){
if (charmatch("hours since", units, nomatch=0) |
charmatch("minutes since", units, nomatch=0) |
charmatch("seconds since", units, nomatch=0)) {
mul <- 1
ref.txt <- substr(units, 15,33)
if (charmatch("minutes", units, nomatch=0)) mul <- 60
if (charmatch("hours", units, nomatch=0)) {
mul <- 3600
ref.txt <- substr(units, 13,31)
}
times <- vals * mul
# if (nchar(ref.txt) == 19){
# ref <- as.POSIXct(ref.txt, tz)
# } else {
# ref <- as.POSIXct(paste(ref.txt, "00", sep=":"), tz)
# }
# time <- as.Date(ref + times)
ref <- strptime(ref.txt, format = "%Y-%m-%d %H:%M:%S", tz = tz)
time <- ref + times
}
if (charmatch("days since", units, nomatch=0)){
time <- as.Date(substr(units, 12,21), "%Y-%m-%d") + vals
}
if (charmatch("months since", units, nomatch=0)) {
ref.yr <- as.numeric(substr(units, 14,17))
ref.month <- as.numeric(substr(units, 19,20))
ref.day <- as.numeric(substr(units, 22,23))
if (is.null(ref.day)) ref.day <- 1
month <- floor((vals+ref.yr*12 + ref.month-1) %% 12) + 1
year <- floor((vals+ref.yr*12 + ref.month-1)/12)
time <- as.Date(ISOdate(year, month, ref.day))
}
if (charmatch("years since", units, nomatch=0)) {
unit.tmp <- paste(strsplit(units, " ")[[1]][3:4], collapse=" ")
#ref.yr <- substr(units, 13,16)
#ref.month <- as.numeric(substr(units, 18,19))
ref.yr <- as.numeric(format(as.Date(unit.tmp), "%Y"))
ref.month <- as.numeric(format(as.Date(unit.tmp), "%m"))
if (is.null(ref.month)) ref.month <- 1
#ref.day <- as.numeric(substr(units, 21,22))
ref.day <- as.numeric(format(as.Date(unit.tmp), "%d"))
if (is.null(ref.day)) ref.day <- 1
year <- floor(vals)
month <- floor((vals*12)%%12)
time <- as.Date(ISOdate(ref.yr + year, ref.month + month, ref.day))
}
if (charmatch("day as", units, nomatch=0)) {
date <- floor(vals)
day <- as.numeric(substr(date, nchar(date)-1, nchar(date)))
if (all(day > 28)) date <- as.character(as.numeric(date) - max(day, na.rm=T) + 28)
date <- paste("000",date, sep="")
date <- substr(date, nchar(date)-7, nchar(date))
time <- as.Date(date, "%Y%m%d")
}
} else {
if (.udunitsInclude){
time <- utCalendar(vals, units, calendar=calendar, style="array")
if (force){
tmp <- try(ISOdatetime(time$year, time$month, time$day, time$hour,
time$minute, time$second, tz), silent=T)
if (class(tmp)[1] != "try-error") time <- tmp
}
} else {
stop("Package udunits cannot be loaded or initialized via utInit()")
}
}
time
}
......@@ -5,12 +5,17 @@ function (x, ...){
`plot_colourbar.default` <-
function(levs, cols, side=1, ylab="", labels=NULL,
xlab="", nn=1, center=F, cex.axis=1, sea.col=NULL,...){
xlab="", nn=1, center=F, cex.axis=1, sea.col=NULL, sea.lab=NULL, border=par("fg"), triangle_ends=NULL, ...){
# plots a colour bar given the levs and cols
# centers the axis labelling instead of at the
# boundary when center==TRUE
# if sea.col is set, the colourbar is prolonged by one
# colour with the uncentered label sea underneath
# Check triangle_ends
if (!is.null(triangle_ends) && (!is.logical(triangle_ends) || length(triangle_ends) != 2)) {
stop("Parameter 'triangle_ends' must be a logical vector with two elements.")
}
las<-c(1,3,1,3)
......@@ -21,19 +26,115 @@ function(levs, cols, side=1, ylab="", labels=NULL,
}
ncols <- length(cols)
lev.arr <- array(1:ncols, c(1, ncols))
triangles.allowed<-FALSE
if (sum(triangle_ends) > 0 && !center && sea.add==0) triangles.allowed<-TRUE
#... colorbar at bottom side=1 or top side=3
if (side %in% c(1,3)){
lev.arr <- t(lev.arr)
}
if (side %in% c(1,3)){
image(1:ncols, 1, lev.arr, axes=F, breaks=1:(ncols+1)-0.5, col=cols,
if (triangles.allowed) {
#... triangles at the ends
image(1:ncols, 1, lev.arr, axes=F, breaks=1:(ncols+1)-0.5, col=rep(NA,ncols),
ylab=NA, xlab=NA, ...)
usr <- par( "usr")
triangle_height <- (usr[2]-usr[1])/ncols
left_triangle <- list(x = c(usr[1], usr[1]+triangle_height, usr[1]+triangle_height) ,
y = c(usr[3]+(usr[4]-usr[3])/2, usr[4], usr[3]))
right_triangle <- list(x = c(usr[2], usr[2]-triangle_height, usr[2]-triangle_height) ,
y = c(usr[3]+(usr[4]-usr[3])/2, usr[4], usr[3]))
if (sum(triangle_ends) == 2) {
#... triangles at both ends
col.start <- 2
col.end <- ncols-1
lab.start<-2
lab.end<-ncols
image(1:ncols, 1, lev.arr, axes=F, add=TRUE, breaks=1:(ncols+1)-0.5, col=c(NA,cols[col.start:col.end],NA),
ylab=ylab, xlab=xlab, ...)
abline(v=col.start:col.end + 0.5, col=border)
polygon(left_triangle$x, left_triangle$y,col=cols[1],border=border)
polygon(right_triangle$x, left_triangle$y,col=cols[ncols],border=border)
} else if (triangle_ends[1]) {
#... triangle at the left end
col.start <- 2
col.end <- ncols
lab.start<-2
lab.end<-ncols+1
image(1:ncols, 1, lev.arr, axes=F, add=TRUE, breaks=1:(ncols+1)-0.5, col=c(NA,cols[col.start:col.end]),
ylab=ylab, xlab=xlab, ...)
abline(v=col.start:col.end + 0.5, col=border)
polygon(left_triangle$x, left_triangle$y,col=cols[1],border=border)
} else {
#... triangle at the right end
col.start <- 1
col.end <- ncols-1
lab.start<-1
lab.end<-ncols
image(1:ncols, 1, lev.arr, axes=F, add=TRUE, breaks=1:(ncols+1)-0.5, col=c(cols[col.start:col.end],NA),
ylab=ylab, xlab=xlab, ...)
abline(v=col.start:col.end + 0.5, col=border)
polygon(right_triangle$x, left_triangle$y,col=cols[1],border=border)
}
} else {
#... no triangles at the ends
image(1:ncols, 1, lev.arr, axes=F, breaks=1:(ncols+1)-0.5, col=cols,
ylab=ylab, xlab=xlab, ...)
abline(v=1:(ncols-1) + 0.5)
if (!is.na(border)) abline(v=1:(ncols-1) + 0.5, col=border)
}
} else {
image(1, 1:ncols, lev.arr, axes=F, breaks=1:(ncols+1)-0.5, col=cols,
ylab=ylab, xlab=xlab, ...)
abline(h=1:(ncols-1) + 0.5)
#... colorbar left side=2 or right side=4
if (triangles.allowed) {
#... triangles at the ends
image(1, 1:ncols, lev.arr, axes=F, breaks=1:(ncols+1)-0.5, col=rep(NA,ncols),
ylab=NA, xlab=NA, ...)
usr <- par( "usr")
triangle_height <- (usr[4]-usr[3])/ncols
bottom_triangle <- list(x = c(usr[1]+(usr[2]-usr[1])/2, usr[2],usr[1]) ,
y = c(usr[3], usr[3]+triangle_height, usr[3]+triangle_height))
top_triangle <- list(x = c(usr[1]+(usr[2]-usr[1])/2, usr[2],usr[1]) ,
y = c(usr[4], usr[4]-triangle_height, usr[4]-triangle_height))
if (sum(triangle_ends) == 2) {
#... triangles at both ends
col.start <- 2
col.end <- ncols-1
lab.start<-2
lab.end<-ncols
image(1, 1:ncols, lev.arr, axes=F, add=TRUE, breaks=1:(ncols+1)-0.5, col=c(NA,cols[col.start:col.end],NA),
ylab=ylab, xlab=xlab, ...)
abline(h=col.start:col.end + 0.5, col=border)
polygon(bottom_triangle$x, bottom_triangle$y,col=cols[1],border=border)
polygon(top_triangle$x, top_triangle$y,col=cols[ncols],border=border)
} else if (triangle_ends[1]) {
#... triangle at the bottom end
col.start <- 2
col.end <- ncols
lab.start<-2
lab.end<-ncols+1
image(1, 1:ncols, lev.arr, axes=F, add=TRUE, breaks=1:(ncols+1)-0.5, col=c(NA,cols[col.start:col.end]),
ylab=ylab, xlab=xlab, ...)
abline(h=col.start:col.end + 0.5, col=border)
polygon(bottom_triangle$x, bottom_triangle$y,col=cols[1],border=border)
} else {
#... triangle at the right end
col.start <- 1
col.end <- ncols-1
lab.start<-1
lab.end<-ncols
image(1, 1:ncols, lev.arr, axes=F, add=TRUE, breaks=1:(ncols+1)-0.5, col=c(cols[col.start:col.end],NA),
ylab=ylab, xlab=xlab, ...)
abline(h=col.start:col.end + 0.5, col=border)
polygon(top_triangle$x, top_triangle$y,col=cols[1],border=border)
}
} else {
#... no triangles at the ends
image(1, 1:ncols, lev.arr, axes=F, breaks=1:(ncols+1)-0.5, col=cols,
ylab=ylab, xlab=xlab, ...)
if (!is.na(border)) abline(h=1:(ncols-1) + 0.5, col=border)
}
}
if (center){
at.lev <- seq(1+sea.add,ncols,nn)
if (is.null(labels)){
......@@ -41,21 +142,49 @@ function(levs, cols, side=1, ylab="", labels=NULL,
}
axis(side, at=at.lev, labels=labels, las=1, tick=F, cex.axis=cex.axis)
} else {
if (triangles.allowed) {
at.lev <- seq(lab.start,lab.end,nn)
if (is.null(labels)){
labels <- levs[at.lev]
}
} else {
at.lev <- seq(2+sea.add,ncols,nn)
if (is.null(labels)){
labels <- levs[at.lev - sea.add]
}
axis(side, at=at.lev-0.5,labels=labels, las=1, cex.axis=cex.axis)
}
# if (!is.na(border)) {
# axis(side, at=at.lev-0.5,labels=labels, las=1, cex.axis=cex.axis)
# } else {
# axis(side, at=at.lev-0.5,labels=labels, las=1, col = NA, col.ticks = 1, cex.axis=cex.axis)
axis(side, at=at.lev-0.5,labels=labels, las=1, col = border, col.ticks = 1, cex.axis=cex.axis)
# }
}
if (!is.null(sea.col)){
axis(side, at=1, labels="water", las=las[side], tick=F, cex.axis=cex.axis)
axis(side, at=1, labels=sea.lab, las=las[side], tick=F, cex.axis=cex.axis)
}
box()
}
if (triangles.allowed) {
mai.bak <- par("mai")
mai<-mai.bak
pin <- par("pin")
if (side %in% c(1,3)){
dx<-pin[1]/ncols
if (triangle_ends[1]) mai[2]<-mai.bak[2]+dx
if (triangle_ends[2]) mai[4]<-mai.bak[4]+dx
} else {
dy<-pin[2]/ncols
if (triangle_ends[1]) mai[1]<-mai.bak[1]+dy
if (triangle_ends[2]) mai[3]<-mai.bak[3]+dy
}
par(mai=mai)
}
if (!is.na(border)) box(col=border)
if (triangles.allowed) par(mai=mai.bak)
}
`plot_colourbar.plotmap` <-
function(x, incl.units=T, side=1, cex.axis=1, center=F, labels=NULL, ...){
function(x, incl.units=T, side=1, cex.axis=1, labels=NULL, ...){
if (!is.null(x$flag_values)){
center <- TRUE
labels <- x$flag_values
......@@ -65,8 +194,8 @@ function(x, incl.units=T, side=1, cex.axis=1, center=F, labels=NULL, ...){
}
las<-c(1,3,1,3)
plot_colourbar(x$lev, x$col, sea.col=x$sea.col, side=side,
cex.axis=cex.axis, center=center, labels=labels, ...)
plot_colourbar(x$lev, x$col, sea.col=x$sea.col, sea.lab=x$sea.lab, side=side,
cex.axis=cex.axis, labels=labels, ...)
# plot unit-string in axis
if (incl.units){
......
plotmap <-
function(file, sponge=8,
varname = NULL, lsm.file=NULL,
units=NULL, longname=NULL, add.offset=0, scale_factor=1,
col=NULL, levels=NULL, sea.col=NULL, na.col=NULL, rivers=T,
units=NULL, longname=NULL, add.offset=0, scale.factor=1,
col=NULL, levels=NULL, sea.col=NULL, sea.lab="water", na.col=NULL, rivers=T, lakes=T,
cities=T, label=TRUE, minpop=NULL, ncities=10, city.pch=19,
alt.contour=F, alt.lev=NULL,
grid=TRUE, grid.txt=c(T,T,T,T), grid.lty=2, grid.lwd=1,
......@@ -128,7 +128,7 @@ function(file, sponge=8,
lsm.name = NA
}
data <- ncvar_get(nc, varname) * scale_factor + add.offset
data <- ncvar_get(nc, varname) * scale.factor + add.offset
if (length(dim(data)) == 3){
if (length(nc$var[[varname]]$size) == 3){
if (nc$var[[varname]]$dim[[3]]$name == 'time') {
......@@ -287,6 +287,12 @@ function(file, sponge=8,
}
if (lakes){
lak.dat <- map('lakes', plot=F)
lakes.rot <- geo2mod(lon.ext, lat.ext, lak.dat$x, lak.dat$y ,gridmapping)
}
data.tmp <- data
data.tmp[!lsm] <- NA
if (smoothed) {
......@@ -327,6 +333,11 @@ function(file, sponge=8,
if (is.null(riv.col)) riv.col <- "white"
lines(rivers.rot$x, rivers.rot$y, col=riv.col)
}
if (lakes){
lak.col <- "black"
# if (is.null(lak.col)) lak.col <- "white"
lines(lakes.rot$x, lakes.rot$y, col=lak.col)
}
if