Raster datasets defined as S2 cells?

[florisvdh]

Just an idea, maybe it's nuts or maybe it is possible already - having almost no experience with S2 at this stage (I've read the sf vignette and some of the documentation at s2geometry.io).

Is there an object format (class) in R (or would one be planned) that is appropriate to contain an S2 grid dataset in R with origin, resolution and extent? I guess resolution would mean choosing a particular level of the S2 cell hierarchy. And it would accommodate at least one data attribute (at cell level).

[paleolimbot]

The short answer is that it's outside the scope of S2 and that you probably don't need S2 to implement it: each face of the S2 cube is (I think) an azimuthal projection assuming a spherical earth (PROJ can do this, I think). I'm working on a primitive in wk that is just an array + an extent that is a potential structure for hosting it. We've got some other development hopes and dreams for now, but look forward to seeing what you come up with!

[florisvdh]

you probably don't need S2 to implement it: each face of the S2 cube is (I think) an azimuthal projection assuming a spherical earth (PROJ can do this, I think)

Since all edges in S2 are spherical geodesics, yes, this should mean that an azimuthal projection with viewpoint at the Earth's centre (a gnomonic projection) is the operation needed to link a cube face with a corresponding face cell, and a straight line on the cube face with a corresponding spherical geodesic. Thanks!

Perhaps the challenge of generating a spherical coordinate system (CS) of equal-sized 'sphere cells' in which a spherical raster could be represented, could be approached in steps as follows:

1. Consider a collection of planar CSs (one for each cube face), in which we consider only the area within a square cube face.
2. Define equal-sized, square cells in each rectangular CS.
3. Now we need a projection pipeline to map equal-sized cells from the planar CS to equal sized sphere cells (geodesic quadrilaterals if that's the correct term), in order to link planar cells to locations in a geodetic CS (and eventually real locations, by adding a datum). In my imagination that might involve the inverse of the following pipeline (i.e. inverse of projecting sphere to plane): a. gnomonic projection from sphere to plane. It distorts areas (growing larger towards the CS edge. b. rescale distances along direction towards the origin (centre of the cube face) in order to achieve an equal-area projection (much like the equal-area azimuthal projection does).

I'm not aware of existing conversion formulae to create an equal-area gnomonic projection (and its inverse), but I'll have a look, perhaps will try some things in PROJ, it may be possible indeed.

Let me know if you saw it differently - this is what I make of it :wink:. I'm not so familiar with these subjects!

[paleolimbot]

The way I always envisoned it was that a raster could be represented by (1) An S2Cell representing an extent and (2) an S2Cell level representing the pixel. This means the pixel is always a cell (so can be queried using S2). This means you don't get to choose your pixel size but could use the machinery of S2 a bit more efficiently. (All speculation, here)

[florisvdh]

Sounds good! Indeed, not needing projections, but using S2 to do that would feel more elegant. Interfacing with S2 goes beyond my technical background though.

There is a (minor?) pitfall with the approach of an equal-area gnomonic projection: the scaling to obtain equal area will curve straight lines of the gnomonic projection - except those going in the direction of the origin, or perpendicular to that. But the cell edges would still remain straight, and that's probably all that's needed. Just to remind that oblique lines in the plane won't match spherical geodesics anymore (nor are their planar distances correct, as with all projections).

[edzer]

s2 might work slightly different from what you think, see https://docs.google.com/presentation/d/1Hl4KapfAENAOf4gv-pSngKwvS_jwNVHRPZTTDzXXn6Q/view#slide=id.i28 and that whole presentation.

[paleolimbot]

That's a good point...point->cell is unlikely to be repeatable using an existing PROJ. I pictured it being useful as a slippy tile-like system for distributing global data sets rather than something anybody would use for processing. GDAL probably can't do it but wk might be able to soon (next version will probably have a version of warp that can use an arbitrary C callable as a transform).

[florisvdh]

Thank you for referring these slides @edzer.

This confirms that a gnomonic projection is used for each cube face ((x,y,z) -> (face,u,v)), followed by an operation to ± equalize area ((face,u,v) -> (face,s,t)). The result of the applied algorithm is not perfect though: largest S2 cells are about twice the size of the smallest, cf. slide 14 - and also seen in https://s2geometry.io/resources/s2cell_statistics. So a raster S2 implementation will have (implicit) non-constant cell size (with a size range depending on the raster extent). Should be taken into account when e.g. summarizing area of raster subsets.

[paleolimbot]

Closing since this isn't an s2 development priority at the moment (the formulas are all that's needed to implement this elsewhere)!