Skip to content
Ben Tupper edited this page Aug 12, 2015 · 1 revision

library(raster) library(gridExtra) library(spnc)

We start by defining a bounded box for the region of interest.

bb <- c(-72,-63,39,46)
# opendap
X <- SPNC('http://oceandata.sci.gsfc.nasa.gov:80/opendap/MODISA/L3SMI/2015/188/A2015188.L3m_DAY_CHL_chlor_a_4km.nc', bb = bb)
X
# Reference Class: "L3SMIRefClass" 
#   flavor: source=L3SMI type=raster local=FALSE 
#   state: opened 
#   bounding box: -72 -63 39 46 
#   VARS: palette chlor_a 
#   DIMS: eightbitcolor=256 lat=4320 lon=8640 rgb=3 
#   LON: [ -179.979172, 179.979172] 
#   LAT: [ 89.979164, -89.979179] 
#   TIME: [ 2015-07-07, 2015-07-07] 

# local downloaded via
# http://oceandata.sci.gsfc.nasa.gov/cgi/getfile/A2015188.L3m_DAY_CHL_chlor_a_4km.nc
x <- SPNC("~/Downloads/A2015188.L3m_DAY_CHL_chlor_a_4km.nc", bb=bb)
x
# Reference Class: "L3SMIRefClass" 
#   flavor: source=L3SMI type=raster local=TRUE 
#   state: opened 
#   bounding box: -72 -63 39 46 
#   VARS: chlor_a palette 
#   DIMS: lat=4320 lon=8640 rgb=3 eightbitcolor=256 
#   LON: [ -179.979172, 179.979172] 
#   LAT: [ 89.979164, -89.979179] 
#   TIME: [ 2015-07-07, 2015-07-07] 

Other than the order in which thew variables are stored, there is not a lot to distinguish the two. That's good! Now let's do the extraction and transformation into a SpatialRaster object for each. Note that we must specify the variable to retrieve.

# get the rasters
R <- X$get_raster("chlor_a")
r <- x$get_raster("chor_a")

We can show these are the same with simple subtraction. Technically, we should do something more vigorous as what R reports as 0 may in fact be some small number, but this is fine for our purposes.

# are they the same?
summary(R-r)
#         layer
# Min.        0
# 1st Qu.     0
# Median      0
# 3rd Qu.     0
# Max.        0
# NA's    27787

Or we can plot them side-by-side for visual confirmation. (Feel better now?)

g <- gridExtra::grid.arrange(spplot(log10(R)), spplot(log10(r)), ncol = 2)

Clone this wiki locally