geobox_dsn.md

GeoBox from a DSN

def extent(dimension, transform): 
  return [transform[0], transform[0] + dimension[0] * transform[1], 
  transform[3] + dimension[1] * transform[5], transform[3]]
  
  
def bbox(dimension, transform): 
  return [transform[0], transform[3] + dimension[1] * transform[5],
          transform[0] + dimension[0] * transform[1], transform[3]]
  

def geobox_dsn(dsn, dimension = None, crs = None):
  ds = gdal.Open(dsn)
  ## let user set this (or resolution)
  if dimension is None: 
    dimension = [ds.RasterXSize,ds.RasterYSize]
  dsn_crs = ds.GetSpatialRef().ExportToWkt()
  
  bbox_src = bbox([ds.RasterXSize,ds.RasterYSize], ds.GetGeoTransform())
  ## allow user override (not transform)
  if crs is None:
    crs = dsn_crs
  ds = None
  return GeoBox.from_bbox(bbox_src, shape = dimension, crs = crs)
  
from odc.geo.geobox import GeoBox
from osgeo import gdal

dsn = "/vsicurl/https://gebco2023.s3.valeria.science/gebco_2023_land_cog.tif"

geobox_dsn(dsn)

geobox_dsn("/vsicurl/https://raw.githubusercontent.com/mdsumner/rema-ovr/main/REMA-2m_dem_ovr.vrt", dimension = [2725, 2921])

def extent(dimension, transform): 
  return [transform[0], transform[0] + dimension[0] * transform[1], 
  transform[3] + dimension[1] * transform[5], transform[3]]
  
  
def bbox(dimension, transform): 
  return [transform[0], transform[3] + dimension[1] * transform[5],
          transform[0] + dimension[0] * transform[1], transform[3]]
  

def geobox_dsn(dsn, dimension = None, crs = None):
  ds = gdal.Open(dsn)
  ## let user set this (or resolution)
  if dimension is None: 
    dimension = [ds.RasterXSize,ds.RasterYSize]
  dsn_crs = ds.GetSpatialRef().ExportToWkt()
  
  bbox_src = bbox([ds.RasterXSize,ds.RasterYSize], ds.GetGeoTransform())
  ## allow user override (not transform)
  if crs is None:
    crs = dsn_crs
  ds = None
  return GeoBox.from_bbox(bbox_src, shape = dimension, crs = crs)
  
from odc.geo.geobox import GeoBox
from osgeo import gdal
gdal.UseExceptions()
import odc.stac
dsn = "/vsicurl/https://raw.githubusercontent.com/mdsumner/rema-ovr/main/REMA-2m_dem_ovr.vrt"
geo = geobox_dsn(dsn, dimension = [2725, 2921])


# import os
# import json
import rasterio
# import urllib.request
import pystac
# 
from datetime import datetime, timezone
from shapely.geometry import Polygon, mapping
# from tempfile import TemporaryDirectory

datetime_utc = datetime.now(tz=timezone.utc)

      
ds = gdal.Open(dsn)
bb = bbox([ds.RasterXSize, ds.RasterYSize], ds.GetGeoTransform())## needs the footprint doh
        

item = pystac.Item(id='remav2',
                 geometry=geo.footprint("4326").geom,
                 bbox=bb,
                 datetime=datetime_utc,
                 properties={})
                 
item.add_asset(
    key='dem',
    asset=pystac.Asset(
        href=dsn,
        ## it's VRT, but all COGS underneath
        media_type=pystac.MediaType.GEOTIFF
    )
)

rema = odc.stac.load([item], geobox = geobox_dsn(dsn, dimension = [2500, 2500]), chunks = {})
rema.load()
 
<xarray.Dataset> Size: 25MB
Dimensions:      (y: 2500, x: 2500, time: 1)
Coordinates:
  * y            (y) float64 20kB 3.343e+06 3.341e+06 ... -2.495e+06 -2.497e+06
  * x            (x) float64 20kB -2.7e+06 -2.698e+06 ... 2.746e+06 2.748e+06
    spatial_ref  int32 4B 3031
  * time         (time) datetime64[ns] 8B 2024-04-19T01:45:58.522406
Data variables:
    dem          (time, y, x) float32 25MB nan nan nan ... -26.85 -26.84 -26.82