Plotting_LPJmL_output_in_R - PIK-LPJmL/LPJmL GitHub Wiki
Plotting LPJmL output in R
You find the output files of LPJmL in the output directory that you
defined in the lpjmlMexicoLanduse.conf. In this case, the directory is
$LPJROOT/outputMexicoLanduse. Here, you will see a number of binary
files that your can read into R and view the output, e.g. as scalar
fields projected on a map.
Output files starting with an “m” are monthly output, e.g. the mnpp.bin file stands for “monthly net primary production”.
Start R either on your local machine, or on the cluster. On the cluster, you need to load the module for R.
The advantage of using your local machine is that you can easily debug your R script within the R GUI or RStudio. On the cluster, you can only use the command line input. The disadvantage of using your local machine is that you cannot read in very large input files (e.g. global ones). For this tutorial, however, it is better to use your local machine so you can debug and change your code according to your needs.
Let’s assume you want to plot the average monthly primary production (NPP) in the simulation years 1980-2005 and the average across the entire period. So the output file of interest is mnpp.bin. As this file contains monthly data, you will have to do some averaging in R.
The R script starts with defining some general variables. Each block of code is commented here to lead you through the code.
settings
rm(list=ls(all=TRUE)) # deletes all variables
# load required libaries (install before usage!)
require(fields)
require(maps)
require(plotrix)
################################################### define custom variables ##################################
res <- 0.5 #grid resolution in degrees, 0.5 is LPJmLs default resolution
#dpath <- "/iplex/01/biodiversity/alice/LPJ/" #path for running on the cluster
dpath <- "W://biodiversity/alice/LPJ/" #path for debugging on windows with the mapped drive
dir <- "outputMexicoLanduse/" #your output directory
title <- "Net primary production (gC m^-2 month^-1) from file: " #title for the plot
title.nbands <- " " #not needed here
grid <- paste(dpath,dir,"grid.bin",sep="") #the grid is needed to plot the map
ofiles <- c("mnpp.bin") #the output file of interest (NPP in this case)
yearMin <- 1980 #the first year of the evaluation period
yearMax <- 2005 #the last year of the evaluation period
sizeof.data <- 4 #4 bytes for LPJ output, 2 bytes for climate data
sizeof.data.grid <- 2 # 2 bytes for lat/long pairs
data.type <- numeric() # numeric() for LPJ output, integer() for climate data
noutput_years = yearMax - yearMin #number of output years that should be plotted
#define the data range of the axes:
#x-, y-axis are for the grid, and the z-axis is for the output variable, e.g. NPP
set.zlim <- 1 #0=automatic, 1=disp.zlim as limit: if you want to define min/max of the axis yourself, use 1 here
disp.zlim=c(-10,200) #define custom limits for NPP
set.ylim <- 0 #keep automatically set limits
disp.ylim=c(-180,180) #unused in this case
set.xlim <- 0 #keep automatically set limits
disp.xlim=c(-90,90) #unused in this case
npixel <- file.info(grid)$size / sizeof.data.grid / 2 #number of pixels in the image. Divide by 2 because there are 2 values in the grid (lat/lon)
npixel
nyears <- file.info(paste(dpath,dir,"vegc.bin",sep=""))$size/sizeof.data/npixel #total number of years in the data set
nyears
colors <- matrix(0,nrow=length(ofiles),ncol=100) #set color range
for (i in 1:length(ofiles)) { colors[i,] <- tim.colors(100)} #default blue to red
reading the grid
Next, you have to read in the underlying grid of Mexico:
################################################### read in the underlying grid ##################################
#read underlying grid (z-values of NPP are not read yet!)
grid.fn <- file(grid,"rb") #open binary file of the grid for reading (flag: 'rb')
seek(grid.fn, where=0, origin="start") #go to the starting pointer of the grid
grid.data <- matrix(readBin(grid.fn, integer(), n=npixel*2, size=2),ncol=2,byrow=T) / 100 #read in the grid
#define extent of the grid
ext.lon <- c(min(grid.data[,1]),max(grid.data[,1]))
ext.lat <- c(min(grid.data[,2]),max(grid.data[,2]))
nrow <- (ext.lat[2]-ext.lat[1])/res+1
ncol <- (ext.lon[2]-ext.lon[1])/res+1
#the image that is going to be created has pixels with i,j integer indices. We have to define
#the position in the image matrix for each coordinate to convert the lat/lon info to i,j position in the image matrix
grid.ilon <- as.integer((grid.data[,1]-ext.lon[1])/res + 1.01)
grid.ilat <- as.integer((grid.data[,2]-ext.lat[1])/res + 1.01)
close(grid.fn) #close the file
#create the image matrix for a single band (here, a month)
img <- array(NA, dim=c(nrow,ncol)) # create matrix with established dimensions filled with NA
if(!set.xlim) disp.xlim <- ext.lon #just making sure that the axis limits are set
if(!set.ylim) disp.ylim <- ext.lat
reading in data
Now you are ready to read in your z-values from the LPJmL output file, NPP in this case. Note that the following block of code opens a loop over all output files, so you may have several output files processed, if you wish. In this case, we only process one at a time (mnpp.bin).
################################################### read in NPP (z-values) ##################################
for (i in 1:length(ofiles)){ #open the loop over the output files (1 in this case)
#z-values are read in here
file <- ofiles[i] #file handler
col <- colors[i,] #colors
nbands <- file.info(paste(dpath,dir,file,sep=""))$size/sizeof.data/npixel/nyears #number of bands is 12 because of monthly data in this case
nbands
offset <- nbands * noutput_years * sizeof.data * npixel #offset in the data grid to only plot a certain range of years, e.g. 1980-2005
offset
file.f <- file(paste(dpath,dir,file,sep=""),"rb") #open binary file for reading (flag: 'rb')
seek(file.f, where=offset, origin="start") #seek for the desired years with the offset defined above
#read in file with z-values
# all pixels of one band in a row
# time series for a pixel in a column
file.data <- matrix(readBin(file.f, data.type, n=npixel*noutput_years*nbands, size=sizeof.data),ncol=npixel,nrow=noutput_years*nbands,byrow=T)
close(file.f)
if(!set.zlim) disp.zlim=c(min(file.data),max(file.data)) #make sure the z-axis range is set
slicing data
Next, a 3D plot matrix is composed of all the 2D images for each month which are stacked upon each other. That is, the 3D plot matrix contains all (still monthly) time slices that we want to average later. Afterwards, the loop over the output files is closed.
################################################### create a 3D temporal matrix of 2D images ##################################
#construct the 3D matrix
plot_matrix_all <- array(NA,dim=c(nrow,ncol,nbands*noutput_years)) # this contains all available time slices
for (b in seq(1,(nbands*noutput_years)))#for each band = each month
{
#for each pixel in the plot, write z-value (NPP) into 2D image at the matching i,j (ilon/ilat) position)
for(p in seq(1,npixel)) { img[grid.ilat[p],grid.ilon[p]] <- file.data[b,p] }
plot_matrix_all[,,b] <- img #write the entire img into 3D plot matrix (stack on top)
if(set.zlim){
img <- ifelse(img<disp.zlim[1],disp.zlim[1],img)
img <- ifelse(img>disp.zlim[2],disp.zlim[2],img)
}
} #endfor each band
}#endfor (outfiles)
averaging over time
Then, the averaging across the years is done.
################################################### averaging data across years #################################################################
#aggregate data to average NPP over each year and over the entire period
plot_matrix_ann_avg <- array(NA,dim=c(nrow,ncol,noutput_years)) #final output matrix with annual averages of NPP
plot_matrix_avg <- array(0,dim=c(nrow,ncol,1)) #final output matrix with average of NPP across the entire simulation period
nmonths = dim(plot_matrix_all)[3] #total number of nbands = months in the simulation
for (i in seq(1:noutput_years))#loop over the output years you want to plot
{
NPP_ann_avg <- array(0,dim=c(nrow,ncol,noutput_years)) #temporary storage for summation and averaging
#remember that nrow = dim(plot_matrix_all)[1] and ncol = dim(plot_matrix_all)[2]
for ( latitude in seq(1:nrow) )#loop over matrix dimension 1 = latitude
{
for ( longitude in seq(1:ncol) )#loop over matrix dimension 2 = longitude
{
for (j in seq(1:nbands))#loop over months
{
#ID of the respective month in the matrix per year
#plot_matrix_all has dim 36x61x300, 300 because of 25 (years) * 12 (months)
#sum up over each year: (i-1)*nbands + j is e.g. 0+j, j=1,...12 for year 1 and 12+1=13 for year 2
monthID = (i-1)*nbands + j
NPP_ann_avg[latitude,longitude,i] = NPP_ann_avg[latitude,longitude,i] + 1/nbands * plot_matrix_all[latitude,longitude,monthID]#sum up over the months of year i and build arithm. average
plot_matrix_avg[latitude,longitude,1] = plot_matrix_avg[latitude,longitude,1] + 1/nmonths * plot_matrix_all[latitude,longitude,monthID]#sum up over all months and build average
#cat("MonthID: ", monthID, "\n", sep="") #debug
#cat("Matrix value: ", plot_matrix_all[latitude,longitude,monthID], "\n", sep="") #debug
#cat("Summed value: ", NPP_ann_avg[latitude,longitude,i], "\n", sep="") #debug
}
}#loop over lon
}#loop over lat
#store permanently
plot_matrix_ann_avg[,,i] = NPP_ann_avg[,,i]
#plot annual averages on the fly
image.plot(x=seq(ext.lon[1],ext.lon[2],by=res), y=seq(ext.lat[1],ext.lat[2],by=res),z=t(plot_matrix_ann_avg[,,i]), zlim=disp.zlim, xlim=disp.xlim, ylim=disp.ylim,xlab="",ylab="",asp=1,col=rev(col),tcl=-0.2,yaxp=c(disp.ylim,4),xaxp=c(disp.xlim,8),axes=T)
ititle <- paste("NPP (gC m^-2 month^-1), annual avg. of ", (yearMin + i - 1), sep = " ")
title(paste(title,file),ititle,cex.sub=1.5)
map("world",add=TRUE,fill=FALSE)
}#end (noutput_years)
#plot average across all years #
#plot total average on the fly
image.plot(x=seq(ext.lon[1],ext.lon[2],by=res), y=seq(ext.lat[1],ext.lat[2],by=res),z=t(plot_matrix_avg[,,1]), zlim=disp.zlim, xlim=disp.xlim, ylim=disp.ylim,xlab="",ylab="",asp=1,col=rev(col),tcl=-0.2,yaxp=c(disp.ylim,4),xaxp=c(disp.xlim,8),axes=T)
ititle <- paste("NPP (gC m^-2 month^-1), avg. ", yearMin, " - ", yearMax, sep = " ")
title(paste(title,file),ititle,cex.sub=1.5)
map("world",add=TRUE,fill=FALSE)
write to PDF
After gathering all the averages you need, write them out to PDF for selected years.
################################################### write out to PDF ##################################
#output annual averages in pdf file
pdf.file <- paste(dpath,dir,"NPP_Mexico.pdf",sep='') #choose file name
pdf(pdf.file, width=600,height=600,paper="default") #open the file
par(mfrow=c(3,2)) #page layout
yearIDs = c(1,6,12,18,25) #year IDs starting from yearMin = 1980 you like to plot
for (i in 1:length(yearIDs))#for each year you like to plot
{
yearID = yearIDs[i]
#plot the image (it will be written out as PDF due to previous 'pdf' command
image.plot(x=seq(ext.lon[1],ext.lon[2],by=res), y=seq(ext.lat[1],ext.lat[2],by=res),z=t(plot_matrix_ann_avg[,,yearID]), zlim=disp.zlim, xlim=disp.xlim, ylim=disp.ylim,xlab="lon",ylab="lat",asp=1,col=rev(col),tcl=-0.2,yaxp=c(disp.ylim,4),xaxp=c(disp.xlim,8),axes=T)
par(new=TRUE)
map("world",add=TRUE,fill=FALSE)
ititle <- paste("NPP (gC m^-2 month^-1) in ", yearMin + yearID - 1, sep = " ")
title(ititle,cex.sub=1.5)
}
#add the total avg. over all years from yearMin to yearMax
image.plot(x=seq(ext.lon[1],ext.lon[2],by=res), y=seq(ext.lat[1],ext.lat[2],by=res),z=t(plot_matrix_avg[,,1]), zlim=disp.zlim, xlim=disp.xlim, ylim=disp.ylim,xlab="lon",ylab="lat",asp=1,col=rev(col),tcl=-0.2,yaxp=c(disp.ylim,4),xaxp=c(disp.xlim,8),axes=T)
ititle <- paste("NPP (gC m^-2 month^-1), avg. ", yearMin, " - ", yearMax, sep = " ")
title(ititle,cex.sub=1.5)
map("world",add=TRUE,fill=FALSE)
dev.off() #close output to PDF
Finally, you can also write out the results as ASCII grids, e.g. for reading them into ArcGIS.
#######################################write out to ASCII grid ################################
###the avg over each specific year
yearIDs = c(1,6,12,18,25)
for (i in 1:length(yearIDs))
{
yearID = yearIDs[i]
matrixRotated <- rot90(plot_matrix_ann_avg[,,yearID], -1) #matrix must be rotated anti-clockwise to be displayed correctly in ArcGIS
v_lpjoutput <- as.vector(matrixRotated) #convert to a vector
#establish a grid topology of correct size
#arguments are: lower left corner of the country grid (Mexico in this case), grid resolution (half degree in this case), and dimensions of the read-in plot matrix
#you can get the first argument for the grid from variable ext.lon[1] and ext.lat[1] (this is lower left corner on grid))
xy <- GridTopology(c(ext.lon[1], ext.lat[1]), c(res, res), c(dim(plot_matrix_avg)[2], dim(plot_matrix_avg)[1])) #turn dimension of plot matrix around because it is rotated anti-clockwise by 90 degrees
#do projections for the spatial object of type spatial grid data frame..crs is the coordinate reference system, achieves conversion between coordinate systems
sgdf_lpjoutput <- SpatialGridDataFrame(xy,proj4string=CRS("+proj=lonlat"),data=data.frame(v_lpjoutput))
names(sgdf_lpjoutput) <- c('NPP')
#write out the ASCII grid
# na.value is the no-data value to prevent the grid from crashing at undefined values
write.asciigrid(sgdf_lpjoutput,paste(dpath,dir,"NPP_avg_", (yearMin + yearID - 1) , "_mexico.txt",sep=''),na.value= -9999)
}
#If you don't have ArcGIS, you can read an ASCII grid back in for testing, e.g. by
#fname <- paste(dpath,dir,"NPP_avg_1981_mexico.txt",sep='')
#x <- read.asciigrid(fname, as.image = FALSE, plot.image = TRUE, colname = fname, proj4string = CRS(as.character(NA)))
#image(x)
###the avg over all years
matrixRotated <- rot90(plot_matrix_avg[,,1], -1)
v_lpjoutput <- as.vector(matrixRotated)
xy <- GridTopology(c(ext.lon[1], ext.lat[1]), c(res, res), c(dim(plot_matrix_avg)[2], dim(plot_matrix_avg)[1])) #when rotated by 90 degrees anti-clockwise
sgdf_lpjoutput <- SpatialGridDataFrame(xy,proj4string=CRS("+proj=lonlat"),data=data.frame(v_lpjoutput))
names(sgdf_lpjoutput) <- c('NPP')
write.asciigrid(sgdf_lpjoutput,paste(dpath,dir,"NPP_avg_1980-2005_mexico.txt",sep=''),na.value= -9999)
An example of the PDF output from this R script is attached to this page below.
plotting raster maps in alternative projections
This can be a real pain in the a.., er, neck. Here’s an example of how it works.
rm(list=ls(all=TRUE))
gc()
#require(maptools)
#require(ncdf4)
require(raster)
require(rgdal)
require(ggplot2)
require(gridExtra)
#### crop names
cr.names <- c("mai","whe","ric","soy")
#### prepare CWP data from NC file
f <- "D:/data/GGCMI/median_delta_cwp_2050.nc" # NC file with 4 layers for maize, wheat, rice and soy, global extent
cwp <- raster(f)
#plot(cwp)
#read individual layers to separate variables
for( i in 1:cwp@file@nbands ){
assign( paste("cwp" , i , sep = ".") , raster(f, band=i))
}
# function to project raster data and convert to data frame as needed by ggplot
project.raster <- function(data,pr="+proj=robin",str="CWP",newex=c(-180,180,-65,85)){
newext <- newex
data <- crop(data,newext)
data.p <- projectRaster(data,crs=pr)
# inspired by http://nrelscience.org/2013/05/30/this-is-how-i-did-it-mapping-in-r-with-ggplot2/
#convert the raster to points for plotting
map.p <- rasterToPoints(data.p)
#Make the points a dataframe for ggplot
df.data <- data.frame(map.p)
#Make appropriate column headings
colnames(df.data) <- c("long", "lat", str)
df.data
}
#### read spatial data for nice look
# inspired by http://rpsychologist.com/working-with-shapefiles-projections-and-world-maps-in-ggplot
setwd("D:/data/GGCMI/")
# downloaded from http://www.naturalearthdata.com/downloads/110m-cultural-vectors/110m-admin-0-countries/
# add country borders
countries <- readOGR("ne_110m_admin_0_countries", layer="ne_110m_admin_0_countries")
countries_robin <- spTransform(countries, CRS("+init=ESRI:54030"))
countries_robin_df <- fortify(countries_robin)
# downloaded from http://www.naturalearthdata.com/downloads/110m-physical-vectors/110m-graticules/
grat <- readOGR("ne_110m_graticules_all", layer="ne_110m_graticules_15") # this is the 15° grat. There are finer and coarser data available, too
bbox <- readOGR("ne_110m_graticules_all", layer="ne_110m_wgs84_bounding_box")
grat_robin <- spTransform(grat, CRS("+proj=robin")) # reproject graticule
grat_robin_df <- fortify(grat_robin)
bbox_robin <- spTransform(bbox, CRS("+proj=robin")) # reproject bounding box
bbox_robin_df <- fortify(bbox_robin)
#### plot things using ggplot
# taken from https://github.com/hadley/ggplot2/wiki/Share-a-legend-between-two-ggplot2-graphs
grid_arrange_shared_legend <- function(...) {
plots <- list(...)
g <- ggplotGrob(plots[1](/PIK-LPJmL/LPJmL/wiki/1) + theme(legend.position="bottom"))$grobs
legend <- g[which(sapply(g, function(x) x$name) == "guide-box")](/PIK-LPJmL/LPJmL/wiki/which(sapply(g,-function(x)-x$name)-==-"guide-box"))
lheight <- sum(legend$height)
grid.arrange(
do.call(arrangeGrob, lapply(plots, function(x)
x + theme(legend.position="none"))),
legend,
ncol = 1,
heights = unit.c(unit(1, "npc") - lheight, lheight))
}
plot.proj.raster <- function(bbox,rast.df,countries.df,grat.df,tit){
# inspired by http://rpsychologist.com/working-with-shapefiles-projections-and-world-maps-in-ggplot
# some settings, set most thing to empty, background to grey
theme_opts <- list(theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_rect(fill="#e6e8ed"),
panel.border = element_blank(),
axis.line = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
plot.title = element_text(size=10)))
# conversion of raster to points seems to be the only workable method to combine grid lines with raster data in specific projections....
ggplot(bbox, aes(long,lat)) +
geom_polygon(fill="white") +
#geom_polygon(data=df.cwp.robin, aes(long,lat, group=group)) +
geom_raster(aes(long,lat,fill=CWP),data=rast.df) +
geom_path(data=countries.df, aes(long,lat, group=group), color="black", size=0.3) +
geom_path(data=grat.df, aes(long, lat, group=group, fill=NULL), linetype="dashed", color="grey50",size=0.3) +
labs(title=tit) +
coord_equal() +
theme_opts +
scale_fill_distiller(palette="Spectral")
}
#### for one layer tryout
df.cwp <- project.raster(cwp.1)
plot.proj.raster(bbox_robin_df,df.cwp,countries_robin_df,grat_robin_df,"CWP")
#### plot to 2x2 png file with one legend
png("CWP_robinson_R_tryout.png",width=5*300,height=4*300,res=300,pointsize=6)
for(i in 1:4){
cwp.act <- get(paste("cwp",i,sep="."))
df.cwp.act <- project.raster(cwp.act)
assign(paste("p",i,sep="."),plot.proj.raster(bbox_robin_df,df.cwp.act,countries_robin_df,grat_robin_df,cr.names[i]))
}
#grid.arrange(p.1,p.2,p.3,p.4)
grid_arrange_shared_legend(p.1, p.2, p.3, p.4)
dev.off()