Illustrate "non-catalog" request

[mikejohnson51]

library(opendap.catalog)
library(AOI)

# Define Area of interest as Colorado
AOI = aoi_get(state = "CO")

# Build URL (could do w/ glue and a date)
url = 'https://www.ncei.noaa.gov/thredds-ocean/dodsC/ncei/archive/data/0129374/daily/livneh_NAmerExt_15Oct2014.195001.nc'

# All data!
system.time({ out = dap(url,AOI = AOI) })
#> source:   https://www.ncei.noaa.gov/thredds-ocean/dodsC/ncei/archive/d... 
#> varname(s):
#>    > Prec [mm] (NA)
#>    > Tmax [C] (NA)
#>    > Tmin [C] (NA)
#>    > wind [m/s] (mean daily wind speed)
#> ==================================================
#> diminsions:  113, 66, 31 (names: lon,lat,time)
#> resolution:  0.062, 0.062, 1 days
#> extent:      -109.06, -102, 36.94, 41.06 (xmin, xmax, ymin, ymax)
#> crs:         +proj=longlat +a=6378137 +f=0.00335281066474748 +p...
#> time:        1950-01-01 to 1950-01-31
#> ==================================================
#> values: 924,792 (vars*X*Y*T)
#>    user  system elapsed 
#>   0.944   0.146  10.100

str(out, max.level = 1)
#> List of 4
#>  $ Prec:Formal class 'SpatRaster' [package "terra"] with 1 slot
#>  $ Tmax:Formal class 'SpatRaster' [package "terra"] with 1 slot
#>  $ Tmin:Formal class 'SpatRaster' [package "terra"] with 1 slot
#>  $ wind:Formal class 'SpatRaster' [package "terra"] with 1 slot

# Just Prec
system.time({  out2 = dap(url, AOI = AOI, varname = "Prec") })
#> source:   https://www.ncei.noaa.gov/thredds-ocean/dodsC/ncei/archive/d... 
#> varname(s):
#>    > Prec [mm] (NA)
#> ==================================================
#> diminsions:  113, 66, 31 (names: lon,lat,time)
#> resolution:  0.062, 0.062, 1 days
#> extent:      -109.06, -102, 36.94, 41.06 (xmin, xmax, ymin, ymax)
#> crs:         +proj=longlat +a=6378137 +f=0.00335281066474748 +p...
#> time:        1950-01-01 to 1950-01-31
#> ==================================================
#> values: 231,198 (vars*X*Y*T)
#>    user  system elapsed 
#>   0.561   0.052   3.033

str(out2, max.level = 1)
#> List of 1
#>  $ Prec:Formal class 'SpatRaster' [package "terra"] with 1 slot

^{Created on 2022-04-13 by the reprex package (v2.0.1)}

[mikejohnson51]

library(terra)
library(opendap.catalog)
library(AOI)
library(zonal)
library(rvest)

# Find Files
base = 'https://psl.noaa.gov/thredds/catalog/Datasets/livneh/metvars/'

dat = read_html(paste0(base, 'catalog.html')) |> 
   html_nodes("a")

files = paste0(gsub("catalog", "dodsC", base), 
                basename(html_attr(dat, "href")))

length(files)
#>[1] 376

head(files)
#>[1] "https://psl.noaa.gov/thredds/dodsC/Datasets/livneh/metvars/metvars"        
#>[2] "https://psl.noaa.gov/thredds/dodsC/Datasets/livneh/metvars/cdd.mon.ltm.nc" 
#>[3] "https://psl.noaa.gov/thredds/dodsC/Datasets/livneh/metvars/cdd.mon.mean.nc"
#>[4] "https://psl.noaa.gov/thredds/dodsC/Datasets/livneh/metvars/hdd.mon.ltm.nc" 
#>[5] "https://psl.noaa.gov/thredds/dodsC/Datasets/livneh/metvars/hdd.mon.mean.nc"
#>[6] "https://psl.noaa.gov/thredds/dodsC/Datasets/livneh/metvars/prec.1915.nc"   

# Define AOI
AOI = aoi_get(state = "CO", county = "all")

# Get 'prec' data from file '100' 
out = dap(files[100], AOI = AOI, varname = "prec")
#> source:   https://psl.noaa.gov/thredds/dodsC/Datasets/livneh/metvars/p... 
#> varname(s):
#>    > prec [mm] (Daily Total Precipitation)
#> ==================================================
#> diminsions:  113, 66, 365 (names: lon,lat,time)
#> resolution:  0.062, 0.062, 1 days
#> extent:      -109.06, -102, 36.94, 41.06 (xmin, xmax, ymin, ymax)
#> crs:         +proj=longlat +a=6378137 +f=0.00335281066474748 +p...
#> time:        2009-01-01 to 2009-12-31
#> ==================================================
#> values: 2,722,170 (vars*X*Y*T)

# Find zonal "mean"
out2 = execute_zonal(out[[1]], geom = AOI, ID = "fip_code", FUN = "mean")

# Plot
plot(out2[paste0("V", 1:16)])