# Example which highlights the limitations of NetCDF-style coordinates for large geospatial rasters

**URL:** <https://discourse.pangeo.io/t/example-which-highlights-the-limitations-of-netcdf-style-coordinates-for-large-geospatial-rasters/4140>\
**Category:** Uncategorized\
**Created:** [April 3, 2024, 4:05pm UTC](https://discourse.pangeo.io/t/example-which-highlights-the-limitations-of-netcdf-style-coordinates-for-large-geospatial-rasters/4140 "2024-04-03T16:05:39Z")\
**Posts on this page:** 1\
**Showing post:** 6

<div class="post-metadata">

**Author:** ![kirill.kzb](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/kirill.kzb/32/1276_2.png) [@kirill.kzb](https://discourse.pangeo.io/u/kirill.kzb)\
**Post date:** [April 3, 2024, 11:47pm UTC](https://discourse.pangeo.io/t/example-which-highlights-the-limitations-of-netcdf-style-coordinates-for-large-geospatial-rasters/4140/6 "2024-04-03T23:47:53Z")

</div>

This issue comes up every now and then:

> <https://github.com/opendatacube/odc-geo/issues/127>
>
> I ran into this issue while reprojecting some data but I think the underlying ca…use is that \`GeoBox\` objects created by different means for the same grid location will have affine transforms that are slightly different. It may be specific to EPSG 4326 since I haven't run into this issue with other projections.
> 
> When reprojecting to EPSG 4326, the geobox for the result does not match what I input. A short example:
> 
> \`\`\`py
> import numpy as np
> import xarray as xr
> import rioxarray as riox
> 
> src = xr.DataArray(
> np.ones((1, 10, 10)),
> dims=("band", "y", "x"),
> coords=(\[1\], np.arange(10), np.arange(10)\[::-1\]),
> ).rio.write\_crs("EPSG:3310")
> 
> dst\_gb = src.odc.geobox.to\_crs(4326)
> dst = src.odc.reproject(dst\_gb)
> print("GeoBox comparison:", dst.odc.geobox == dst\_gb)
> print("Shape comparison:", dst.odc.geobox.shape == dst\_gb.shape)
> print("CRS comparison:", dst.odc.geobox.crs == dst\_gb.crs)
> print("Affine comparison:", dst.odc.geobox.affine == dst\_gb.affine)
> print()
> print("dst\_gb affine")
> print(repr(dst\_gb.affine))
> print("result affine")
> print(repr(dst.odc.geobox.affine))
> print()
> print("Affine comparison with tolerance:", np.allclose(list(dst\_gb.affine), list(dst.odc.geobox.affine)))
> \`\`\`
> \`\`\`
> GeoBox comparison: False
> Shape comparison: True
> CRS comparison: True
> Affine comparison: False
> 
> dst\_gb affine
> Affine(1.0200436314046532e-05, 0.0, -120.00002388788778,
> 0.0, -1.0200436314046532e-05, 38.01647531889047)
> result affine
> Affine(1.0200435637076365e-05, 0.0, -120.00001592863761,
> 0.0, -1.0200435637131023e-05, 38.01647279736898)
> 
> Affine comparison with tolerance: True
> \`\`\`
> 
> I'm using version 0.4.2.

It’s important to understand that GeoBox is recomputed from coordinates, that’s needed to support slicing into geo-referenced data, and also to support data constructed by other libraries. BUT there is no guarantee that this recomputed GeoBox will produce exactly the same coordinates when used to create a new array from it. Essentially `GeoBox -> coords -> GeoBox` is not guaranteed to be lossless. It can only be lossless when both resolution and translation components of the Affine matrix are basically integers. Sentinel-2 has that property for example, `scale=+/-10 tx,ty=10*N` where `N` is an integer.

Not sure what a proper solution for this should be. We can keep track of the original resolution in an attribute of the coordinate, and use exactly that value when extracting GeoBox from coords (with a check to deal with `xx[::10, ::10]` type of slicing. I guess we can also keep track of original translation, and only recompute that if array has been sliced. I’ll probably implement that for the next version of `odc-geo`, actually.

The problem of creating sub-geobox that will be able to produce exactly the same coords as the original geobox in the sliced section is much harder to address. We want this invariant:

`gbox[roi].coords == gbox.coords[roi]`

That’s not possible unless `gbox[roi]` retains parent and slice, and essentially return `self.parent.coords[self.roi]`. Or we implement “rounding to some fraction of a pixel” kinda logic.

---

_[View the full topic](https://discourse.pangeo.io/t/example-which-highlights-the-limitations-of-netcdf-style-coordinates-for-large-geospatial-rasters/4140)._