Support multiband GeoTIFFs such that all bands are processed

[mdales]

Currently Yirgacheffe works with a single band within any raster file. It defaults to picking the first band in a GeoTIFF, but you can select other bands if you wish by providing the band=n argument when you load a raster.

What you can't do is process all the bands within the raster. This would be useful if you wanted to say crop a multi band dataset to a GeoJSON (i.e., working with data like the 128-band Tessera tiles).

I can imagine we want a way to mark a raster as being multi-band at load time, meaning that you apply to it all the operations in parallel. For example, something perhaps like:

import yirgacheffe as yg

with yg.read_raster("my_multiband_tile.tif", bands=yg.ALL) as raster_layer:
  with yg.read_shape("my_clipping_data.geojson") as clipping:
    clipped = raster_layer * clipping
    clipped.to_geotiff("result.tif")

The result here would be a multiband raster with all the bands clipped.

The challenge then comes when a multiband layer is applied to other things, what does it even mean if someone tries to write:

import yirgacheffe as yg

with yg.read_raster("my_multiband_tile.tif", bands=yg.ALL) as raster_layer:
  with yg.read_shape("my_clipping_data.geojson") as clipping:
    wtf = clipping / raster_layer
    wtf.to_geotiff("result.tif")

How can a multiband layer sit on the RHS of an expression like that? Do we just error? Do we pick the first layer? Do we expand clipping internally to become multiband? That last one is very tempting until someone tries to do:

import yirgacheffe as yg

with yg.read_raster("my_multiband_tile.tif", bands=yg.ALL) as raster_layer:
  wtf = raster_layer * raster_layer
  wtf.to_geotiff("result.tif")

Are we going to generate a result that has the number of bands squared? I feel that is a dangerous thing to make easy to do with the amount of data in a geospatial pipeline.

---

One thing I do like about this though is that in theory it's another way to start to express parallelism without having to have the ecologist worry about parallelism.

For example, in the classic AOH example we take a set of species ranges as polygon data, and multiply them by habitat and elevation maps. Could we start to use this to do something like:

import yirgacheffe as yg

with yg.read_raster("habitat.tif") as habitat_map:
  with yg.read_raster("elevation.tif") as elevation_map:
    with yg.read_shape("species.gpkg", bands=yg.ALL) as species:
      aohs = habitat_map * elevation_map * species
      aohs.to_geotiff("aohs.tif") 

That doesn't really work, as the habitat and elevation layers in AOH maps need to be adjusted per species, but I feel there's something in this idea that would let yirgacheffe unlock parallelism without the ecologist needing to think in computer terms.

[mdales]

If we have an expression with two multiband layers, do we do the multiplicative layer results as the output, or do we just apply layer to layer? The later feels more fragile (we have to ensure all layers match) but more likely to be what people want.

[mdales]

Expanding on that AOH example from above:

import yirgacheffe as yg

with yg.read_raster("habitat.tif") as habitat_map:
  with yg.read_raster("elevation.tif") as elevation_map:
    with yg.read_shape("species.gpkg", bands=yg.ALL) as species:
      aohs = habitat_map * elevation_map * species
      aohs.to_geotiff("aohs.tif") 

That, as I mentioned, is a gross simplification, as we need to filter the layers by species specific habitat and elevation preferences. However, in a GPKG you can have other metadata beyond the polygon. If we had a way to expose the other properties on the layer (which we should do anyway, as in my current AOH code I load the GeoJSON twice, once as a yirgacheffe layer and once as a JSON file to get the metadata, which now I see is silly - filed as #64), then we could also do that per band...

import yirgacheffe as yg

with yg.read_raster("habitat.tif") as habitat_map:
  with yg.read_raster("elevation.tif") as elevation_map:
    with yg.read_shape("species.gpkg", bands=yg.ALL) as species:
      filtered_elevation = (elevation_map > species.lower_elevation) && (elevation_map < species.upper_elevation)
      filtered_habitat = habitat_map.isin(species.habitat_list)
      aohs = filtered_habitat * filtered_elevation * species
      aohs.to_geotiff("aohs.tif") 

So long as we use the same "band" from the GPKG for all bits of the internal s-expression, that would work quite nicely I think.

[mdales]

In terms of not making perfect the enemy of good enough, to cater to the Tessera case we can just do a pass over the s-expression and limit multiband layers to:

• Just one per s-expression
• Can only appear in certain places

And I think that'd be a nice way to start. But I do want to explore the API a bit first so we don't end up making the data-parallelism approach harder to do.

[mdales]

Just one limitation here is that the TIFF image format, one which GeoTIFF is based has a limit to the number of bands due to only storing the number in a SHORT value, so 16 bit unsigned value, limiting it to 64k layers. Whilst on one had that is a lot (more than the number of terrestrial species processed in the LIFE pipeline) it is quite low for other domains (there are 350k species of plant).

Still, I think for many domains this is a workable limit.

[mdales]

Another operation on multi band data would be a flatten option. For example, if you had a GeoTIFF or GPKG with species AOH or ranges separated into bands you could "flatten" them with an '+' operator and get the species richness result.

The question is just how to express that without resorting to lambdas etc. which are not very user friendly.

[mdales]

I think there needs to be an operator that combines multiple calculators into a set of layers. I wonder if we can somehow override the array operator? Probably not, but something where you can do:

with yg.read_raster("foo.tif") as input:
  scaled_1 = input * 2
  scaled_2 = input * 3

  all_scales = yg.assemble_bands(scaled_1, scaled_2)
  all_scales.to_geotiff("bar.tif")

all_scales is then a multiband layer, the same as if I'd opened a multi band raster