Introduce special state _epsg=0, that implies "epsg code is not known yet", _epsg=None means "epsg code is known to be absent".

Force all CRS construction code paths through cache.

UTM convenience methods

1. CRS.utm(some_geom) picks appropriate utm zone for a given geometry and returns corresponding CRS object
2. When transforming geometry or geobox accept crs="utm|utm-n|utm-s" as a valid destination crs, meaning pick "appropriate" utm CRS for an object being transformed.
3. ~"appropriate UTM zone" is currently defined by the lon,lat of the centroid of the geometry~
4. UTM zones that overlap with some_geom are located using pyproj, then zone with biggest overlap is chosen

Snap to

Adds GeoBox.snap_to(other), this method adjust source geobox slightly by shifting footprint by a sub-pixel amount such that pixel edges align exactly with pixel edges of the other geobox. Method returns adjusted copy of self. Method only works when it is possible to align pixels by pure translation only - same projection, same resolution, same orientation. ValueError is raised when it's not possible.

Closes #59

Enclosing

Related to snap_to somewhat. GeoBox construction method that creates a new "pixel grid compatible" geobox that minimally covers some area.

Generalized Slicing into GeoBox

Overload geobox[roi] to support geometry types on input, footprint in the world can be converted to pixel space, from which pixel slice is computed.

Closes #31

• Adding library description to readme and docs

• Adding badges to readme on github

• Adding readthedocs integration

• Bumping version and fixing description in the wheel

• Fixing up add some more __repr__ methods

• Reverting changes to multipolygon constructor function

• rename odc.geo._overlap to just odc.geo.overlap

• fix in compute_reproject_roi when dst is an integer shrink of source

There were a bunch of places where it was not clear whether some tuple contains X,Y or Y,X data. It's confusing on the interface and confusing when reading code. This PR adds an XY[T] type and a bunch of constructor methods that make it easy to understand what order x,y parameters are added.

Instead of using plain tuples those new types are used instead. On the input shape= parameters are still allowed to be supplied as plain (int, int) tuples and are equivalent to iyx_(y, x). Internally those properties are kept in XY[int] form. Similarly a 2d index is also allowed to be supplied as (int, int), as this is what allows some_thing[ix, iy] or some_thing[row, col] syntax, order depends on use case, but should be clearly stated in the docs.

Resolution input parameters are allowed to be a single float, but plain tuples are not allowed as this causes order confusion. Single value resolution is equivalent to resyx_(-value, value).

Change GeoBox constructor to accept shape

Geobox(width, height, ...) interface is problematic for several reasons:

1. not compatible with ndarray, where shape is used instead
2. no type-safety, easy to get order wrong
3. two parameters to configure single property

Replacing width, height with a single parameter shape. This could be a simple (int, int) tuple for easier interop with ndarray. But one can choose to use more typesafe XY[int].

Example of the change at call site

   # before change:
   GeoBox(512, 256, A, "espg:4326")

   # becomes
   ## if array order is "native" at call-site
   GeoBox((256, 512), A, "epsg:4326")
   ## if width/height orded is preferred
   GeoBox(ixy_(512, 256), A, "epsg:4326")

this failed

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
/tmp/ipykernel_19782/44836221.py in <cell line: 1>()
----> 1 newcheckdask = odc.geo.xr.colorize(newcheck['test'], cmap=newcmp,vmin=1,vmax=240)

~/exploracornzarr/lib/python3.8/site-packages/odc/geo/_rgba.py in colorize(x, cmap, attrs, clip, vmin, vmax)
    240 
    241         assert x.chunks is not None
