Net primary production 2005

## =============================================================================
## Reads the monthly primary production data,
## converts them to cross tabular format and
## makes yearly summed primary production.
## implemented by Karline Soetaert
## =============================================================================

require (OceanView)
url <- "http://www.science.oregonstate.edu/ocean.productivity"

# These are the names of the xyz files, which have been ftp-d from the 
# primary production website; 
# the first line has been manually edited so that  
# it contains only the column headings

#jan <- read.table("2005/vgpm.2005001.all.xyz", header = TRUE)
#feb <- read.table("2005/vgpm.2005032.all.xyz", header = TRUE)
#mar <- read.table("2005/vgpm.2005060.all.xyz", header = TRUE)
#apr <- read.table("2005/vgpm.2005091.all.xyz", header = TRUE)
#may <- read.table("2005/vgpm.2005121.all.xyz", header = TRUE)
#jun <- read.table("2005/vgpm.2005152.all.xyz", header = TRUE)
#jul <- read.table("2005/vgpm.2005182.all.xyz", header = TRUE)
#aug <- read.table("2005/vgpm.2005213.all.xyz", header = TRUE)
#sep <- read.table("2005/vgpm.2005244.all.xyz", header = TRUE)
#oct <- read.table("2005/vgpm.2005274.all.xyz", header = TRUE)
#nov <- read.table("2005/vgpm.2005305.all.xyz", header = TRUE)
#dec <- read.table("2005/vgpm.2005335.all.xyz", header = TRUE)

#makecrosstab <- function(input) {
#  Dat <- db2cross(input)
#  names(Dat) <- c("lon", "lat", "value")
#  Dat$value[Dat$value == - 9999] <- NA
#  Dat$name = "Net Primary Production"
#  Dat$units = "mgC/m2/day"
#  Dat$month = deparse(substitute(input))
#  Dat
#}

# make a cross table of each data set
#pp.jan <- makecrosstab(jan)
#pp.feb <- makecrosstab(feb)
#pp.mar <- makecrosstab(mar)
#pp.apr <- makecrosstab(apr)
#pp.may <- makecrosstab(may)
#pp.jun <- makecrosstab(jun)
#pp.jul <- makecrosstab(jul)
#pp.aug <- makecrosstab(aug)
#pp.sep <- makecrosstab(sep)
#pp.oct <- makecrosstab(oct)
#pp.nov <- makecrosstab(nov)
#pp.dec <- makecrosstab(dec)

#save(file = "ppmonths.rda" , pp.jan, pp.feb, pp.mar, pp.apr,
#     pp.may, pp.jun, pp.jul, pp.aug, pp.sep, pp.oct, pp.nov, pp.dec)

# if, later we want to work on these data, we simply load the file
load(file = "ppmonths.rda")
par (mfrow = c(2, 2), oma = c(2, 0, 0, 0))
image2D(pp.jan$value, rasterImage = TRUE, log = "z", main = "januari", 
  zlim = c(10, 20000), NAcol = "black", clab = "mgC/m2/d")
image2D(pp.apr$value, rasterImage = TRUE, log = "z", main = "april", 
  zlim = c(10, 20000), NAcol = "black", clab = "mgC/m2/d")
image2D(pp.jul$value, rasterImage = TRUE, log = "z", main = "july", 
  zlim = c(10, 20000), NAcol = "black", clab = "mgC/m2/d")
image2D(pp.oct$value, rasterImage = TRUE, log = "z", main = "october", 
  zlim = c(10, 20000), NAcol = "black", clab = "mgC/m2/d")
mtext(outer = TRUE, line = -1, cex = 1.2, "net primary production 2005")