## ----initialize, echo = FALSE------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ## ----setup, include = FALSE--------------------------------------------------- ### needed libraries library(cmocean) library(mapdata) library(plotdap) library(rerddap) library(rerddapXtracto) ## ----install, eval = FALSE---------------------------------------------------- # install.packages("ncdf4", dependencies = TRUE) # install.packages("parsedate", dependencies = TRUE) # install.packages("plotdap", dependencies = TRUE) # install.packages("rerddap", dependencies = TRUE) # install.packages("sp", dependencies = TRUE) ## ----install_package, eval = FALSE-------------------------------------------- # install.packages("rerddapXtracto", dependencies = TRUE) ## ----installGit, eval = FALSE------------------------------------------------- # install.packages("devtools") # remotes::install_github("rmendels/rerddapXtracto", subdir = 'development') ## ----eval = FALSE------------------------------------------------------------- # library("rerddapXtracto") ## ----eval = FALSE------------------------------------------------------------- # library("gganimate") # library("ggplot2") # library("plotdap") ## ----info, echo = TRUE, eval = FALSE------------------------------------------ # require("rerddap") # ## base URL does not need to given because it is the default one # dataInfo <- info('erdMBchla1day') # dataInfo ## ----dataInfo, echo = TRUE, eval = FALSE-------------------------------------- # erdMBchla1day # Base URL: https://upwell.pfeg.noaa.gov/erddap/ # Dimensions (range): # time: (2006-01-01T12:00:00Z, 2020-10-27T12:00:00Z) # altitude: (0.0, 0.0) # latitude: (-45.0, 65.0) # longitude: (120.0, 320.0) # Variables: # chlorophyll: # Units: mg m-3 # ## ----getMarlinChl, eval = FALSE, echo = TRUE---------------------------------- # require("rerddap") # require("rerddapXtracto") # # # First we will copy the Marlintag38606 data into a variable # # called tagData so that subsequent code will be more generic. # # tagData <- Marlintag38606 # xpos <- tagData$lon # ypos <- tagData$lat # tpos <- tagData$date # zpos <- rep(0., length(xpos)) # swchlInfo <- rerddap::info('erdSWchla8day') # swchl1 <- rxtracto(swchlInfo, parameter = 'chlorophyll', # xcoord = xpos, ycoord = ypos, tcoord = tpos, zcoord = zpos, # xlen = .2, ylen = .2, progress_bar = TRUE) ## ----meantrackPlot, eval = FALSE, echo = TRUE--------------------------------- # require("ggplot2") # require("plotdap") # # myPlot <- plotTrack(swchl1, xpos, ypos, tpos, plotColor = 'algae') # myPlot ## ----animateTrack, echo = TRUE, eval = FALSE---------------------------------- # myPlot <- plotTrack(swchl1, xpos, ypos, tpos, plotColor = 'algae', # animate = TRUE, cumulative = TRUE) # ## ----tidy_track, echo = TRUE, eval = FALSE------------------------------------ # # swchl1_tidy <- as.data.frame(swchl1) # ## ----topotag, eval = FALSE, echo = TRUE--------------------------------------- # require("ggplot2") # require("plotdap") # require("rerddap") # require("rerddapXtracto") # ylim <- c(15, 30) # xlim <- c(-160, -105) # topoInfo <- rerddap::info('etopo360') # topo <- rxtracto(topoInfo, parameter = 'altitude', xcoord = xpos, ycoord = ypos, xlen = .1, ylen = .1) # myFunc = function(x) -x # topoPlot <- plotTrack(topo, xpos, ypos, NA, plotColor = 'dense', # name = 'Depth', myFunc = myFunc) # topoPlot ## ----extract3D, echo = TRUE, eval = FALSE------------------------------------- # require("rerddap") # urlBase <- "https://erddap.marine.ie/erddap/" # parameter <- "Sea_water_temperature" # dataInfo <- rerddap::info("IMI_CONN_3D", url = urlBase) # #get the actual last 3 times, and extract from data frame # dataInfo1 <- read.csv("https://erddap.marine.ie/erddap/griddap/IMI_CONN_3D.csv0?time[last-2:1:last]", stringsAsFactors = FALSE, header = FALSE, row.names = NULL) # sstTimes <- dataInfo1[[1]] # sstLats <- c(53.505758092414446, 53.509303546859805, 53.51284900130517) # sstLons <- c(-10.25975390624996, -10.247847656249961, -10.23594140624996) # sstDepths <- c(2, 6, 10) # sstTrack <- rxtracto(dataInfo, parameter = parameter, xcoord = sstLons, ycoord = sstLats, tcoord = sstTimes, zcoord = sstDepths, xlen = .05, ylen = .05, zlen = 0., zName = 'altitude') # str(sstTrack) ## ----extract3D_struct, echo = TRUE, eval = FALSE------------------------------ # List of 13 # $ mean Sea_water_temperature : num [1:3] 11.7 11.7 11.6 # $ stdev Sea_water_temperature : num [1:3] 0.0407 0.0555 0.0859 # $ n : int [1:3] 493 491 484 # $ satellite date : chr [1:3] "2020-10-31T22:00:00Z" "2020-10-31T23:00:00Z" "2020-11-01T00:00:00Z" # $ requested lon min : num [1:3] -10.3 -10.3 -10.3 # $ requested lon max : num [1:3] -10.2 -10.2 -10.2 # $ requested lat min : num [1:3] 53.5 53.5 53.5 # $ requested lat max : num [1:3] 53.5 53.5 53.5 # $ requested z min : num [1:3] 2 6 10 # $ requested z max : num [1:3] 2 6 10 # $ requested date : chr [1:3] "2020-10-31T22:00:00Z" "2020-10-31T23:00:00Z" "2020-11-01T00:00:00Z" # $ median Sea_water_temperature: num [1:3] 11.7 11.7 11.6 # $ mad Sea_water_temperature : num [1:3] 0.0208 0.0361 0.0986 # - attr(*, "row.names")= chr [1:3] "1" "2" "3" # - attr(*, "class")= chr [1:2] "list" "rxtractoTrack" # ## ----dateline_track, echo = TRUE, eval = FALSE------------------------------- # dataInfo <- rerddap::info('jplMURSST41mday') # parameter <- 'sst' # xcoord <- c(179.7, 179.8, 179.9, 180., 180.1, 180.2, 180.3, 180.4) # ycoord <- c(40, 40, 40, 40, 40, 40, 40, 40) # tcoord <- c('2018-03-16', '2018-03-16', '2018-03-16','2018-03-16','2018-03-16', # '2018-03-16','2018-03-16','2018-03-16') # xlen <- .05 # ylen <- .05 # extract <- rxtracto(dataInfo, parameter = parameter, xcoord = xcoord, # ycoord = ycoord, tcoord = tcoord, # xlen = xlen, ylen = ylen) # str(extract) ## ----dateline_track_struct, echo = TRUE, eval = FALSE------------------------ # List of 13 # $ mean sst : num [1:8] 11.1 11.1 11.1 11.2 11.1 ... # $ stdev sst : num [1:8] 0.01192 0.00602 0.01025 0.00876 0.01446 ... # $ n : int [1:8] 30 30 35 25 30 30 30 35 # $ satellite date : chr [1:8] "2018-03-16T00:00:00Z" "2018-03-16T00:00:00Z" "2018-03-16T00:00:00Z" "2018-03-16T00:00:00Z" ... # $ requested lon min: num [1:8] 180 180 180 180 180 ... # $ requested lon max: num [1:8] 180 180 180 180 180 ... # $ requested lat min: num [1:8] 40 40 40 40 40 ... # $ requested lat max: num [1:8] 40 40 40 40 40 ... # $ requested z min : logi [1:8] NA NA NA NA NA NA ... # $ requested z max : logi [1:8] NA NA NA NA NA NA ... # $ requested date : chr [1:8] "2018-03-16" "2018-03-16" "2018-03-16" "2018-03-16" ... # $ median sst : num [1:8] 11.1 11.1 11.1 11.2 11.1 ... # $ mad sst : num [1:8] 0.01416 0.0052 0.01149 0.00887 0.01744 ... # - attr(*, "row.names")= chr [1:8] "1" "2" "3" "4" ... # - attr(*, "class")= chr [1:2] "list" "rxtractoTrack" # ## ----echo = TRUE, eval = FALSE----------------------------------------------- # interp <- c('interpolation method', 'number of neighbors') ## ----echo = TRUE, eval = FALSE----------------------------------------------- # interp <- c('Mean', '16') ## ----echo = TRUE, eval = FALSE----------------------------------------------- # my_interp <- c('Mean', '16') # # # # use the coastwatch ERDDAP as it is a more recent version # # # swchlInfo <- rerddap::info('erdSWchla8day', # url = 'https://coastwatch.pfeg.noaa.gov/erddap/') # swchl1 <- rxtracto(swchlInfo, parameter = 'chlorophyll', # xcoord = xpos, ycoord = ypos, tcoord = tpos, zcoord = zpos, # interp = my_interp, progress_bar = TRUE) # ## ----echo = TRUE, eval = FALSE----------------------------------------------- # PB_Argos_subset <- PB_Argos[1:100, ] # # datetime is not in format for ERDDAP # PB_Argos_subset$DateTime <- lubridate::as_datetime(PB_Argos_subset$DateTime) # head(PB_Argos_subset, 4) # # # A tibble: 4 × 4 # DateTime QualClass Lat Lon # # 1 2009-04-20 17:01:39 B 70.4 -132. # 2 2009-04-20 17:23:00 A 71.0 -131. # 3 2009-04-20 18:12:15 A 71.0 -131. # 4 2009-04-20 20:43:17 A 70.9 -131. # ## ----echo = TRUE, eval = FALSE----------------------------------------------- # url= 'https://polarwatch.noaa.gov/erddap/' # datasetid = 'nsidcG02202v4nh1day' # dataInfo <- rerddap::info(datasetid, url) # proj_crs_code_index <- which(dataInfo$alldata$NC_GLOBAL$attribute_name == "proj_crs_code" ) # proj_crs_code <- dataInfo$alldata$NC_GLOBAL$value[proj_crs_code_index] # proj_crs_code # [1] "EPSG:3411" ## ----echo = TRUE, eval = FALSE----------------------------------------------- # temp_df <- data.frame(Lat = PB_Argos_subset$Lat, Lon = PB_Argos_subset$Lon) # # transform PB_Argos_subset to sf object # # EPSG:4326 is basic lat-lon coordinates # temp_df <- sf::st_as_sf(temp_df, coords = c("Lon", "Lat"), crs = 4326) # # project data # temp_df <- sf::st_transform(temp_df, crs = 3411) # # get projection coordinates # coordinates <- sf::st_coordinates(temp_df) ## ----echo = TRUE, eval = FALSE----------------------------------------------- # xcoord <- coordinates[, 1] # ycoord <- coordinates[, 2] # # R datetimes are passed as number, require ISO character string # tcoord <- as.character(PB_Argos_subset$DateTime) # parameter <- 'cdr_seaice_conc' # extract <- rxtracto(dataInfo, # xName="xgrid", # yName="ygrid", # tName="time", # parameter=parameter, # xcoord = xcoord, # ycoord = ycoord, # tcoord = tcoord # ) # str(extract) # List of 13 # $ mean cdr_seaice_conc : num [1:100] 1 1 1 1 1 1 1 1 1 1 ... # $ stdev cdr_seaice_conc : num [1:100] NA NA NA NA NA NA NA NA NA NA ... # $ n : int [1:100] 1 1 1 1 1 1 1 1 1 1 ... # $ satellite date : chr [1:100] "2009-04-21T00:00:00Z" "2009-04-21T00:00:00Z" "2009-04-21T00:00:00Z" "2009-04-21T00:00:00Z" ... # $ requested x min : num [1:100] -2143958 -2077605 -2077124 -2082576 -2036149 ... # $ requested x max : num [1:100] -2143958 -2077605 -2077124 -2082576 -2036149 ... # $ requested y min : num [1:100] -118477 -131404 -131955 -130368 -120720 ... # $ requested y max : num [1:100] -118477 -131404 -131955 -130368 -120720 ... # $ requested z min : logi [1:100] NA NA NA NA NA NA ... # $ requested z max : logi [1:100] NA NA NA NA NA NA ... # $ requested date : chr [1:100] "2009-04-20 17:01:39" "2009-04-20 17:23:00" "2009-04-20 18:12:15" "2009-04-20 20:43:17" ... # $ median cdr_seaice_conc: num [1:100] 1 1 1 1 1 1 1 1 1 1 ... # $ mad cdr_seaice_conc : num [1:100] 0 0 0 0 0 0 0 0 0 0 ... # - attr(*, "row.names")= chr [1:100] "1" "2" "3" "4" ... # - attr(*, "class")= chr [1:2] "list" "rxtractoTrack" # ## ----echo = TRUE, eval = FALSE----------------------------------------------- # xgrid <- 'requested x min' # ygrid <- 'requested y min' # temp_df <- data.frame(xgrid = extract[[xgrid]], ygrid = extract[[ygrid]]) # temp_df <- sf::st_as_sf(temp_df, coords = c("xgrid", "ygrid"), crs = 3411) # temp_df <- sf::st_transform(temp_df, crs = 4326) # coordinates <- sf::st_coordinates(temp_df) # ## ----VIIRSchla, echo = TRUE, eval = FALSE------------------------------------- # require("rerddap") # require("rerddapXtracto") # # xpos <- c(-125, -120) # ypos <- c(39, 36) # tpos <- c("last", "last") # tpos <- c("2017-04-15", "2017-04-15") # VIIRSInfo <- rerddap::info('erdVH3chlamday') # VIIRS <- rxtracto_3D(VIIRSInfo, parameter = 'chla', xcoord = xpos, ycoord = ypos, tcoord = tpos) ## ----VIIRSLogPlot, echo = TRUE, eval = FALSE---------------------------------- # require("ggplot2") # require("plotdap") # myFunc <- function(x) log(x) # chlalogPlot <- plotBBox(VIIRS, plotColor = 'algae', myFunc = myFunc) # chlalogPlot ## ----viirs_tidy, echo = TRUE, eval = FALSE----------------------------------- # VIIRS_tidy <- tidy_grid(VIIRS) ## ----viirs_tidy_str, echo = TRUE, eval = FALSE------------------------------- # 'data.frame': 8833 obs. of 4 variables: # $ time : POSIXlt, format: "2017-04-15" "2017-04-15" ... # $ latitude : num [1:8833(1d)] 36 36.1 36.1 36.1 36.2 ... # $ longitude: num [1:8833(1d)] -125 -125 -125 -125 -125 ... # $ chla : num 0.194 0.193 0.186 0.18 0.163 ... # ## ----dateline_3D, echo = TRUE, eval = FALSE---------------------------------- # dataInfo <- rerddap::info('jplMURSST41mday') # parameter <- 'sst' # xcoord <- c(175, 185) # ycoord <- c(40, 50) # tcoord <- c('2019-03-16', '2019-03-16') # mur_dateline <- rxtracto_3D(dataInfo, parameter, xcoord = xcoord, ycoord = ycoord, # tcoord = tcoord) ## ----world2hires, echo = TRUE, eval = FALSE----------------------------------- # xlim <- c(170, 190) # ylim <- c(40, 55) # remove <- c("UK:Great Britain", "France", "Spain", "Algeria", "Mali", "Burkina Faso", "Ghana", "Togo") # w <- map("world2Hires", xlim = xlim, ylim = ylim, fill = TRUE, plot = FALSE) # w <- map("mapdata::world2Hires", regions = w$names[!(w$names %in% remove)], plot = FALSE, fill = TRUE, ylim = ylim, xlim = xlim) # ## ----world2hires_map, echo = TRUE, eval = FALSE------------------------------- # mapFrame <- function(longitude, latitude, my_data) { # my_data_name <- names(my_data) # temp_data <- drop(my_data[[1]]) # dims <- dim(temp_data) # temp_data <- array(temp_data, dims[1] * dims[2]) # my_frame <- expand.grid(x = longitude, y = latitude) # my_frame[my_data_name] <- temp_data # return(my_frame) # } # mur_frame <- mapFrame(mur_dateline$longitude, mur_dateline$latitude, mur_dateline['sst']) # mycolor <- cmocean::cmocean('thermal')(256) # myplot <- ggplot(data = mur_frame, aes(x = x, y = y, fill = sst)) + # geom_polygon(data = w, aes(x = long, y = lat, group = group), fill = "grey80") + geom_tile(interpolate = FALSE) + # scale_fill_gradientn(colours = mycolor, na.value = NA) + # theme_bw() + ylab("latitude") + xlab("longitude") + # coord_fixed(1.3, xlim = xlim, ylim = ylim) # myplot ## ----mbnmsChla, echo = TRUE, eval = FALSE------------------------------------- # require("rerddapXtracto") # dataInfo <- rerddap::info('erdVH3chlamday') # parameter = 'chla' # tpos <- c("2014-09-01", "2014-10-01") # #tpos <-as.Date(tpos) # xpos <- mbnms$Longitude # ypos <- mbnms$Latitude # sanctchl <- rxtractogon(dataInfo, parameter = parameter, xcoord = xpos, ycoord = ypos, tcoord = tpos) # str(sanctchl) ## ----mbnmsChla_struct, echo = TRUE, eval = FALSE------------------------------ # List of 6 # $ chla : num [1:50, 1:57, 1:2] NA NA NA NA NA NA NA NA NA NA ... # $ datasetname: chr "erdVH3chlamday" # $ longitude : num [1:50(1d)] -123 -123 -123 -123 -123 ... # $ latitude : num [1:57(1d)] 35.6 35.6 35.6 35.7 35.7 ... # $ altitude : logi NA # $ time : POSIXlt[1:2], format: "2014-09-15" "2014-10-15" # - attr(*, "class")= chr [1:2] "list" "rxtracto3D" # ## ----mbnmsChlaPlot, echo = TRUE, eval = FALSE-------------------------------- # require("ggplot2") # require("plotdap") # myFunc <- function(x) log(x) # sanctchl1 <- sanctchl # sanctchl1$chla <- sanctchl1$chla[, , 2] # sanctchl1$time <- sanctchl1$time[2] # sanctchlPlot <- plotBBox(sanctchl1, plotColor = 'algae', myFunc = myFunc) # sanctchlPlot ## ----animate, eval = FALSE---------------------------------------------------- # require("gganimate") # #> Loading required package: gganimate # require("ggplot2") # require("plotdap") # myFunc <- function(x) log(x) # sanctchlPlot <- plotBBox(sanctchl, plotColor = 'algae', myFunc = myFunc, time = identity, animate = TRUE) ## ----mbnmsBathy, echo = TRUE, eval = FALSE----------------------------------- # require("rerddap") # dataInfo <- rerddap::info('etopo180') # xpos <- mbnms$Longitude # ypos <- mbnms$Latitude # bathy <- rxtractogon(dataInfo, parameter = 'altitude', xcoord = xpos, ycoord = ypos) # str(bathy) ## ----mbnmsBathy_struct, echo = TRUE, eval = FALSE---------------------------- # List of 6 # $ depth : num [1:123, 1:141, 1] NA NA NA NA NA NA NA NA NA NA ... # $ datasetname: chr "etopo180" # $ longitude : num [1:123(1d)] -123 -123 -123 -123 -123 ... # $ latitude : num [1:141(1d)] 35.5 35.6 35.6 35.6 35.6 ... # $ altitude : logi NA # $ time : logi NA # - attr(*, "class")= chr [1:2] "list" "rxtracto3D" # ## ----mbnmsBathyPlot, echo = TRUE, eval = FALSE-------------------------------- # require("ggplot2") # require("mapdata") # myFunc = function(x) -x # bathyPlot <- suppressMessages((plotBBox(bathy, plotColor = 'dense', myFunc = myFunc, name = 'Depth'))) # bathyPlot ## ----soda70, echo = TRUE, eval = FALSE--------------------------------------- # require("rerddap") # dataInfo <- rerddap::info('erdSoda331oceanmday') # xpos <- c(185.25, 240.25) # ypos <- c(20.25, 60.25) # zpos <- c(76.80285, 76.80285) # tpos <- c('2010-12-15', '2010-12-15') # soda70 <- rxtracto_3D(dataInfo, parameter = 'temp', xcoord = xpos, ycoord = ypos, tcoord = tpos, zcoord = zpos, zName = 'depth') ## ----soda70Plot, echo = TRUE, eval = FALSE------------------------------------ # require("ggplot2") # require("plotdap") # sodaPlot <- plotBBox(soda70, plotColor = 'thermal', name = 'temp_at_70m', maxpixels = 30000) # sodaPlot # ## ----NAtlSSS, eval = FALSE, echo = TRUE--------------------------------------- # require("rerddap") # urlBase <- "https://erddap.marine.ie/erddap/" # parameter <- "sea_surface_salinity" # sssTimes <- c("last", "last") # sssLats <- c(48.00625, 57.50625) # sssLons <- c(-17.99375, -1.00625) # dataInfo <- rerddap::info("IMI_NEATL", url = urlBase) # NAtlSSS <- rxtracto_3D(dataInfo, parameter = parameter, xcoord = sssLons, ycoord = sssLats, tcoord = sssTimes) # ## ----NAtlSSSplot, eval = FALSE, echo = TRUE----------------------------------- # require("ggplot2") # require("plotdap") # NAtlSSSPlot <- plotBBox(NAtlSSS, plotColor = 'haline', name = "salinity", maxpixels = 30000) # NAtlSSSPlot ## ----NAtlSSSplot1, eval = FALSE, echo = TRUE---------------------------------- # require("ggplot2") # require("plotdap") # haline = cmocean::cmocean('haline')(256) # add_ggplot(NAtlSSSPlot, scale_colour_gradientn(colours = haline, na.value = NA, limits = c(32, 36)), scale_fill_gradientn(colours = haline, na.value = NA, limits = c(32, 36))) ## ----nearGrid, eval = FALSE--------------------------------------------------- # latitude[which.min(abs(latitude - ypos[1]))] # minimum latitude # latitude[which.min(abs(latitude - ypos[2]))] # maximum latitude # longitude[which.min(abs(longitude- xpos[1]))] # minimum longitude # longitude[which.min(abs(longitude - xpos[2]))] # maximum longitude # isotime[which.min(abs(time - tpos[1]))] # minimum time # isotime[which.min(abs(time - tpos[2]))] # maximum time