mdsumner

this query

https://earth-search.aws.element84.com/v1/search?collections=sentinel-2-c1-l2a&bbox=110.477,-66.332,110.577,-66.232&datetime=2021-11-22T00:00:00Z/2021-11-22T23:59:59Z

has only assets coastal, nir09, preview - other bands are missing.

mdsumner

library(shiny)
library(shinyglide)

ui <- fluidPage(
  glide(
    screen(
      h2("Image 1"),
      img(src = "image1.png", width = "100%")
    ),
    screen(

mdsumner

this appears to be working (on setonix I just had to uv venv .venv and install virtual-tiff and it runs)

import obstore
from virtualizarr.registry import ObjectStoreRegistry
from virtual_tiff import VirtualTIFF
import xarray as xr

# Configuration
endpoint = "https://projects.pawsey.org.au/"

mdsumner

You should update your understanding

don't need anything client side to read a single polygon or field from within a shape file at a URL (or in a zip at a URL). Same goes for parquet

I think you're confused about tools you use.

accuracy matters as does not spreading bullshit

https://bsky.app/profile/postholer.com/post/3lmb6nxh7fs2x

mdsumner

minute check on this PR in GDAL for GetRawBlockInfo

OSGeo/gdal#13260

from osgeo import gdal
gdal.UseExceptions()
ds = gdal.OpenEx("/vsicurl/https://thredds.nci.org.au/thredds/fileServer/gb6/BRAN/BRAN2023/daily/ocean_temp_2010_01.nc", gdal.OF_MULTIDIM_RASTER)
info  = temp.GetRawBlockInfo([0,0,0,0])
info.GetOffset()

mdsumner

<VRTDataset>
  <Group name="/">
    <Dimension name="Time" type="TEMPORAL" size="5479" indexingVariable="Time" />
    <Dimension name="st_ocean" type="VERTICAL" direction="DOWN" size="51" indexingVariable="st_ocean" />
    <Dimension name="xt_ocean" type="HORIZONTAL_X" size="3600" indexingVariable="xt_ocean" />
    <Dimension name="yt_ocean" type="HORIZONTAL_Y" size="1500" indexingVariable="yt_ocean" />
    <Array name="Time">
      <DataType>Float64</DataType>
      <DimensionRef ref="Time" />
      <RegularlySpacedValues start="11323.5" increment="1" />

mdsumner

#that for me needs env vars set: 
#COPERNICUSMARINE_SERVICE_USERNAME 
#COPERNICUSMARINE_SERVICE_PASSWORD

#that is able to obtain the "Authorization Bearer: <token>", which fwiw can be used directly with any GDAL tool 
## via GDAL_HTTP_HEADERS or GDAL_HTTP_HEADER_FILE (and we can use multidim, coming soon to R)
## so you don't have to use the packages stars' interface (lots of options): 
#https://github.com/pepijn-devries/CopernicusMarine/blob/a006052d343b09a52e0b9ff8e7aadd998ad9000f/R/cms_login.r#L25-L33

mdsumner

new gdal mdim mosaic, with 180 Bluelink daily ocean_temp netcdf files from Thredds

## obtain a list of Thredds urls (we previously cached, we grep out daily_ocean_temp, 180 NetCDFs)
files <- arrow::read_parquet("https://github.com/mdsumner/dryrun/releases/download/latest/BRAN-netcdf-2023-netcdf.parquet") |> 
  dplyr::filter(stringr::str_detect(url, ".*daily.*temp.*\\.nc$"))

## generate gdal mdim mosaic command and run, a few minutes

mdsumner

maybe fix for terra::stretch(histeq = TRUE)

stretch_hist <- function(x, ...) {
  ## stretch as if all the pixels were in the same band (not memory safe)
  rv <- terra::stretch(terra::rast(matrix(terra::values(x))), histeq = TRUE)
  ## set the values to the input, then stretch to 0,255
  terra::stretch(terra::setValues(x, c(terra::values(rv))), histeq = FALSE)
}

mdsumner

Just an example that was mentioned, I've been playing with vrtility today looking at other parts of the world too and learning how it works.

A lot to say, this is a pretty big area for Sentinel (~180x220km), which is 5561, 4483 pixels at 40m resolution (so reduce the size of extend_x, extend_y (in degrees), and move the POINT for where you want the centre. Can go as fine as 10m resolution.

You don't really need median calc done if there isn't a lot of cloud, and we don't necessarily need mask either, date can be as fine as a two-day window (but there's only imagery every few days, depending).

library(vrtility)
# La Pampa
x0 <- c(xmin = -65.5, ymin = -39.5, xmax = -62, ymax = -34)