fold in boundary to just give the cell index from dimension

[mdsumner]

i.e.

https://forum.earthdata.nasa.gov/viewtopic.php?t=6462

https://n5eil01u.ecs.nsidc.org/VIIRS/VNP10A1.002/2025.05.06/

bound <- function(dm) {

    c(vaster::cell_from_row(dm,  1),
         vaster::cell_from_col(dm, dm[1]),
         rev(vaster::cell_from_row(dm, dm[2])),
         rev(vaster::cell_from_col(dm, 1)))

}

f <- "~/S3B_OLCI_ERRNT.20250101T164657.L2.IOP.NRT.nc"
library(terra)
ll <- flip(c(rast(f, "/navigation_data/longitude"), 
  rast(f, "/navigation_data/latitude")), "horizontal")

plot(ll[bound(dim(ll)[2:1])], pch = ".")
maps::map(add = T)
r0 <- rast(f, "/geophysical_data/bbp_443")
set.crs(r0, "EPSG:4326")
plot(project(r0, rast(ext(160, 180, -90, -60), res = .2, crs = "EPSG:4326"), by_util = TRUE, method = "sum"))
plot(project(r0, rast(ext(-180, -140, -80, 0), res = .2, crs = "EPSG:4326"), by_util = TRUE, method = "sum"))

dsn <- sprintf("NetCDF:\"%s\":/geophysical_data/bbp_443", normalizePath(f))
library(gdalraster)
s <- 1.2e7
lon <- seq(-160, -40, by = 5)
#i <- 0
i <- 12
prj <- sprintf("+proj=laea +lon_0=%i +lat_0=-10", lon[i])
w <- warp(dsn, tf <- tempfile(tmpdir = "/vsimem", fileext = ".tif"), prj, cl_arg = c("-te", -s, -s, s, s, "-ts",  1024, 0))
plot_raster(read_ds(new(GDALRaster, tf)))
# plot(reproj::reproj_xy(ll[bound(dim(ll)[2:1])], prj), pch = ".", col = "firebrick")
points(reproj::reproj_xy(ll[bound(dim(ll)[2:1])], prj), pch = ".", col = "firebrick")

m <- reproj::reproj_xy(do.call(cbind, maps::map(plot = F)[1:2]), prj)
lines(m)

[mdsumner]

a few helpers

area <- function(x) {
    area::polygon_area(x)
}

mid <- function(x) {
    matrix(apply(x, 2, mean), ncol = 2L)
}
proj_xyz <- function(x) {
    reproj::reproj_xyz(x, "+proj=geocent", source = "EPSG:4326")
}
proj_llz <- function(x) {
    reproj::reproj_xyz(x, "EPSG:4326", source = "+proj=geocent")
}
midll <- function(p) {
    proj_llz(matrix(apply(proj_xyz(p), 2L, mean), ncol = 3L))[,1:2, drop = F]
}
## treat these as a closure for now
params <- function(x, secant = FALSE) {
    lon <- as.integer(midll(x)[1])
    lat <- as.integer(midll(x)[2])
    if (secant) {
        sec <-  as.integer(range(x[,2]))
        return(list(lon_0 = lon, lat_0 = lat, lat_1 = sec[1], lat_2 = sec[2]))
    } 
    list(lon_0 = lon, lat_0 = lat)
}
params_fold <- function(x) {
    unname(apply(cbind(names(x), unlist(x)), 1, paste0, collapse = "="))
}

proj <- function(params, family   = NULL) {
    if (is.null(family)) stop("family must be named (i.e. laea, stere, lcc")
    params$proj <- family
    nms <- setdiff(names(params), "proj")
    params <- params[c("proj", nms)]
 
 paste0( sprintf("+%s", params_fold(params)), collapse = " ")
}
proj(params(m), "stere")
(prj <- proj(params(m, secant = TRUE), "lcc"))
(prj <- proj(params(m), "ortho"))
# alpha <- seq(-50, 50, by = 5)
# i <- 0
# i <- i + 1
prj <- sprintf("+proj=omerc +lonc=-113 +lat_0=-20 +alpha=%f +gamma=1", 15)
proj_xy <- function(x) {
    reproj::reproj_xy(x, prj, source = "EPSG:4326")
}
proj_ll <- function(x) {
    reproj::reproj_xy(x, "EPSG:4326", source = prj)
}
plot(proj_xy(m), pch = ".", col = "firebrick", asp = 1)
lines(proj_xy(m), col = "firebrick")
points(proj_xy(midll(m)))
mp <- do.call(cbind, maps::map(plot = F)[1:2])
points(proj_xy(mp), pch = ".")


ex <- c(apply(proj_xy(m), 2, range))
dsn <- sprintf("NetCDF:\"%s\":/geophysical_data/bbp_443", normalizePath(f))
library(gdalraster)
w <- warp(dsn, tf <- tempfile(tmpdir = "/vsimem", fileext = ".tif"), prj, cl_arg = c("-te", ex[1], ex[3], ex[2], ex[4], "-ts",  1024, 0))
plot_raster(read_ds(new(GDALRaster, tf)))
image(sqrt(rast(tf)), col = hcl.colors(256))