Refactor Geography and add (scalar) creation functions

[benbovy]

This PR refactors the Geography (sub)classes in spherely and adds geography creation scalar functions.

Supersedes #26.
Closes #41.
Closes #27.

• Goal is to remove the Geography subclasses (Point, Polygon, etc.): only one Geography wrapper class remains
• Geography type (and also whether the geography is empty or not) is inferred at object creation
• Geography objects are created via the top level creation functions (only scalar for now except for points)

TODO:

• fix segfault when creating a geography collection itself including collection(s)
• add multipolygon creation function
• add EMPTY geographies creation function overloads (+ tests errors when trying to build MULTI geographies from EMPTY geographies)
• move geography creation tests into test_creation.py
• raise when trying to create linestring with one point only

[jorisvandenbossche]

I know this is not yet ready, was just already taking a quick look.

One thought regarding the scope of this PR, since it is a big one: maybe you could limit this PR to only cover what essentially already existed (to get the current tests (updated for new function names) passing). For example just provide scalar creation functions for linestring/polygon, and leave the vectorized creation functions (for everything other than point) for a follow-up.

[benbovy]

Yes I've been working further on this (last updates not uploaded yet) but that might be too much of a change indeed. I'll restrict this PR to the refactoring of Geography wrapper classes.

Also I'm not exactly sure whether we should closely reuse Shapely's API for creation functions. I find that API convoluted with the mixing of a unique coords argument vs. separate x, y (and z) arguments as well as the different behavior depending on indices. We should probably reuse the same API for compatibility although it makes things more cumbersome to implement.

[jorisvandenbossche]

Yeah, the vectorized creation functions in shapely are quite complex, and I am also not sure if they are always that useful (like creating polygons from a single 3D array, which I think in practice is only useful to generate some random triangles or to create bboxes), and agreed we should carefully think about which API we exactly want here.
That's also the reason I mentioned leaving out the vectorized creation functions for now, and only providing the scalars ones (that already exist and that are used in the tests). For points, I would keep the vectorized function you already have, since in this case it is quite simple, I think? (and you already have that)

[benbovy]

For points, I would keep the vectorized function you already have, since in this case it is quite simple, I think? (and you already have that)

