## devtools::install_github("hypertidy/tidync")
## example thanks to Ben Tupper
uri <- 'https://oceandata.sci.gsfc.nasa.gov:443/opendap/MODISA/L3SMI/2016/001/A20160012016032.L3m_R32_SST_sst_9km.nc'
library(tidync)
## Warning: Installed Rcpp (0.12.12) different from Rcpp used to build dplyr (0.12.11).
## Please reinstall dplyr to avoid random crashes or undefined behavior.
(nc <- tidync(uri))
## not a file:
## ' https://oceandata.sci.gsfc.nasa.gov:443/opendap/MODISA/L3SMI/2016/001/A20160012016032.L3m_R32_SST_sst_9km.nc '
##
## ... attempting remote connection
## Connection succeeded.
##
## Data Source (1): A20160012016032.L3m_R32_SST_sst_9km.nc ...
##
## Grids (4) <dimension family> : <associated variables>
##
## [1] D0,D3 : palette **ACTIVE GRID** ( 768 values per variable)
## [2] D2,D1 : qual_sst, sst
## [3] D1 : lat
## [4] D2 : lon
##
## Dimensions (4):
##
## dimension id name length unlim coord_dim
## <chr> <dbl> <chr> <dbl> <lgl> <lgl>
## 1 D0 0 eightbitcolor 256 FALSE FALSE
## 2 D1 1 lat 2160 FALSE TRUE
## 3 D2 2 lon 4320 FALSE TRUE
## 4 D3 3 rgb 3 FALSE FALSE
## we have to activate the virtual table for dimensions [2,1]
(nc <- activate(nc, "D2,D1"))
##
## Data Source (1): A20160012016032.L3m_R32_SST_sst_9km.nc ...
##
## Grids (4) <dimension family> : <associated variables>
##
## [1] D0,D3 : palette
## [2] D2,D1 : qual_sst, sst **ACTIVE GRID** ( 9331200 values per variable)
## [3] D1 : lat
## [4] D2 : lon
##
## Dimensions (4):
##
## dimension id name length unlim coord_dim
## <chr> <dbl> <chr> <dbl> <lgl> <lgl>
## 1 D0 0 eightbitcolor 256 FALSE FALSE
## 2 D1 1 lat 2160 FALSE TRUE
## 3 D2 2 lon 4320 FALSE TRUE
## 4 D3 3 rgb 3 FALSE FALSE
## now let's explore
nc %>% hyper_filter()
## # A tibble: 1 x 2
## access
## <dttm>
## 1 2017-07-20 22:13:58
## # ... with 1 more variables: source <chr>
## filtered dimension summary:
## # A tibble: 2 x 5
## name coord_dim min max length
## <chr> <lgl> <dbl> <dbl> <int>
## 1 lon TRUE -179.95833 179.95836 4320
## 2 lat TRUE -89.95834 89.95834 2160
## we'll have to ask for two pieces
## I am a bit frightened of processing these two requests automatically
## but will give it some thought
west <- nc %>% hyper_filter(lon = lon <= -90)
east <- nc %>% hyper_filter(lon = lon >= 100)
## this is very DIY, but could be wrapped with an R function
vname <- "sst"
west_slab <- hyper_slice(west, select_var = vname)
east_slab <- hyper_slice(east, select_var = vname) ## use select_var, else we get qual_sst as well
## these are lists because the virtual tables of D-imension combinations
## can have multiple variables
str(west_slab)
## List of 1
## $ sst: num [1:1080, 1:2160] NA NA NA NA NA NA NA NA NA NA ...
str(east_slab)
## List of 1
## $ sst: num [1:960, 1:2160] NA NA NA NA NA NA NA NA NA NA ...
library(raster)
## Loading required package: sp
## we need this transpose, so rbind
r <- setExtent(raster(t(rbind(east_slab[[vname]], west_slab[[vname]]))),
extent(100, 180 + 90, -90, 90))
## note that the extent there is a lucky guess, this needs more careful handling
## to ensure we got the right cell-extent from the hyper_filter info
## but, this is "good progress" from my perspective
## and, I maintain, this is the most efficient possible opendap extraction
## and will work with more complex examples that won't work with raster in the
## usual workflow
plot(r)
library(maps)
map("world2", add = TRUE)
## extra illustration, let's go ggplot2 route (you wouldn't do this for
## a long time series extraction, hyper_tibble uses hyper_slice under the hood and
## we could save the coord expansion every time)
## also, we are re-reading the source here (which a future workflow will avoid)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
##
## intersect, select, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
tib <- dplyr::bind_rows(
hyper_tibble(west) %>% dplyr::mutate(lon = lon + 360), hyper_tibble(east)
) %>% dplyr::filter(!is.na(sst))
library(ggplot2)
ggplot(tib, aes(x = lon, y = lat, fill = sst, border = NA)) + geom_raster()
