--- title: "Data Access in R - OceanHackWeek 2022" author: "Johnathan Evanilla" date: "8/10/2022" output: html_document --- ```{r} suppressPackageStartupMessages({ # from CRAN library(dplyr) library(sf) library(stars) library(readr) library(rerddap) library(ncdf4) library(leaflet) library(ggplot2) # from Github library(xyzt) library(ghrsst) }) ``` ## 1) ERDDAP [rerddap](https://docs.ropensci.org/rerddap/index.html) ### Get a listing of servers that rerddap knows about ```{r} s <- servers() s ``` ```{r} s |> filter(short_name == "OOI") ``` ### What services are available in the server? ```{r} ds <- ed_datasets(url="http://erddap.dataexplorer.oceanobservatories.org/erddap/") paste("We have", sum(ds$tabledap!=""), "tabledap", sum(ds$griddap!=""), "griddap, and", sum(ds$wms!=""), "wms.", sep=" ") ``` ### Let's query all the datasets that have the _standard_name_ of _sea_water_practical_salinity_. ```{r} out <- ed_search(query = "sea_water_practical_salinity", which = "tabledap", url="http://erddap.dataexplorer.oceanobservatories.org/erddap/") out ``` ```{r} paste("Found", length(out$alldata), "datasets", sep=" ") ``` ### Let us narrow our search to deployments that are within a lon/lat/time extent. [leaflet](https://rstudio.github.io/leaflet/) ```{r} lon <- c(-72, -69) lat <- c(38, 41) leaflet() |> addTiles() |> setView(lng=sum(lon)/2, lat=sum(lat)/2, zoom=6) |> addRectangles(lng1 = lon[1], lng2 = lon[2], lat1 = lat[1], lat2 = lat[2]) ``` ```{r} out_bb <- ed_search_adv(standard_name= "sea_water_practical_salinity", which = "tabledap", url="http://erddap.dataexplorer.oceanobservatories.org/erddap/", minLat=38, maxLat=41, minLon=-72, maxLon=-69, minTime = "2016-07-10T00:00:00Z", maxTime = "2017-02-10T00:00:00Z") dataset_ids <- out_bb$info$dataset_id out_bb ``` ```{r} paste("Found", length(out_bb$alldata), "Datasets:", sep=" ") dataset_ids ``` ```{r} dataset_ids[1] ``` ```{r} info(dataset_ids[1], url="http://erddap.dataexplorer.oceanobservatories.org/erddap/") ``` ```{r} zz <- tabledap(x = dataset_ids[1], fields = c("z", "latitude", "longitude", "sea_water_practical_salinity_profiler_depth_enabled", "sea_water_temperature_profiler_depth_enabled", "time"), url = "http://erddap.dataexplorer.oceanobservatories.org/erddap") ``` ```{r} head(zz) ``` ```{r} summary(zz) ``` ```{r} zz_sub <- zz |> tibble() |> mutate(time = as.POSIXct(time), z = as.numeric(z), temperature = as.numeric(sea_water_temperature_profiler_depth_enabled), salinity = as.numeric(sea_water_practical_salinity_profiler_depth_enabled)) |> select(-sea_water_temperature_profiler_depth_enabled, -sea_water_practical_salinity_profiler_depth_enabled) |> filter(between(time, as.POSIXct("2016-01-01 00:00:00"), as.POSIXct("2016-12-31 00:00:00")), z == min(z)) ``` ```{r} summary(zz_sub) ``` ```{r} plot(x=zz_sub$time, y=zz_sub$salinity) ``` ```{r} plot(x=zz_sub$salinity, y=zz_sub$temperature) ``` ```{r} ggplot(data=zz_sub, aes(x=salinity, y=temperature)) + geom_point() ``` ## 2) OPeNDAP [NetCDF4 Cheat Sheet](https://www.r-bloggers.com/2016/08/a-netcdf-4-in-r-cheatsheet/) ### Bigelow Homegrown R Packages [Intro Video](https://youtu.be/2ZWFYWMwq1c) [xyzt](https://github.com/BigelowLab/xyzt) [ghrsst](https://github.com/BigelowLab/ghrsst) ```{r} url <- "https://opendap.jpl.nasa.gov/opendap/OceanTemperature/ghrsst/data/GDS2/L4/GLOB/JPL/MUR/v4.1/2022/221/20220809090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc" X <- nc_open(url) ``` ### Retrieve the variables contained in the netcdf object ```{r} mur_vars(X) ``` ### Retrieve the spatial resolution stated in the global metadata ```{r} mur_res(X) ``` ## Read in CalCOFI Station Data ```{r} x <- xyzt::read_calcofi() x ``` ## Change x to be an sf object containing points ```{r} points <- x |> as_POINT() ``` ```{r} covars <- ghrsst::extract(points, X, varname = "analysed_sst") (y <- dplyr::bind_cols(x, covars)) y ``` ```{r} summary(y) ``` ## Change x to be an sf object containing a bounding box ```{r} bbox <- x |> as_BBOX() ``` ### Extract the bounding box ```{r} layer <- ghrsst::extract(bbox, X, varname = "analysed_sst") layer ``` ### Plot the layer we extracted and the points on top ```{r} par(mfrow = c(1,2)) plot(layer, attr = 'analysed_sst', col = sf.colors(n=6), axes = TRUE, reset = FALSE) plot(sf::st_geometry(points), add = TRUE, col = "black", pch = 19, cex = 1) ``` ## 3) WMS ### Learning Objectives: - Add a WMS layer to an interactive map. ("GEBCO - General Bathymetric Chart of the Oceans") ```{r} leaflet() |> addTiles() |> setView(lng = -88, lat = 30, zoom = 5) |> addWMSTiles( "https://www.gebco.net/data_and_products/gebco_web_services/web_map_service/mapserv?", layers = "GEBCO_latest", options = WMSTileOptions(format = "image/png", transparent = TRUE), attribution="Weather data © 2012 IEM Nexrad") ```