Yes. For a single array this is currently limited to 2D (i.e., not a true ufunc like in shapely) but that's fine for now. I don't think that we can use pybind11.vectorize to extend support to a single nD input array (pybind/pybind11#1294).

I wouldn't be too hard I think to mimic the behavior of a ufunc by relying on pybind11's array view & indexing capabilities (+ some manual computation of nD indices), though. That would be even easier to let xtensor(-python) takes care of all of that (also for creating temporary arrays of geographies on the C++ side while releasing the GIL), but I'm not sure if adding such build dependencies is worth it.

[benbovy]

This is getting close. I'd just like to investigate the segfault when creating nested geography collections before making this PR ready (EDIT: fixed segfault).

A few other comments that we can address later:

• Empty geographies: there is no overload function to create MULTI geographies without argument (or with None) ; also when empty or 1-length sequences are given the multi-types are demoted (see Singleton Multi-types automatically demoted paleolimbot/s2geography#17)
  • I guess the issue is the same in ENH: add from_wkt/to_wkt functions #50 and maybe Refactor Geography and add (scalar) creation functions #51 ?
  • We could address this in Spherely by forcing the geography type in the spherely.Geography wrapper, but it depends on what we want to do.
• The collection type is named GeographyType::GeographyCollection but both the Geography repr and WKT have the GEOMETRYCOLLECTION name. Do we want to uniformize things and rename GeographyType::GeographyCollection to GeographyType::GeometryCollection?
• Full polygon not supported (yet)
• It is not currently possible to create a MULTIPOLYGON with duplicate rings, due to the fact that s2geography::PolygonGeography reuse the same underlying S2Polygon for multi-polygons. Maybe if we deactivate s2 normalization and/or validity check it will pass but will it work ?
• No test added for vectorised point creation functions, but the API may not be final.
• Scalar creation function naming: still no strong opinion, although I would lean towards either point(), linestring(), etc. or create_point(), create_linestring(), etc.

[jorisvandenbossche]

Nice!

[jorisvandenbossche]

I don't know what the performance overhead of extract_geog_properties would be, but should we allow to create Geography() with a given geography type (like for clone_s2geography you have a variant that infers vs receives the type, we could have the same for extract_geog_properties). For certain functions (like creation functions), you always know the type in advance.

[jorisvandenbossche]

Actually, I should have commented that on the case below (creating from s2geography), because for the line above in theory we could propagate the attributes from the passed Geography object? (like is done in the move operator implementation. Or if you want to avoid this overhead, just use move instead?)

[benbovy]

Actually, I should have commented that on the case below (creating from s2geography), because for the line above in theory we could propagate the attributes from the passed Geography object?

Yes indeed we should simply propagate the attributes in this case.

In the case of creating from an s2geography object when we know the geography type, we could indeed provide the type explicitly instead of inferring it (and therefore avoid the inference overhead). We certainly want to avoid passing inconsistent s2geography object type vs. geography type, though, in order to prevent odd behavior or segfaults happening later. Checking for this consistency at runtime would be pretty similar to extract_geog_properties in terms of performance overhead I think. If we really want to optimize this, perhaps we could use template constructors + SFINAE or if constexpr...

[jorisvandenbossche]

This is getting close.

A few other comments that we can address later:

Agreed! Let's indeed get this in (to have the basics working, and will also be easier for other PRs adding features), and the things below can be handled in follow-ups.

• Empty geographies: there is no overload function to create MULTI geographies without argument (or with None) ; also when empty or 1-length sequences are given the multi-types are demoted (see Singleton Multi-types automatically demoted paleolimbot/s2geography#17)

Personally I am fine with that limitation for now, especially given that s2geography essentially does not make the distinction between single and multi geometries anyway (as you linked), and that it is only for output (like WKT, or asking for the type) that the distinction is introduced, based on the number of elements stored in the object.

• The collection type is named GeographyType::GeographyCollection but both the Geography repr and WKT have the GEOMETRYCOLLECTION name. Do we want to uniformize things and rename GeographyType::GeographyCollection to GeographyType::GeometryCollection?

I would maybe lean towards GeographyType::GeometryCollection, although that is certainly odd as well ..

Another case where this occurs is in the spherely.geography_collection(..) constructor. For that case, we could also shorten it to just collection() to avoid that issue? (or at least when we add a prefix like create_.., that might be sufficient)

• It is not currently possible to create a MULTIPOLYGON with duplicate rings, due to the fact that s2geography::PolygonGeography reuse the same underlying S2Polygon for multi-polygons. Maybe if we deactivate s2 normalization and/or validity check it will pass but will it work ?

A multipolygon with duplicate rings also sounds like "wrong" .. (in shapely / GEOS, that would be an invalid geometry anyway, and thus give unexpected results), so I am personally fine with not allowing to create it.

Although I suppose at some point we should allow to create invalid geographies, end then have a method to check if they are valid / what the reason is for the invalidity (for example s2_make_polygon has a check keyword to disable validation when constructing a polygon)

This is also another case (together with the geometrycollection / geographycollection above) where it points the the mismatch between the "simple features" model we try to emulate and how s2 works.
For the constructor functions you implemented here, it is nice to make the distinction between all the types, I think (as it gives a more user friendly interface), but for accessors like the geography type, we could also stay more closely to s2geography and have a function that just returns point/line/polygon/collection

• Scalar creation function naming: still no strong opinion, although I would lean towards either point(), linestring(), etc. or create_point(), create_linestring(), etc.

Let's open a separate issue for this and see if we can get some more opinions

[benbovy]

I would maybe lean towards GeographyType::GeometryCollection, although that is certainly odd as well ..

Another case where this occurs is in the spherely.geography_collection(..) constructor. For that case, we could also shorten it to just collection() to avoid that issue?

I renamed GeographyType::GeographyCollection -> GeographyType::GeometryCollection and geography_collection() to collection(), I also think it looks nicer for now.

[jorisvandenbossche]

OK, merged this, so we can move forward with other PRs and follow-ups. Thanks @benbovy!

[benbovy]

👍 thanks for the reviews @jorisvandenbossche !