stac_diy.md

x <- readLines("https://earth-search.aws.element84.com/v1/search?limit=500&collections=sentinel-2-l2a&datetime=2023-10-20T00:00:00Z%2F..&bbox=100,-67,112,-62")
js <- jsonlite::fromJSON(x)
imred <- vapour::gdal_raster_data(js$features$assets$red$href[1], target_dim = c(1024, 0))
imblue <- vapour::gdal_raster_data(js$features$assets$blue$href[1], target_dim = c(1024, 0))
imgreen <- vapour::gdal_raster_data(js$features$assets$green$href[1], target_dim = c(1024, 0))

library(terra)
## me mangling the metadata from my vapour output, but it's just extent/dimension/crs
r <- rast(ext(attr(imred, "extent")), nrows = attr(imred, "dimension")[2], ncols = attr(imred, "dimension")[1], nlyrs  = 3, vals = cbind(imred[[1]], imgreen[[1]], imblue[[1]]))
plot(terra::project(terra::vect(x), attr(imred, "projection")))
## terra automatically stretches the Int16 but I'm not sure how yet (how to decide the max?)
plotRGB(r, add = T)

using new osgeo rocker:

sys.path.append('/usr/local/lib/python3/dist-packages')
from osgeo import gdal
gdal.UseExceptions()
urls = ['/vsicurl/https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/10/S/EG/2022/7/S2B_10SEG_20220728_0_L2A/B04.tif', '/vsicurl/https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/10/S/EG/2022/7/S2B_10SEG_20220728_0_L2A/B03.tif',  '/vsicurl/https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/10/S/EG/2022/7/S2B_10SEG_20220728_0_L2A/B02.tif']

urls = [f'vrt://{url}?scale=0,16383&ot=Byte' for url in urls ]

vrt = gdal.BuildVRT('/vsimem/built.vrt', urls, separate = True)
vrt.RasterCount()
3

## now we can gdal.Warp() that same as gdal_raster_data above
w = gdal.Warp('/vsimem/out.tif', vrt, width = 1024, height = 0)
w.RasterXSize ## 1024
w.RasterYSize ## 1024
w.RasterCount ## 3