# Generating COGs and STAC items from DataArrays

**URL:** https://discourse.pangeo.io/t/generating-cogs-and-stac-items-from-dataarrays/1913
**Category:** Data
**Created:** [November 10, 2021, 3:45pm UTC](https://discourse.pangeo.io/t/generating-cogs-and-stac-items-from-dataarrays/1913 "2021-11-10T15:45:51Z")
**Posts on this page:** 5
**Page:** 1

<div class="post-metadata">

### Author: ![TomAugspurger](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/tomaugspurger/32/21_2.png) [@TomAugspurger](https://discourse.pangeo.io/u/TomAugspurger)
#### Post date: [November 10, 2021, 3:45pm UTC](https://discourse.pangeo.io/t/generating-cogs-and-stac-items-from-dataarrays/1913/1 "2021-11-10T15:45:51Z")

</div>

Hi all,

I’m looking for guidance / best practices on writing an xarray object to (a collection of) COGs. Let’s start with a common case of a DataArray that’s indexed by `(time, band, y, x)`. Let’s also assume that it’s a chunked DataArray, with a chunksize of 1 for `time` and `band`, and it might be chunked along `y` and `x` as well.

My high-level questions:

1. Does rioxarray’s `.rio.to_raster(path, driver="COG")` have the right defaults? Is there anything special we should do to make sure we write “good” COGs for a single chunk?
2. Is there an established convention for organizing a directory of COG files that represent a 4-d datacube? Maybe some Open Data Cube folks have run into this?

I’m particularly interested in item 2, since I’d like to easily generate STAC Items for this collection of COGs. With a well-structured directory convention, we can easily group multiple STAC Assets (one per COG) that cover the same area in space and time into a STAC item. My proposed naming convention is

```auto
<prefix>/time=<time>/band=<band>-y=<y-coord>-x=<x-coord>.tif

```

This works well for xarray: we have coordinate information available when writing the chunk, so we can safely generate a unique name for a chunk using the `(time, band, y, x)` coordinates of, say, the top-left value in the chunk.

To make things concrete, here’s a small proof of concept. Full code is at [https://nbviewer.org/gist/TomAugspurger/06ea84bd63b19f86cd5560d9c6a8ff9d](https://nbviewer.org/gist/TomAugspurger/06ea84bd63b19f86cd5560d9c6a8ff9d).

Suppose we have some DataArray with dimensions `(time, band, y, x)` and enough metadata for rioxarray to know the CRS, extent, etc.

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/pangeo/original/2X/6/6489b4ac8ff67ec05c1c389ba99f7904853a8416.png)

`xcog` has a little utility for writing that DataArray out to a collection of COGs (one per chunk):

```auto
from pathlib import Path
import xcog

dst = Path("/tmp/cogs/")
dst.mkdir(parents=True, exist_ok=True)

items = xcog.to_cog_and_stac(data, prefix=str(dst), storage_options=dict(auto_mkdir=True))

```

That writes out these files:

```auto
$ tree /tmp/cogs/
/tmp/cogs/
├── time=2021-01-01T17:07:19.024000
│ ├── band=B02-y=4745100.0-x=399960.0.tif
│ ├── band=B02-y=4745100.0-x=454860.0.tif
│ ├── band=B02-y=4800000.0-x=399960.0.tif
│ ├── band=B02-y=4800000.0-x=454860.0.tif
│ ├── band=B03-y=4745100.0-x=399960.0.tif
│ ├── band=B03-y=4745100.0-x=454860.0.tif
│ ├── band=B03-y=4800000.0-x=399960.0.tif
│ ├── band=B03-y=4800000.0-x=454860.0.tif
│ ├── band=B04-y=4745100.0-x=399960.0.tif
│ ├── band=B04-y=4745100.0-x=454860.0.tif
│ ├── band=B04-y=4800000.0-x=399960.0.tif
│ └── band=B04-y=4800000.0-x=454860.0.tif
└── time=2021-01-04T17:17:19.024000
    ├── band=B02-y=4745100.0-x=399960.0.tif
    ├── band=B02-y=4745100.0-x=454860.0.tif
    ├── band=B02-y=4800000.0-x=399960.0.tif
    ├── band=B02-y=4800000.0-x=454860.0.tif
    ├── band=B03-y=4745100.0-x=399960.0.tif
    ├── band=B03-y=4745100.0-x=454860.0.tif
    ├── band=B03-y=4800000.0-x=399960.0.tif
    ├── band=B03-y=4800000.0-x=454860.0.tif
    ├── band=B04-y=4745100.0-x=399960.0.tif
    ├── band=B04-y=4745100.0-x=454860.0.tif
    ├── band=B04-y=4800000.0-x=399960.0.tif
    └── band=B04-y=4800000.0-x=454860.0.tif

2 directories, 24 files

```

It also returned a list of STAC items with one STAC item per `(time, y, x)` chunk. Notably, the separate bands have been merged into a single STAC item.

```auto
>>> items[0].assets
{'B02': <Asset href=/tmp/cogs/time=2021-01-01T17:07:19.024000/band=B02-y=4745100.0-x=399960.0.tif>,
 'B03': <Asset href=/tmp/cogs/time=2021-01-01T17:07:19.024000/band=B03-y=4745100.0-x=399960.0.tif>,
 'B04': <Asset href=/tmp/cogs/time=2021-01-01T17:07:19.024000/band=B04-y=4745100.0-x=399960.0.tif>}

```

Thanks to rioxarray and xarray’s coordinates, this has all the `proj` information necessary for tools like stackstac, so kind of a round-trip from DataArray → {STAC + COG} → DataArray is possible:

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/pangeo/original/2X/6/654ad7815a2cec5c02f568ab739b9a907c9afca3.png)

