new stuff

[mdsumner]

lot in progress

and somewhat relevant to: https://github.com/AustralianAntarcticDivision/blueant/issues/12

## full L1B modis workflow

#https://ladsweb.modaps.eosdis.nasa.gov/archive/allData/61/MYD021KM/2023/243
# wget -e robots=off -m -np -R .html,.tmp -nH --cut-dirs=3 "https://ladsweb.modaps.eosdis.nasa.gov/archive/allData/61/MYD021KM/2023/243/" --header "Authorization: Bearer  <see wget from page>" -P .

prj <- "OGC:CRS84"
ex <- c(-180, 180, -90, 90)

sds_l1b <- function(x) {
  unname(vapply(x, function(.x) vapour::vapour_sds_names(.x)[1L], ""))
}

f <- unname(fs::dir_ls("/rdsi/PRIVATE/raad/data_staging/ladsweb.modaps.eosdis.nasa.gov/archive/allData/61/MYD021KM/2023/243", recurse = TRUE, regexp = ".*hdf$"))
sds <- sds_l1b(f)

library(vapour)

im <- gdal_raster_data(sds, target_crs = prj, target_ext = ex, target_dim = c(1256, 0), options = c("-multi", "-wo", "NUM_THREADS=ALL_CPUS"))

prj <- "+proj=laea +lon_0=109 +lat_0=-61"
ex <- c(-1, 1, -1, 1) * 50000
imloc <- gdal_raster_data(sds, target_crs = prj, target_ext = ex, target_dim = c(1256, 0), options = c("-multi", "-wo", "NUM_THREADS=ALL_CPUS"))

dsn::

hdfL1B <- function(x, varname) {
  sprint("HDF4_EOS:EOS_SWATH_GEOL:\"%s\":MODIS_SWATH_Type_L1B:%s", x, varname)
}

HDF4_EOS:EOS_SWATH_GEOL:"/rdsi/PRIVATE/raad/data_staging/ladsweb.modaps.eosdis.nasa.gov/archive/allData/61/MYD021KM/2023/243/MYD021KM.A2023243.0000.061.2023243154235.hdf":MODIS_SWATH_Type_L1B:Longitude
HDF4_EOS:EOS_SWATH_GEOL:"/rdsi/PRIVATE/raad/data_staging/ladsweb.modaps.eosdis.nasa.gov/archive/allData/61/MYD021KM/2023/243/MYD021KM.A2023243.0000.061.2023243154235.hdf":MODIS_SWATH_Type_L1B:Latitude
HDF4_EOS:EOS_SWATH:"/rdsi/PRIVATE/raad/data_staging/ladsweb.modaps.eosdis.nasa.gov/archive/allData/61/MYD021KM/2023/243/MYD021KM.A2023243.0000.061.2023243154235.hdf":MODIS_SWATH_Type_L1B:EV_1KM_RefSB

[mdsumner]

get the boundary something like this, so we can trim which granules to make an image from

get_geol_boundary <- function(file) {
lon <- vapour::vapour_read_raster(hdfL1B(file, "Longitude"), native = T)
lon <- matrix(lon[[1]], 406, byrow = TRUE)
ll <- c(lon[1, ], lon[-1,ncol(lon)], rev(lon[nrow(lon), ]), rev(lon[-c(1, nrow(lon)),1]))
lat <- vapour::vapour_read_raster(hdfL1B(file, "Latitude"), native = T)
lat <- matrix(lat[[1]], 406, byrow = TRUE)
lt <- c(lat[1, ], lat[-1,ncol(lat)], rev(lat[nrow(lat), ]), rev(lat[-c(1, nrow(lat)),1]))
cbind(ll, lt)
}

determine intersections with region of interest thus:

l <- list()
for (i in seq_along(f)) {
 bb <- reproj::reproj_xy(get_geol_boundary(f[i]), prj, source = "OGC:CRS84")
 wk <- wk::wk_polygon(wk::xy(bb[,1], bb[,2]))
 int <- geos::geos_intersection(wk, p)

 if (!geos::geos_is_empty(int)) {
 l <- c(l, list(int))
}
}