landsat8-c2-l2 proj:bbox is off by a half-pixel

[TomAugspurger]

Hi all, looking for some debugging help here.

We have an item, https://planetarycomputer.microsoft.com/api/stac/v1/collections/landsat-8-c2-l2/items/LC08_L2SP_231060_20200112_20200823_02_T1, where the proj:bbox doesn't match what's given by rasterio:

In [1]: import pystac, planetary_computer, rasterio

In [2]: item = pystac.Item.from_file("https://planetarycomputer.microsoft.com/api/stac/v1/collections/landsat-8-c2-l2/items/LC08_L2SP_231060_20200112_20200823_02_T1")

In [3]: ds = rasterio.open(planetary_computer.sign(item.assets["SR_B3"].href))

In [4]: item.properties["proj:bbox"]
Out[4]: [60900.0, -116400.0, 288900.0, 115800.0]

In [5]: list(ds.bounds)
Out[5]: [60885.0, -116415.0, 288915.0, 115815.0]

The bounds seem to be off by 15, half a pixel.

---

https://gist.github.com/TomAugspurger/cf5cdc4d6e94ec3757411dadae6e8558 has the XML MTL file for this item. IIUC, we get proj:bbox at

stactools/stactools_landsat/stactools/landsat/mtl_metadata.py

Lines 96 to 120 in 0e7509b

	@property
	def proj_bbox(self) -> List[float]:
	xs = [
	self._get_float(
	"PROJECTION_ATTRIBUTES/CORNER_UL_PROJECTION_X_PRODUCT"),
	self._get_float(
	"PROJECTION_ATTRIBUTES/CORNER_UR_PROJECTION_X_PRODUCT"),
	self._get_float(
	"PROJECTION_ATTRIBUTES/CORNER_LL_PROJECTION_X_PRODUCT"),
	self._get_float(
	"PROJECTION_ATTRIBUTES/CORNER_LR_PROJECTION_X_PRODUCT")
	]
	ys = [
	self._get_float(
	"PROJECTION_ATTRIBUTES/CORNER_UL_PROJECTION_Y_PRODUCT"),
	self._get_float(
	"PROJECTION_ATTRIBUTES/CORNER_UR_PROJECTION_Y_PRODUCT"),
	self._get_float(
	"PROJECTION_ATTRIBUTES/CORNER_LL_PROJECTION_Y_PRODUCT"),
	self._get_float(
	"PROJECTION_ATTRIBUTES/CORNER_LR_PROJECTION_Y_PRODUCT")
	]
	return [min(xs), min(ys), max(xs), max(ys)]

. Looking at the file, we do indeed have 60900.000 for the left-most bound. Is this pretty clearly a case of the metadata file being wrong, which should be reported to USGS? Or might rasterio / GDAL be doing something wrong? Or might stactool's derivation of proj:bbox failing to take something into account? Any help would be appreciated.

[gadomski]

The USGS docs (pdf) specify that those corner coordinates are center of pixel:

But GDAL geotransforms are measured to the corner of the pixel. So stactools is wrong, it needs to subtract/add half a pixel to get the correct bounding box.

[TomAugspurger]

Thanks for digging up those docs!