Skip to content

Instantly share code, notes, and snippets.

@mdsumner
Last active June 26, 2025 04:03
Show Gist options
  • Save mdsumner/2af156a8250dd825221c634f1b6fdeba to your computer and use it in GitHub Desktop.
Save mdsumner/2af156a8250dd825221c634f1b6fdeba to your computer and use it in GitHub Desktop.

Here's my take on defining a "geobox" at a longlat point, where I only need a longitude, latitude, distance, and n-pixels.

I use a local projection to ease the geographic extent definition, then transform that to (auto) UTM, then get the longlat bbox (for STAC query), and the UTM GeoBox (for odc to render to).

from odc.geo import BoundingBox
from odc.geo.geobox import GeoBox
import math
ltln = (-54.50,158.93)
bf = 1000  ## radius around our ltln point
dm = 1024  ## shape of the grid (n-pixels) if we moke x,y for bf we need to change aspect ratio here too (like gdal_translate -outsize 1024 0)

gb0 = GeoBox.from_bbox(BoundingBox(left=-bf, bottom=-bf, right=bf, top=bf,   
                crs=f"+proj=laea +lon_0={ltln[1]} +lat_0={ltln[0]} +datum=WGS84"), shape=(dm, dm))

gb = gb0.to_crs(f'EPSG:32{6 if ltln[0] < 0 else 7 }{str(math.ceil((ltln[1] + 180) / 6))}')
bbox = gb.to_crs("EPSG:4326").boundingbox
gb.explore()

image

https://rstats.me/@mdsumner/114747599744750602

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment