[GIS] How to open and plot Globcolour chlorophyll image in R

Question

Level 3 globcolour chlorophyll images (download from hermes) are in netCDF format. Those files can be imported via ncdf package to R 3.1.

ch.nc = open.ncdf("E:/Data/globclour/L3b_20140101-20140131__GLOB_4_AV-MOD_CHL1_MO_00.nc")

How can I plot CHL1_mean? I follow guide by Mark but not getting any success.

Answer 1

Here is an example with the 25-km dataset, for January 2014, followed by an example for the 4-km dataset, January 2014 (improvement are welcome). The 4-km version is read using the example given here. You need to download functions from the same website: val2col, image.scale and earth.dist.

# Required librairies
library(rasterVis)
library(rgdal)
library(grid)

# Land mask using Natual Earth data (http://www.naturalearthdata.com/downloads/10m-physical-vectors/)
msk <- readOGR(".", "ne_10m_admin_0_countries_lakes")

# 25 km
CHL1 <- raster('L3m_20140101-20140131__GLOB_25_AV-MOD_CHL1_MO_00.nc', varname='CHL1_mean')
names(CHL1) <- 'CHL1_mean'

# Get values from the RasterLayer
z <- getValues(CHL1)

# Create breaks
brks <- c(0,0.01,0.02,0.03,0.05,0.1,0.2,0.3,0.5,1,2,3,5,10,20,30,50,100)

# Substitute chlorophyll values to classes
new_z <- findInterval(z, brks)

# Write new values in RasterLayer
CHL1 <- setValues(CHL1, new_z)

# Create color palette
pal <- rev(colorRampPalette(brewer.pal(11, 'Spectral'))(100))

# Customize legend
ckey <- list(at=seq_along(brks),
             labels=list(at=seq_along(brks),
                         labels=brks,
                         cex=1))
 # Plot chlorophyll
Obj <- 
  levelplot(CHL1, margin=FALSE, col.regions=pal, at=seq_along(brks), colorkey=ckey) +
  layer(sp.polygons(msk, fill='grey', col='transparent'))

print(Obj)
trellis.focus("legend", side="right", clipp.off=TRUE, highlight=FALSE)
grid.text(expression(paste("[",mg/m^3,"]",sep="")),1,0.50,rot=90)
trellis.unfocus()

Another approach for the 4-km version, adapted from here. This example is plotting chlorophyll around Florida.

# Required libraries
library(ncdf4)
library(classInt)
library(RColorBrewer)
library(rgdal)
library(rgeos)
library(raster)

# Create "extent" object for cropping, here focus on Florida
e <- extent(-85,-80,20,30)

# Land mask using Natual Earth data (http://www.naturalearthdata.com/downloads/10m-physical-vectors/)
msk <- readOGR(".", "ne_10m_admin_0_countries_lakes")
CP <- as(e, "SpatialPolygons")
proj4string(CP) <- CRS(proj4string(msk))
ne_10m_land <- gIntersection(msk, CP, byid=TRUE)

# Open the NetCDF file
nc <- nc_open('L3b_20140101-20140131__GLOB_4_AV-MOD_CHL1_MO_00.nc')
idxrow <- ncvar_get(nc, "row")+1
idxcol <- ncvar_get(nc, "col")+1
val <- ncvar_get(nc, "CHL1_mean")
nc_close(nc)

# From http://menugget.blogspot.com/2012/04/working-with-globcolour-data.html
Re <- 6378.145
Nlat <- 4320

dr <- (pi*Re)/Nlat
dPHI <- pi/Nlat

PHI <- -(pi/2)+(1:Nlat)*dPHI-(dPHI/2)

p <- 2*pi*Re*cos(PHI)
Nlon <- round(p/dr, 0)
dlon <- p/Nlon
dphi <- (2*pi)/Nlon
Ntot <- sum(Nlon)
lat <- PHI*180/pi