--> 242         data = da.map_blocks(
    243             _impl,
    244             x.data,

~/exploracornzarr/lib/python3.8/site-packages/dask/array/core.py in map_blocks(func, name, token, dtype, chunks, drop_axis, new_axis, meta, *args, **kwargs)
    732         adjust_chunks = None
    733 
--> 734     out = blockwise(
    735         func,
    736         out_ind,

~/exploracornzarr/lib/python3.8/site-packages/dask/array/blockwise.py in blockwise(func, out_ind, name, token, dtype, adjust_chunks, new_axes, align_arrays, concatenate, meta, *args, **kwargs)
    274         from .utils import compute_meta
    275 
--> 276         meta = compute_meta(func, dtype, *args[::2], **kwargs)
    277     return new_da_object(graph, out, chunks, meta=meta, dtype=dtype)
    278 

~/exploracornzarr/lib/python3.8/site-packages/dask/array/utils.py in compute_meta(func, _dtype, *args, **kwargs)
    158                 return None
    159 
--> 160         if _dtype and getattr(meta, "dtype", None) != _dtype:
    161             with contextlib.suppress(AttributeError):
    162                 meta = meta.astype(_dtype)

~/exploracornzarr/lib/python3.8/site-packages/dask/delayed.py in __bool__(self)
    588 
    589     def __bool__(self):
--> 590         raise TypeError("Truth of Delayed objects is not supported")
    591 
    592     __nonzero__ = __bool__

TypeError: Truth of Delayed objects is not supported

this works - but as your docs suggest, gives strange results

newcheckdask = odc.geo.xr.colorize(newcheck['test'], cmap=newcmp,vmin=1,vmax=240)

Have I done something wrong, or is there a bug?

• Closes #52
• maybe addresses #51
• Add GeoBox.compat to convert to datacube version of GeoBox
• In reproject use nan as fill value for float32 outputs if output nodata is not configured
  • fixes boundary effects when plotting float32 input on a map
• In colorize handle nodata values properly for non float inputs
  • fixes boundary effects when plotting non-float inputs on a map
• In colorize drop nodata attribute on the output array
• In geobox display support non-epsg authority codes like ESRI
• Better Dask support for GeoBox, CRS and GeoBoxTiles classes (fixes hashing)
• Add GCP extraction utilities
• Add Geometry[Point] -> XY[float] conversion
• Deprecated Geometry.__iter__ interface to follow shapely deprecation of that
• Fixes .odc.add_to for maps with custom projections

Motivation

Code in xx.odc.geobox assumes that xarray raster data is axis aligned. That is the off-diagonal terms of the affine matrix are zero and so entire row of pixels will have the same Y coordinate and similarly entire column same X coordinate:

|y\x|0|1|2|3| |-|-|-|-|-| |0|0,0|1,0|2,0|3,0| |1|0,1|1,1|2,1|3,1|

This limitation means one can not use .odc.geobox on rotated images. For a lot of satellite data one swatch of observation happens on an angle, when this gets rendered into axis aligned image you end up with a large proportion of the mosaic image being empty. Supporting arbitrary affine transform could make it much easier to work with large scale mosaics for such data sources.

Possible Implementations

2-d axis

Xarray supports 2d coordinate, so one could simply compute X,Y coordinate for every pixel and store that in two 2d arrays, one for X and one for Y.

pix = img[r,c]
coord = (img.x[c], img.y[r])  # axis aligned representation
...
coord = (img.x[r,c], img.y[r,c]) # 2d-axis representation

.odc.geobox could then detect that 2d axis are used and reconstruct affine matrix from coordinates similar to current implementation but it will need to be computed from 3 points (Affine matrix has 6 degrees of freedom).

• :-1: Significant memory overhead (2 float32 values per pixel, so 8 bytes per pixel)
• :-1: This representation demands axis being Dask arrays as well.
• :+1: Same representation can be used for completely non-linear pixel to world mapping as well
• :+1: Direct support for plotting from xarray itself

Pixel coords + Affine matrix in attributes

In this scheme we keep 1-d X,Y axis, but they contain pixel coordinates and not real world coordinates, so x = [0.5, 1.5, 2.5... W-0.5], that sort of thing. We then also need to store actual affine matrix in some attributes attached to x,y axis.

In this scenario .odc.geobox needs to detect the fact that X,Y axis contain pixel coordinates rather than real world coordinates and instead extract affine matrix from attributes/encoding. Slicing should still work as we will have access to the original image coordinates which is what we need to compute world coordinates from. Need to experiment with the location and representation for affine matrix data and how that works with things like .to_netcdf|to_zarr.

• affine format: 6 float values (string vs array vs dict)

• attachment point: x,y axis vs crs axis vs band

• :+1: Same memory requirements as axis aligned representation

• :-1: xarray plotting won't display more than one dataset on the same axis properly

Currently testing transitioning from datacube to odc-geo in geocube: https://github.com/corteva/geocube/pull/95

Have some scenarios in the tests where the shape is off by one: odc-geo: (y: 100001, x: 11) datacube: (y: 100001, x: 12)

And the boundary differs: odc-geo: [1665945., 7018509., 1665478., 7018306.] datacube: [1665478., 7018306., 1665945., 7018509.]

Are these changes expected?

GCP

When fitting polynomial use lower rank if only few points available.

Closes #67

New features

• BoundingBox.boundary() generate points along the perimiter
• {BoundingBox|GeoBox}.qr2sample quasi-random 2d sampling using R2 sequence
• Geometry.assign_crs use different CRS without changing geometry itself
• Geometry.transform(..., crs=) can now optionally change CRS of the transformed geometry
• odc.geo.geom.triangulate wrapper for for shapely method
• MaybeCRS now includes special Unset case to differentiate from None
• More robust Geometry.geojson (support composite and crs=None geometries)
• rio_geobox helper method in converters

• Support matplotlib colormaps in colorize/add_to
• Change to_grba to use vmin=, vmax= instead of non-standard clamp=
• Re-order API docs
• Add documentation section for xarray extension
• Add documentation section for working with rasters
• Extended sample data to include country geometries from geopandas

1. BoundingBox was an odd one out without CRS info attached, now it has that

   • This enables other conveniences, two way conversion polygon<->bbox

2. Changes to GeoBox.from_bbox

   • Can get CRS info directly from bounding box now
   • Replace confusing "align" that is pretty hard to explain with "anchor"

Essentially if you take 0,0 in the real world, bring it to pixel space and locate pixel that contains it, then you want location within that pixel to be a certain point in [0, 1), [0, 1), so 0,0 is edge aligned, same as align=(0,0) used to be, and to get center aligned pixel you pick (0.5, 0.5). Also you can just not align at all, in which case resulted geobox will start at left,top and might go a little bit over right, bottom.

Also adds some more convenience methods to bounding box and geobox.

Closes #25

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?

Currently odc-geo is tested in one single configuration: python 3.8 with optional dependencies installed to maximize test coverage. We should be testing things with minimal dependencies and with other versions of Python.

Related to: opendatacube/datacube-core#752

Currently the resolution inside of the coordinate attributes are used:

https://github.com/opendatacube/datacube-core/blob/05093b75a8f15b643c7047502ae27b662af2d9b4/datacube/utils/xarray_geoextensions.py#L133-L139

GeoTransform is a property used by GDAL to store this information (https://gdal.org/drivers/raster/netcdf.html#georeference)

Here is an example of reading in this property: https://github.com/corteva/rioxarray/blob/f658d5f829a88204b2ec7b239735aba31694ed0a/rioxarray/rioxarray.py#L540-L545

Thoughts about supporting this as well?

• [DONE] ~Output folium or ipyleaflet maps when libraries are available when displaying geometries/geoboxes/bounding boxes.~
• Provide GeoJSON outputs for more types
• Conversion to geopandas when available

Most of the low-level utilities needed for Dask-backed reprojection are already in odc-geo. This is mostly

https://github.com/opendatacube/odc-geo/blob/11c59371d306e010e8cbc4b4b835a80005e3ec4e/odc/geo/geobox.py#L909

and also: https://github.com/opendatacube/odc-geo/blob/11c59371d306e010e8cbc4b4b835a80005e3ec4e/odc/geo/overlap.py#L130

Expected interface:

xx = dc.load(.., dask_chunks= {}) # or any other supported load backend

# automatically choose resolution and bounding box, align pixel edges to 0
# automatically choose chunk size
yy = xx.odc.to_crs("epsg:3857")

# fully defined destination pixel plane
# configurable destination chunking
yy = xx.odc.reproject(GeoBox.from_bbox(..), chunks={'x': 2048, 'y': 4096})

• Support Dask and non-Dask inputs
• Support xarray and dask.dataarray
• Support xarray Datasets as well as DataArray
• Support reasonable automatic resolution, bbox and chunking determination
• Support any number of leading dimensions as well as optional interleaved band dimension: ..., y, x[,band]

We are caching pyproj.CRS objects here:

https://github.com/opendatacube/odc-geo/blob/202cd8aeb3b03a54d9e484c3514c19a8b57b1a7b/odc/geo/_crs.py#L31-L32

And pyproj transformers here:

https://github.com/opendatacube/odc-geo/blob/202cd8aeb3b03a54d9e484c3514c19a8b57b1a7b/odc/geo/_crs.py#L47-L48

• What that means for multi-threaded access?
• Should we use lock when populating cache?
• Should we use thread local cache instead of locking?
• What are the constraints of multi-threaded access in pyproj/PROJ itself
  • currently we are assuming that sharing CRS objects across threads is fine?
• We should also add purging rules to those caches, or at least expose manual purge option.
• We should understand cost of caching in terms of RAM, especially for transformers cache.