(I noticed that I’m off by one on the size there).

* * *

Hopefully that example makes the use-case clear. Looking forward to any thoughts from people who have dealt with this before. Thanks!

---

<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: [November 22, 2021, 11:16pm UTC](https://discourse.pangeo.io/t/generating-cogs-and-stac-items-from-dataarrays/1913/2 "2021-11-22T23:16:04Z")

</div>

Naming conventions are hard to design, especially for generic case. On one hand you want readability and compactness on the other hand you want to include enough information to have unique paths and possibly filenames too, you also want that information to be “glance-able by a human” (hence compactness is a plus).

I think for library/tool design you are better off starting with some sort of “path construction DSL”, then come up with reasonable default for that. This DSL could be as simple as Python format string with a bunch of properties made available to it, things like band name and various kinds of date/time representations and who knows what else, like keys from STAC properties.

Regarding specific examples in your post, I must say I’m not a fan ☹, sorry.

1. I doubt a lot of people have a good idea what UTM coordinate `-y=4800000.0-x=454860.0` corresponds to, also without knowing which UTM zone that is it’s not even possible to tell.
2. Other problem with raw coordinates is that they can be negative, and in this case you end up with two kinds of meaning for `-` character: a separator and a negative sign. This also ruins column alignment as negative coordinates have extra character to them, you might want to balance that with a `+` for positive values, but `+` is problematic for URLs so is not a good idea to include that. And then there are lon/lat which require a lot of precision. I’d say coordinates in a file path are just not practical.
3. Band is a primary property, so maybe replace `band={band}` with just `{band}`
4. Using correct amount of precision for timestamps would be good, but might be hard without out-of-band information provided by an operator.

I personally prefer notation where file name is just `{band}.tif` and all other identifying information is in the directory path. The path is then either `{time}/{space}` or `{space}/{time}` with various levels of nesting for both `{time}` and `{space}` components. Problem with that notation however is that file paths are not unique so you can’t easily just download a bunch into a folder for debugging without renaming as you go.

Regarding “off by one” error: `stackstac` defaults to using top-left corner for your `x/y` coordinates, but `rio-xarray` probably expects pixel centers. So you might end up with off by half a pixel transform written to TIFF, then when `stackstac` loads it back it does “pixel grid snapping” and so wants to offset everything by half a pixel. I wish `stackstac` had `xy_coords="center"` as a default instead. I haven’t actually ran any of that it’s just a guess…

---

<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: [November 23, 2021, 12:55am UTC](https://discourse.pangeo.io/t/generating-cogs-and-stac-items-from-dataarrays/1913/3 "2021-11-23T00:55:31Z")

</div>

Also see this issue in stackstac

> <https://github.com/gjoseph92/stackstac/issues/68>
>
> Hi all! 
> Either the shift in pixels that I have observed is actually unexpected… and needs to somehow be adjusted in the \`to\_coords\` function or I'm simply missing something or misunderstanding the concept behind geotransforms and it's only unexpected for me at the moment. Hope someone can help me out either way 🙂 
> 
> In the following screenshots, the colored image is always the original raster, while the greyscale image is always the one created with these lines of code:
> \`\`\`
> scene\_stack = stackstac.stack(items=stac\_obj, xy\_coords=xy\_coords, dtype='float32', rescale=False)
> scene\_stack.rio.to\_raster('./scene\_stack.tif')
> \`\`\` 
> Transform and projection of the raster: \`'proj:transform': \[10, 0, 399960, 0, -10, 52000020\]\` & \`'proj:epsg': 32632\` (UTM 32N)
> 
> 1. \`xy\_coords=False\`
> !\[xycoords\_false\](https://user-images.githubusercontent.com/56583917/130591389-14c6af0e-5584-4e6c-b650-8d8a1b58ad78.png)
> 
> 2. \`xy\_coords='center'\` | pixel shift: 0, 1 (expected: 0.5, -0.5)
> !\[xycoords\_center\](https://user-images.githubusercontent.com/56583917/130591410-7bf632cb-c978-4742-9ad8-03404e1e9530.png)
> 
> 3. \`xy\_coords='topleft'\` | pixel shift: -0.5, 0.5 (expected: 0, 0)
> !\[xycoords\_topleft\](https://user-images.githubusercontent.com/56583917/130591422-121d3ea7-37b4-43cc-a8ea-b07fc1fdebd4.png)
> 
> As expected, no shift is observed in (1). Now with (2) and (3), I'm confused. My original raster is a COG following the usual raster conventions, e.g. being top left aligned. Therefore I'd assume that using the default as in (3) results in no pixel shift at all. I'd also assume that the pixels are shifted by 0.5 and -0.5 in x and y respectively when using 'center' as in (2). 
> 
> Looking forward to hearing some other thoughts on this! 
> 
> Cheers, Marco

---

<div class="post-metadata">

### Author: ![TomAugspurger](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/tomaugspurger/32/21_2.png) [@TomAugspurger](https://discourse.pangeo.io/u/TomAugspurger)
#### Post date: [November 24, 2021, 9:37pm UTC](https://discourse.pangeo.io/t/generating-cogs-and-stac-items-from-dataarrays/1913/4 "2021-11-24T21:37:19Z")

</div>

Thanks for the feedback!

> On one hand you want readability and compactness on the other hand you want to include enough information to have unique paths and possibly filenames too, you also want that information to be “glance-able by a human”

This is a really important point. One additional thing to consider: now that we have STAC, maybe the balance should swing more towards uniqueness? In my mind, this pile of COGs would be easiest to explore / debug if you had STAC items for them. An initial prototype of this just used `dask.base.tokenize` on each chunk, which generated a unique string but was completely unreadable by humans. So the combination of `time/band, y, x` was a way to try to get uniqueness in a more human-friendly way.

My workflow was assuming blob storage, so I didn’t put much thought into nested directories vs. a flat namespace.

---

<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: [November 24, 2021, 11:01pm UTC](https://discourse.pangeo.io/t/generating-cogs-and-stac-items-from-dataarrays/1913/5 "2021-11-24T23:01:03Z")

</div>

> [@TomAugspurger](#):
>
> This is a really important point. One additional thing to consider: now that we have STAC, maybe the balance should swing more towards uniqueness? In my mind, this pile of COGs would be easiest to explore / debug if you had STAC items for them. An initial prototype of this just used `dask.base.tokenize` on each chunk, which generated a unique string but was completely unreadable by humans. So the combination of `time/band, y, x` was a way to try to get uniqueness in a more human-friendly way.

What is “human-friendly” depends on a workflow I suppose. One extreme is to just use UUID or hash and nothing else, this makes it super clear that filename is not to be interpreted directly. Other extreme is to include almost everything in the filename in a somewhat readable form, this only works for specific collections, see Landsat or Sentinel paths for examples of that. Mid point would be something like `{band}-{place}-{time}-{hash}` or `{band}-{time}-{place}-{hash}`, depending on what you value more “grouping by time” or “grouping by place”.

`{time}` is sort of easy, just need to figure out right level of precision, and support ranges and not just points in time. But `{place}` is much harder to do generically, it really only works for tiled data, and even then you need a naming scheme for those tiles. I’d say `{place}` will have to come from the user, maybe an attribute you can set on an xarray? And if it’s not provided by the user just omit it altogether and rely on `{hash}` for uniqueness. I guess you could encode a bounding box of the image into a filename but this will be even less readable and longer than a hash and harder to tab-complete too. If some sort of `{hash}` is too “unhuman”, one can substitute with some random words, like what docker does for container names.

I guess one could use one of those fancy things like [H3](https://eng.uber.com/h3/) to encode `{place}` in a generic way, but that’s getting complicated quick. I’d say better option is to allow user to override those things via attributes or encoding properties of the xarray.

> [@TomAugspurger](#):
>
> My workflow was assuming blob storage, so I didn’t put much thought into nested directories vs. a flat namespace.

In my experience there is still benefit in avoiding “directories” with too many files even when storage medium is some sort of blob store. For example reading COG from a blob store with thousands of “sibling” blobs can be done quickly in software if you configure GDAL right, but sometimes you want to open it in QGIS and there is no obvious way to tell QGIS not to bother looking for side-car files, so QGIS will struggle at open time even though image is COG.
