Skip to content

Instantly share code, notes, and snippets.

@mdsumner
Created July 21, 2025 23:26
Show Gist options
  • Save mdsumner/7b884d969319acfa71124a5e3c5211cb to your computer and use it in GitHub Desktop.
Save mdsumner/7b884d969319acfa71124a5e3c5211cb to your computer and use it in GitHub Desktop.

https://bsky.app/profile/mdsumner.bsky.social/post/3luj4y4apos2k

  library(sooty)  ## remotes::install_github("mdsumner/sooty") ## for (copies of) the NetCDF sources
library(sds)    ## remotes::install_github("hypertidy/sds")  ## for palettized image urls
library(stringr)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(terra)
#> terra 1.8.61
library(purrr)
july25 <- sooty_files(F) |> dplyr::filter(str_detect(Key, "SEAICE_PS_S25km_202507")) |> 
  transmute(source = file.path(Protocol, Host, Bucket, Key))
july25$date <- as.Date(str_extract(basename(july25$source), "[0-9]{8}"), "%Y%m%d")


## get palettized images by date, we'll safely try/catch and cull to those that work
date <- seq(as.Date("2025-07-01"), Sys.Date(), by = "1 day")
files <- map_chr(date, \(.x) sds::nsidc_seaice(.x))

rs <- map(files, safely(terra::rast))
#> Warning: HTTP response code: 404 (GDAL error 11)
#> Warning: HTTP response code: 404 (GDAL error 11)
ok <- map_lgl(rs, \(.x) is.null(.x$error))

r <- rast(map(rs[ok], \(.x) .x$result))
r
#> class       : SpatRaster 
#> size        : 332, 316, 19  (nrow, ncol, nlyr)
#> resolution  : 25000, 25000  (x, y)
#> extent      : -3950000, 3950000, -3950000, 4350000  (xmin, xmax, ymin, ymax)
#> coord. ref. : NSIDC Sea Ice Polar Stereographic South (EPSG:3412) 
#> sources     : S_20250701_concentration_v3.0.tif  
#>               S_20250702_concentration_v3.0.tif  
#>               S_20250703_concentration_v3.0.tif  
#>               ... and 16 more sources
#> color table : 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19 
#> names       : S_202~_v3.0, S_202~_v3.0, S_202~_v3.0, S_202~_v3.0, S_202~_v3.0, S_202~_v3.0, ...
plot(r, maxnl = nlyr(r))

## get the netcdf sources, take the mean
r2 <- rast(map(july25$source, \(.x) {x <- terra::rast(.x); mean(x, na.rm = TRUE)}))
r2 <- setNames(r2, july25$date)
plot(r2, maxnl = nlyr(r2))

Created on 2025-07-21 with reprex v2.0.2

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment