Primary production august 2009

## ==============================================================================
## Map of primary production around Europe
## implemented by Karline Soetaert
## ==============================================================================
# Pimary production data have been downloaded from 
#  \url{http://www.science.oregonstate.edu/ocean.productivity}

#Behrenfeld, MJ, PG Falkowski, 1997.
#Photosynthetic rates derived from satellite-based chlorophyll concentration .
#Limnology and Oceanography (42): 1-20
#\url{http://www.science.oregonstate.edu/ocean.productivity/references/L&O\%201997a.pdf}
url <- "http://www.science.oregonstate.edu/ocean.productivity"

# They are in database format
  require(OceanViewData)
  data(pp.aug2009.db)

# convert them to cross table
  Aug2009 <- db2cross(pp.aug2009.db)

# dimensionality of the problem
  dim(Aug2009$z)

  range(Aug2009$x)
  range(Aug2009$y)

# Clean up and make image plot
# -9999 -> 0; max = 1000
  PP    <- Aug2009$z
  PP[PP < 0] <- NA

# Primary production near europe
  lon_eur  <- seq(-15, 30, by = 0.1)
  lat_eur  <- seq(30, 70, by = 0.1)
  PP_eur   <- remap(var = PP, x = Aug2009$x, y = Aug2009$y,
                  xto = lon_eur, yto = lat_eur)$var
par(oma = c( 2, 0, 0, 0))
image2D(x = lon_eur, y = lat_eur, z = PP_eur, log = "c", NAcol = "black", 
   rasterImage = TRUE, shade = 0.05, main = "PP aug 2009", clab = c("mgC/m2/d"), 
   colkey = list(width = 0.5, cex.axis = 0.8, length = 0.5, 
   dist = 0.02, cex.clab = 0.8, side.clab = 2, line.clab = 0.5), 
   las = 1, xlab = expression(paste("longitude ",degree*E)),
    ylab = expression(paste("latitude ",degree*N))
)