g.row <- idxrow
g.col <- idxcol
g.lat  <-  lat[g.row]
Nlon_row <- Nlon[g.row]
g.lon <- (360*(g.col-0.5)/Nlon_row) - 180

Nlon_row <- Nlon[g.row]
g.width <- 360/Nlon_row
g.height <- 180/Nlat


# Crating SpatialPointsDataFrame
CHL1 <- data.frame(x=g.lon, y=g.lat, g.width=g.width, g.height=g.height, z=val)
coordinates(CHL1) <- ~x+y
proj4string(CHL1) <- CRS('+init=epsg:4326')

# Cropping the SpatialPointsDataFrame to the area of interest
CHL1 <- crop(CHL1, e)

# Get values from the SpatialPointsDataFrame 
z <- CHL1$z

# Create breaks
brks <- c(0,0.01,0.02,0.03,0.05,0.1,0.2,0.3,0.5,1,2,3,5,10,20,30,50,100)

# Substitute chlorophyll values to classes
new_z <- findInterval(z, brks, rightmost.closed = TRUE)

# Write new values in RasterLayer
CHL1$z <- new_z

# Creating polygons
out <- split(cbind(coordinates(CHL1), CHL1$g.width, CHL1$g.height, 1:length(CHL1)), 1:length(CHL1))
foo <- function(x){
  xx <- c(x[1]-x[3]/2, x[1]-x[3]/2, x[1]+x[3]/2, x[1]+x[3]/2, x[1]-x[3]/2)
  yy <- c(x[2]-x[4]/2, x[2]+x[4]/2, x[2]+x[4]/2, x[2]-x[4]/2, x[2]-x[4]/2)

  Sr <- Polygon(cbind(xx,yy))
  Srs <- Polygons(list(Sr), paste0('s',x[5]))
}
polys <- lapply(out, foo)
polys <- SpatialPolygons(polys, 1:length(polys))
projection(polys) <- CRS('+init=epsg:4326')

# Creating colors
source("val2col.R")
pal0 <- rev(colorRampPalette(brewer.pal(11, 'Spectral'))(100))
colcode <- val2col(new_z, col=pal0)

# Mapping
source("earth.dist.R")
source("image.scale.R")
XLIM <- c(e@xmin, e@xmax) 
YLIM <- c(e@ymin, e@ymax)
x.dist <- earth.dist(XLIM[1], mean(YLIM), XLIM[2], mean(YLIM))
y.dist <- earth.dist(mean(XLIM), YLIM[1], mean(XLIM), YLIM[2]) 
heights <- c(6)
widths <- c(heights[1]*x.dist/y.dist, 1.5)

par(omi=c(0.1, 0.1, 0.1, 0.1), ps=12)
layout(matrix(c(1,2),nrow=1,ncol=2,byrow=TRUE), widths = widths, heights = heights, respect=TRUE)
layout.show(2)
par(mai=c(0.5, 0.5, 0.1, 0.1), las=1)
plot(CHL1$x, CHL1$y, t='n', xlim=XLIM, ylim=YLIM, xaxs="i", yaxs="i", xaxt="n", yaxt="n")
plot(polys, col=colcode, border='transparent', add=TRUE)
plot(ne_10m_land, border='transparent', col='grey', add=TRUE)
axis(1, seq(-85,-80,2), parse(text=degreeLabelsEW(seq(-85,-80,2))))
axis(2, seq(20,30,2), parse(text=degreeLabelsNS(seq(20,30,2))))
box(lwd=2)

par(mai=c(0.5, 0.1, 0.1, 1))
image.scale(new_z, col=pal0, horiz=FALSE, yaxt="n", ylab="", xlab="")
lab <- sort(unique(new_z))[2:(length(unique(new_z))-1)]
axis(4, lab, brks[lab])
mtext(expression(paste("[",mg%.%m^-3,"]", sep="")), side=4, line=3)
box()