opendatacube / odc-geo

GeoBox and geometry utilities extracted from datacube-core
https://odc-geo.readthedocs.io/en/latest/
Apache License 2.0
83 stars 15 forks source link

`GeoBox.zoom_to` with a resolution #79

Closed keewis closed 1 year ago

keewis commented 1 year ago

I've been trying to use odc-geo to define the target grid for regridding with xesmf. This works great so far, but one of the features I've been missing so far is creating a Dataset from the GeoBox.

The naive version I've been using so far

def to_dataset(geobox):
    def to_variable(dim, coord):
        data = coord.values
        units = coord.units
        resolution = coord.resolution
        if resolution < 0:
            data = data[::-1]
            resolution *= -1
        return xr.Variable(
            dims=dim, data=data, attrs={"units": units, "resolution": resolution}
        )

    coords = {
        name: to_variable(name, odc_coord)
        for name, odc_coord in geobox.coordinates.items()
    }
    return xr.Dataset(coords=coords)

which obviously lacks support for affine transformations, assumes 1D coordinates, and the resolution trick should probably configurable (and off by default?)

Did I miss anything? Do you know of a better way to do this? If not, would you be open to adding something like this to odc.geo.xr?

Kirill888 commented 1 year ago

I don't think I follow what you mean by "convert GeoBox to Dataset".

do any of ⬇️ do what you need?

Kirill888 commented 1 year ago

those produce DataArray, or coordinates for it, you can then wrap dataset around that if needed

keewis commented 1 year ago

xr_coords is sufficiently close, I totally missed that. I'd still have to change the resolution, but that's easy enough.

Thanks a lot for the quick reply!

Kirill888 commented 1 year ago

I'd still have to change the resolution, but that's easy enough.

like this? https://odc-geo.readthedocs.io/en/latest/_api/odc.geo.geobox.GeoBox.flipy.html

keewis commented 1 year ago

I'm not sure. This helps with doing the actual transformation, but I'd still have to decide when to call it (and switch between flipx and flipy depending on the coordinate). Basically, what I want is to normalize the resolution to always be positive, since that's the way the grids I'm working with are oriented, and xarray's default indexes require the orientation to be identical (although with custom indexes that might change).

keewis commented 1 year ago

actually, since that particular geobox was created by rescaling a existing geobox I could just pass a Resolution object instead of just a float.

Which leads to a different question: is there a simple way to create a geobox from an existing one, but with changed resolution? Basically, what I'm looking for is something similar to zoom_out or zoom_to, but with a absolute step size instead of a zoom factor or a new shape. I've been using

new_geobox = GeoBox.from_bbox(geobox.boundingbox, resolution=0.01)

but that feels clunky, and that's also where the negative resolution is from.

Kirill888 commented 1 year ago

so like geobox.zoom_to(resolution=0.01)? not a thing currently, but makes sense to have

you can use GeoBox.from_bbox(geobox.boundingbox, resolution=odc.geo.resxy_(0.01, 0.01)) to force non-flipped y-axis at construction time.

keewis commented 1 year ago

thanks. I'll probably use odc.geo.Resolution(x=0.01, y=0.01) since that's a bit more explicit, but I'm happy with the result.

Edit: I'll leave this open to track the geobox.zoom_to(resolution=...) suggestion