projmap <- function(x,y,z,rotpol=c(0.,90.), lonlim,latlim, projection="",parameters=NULL,orientation=NULL, mapdat="world",xmaplim,ymaplim,thin.map=0, col=heat.colors(12),breaks,na.col="white", longrds=seq(-180,180,by=20),longr="yes", latgrds=seq(-80,+80,by=20),latgr="yes", ltygrds=3,lwdgrds=1,con=FALSE,...) { # Generate a pixel plot of matrix z on a specified map projection. # This routine deals with sterepgraphic and lambert projections from the mapproj # package. # # Description: # # Returns a pixel plot of matrix z # # Usage: # # projmap(x,y,z,rotpol=c(0.,90.),lonlim,latlim,projection="",parameters=NULL,orientation=NULL, # mapdat="world",xmaplim,ymaplim,thin.map=0,col=heat.colors(12),breaks,na.col="white", # longrds=seq(-180,180,by=20),longr="no",latgrds=seq(-80,+80,by=20),latgr="no", # ltygrds=3,lwdgrds=1,con=T,...) # # Arguments: # # x: A vector of longitudes # # y: A vector of latitudes # # z: A matrix # # rotpol: defines the coordinates of the rotated pole. Default is # geographical north pole (rotpol=c(0.,90.)). # # lonlim: defines the longitude window to be plotted. If # specified the window of the projection is chosen such that the entire # longitude window is comprised in the plot. # # latlim: defines the latitude window to be plotted. If # specified the window of the projection is chosen such that the entire # longitude window is comprised in the plot. # # projection: either lambert or stereographic provided by mapproject(). If # an empty string is used, the projection used in the last application # of mapproject will be used. The last projection settings are stored # in variable .Last.projection # # orientation: is the orientation parameter of mapproject. # # mapdat: is the name of the dataset from which map data is taken. # Default is "world". # # xmaplim: defines a longitude window from which map data shall # be plotted. In default, all map data comprised in the area of x,y is # used for plotting. In some cases this does not comprise all the map # data of the projection window. # # ymaplim: define a latitude window from which map data shall # be plotted. In default, all map data comprised in the area of x,y is # used for plotting. # # thin.map: whether to thin out details of the map (political, continental # outlines). thin.map = 1 means that every second point on the map is # omitted. thin.map = 2 means that this thinning out is applied twice, etc. # By default thin.map = 0, which means no thinning. # # col: breaks the color palette and breaks to be used for the pixels. # length of breaks needs to be one larger than length of colors. If # breaks is not specified equidistant breaks from min to max are used. # # na.col: the color to use for pixels with NA (in image). # Useful values are "white" and "grey80". # # longrds: array of longitude circles for which grid-lines shall be plotted. # # longr: string indicating if the longitude circles should be plotted. Default is "no". # # latgrds: array of latitude circles for which grid-lines shall be plotted. # # longr: string indicating if the latitude circles should be plotted. Default is "no". # # ltygrds: the type of line for longitude and latitude circles. # # lwdgrds: the line thickness for longitude and latitude circles. # # Details: # # Plots a pixel plot with either the lambert or stereographical projection. # This script is used in the function plotmap.r. # # Value: # # A pixel plot of the chosen area and projection. # # See Also: # # `plotmap' # # Author: # # Dag Johan Steinskog 3 November 2005 # Christoph Frei # # Examples: # # lon<-seq(-180,180,2.5) # lat<-seq(60,90,2.5) # data<-matrix(rnorm(length(lon)*length(lat),sd=10),ncol=length(lat),nrow=length(lon)) # projmap(lon,lat,data,projection="stereographic", xmaplim=c(-180,180), ymaplim=c(60,90),longrds=seq(-180,180,by=60),mapdat="world",latgrds=seq(60,80,by=10)) # # lon<-seq(-20,40,2.5) # lat<-seq(50,70,2.5) # data<-matrix(rnorm(length(lon)*length(lat),sd=10),ncol=length(lat),nrow=length(lon)) # projmap(lon,lat,data,projection="lambert",lonlim=c(-20,40),latlim=c(50,70),xmaplim=c(-20,40),ymaplim=c(49,70),longr="no",latgr="no",param=c(-20,60)) # rotated coordinates or not? if (rotpol[2] != 90.) {rot <- TRUE} else {rot <- FALSE} # define generally useful variables nx <- length(x) ny <- length(y) plon <- rotpol[1] plat <- rotpol[2] dx <- diff(x)[1] dy <- diff(y)[1] rgb.miss <- col2rgb(na.col) col.miss <- rgb(rgb.miss[1],rgb.miss[2],rgb.miss[3], maxColorValue = 255) # convert positions/values into single structured arrays xmat <- matrix(rep(x,ny),nrow=nx) ymat <- t(matrix(rep(y,nx),nrow=ny)) xx <- array(xmat) yy <- array(ymat) zz <- array(z) rm(xmat,ymat) # define a lon/lat frame to define the window limits and the # details of the projection # case 1 : if lonlim and latlim are specified if (!missing(lonlim) & !missing(latlim)) { lon <- seq(lonlim[1],lonlim[2],length=nx) lat <- seq(latlim[1],latlim[2],length=ny) fram <- matrix(rep(0,2*2*(nx+ny)),ncol=2) fram[,1] <- c(lon,rep(lon[nx],ny),rev(lon),rep(lon[1],ny)) fram[,2] <- c(rep(lat[1],nx),lat,rep(lat[ny],nx),rev(lat)) } else { # case 2 : frame from field lon <- x; lat <- y fram <- matrix(rep(0,2*2*(nx+ny)),ncol=2) fram[,1] <- c(lon,rep(lon[nx],ny),rev(lon),rep(lon[1],ny)) fram[,2] <- c(rep(lat[1],nx),lat,rep(lat[ny],nx),rev(lat)) if (rot) { fram <- rot2lonlat(fram,pollam=plon,polphi=plat) } } # define the details of the projection and the window limits in # plotting coordinates. all subsequent calls # to mapproject will use the same setting. ll <- mapproject(x=fram[,1],y=fram[,2],projection=projection, parameters=parameters,orientation=orientation) xxlim <- range(ll$x,finite=TRUE); yylim <- range(ll$y,finite=TRUE) rm(ll,fram) # calculate transformed pixel corner positions # a) calculate corner points xxll <- xx-dx/2 xxlr <- xx+dx/2 xxur <- xx+dx/2 xxul <- xx-dx/2 yyll <- yy-dy/2 yylr <- yy-dy/2 yyur <- yy+dy/2 yyul <- yy+dy/2 # b) rotate coordinates if necessary if (rot) { # centers fram <- matrix(rep(0,2*nx*ny),ncol=2) fram[,1] <- xx; fram[,2] <- yy ll <- rot2lonlat(fram,pollam=plon,polphi=plat) xx <- ll[,1]; yy <- ll[,2] # lower left corner fram <- matrix(rep(0,2*nx*ny),ncol=2) fram[,1] <- xxll; fram[,2] <- yyll ll <- rot2lonlat(fram,pollam=plon,polphi=plat) xxll <- ll[,1]; yyll <- ll[,2] # lower right corner fram <- matrix(rep(0,2*nx*ny),ncol=2) fram[,1] <- xxlr; fram[,2] <- yylr ll <- rot2lonlat(fram,pollam=plon,polphi=plat) xxlr <- ll[,1]; yylr <- ll[,2] # upper right corner fram <- matrix(rep(0,2*nx*ny),ncol=2) fram[,1] <- xxur; fram[,2] <- yyur ll <- rot2lonlat(fram,pollam=plon,polphi=plat) xxur <- ll[,1]; yyur <- ll[,2] # upper left corner fram <- matrix(rep(0,2*nx*ny),ncol=2) fram[,1] <- xxul; fram[,2] <- yyul ll <- rot2lonlat(fram,pollam=plon,polphi=plat) xxul <- ll[,1]; yyul <- ll[,2] # clean up rm(fram,ll) } # c) convert to the desired projection # corners xxx <- c(xxll,xxlr,xxur,xxul) yyy <- c(yyll,yylr,yyur,yyul) nnn <- length(xxll) ll <- mapproject(x=xxx,y=yyy) xxll <- ll$x[1:nnn]; yyll <- ll$y[1:nnn] xxlr <- ll$x[(nnn+1):(2*nnn)]; yylr <- ll$y[(nnn+1):(2*nnn)] xxur <- ll$x[(2*nnn+1):(3*nnn)]; yyur <- ll$y[(2*nnn+1):(3*nnn)] xxul <- ll$x[(3*nnn+1):(4*nnn)]; yyul <- ll$y[(3*nnn+1):(4*nnn)] rm(ll,xxx,yyy,nnn) # determine the map window in lon/lat if (!missing(xmaplim) & !missing(ymaplim)) { # case 1 : if xmaplim and ymaplim are specified xlim.map <- xmaplim; ylim.map <- ymaplim } else { # case 2 : general lon <- x; lat <- y fram <- matrix(rep(0,2*2*(nx+ny)),ncol=2) fram[,1] <- c(lon,rep(lon[nx],ny),rev(lon),rep(lon[1],ny)) fram[,2] <- c(rep(lat[1],nx),lat,rep(lat[ny],nx),rev(lat)) if (rot) { fram <- rot2lonlat(fram,pollam=plon,polphi=plat) } xlim.map <- range(fram[,1],finite=TRUE) ylim.map <- range(fram[,2],finite=TRUE) if (!missing(lonlim) & !missing(latlim)) { # if lonlim and latlim are specified xlim.map <- range(c(xlim.map,lonlim),finite=TRUE) ylim.map <- range(c(ylim.map,latlim),finite=TRUE) } } # increments for plotting the grids bylongrd <- diff(xlim.map)/200 bylatgrd <- diff(ylim.map)/200 # get the map coordinates in the appropriate projection if (mapdat != "none") { if (mapdat == "continents") { # continents must be read from external file. Are not available #Êin maps package, data was extracted from package clim.pact load("~/work/R/data/continents.Rdata") mapcors <- mapproject(x=continents) } else { mapcors <- map(mapdat,xlim=xlim.map, ylim=ylim.map,plot=FALSE,interior=FALSE) if(projection == "stereographic" & ylim.map[1]>0) { outside <- mapcors$y < (ylim.map[1]-1) & !is.na(mapcors$y) mapcors$x <- mapcors$x[-which(outside)] mapcors$y <- mapcors$y[-which(outside)] } if(projection == "stereographic" & ylim.map[1]<0) { outside <- mapcors$y > (ylim.map[2]-1) & !is.na(mapcors$y) mapcors$x <- mapcors$x[-which(outside)] mapcors$y <- mapcors$y[-which(outside)] } mapcors <- mapproject(x=mapcors$x,y=mapcors$y) } # thin the map if required if (thin.map > 0) { for (tt in 1:thin.map) { llll <- length(mapcors$x) keep <- rep(c(TRUE,FALSE),length=llll) keep[(is.na(mapcors$x) | is.na(mapcors$y))] <- TRUE keep[llll] <- TRUE keep[c(is.na(mapcors$x[-1]),TRUE)] <- TRUE keep[c(TRUE,is.na(mapcors$x[-llll]))] <- TRUE mapcors <- list(x=mapcors$x[keep],y=mapcors$y[keep]) } rm(keep) } } # determine default breaks zmax <- max(zz); zmin <- min(zz) if (missing(breaks)) { breaks <- seq(zmin,zmax,length=length(col)+1) } else { if (length(breaks) != length(col)+1) { stop("** ERROR ** length of breaks must be length of colors plus 1") } } # determine colors for each pixel val2col <- function(val) { if (is.na(val)) {return(col.miss)} hh <- which(val <= c(breaks,val))[1] if (hh==(length(breaks)+1)) {return(col.miss)} if (hh==1) {return(col.miss)} col[hh-1] } zcols <- sapply(zz,FUN=val2col) if(con){ res<-contourLines(x,y,z) contours_x <- unlist(sapply(res, function(x) c(x$x, NA))) contours_y <- unlist(sapply(res, function(x) c(x$y, NA))) contours <- list(x=contours_x, y=contours_y) contu<-mapproject(contours) } # AND HERE WE GO ---------------------------------------------- # prepare the plot with the frame plot(x=c(1),y=c(1),type="n",axes=FALSE, xlim=xxlim,ylim=yylim,xlab="",ylab="") # plot the colored pixel polygons for (i in seq(1,length(xx))) { polygon(x=c(xxll[i],xxlr[i],xxur[i],xxul[i]), y=c(yyll[i],yylr[i],yyur[i],yyul[i]), col=zcols[i],border=zcols[i],xpd=FALSE) } if (con) { lines(x=contu$x,y=contu$y,xpd=FALSE) } # add the mapoutlines on top if (mapdat != "none") { lines(x=mapcors$x,y=mapcors$y,xpd=FALSE) } # add grid-lines for lon and lat circles if(longr=="yes"){ for (i in longrds) { #ygr <- seq(-90,90,by=bylatgrd) if(ylim.map[1]>0) ygr <- seq(ymaplim[1]-0.5*(diff(y[1:2])),ymaplim[2],by=bylatgrd) if(ylim.map[1]<0) ygr <- seq(ymaplim[2]+0.5*(diff(y[1:2])),ymaplim[1],by=-bylatgrd) grd <- mapproject(x=rep(i,length(ygr)),y=ygr) lines(x=grd$x,y=grd$y,lty=ltygrds,lwd=lwdgrds,xpd=FALSE) }} if(latgr=="yes"){ for (i in latgrds) { xgr <- seq(-180,180,by=bylongrd) grd <- mapproject(x=xgr,y=rep(i,length(xgr))) lines(x=grd$x,y=grd$y,lty=ltygrds,lwd=lwdgrds,xpd=FALSE) }} # add a box # box() }