sst-tinker.R

# from https://thredds.nci.org.au/thredds/catalog/gb6/BRAN/BRAN2020/daily/catalog.html?dataset=gb6/BRAN/BRAN2020/daily/ocean_temp_1995_06.nc
library(tidyverse)
library(terra)
nci_url <- "https://thredds.nci.org.au/thredds/fileServer/gb6/BRAN/BRAN2020/daily/ocean_temp_1995_06.nc"
vsi <- function(url) glue::glue("/vsicurl/{url}")
sf::gdal_utils("info", vsi(nci_url))

ncdf4::nc_open()

# SST data of Antarctica
# can we start at a station
# can we scale that up

# SST from now (ish) until 1st Sept 1981 every day
# template the year, month, actual day
# always same format
sst_url <- "https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/202408/oisst-avhrr-v02r01.20240801.nc"
curl::curl_download(sst_url, basename(sst_url))

sst_file <- "oisst-avhrr-v02r01.20240801.nc"

subdata <- function(filename,
                    varname = NULL,
                    driver = "NetCDF") {
  stopifnot(!is.null(varname))
  glue::glue('{driver}:"{filename}":{varname}')
}

subdata(sst_file, "sst")
sst_aug <- subdata(sst_file, "sst") |> rast()
plot(sst_aug)
sst_aug <- rast(subdata(vsi(sst_url)))
date <- as.Date("2024-01-01")
aug <- rast("~/Downloads/oisst-avhrr-v02r01.20240801.nc")
aug

"20240801"

underway_url <- "https://github.com/mdsumner/nuyina.underway/raw/main/data-raw/nuyina_underway.parquet"

curl::curl_download(underway_url, basename(underway_url))

data_underway <- nanoparquet::read_parquet(
  basename(underway_url)
)

data_underway$datetime - 1
library(tidyverse)
past_week_of_data <- data_underway |> 
  filter(between(as.Date(datetime), date - 7, date + 7 )
  )

data_underway |> 
  filter(between(as.Date(datetime), date - 7, date + 7 )
  )

data_underway$verysouth <- data_underway$latitude < -60
data_underway$segment <- c(0, cumsum(abs(diff(data_underway$verysouth))))
                           

data_underway$segment[!data_underway$verysouth] <- NA
plot(data_underway$datetime, data_underway$latitude, pch = ".", col = sample(palr::d_pal(data_underway$segment)))

plot(data_underway$datetime, data_underway$sea_water_temperature, pch = ".")



library(duckdb)
dbConnect(duckdb(), "nuyina_underway.parquet")

code as start to comparing sat to underway (will be interesting to look at how the pitch/roll and temperature seen by the ship compares minute by minute to how the daily remotes sensing data looks)

underway_url <- "https://github.com/mdsumner/nuyina.underway/raw/main/data-raw/nuyina_underway.parquet"
nuy <-arrow::read_parquet(underway_url) |> dplyr::filter(as.Date(datetime) > "2022-09-01")
 
 
tst <- nuy$latitude < -55
nuy$seg <- c(0, (cumsum(abs(diff(tst)))))
nuy$seg[!tst]  <- NA
nuy$seg <- as.integer(factor(nuy$seg))
library(ggplot2)

ggplot(nuy |> dplyr::filter(!is.na(seg)), aes(datetime, latitude, col = factor(seg))) + geom_point(pch = ".")

(dt <- as.Date(range(nuy$datetime[nuy$seg == 1], na.rm = T)) + c(0, 1))
library(dplyr)
voyage <- nuy |> dplyr::filter(between(as.Date(datetime), dt[1], dt[2]))

## raadtools lookup   (here Mike should cache this lookup for every Nuyina point)
library(raadtools)
voyage$iceconc <- extract(readice, voyage |> dplyr::select(longitude, latitude, datetime), setNA = T)
voyage$oisst <- extract(readsst, voyage |> dplyr::select(longitude, latitude, datetime))
voyage$ghrsst <- extract(readghrsst, voyage |> dplyr::select(longitude, latitude, datetime))

par(mfrow = c(2, 1))

plot(voyage$datetime, voyage$ghrsst)
lines(voyage$datetime, voyage$sea_water_temperature)


plot(voyage$datetime, voyage$oisst)
lines(voyage$datetime, voyage$sea_water_temperature)


ggplot(voyage, aes(datetime, iceconc)) + geom_path()

dateice <- voyage$datetime[min(which(voyage$iceconc > 0))] 
dateice

sst_url <- glue::glue("https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/{format(dateice, '%Y%m')}/oisst-avhrr-v02r01.{format(dateice, '%Y%m%d')}.nc")
sst_file <- basename(sst_url)
if (!file.exists(sst_file)) curl::curl_download(sst_url, sst_file)



subdata <- function(filename,
                    varname = NULL,
                    driver = "NetCDF") {
  stopifnot(!is.null(varname))
  glue::glue('{driver}:"{filename}":{varname}')
}
#vsi <- function(url) glue::glue("/vsicurl/{url}")

sstraster <- subdata(sst_file, "sst") |> rast()
iceraster <- subdata(sst_file, "ice") |> rast()

plot(voyage$longitude, voyage$latitude, asp = 1/cos(mean(voyage$latitude) * pi/180), type = "n")
plot(iceraster, add = TRUE)
lines(voyage$longitude, voyage$latitude)
points(cbind(voyage$longitude, voyage$latitude)[findInterval(as.POSIXct(dateice, tz = "UTC"), voyage$datetime),, drop = F ], pch = "X", col = "hotpink")

plot(voyage$datetime, voyage$sea_water_temperature, pch = ".")
abline(v = dateice, h = -1.8)

plot(voyage$datetime, pmax(voyage$platform_pitch_fore_up, voyage$platform_heave_down, voyage$platform_roll_starboard_down), pch = ".")
abline(v = (dateice) + c(-1, 0, 1) * 24 * 3600)

Thanks, @mdsumner !

I'm getting this error for raadtools, do I need to install something? Or is this perhaps a GDAL thing?

> library(raadtools)
Loading required package: raster
Loading required package: sp

Attaching package: ‘raster’

The following object is masked from ‘package:dplyr’:

    select

> voyage$iceconc <- extract(readice, voyage |> dplyr::select(longitude, latitude, datetime), setNA = T)
no files found in the 'raadfiles.filename.database'
and no root directories found.
Error in .find_files_generic(pattern) : no files found

oh sorry, you can't install raadtools - you have to be on our system/s for now - which I meant to do. I'm just keeping notes here, but I'll make you an account and email deets :)