Skip to content

Instantly share code, notes, and snippets.

@mdsumner
Last active June 2, 2024 23:50
Show Gist options
  • Save mdsumner/6ebd03458caa62080ba14e8a04536e38 to your computer and use it in GitHub Desktop.
Save mdsumner/6ebd03458caa62080ba14e8a04536e38 to your computer and use it in GitHub Desktop.

Crop GHRSST in R terra.

I've ignored all our detailed stuff for our COG creation - this crops into memory, converts to Celsius, writes to file (I've set up a place to set the output directory but it writes to working directory atm) - I tested from compute in California and it was fast. With that you should be able to generate the days you want with seq(as.Date(), by = "1 day") kind of construct and loop. Note that some days are missing, we have reported them so maybe fixed now, but it's worth wrapping the rast(dsn) call in a try in case the read fails

## https://urs.earthdata.nasa.gov/documentation/for_users/user_token
## the token is a long text string, the header is "Authorization: Bearer {token}"

## this is your earthdata token in Authorization Bearer format, I keep mine in a local file (its filename is not relevant)
earthdata <- normalizePath("~/earthdata")
Sys.setenv("GDAL_HTTP_HEADER_FILE" = earthdata)
## could also set GDAL_HTTP_HEADERS to the header itself (don't need a file)
## or use terra::setGDALconfig(), I just like env vars because they work for all GDAL software, and a file
## because the code is purely relative to my local set up and will work


## customize for where you want the files to go
root <- "."
#dir.create(root, recursive = TRUE)

base <- "https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1"
date <- as.Date("2024-01-01")
url  <- glue::glue("{base}/{format(date, '%Y')}{format(date, '%m')}{format(date, '%d')}090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc")

## now we start with the GDAL syntax, vsi for lazy internet read, DRIVER:url:sds for the subdataset syntax
dsn <- glue::glue("NETCDF:/vsicurl/{url}:analysed_sst")

## now pick your weapon
library(terra)
aoi_ext <- ext(-71.3020833333333, -55.3958333333346, 45.1979166666667, 51.9062499999999)
ghrsst <- rast(dsn)
##to make these files more reliable, first (see https://lists.osgeo.org/pipermail/gdal-dev/2024-June/059052.html)
set.ext(ghrsst, ext(-179.995,180.005,-89.995,89.995)) 

## now crop, convert to C, and write to file
loc <- crop(ghrsst, aoi_ext)
loc <- loc - 273.15
terra::units(loc) <- "celsius"
outfile <- file.path(root, gsub("\\.nc$", "_analysed_sst.tif", basename(url)))
writeRaster(loc, outfile)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment