Making the coarsest resolution the base of the pyramid

[felixcremer]

I discussed with @rafaqz how to deal with Base and DD functions in PyramidScheme and we came to the conclusion that the Base function should work as if the Pyramid is a single layer and the different other layers are just add ons.
If we do it that way, it might be better to invert the pyramid, so that the normal functionality would happen on the coarsest resolution because that would mean, that the user would not run into memory issues when they try to run some Base functionality suddenly on a very large array.
We should add some errors or hints to tell the user, that they are operating only on the coarsest resolution and how to get the higher resolution data.

[rafaqz]

That makes sense to me. Then the size is for the smallest level, and you can actually try things in the REPL on the pyramid data.

Maybe for broadcast where we search for regular arrays and find them, we can just broadcast the top layer instead of the error?

[meggart]

I am not fully convinced on this idea yet. Could you share a few examples where this would be beneficial? I mean, I actually really like the philosophy of applying the operations lazily on all layers as you do with broadcast and then whenever user access a bounding box (or even the full array) we return whatever level fits nicely into memory.

[rafaqz]

That already works, but not with another array that's not a pyramid... It's not clear how that should work.

You need to make all arrays pyramids first. So allowing broadcast on the top layer only with regular arrays would be a faster way to explore the data and ideas that will later need full pyramids

[meggart]

Can't we decide the moment we do the broadcast and pick the pyramid level that fits the array best? I mean, if we broadcast between a pyramid and a regular array it will only work if the base level has the same size as the other array anyway,

In general I would really like to follow a philosophy where a pyramid is a thing without a definite size but with an extent.

[meggart]

To be a bit more clear about what I suggest:

I think the best way to interact with a Pyramid is by defining a function extractdata!(ras, pyr) where ras is a raster defining with X and Y defining the grid to be extracted. The implementation of the extractdata function would choose the pyramid resolution which is closest to the target x and y resolution and return the nearest-neighbor interpolated values on the target grid. This is almost identical to the way we do tile extractions for plotting on different resolutions so the functionality for this is already there.

If a user has a pyramid pyr and a Raster ras and does pyr .+ ras then this might be translated to extractdata!(similar(ras), pyr) .+ ras and could just work in-memory on the selected resolution.

If we formulate broadcasting like this, this would allow us to do in-memory operations interactively not only on the coarsest overview but also on small regional subsets in higher resolutions. However, the price we would need to pay is that Pyramids would probably take a special place and not fit in the DimensionalData/Raster type hierarchy as it would work if they were just equivalent to their finest or coarsest resolution from a type-hierarchical point of view.

[rafaqz]

Ok that's better! So any array size that matches a pyramid layer size could work in a broadcast?

Pretty easy to add to the current code where it errors now. It doesn't need the extractdata! step either unless it needs to a resample (to sizes not in the pyramid).

But I'm not sure a resample should be in a broadcast as it makes assumptions about the data (e.g. categorical/intensive/extensive) that we cant change in a broadcast, (unless that's built into the pyramid I guess)

Also what does it mean for the array interface more generally? Flexible size will be difficult

[meggart]

But I'm not sure a resample should be in a broadcast as it makes assumptions about the data (e.g. categorical/intensive/extensive) that we cant change in a broadcast, (unless that's built into the pyramid I guess)

Yes this would be a design decision although I actually like the idea of building resampling methods into the Pyramid, where the default could be to error when the target resolution does not exist, but more flexible might be selected as well and become a type parameter of the Pyramid.

Also what does it mean for the array interface more generally? Flexible size will be difficult

I agree and if we follow this path then I don't think a Pyramid could be an AbstractArray. I guess we could still support many array operations by applying them lazily to the layers, but in general I would say if users want to work with arrays instead they can always extract a layer.

[meggart]

I am still wondering if what we are proposing is actually specific to Pyramids. Another use case that comes to my mind would be lazily computed DimArrays, for example Top-Of-Atmosphere radiation which can be computed through a formula given a lon, lat and time. This can be materialized on all resolutions, but do we have a way to conveniently express this as a DimArray? One might also consider defining an InterpolatedDimArray where the Interpolation algorithm is defined but the target resolution is not realized. All these and Pyramids have in common that they have named spatial dimensions and an Extent, but no concrete grid. So if we define an AbstractLazyGridDimArray all of these could have a common supertype and we would need to define broadcast rules only once for all of them. The interface for the abstract type would be the extractdata! function mentioned above which materializes everything to one grid. Then the following pseudocode might work:

ras = Raster("normalrasterfile.tif") #Has regular grid
toa = MyTOAObject() #Can compute TOA radiation for any time and place on earth
pyr = Pyramid("somepyramid", resample_method="nearest") # Has several resolutions
ras2 = Raster("secondrasterfile.tif")
ras2_reproj = reproject_raster(ras2, method=...) #we define the reporjection method, but not the target grid, nothing happens here

## Everything will materialize only in this final step, which selects the correct pyramid level, computes toa and reprojects ras2 to the same resolution as ras
my_function.(ras,ras2_reproj, pyr, toa) 

Is there already something like this in the DimArray-verse?

[rafaqz]

Not really, we have GDAL based resample in Rasters and it's all eager.

Maybe DimArrayMaterializer to avoid specifying interpolation as it might do something else?

It could materialise an array inside broadcasts as needed based on the lookups of the other DimArrays in the broadcast.

It could even have non-matching extent and crop or fill with missing values as needed. We could have Interpolations.jl, GDAL resample, PyramidScheme etc as backends.