Wednesday, February 17, 2016

Working with Globcolour data (Part 2)

Just a quick note to announce that the makeGlobcolourField and isin.convert functions have been added to the sinkr package. In addition, the makeGlobcolourField function now used the ncdf4 package to read the .nc files. Both functions are only set up to deal with the higher resolution 4 km data based on the ISIN grid ("L3b").
The following script is an example of extracting data for the Philippines, and produces a map of mean Chl1 values:






Example script:

# Required packages -------------------------------------------------------
 
library(sinkr)
library(maps)
library(mapproj)
library(ncdf4)
library(Matrix)
library(mapdata)
 
 
# Create field ------------------------------------------------------------
 
dat_path <- getwd() # or location of Globcolour .nc files
 
F1 <- makeGlobcolourField(
  folderPath=dat_path,
  yearRan=c(1997,2015),
  lonRan=c(116,127),
  latRan=c(4,17),
  greps=c(".nc", "GSM")
)
 
 
# Create grid and polygons ------------------------------------------------
 
grd <- isin.convert(grd=as.numeric(colnames(F1)))
polys <- isin.convert(grd=grd$grd, polygons = TRUE)
 
 
 
# Plot --------------------------------------------------------------------
 
breaks <- c(0.01, 0.05, 0.1, 0.2, 0.5, 1, 1.5, 2, 2.5, 3)
png("meanChla.png", res=400, units="in", width=6, height=5)
layout(matrix(1:2,1,2), widths=c(5,1), heights=5, respect = TRUE)
op <- par(cex=1, ps=10)
par(mar=c(4,4,1,1))
map("worldHires", "philippines", col=NA, ylim=range(grd$lat), xlim=range(grd$lon))
axis(1); axis(2)
box()
incl <- which(colSums(!is.na(F1)) / nrow(F1) > 0.5) # more that 50% of months are non-NAs
col <- val2col(colMeans(F1, na.rm = TRUE)[incl], col=jetPal(length(breaks)-1), breaks=breaks)
for(j in seq(incl)){
  polygon(polys[[incl[j]]], col=col[j], border=col[j], lwd=0.001)
}
map("worldHires", add=TRUE, col=8, resolution = 0, fill=TRUE)
par(mar=c(4,0,1,4))
imageScale(1, breaks=breaks, axis.pos=4, col=jetPal(length(breaks)-1), log="y", add.axis=FALSE)
axis(4, at=breaks, las=2)
par(op)
dev.off()
Created by Pretty R at inside-R.org

2 comments:

  1. Hi! I downloaded your package and tried to reproduce your example. But I discovered the function makeGlobcolourField was not included in the package. Is this a download error and can you maybe provide the script of the function, so I can use it? Thanks!

    ReplyDelete
    Replies
    1. Hi Svenja,
      I removed this function from the package due to some possible issue building the ncdf4 package. I think you should be able to use it though. I have added the function to a gist:
      https://gist.github.com/a2165a5ce38b6e17971c6eddcaafbe74.git

      Delete