thredds.R

u <- "http://coastwatch.pfeg.noaa.gov/erddap/griddap/erdQSwind3day"

library(tidync)
(tnc <- tidync(u))

## not a file: 
## ' http://coastwatch.pfeg.noaa.gov/erddap/griddap/erdQSwind3day '
## 
## ... attempting remote connection
## Connection succeeded.

## 
## Data Source (1): erdQSwind3day ...
## 
## Grids (5) <dimension family> : <associated variables> 
## 
## [1]   D2,D1,D0,D3 : x_wind, y_wind    **ACTIVE GRID** ( 13058345694  values per variable)
## [2]   D0          : altitude
## [3]   D1          : latitude
## [4]   D2          : longitude
## [5]   D3          : time
## 
## Dimensions (4): 
##   
##   dimension    id      name length unlim coord_dim 
##       <chr> <dbl>     <chr>  <dbl> <lgl>     <lgl> 
## 1        D0     0  altitude      1 FALSE      TRUE 
## 2        D1     1  latitude   1201 FALSE      TRUE 
## 3        D2     2 longitude   2881 FALSE      TRUE 
## 4        D3     3      time   3774 FALSE      TRUE

tnc %>% hyper_filter()

## # A tibble: 1 x 2
##                access
##                <dttm>
## 1 2017-07-07 00:00:36
## # ... with 1 more variables: source <chr>
## filtered dimension summary: 
## # A tibble: 4 x 5
##        name coord_dim       min        max length
##       <chr>     <lgl>     <dbl>      <dbl>  <int>
## 1 longitude      TRUE         0        360   2881
## 2  latitude      TRUE       -75         75   1201
## 3  altitude      TRUE        10         10      1
## 4      time      TRUE 932644800 1258632000   3774

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

## now build an an explicit extraction
ht <- tnc %>% hyper_filter(longitude = longitude > 140 & longitude < 152, 
                                 latitude = between(latitude, -47, -40), 
                                 time = index %% 400 == 0) %>% 
  hyper_tibble()




library(ggplot2)
ggplot(ht, aes(longitude, latitude, fill = sqrt(x_wind^2 + y_wind^2))) + 
  geom_raster() + facet_wrap(~time) + coord_fixed(1/cos(44 * pi/180))

ht[c("x", "y")] <- tibble::as_tibble(rgdal::project(as.matrix(dplyr::select(ht, longitude, latitude)), "+proj=laea +lon_0=160 +lat_0=-10 +datum=WGS84"))
const <- 7200

ggplot(ht %>% dplyr::filter(!is.na(x_wind)) %>% sample_frac(0.1), 
       aes(x, y, xend = x + x_wind * const, yend = y + y_wind * const, colour = sqrt(x_wind^2 + y_wind^2))) + 
  geom_segment() + facet_wrap(~time)