Extending longitude to beyond 360 degrees.

[malektronic]

Rationale

When mapping global environmental data, it can be very helpful to be able to visualise large structures which span mutliple ocean basins. Normally one has to choose which basin to cut, but if a map can be repeated then no choice
need be made. This suggestion was made by Edward R. Tufte in The Visual Display of Quantitative Information (Graphic Press, Cheshire, Connecticut, 2nd Ed, 2001; p99). Though he warned against needless repetition of data, he advocated for this approach for long, cyclical structures, and in my own work I have found it to be useful for interpreting data, and in at least one case had I used it it would have helped me avoid losing time from a mistake.

Another use case is for non-repeating data such as non-cyclic trajectories, as was suggested in this issue: #1633

Since in my suggested code modifications Cartopy does not itself extend the data (this must be done manually, though I've written utilities for that), both use cases are catered for.

I realise that there are a lot of changes here, but it wouldn't work properly without all of them!

Implications

For the user not requiring this new functionality, it should not change anything. For those who choose to use it, simply adding over=True when PlateCarree() or Mercator() projections are invoked will enable the functionality. Nightshade also requires the over=True parameter to work properly, and coastlines and tissot have been modified to work. Other functions like adding gridlines work without modification.

[CLAassistant]

All committers have signed the CLA.

[rcomer]

This white space change seems to be responsible for most of the test failures. Is it intentional?

[malektronic]

Thanks for letting me know! It was unintentional. I've fixed this and I'm now looking at the other test failures, which I'll hopefully be able to resolve today.

[malektronic]

I've marked this ready for review again after fixing the problems with the documentation and resolving the pre-commit errors. All the remaining errors appear to be related to fetching resources from nsidc.org. The service appears to be off line at the moment, so it is independent of anything to do with Cartopy or my modifications.

Note that I have incorporated changes made since my own modifications and rebased.

[greglucas]

Mostly general comments about this feeling like we could potentially be playing whack-a-mole with the +over needing to be added / special cased everywhere. I wonder if there is an easier way to handle this somehow and try to remove some of our hard-coded assumptions on (-180, 180) instead.

[greglucas]

Rather than having this as a separate keyword argument here, shouldn't this be a part of the proj4_params that gets passed in? The subclasses of this should update the proj4_params.

[malektronic]

Yes I thought it ought to be, that would be cleaner for sure. But I couldn't work out how to implement it nicely since all the other proj4_params that are passed have an argument, but +over does not.

[greglucas]

Isn't the same as you have added it to the list here already? i.e. something like proj4_params += [("over", None)]

[greglucas]

This is an interesting case. I guess this is us being nice to users and saying if you passed in 0-360 data we will automatically wrap that to -180 to 180 for you, but if we wanted to be more explicit this should be written as
PlateCarree(over=False).transform_points(..., transform=PlateCarree(over=True)) where the source transform has the over=True and the output doesn't.

[malektronic]

I guess in that case there would be no more reason to check, because you could use self.over instead of False in the first invocation of PlateCarree? I.e.

PlateCarree(over=self.over).transform_points(..., transform=PlateCarree(over=True))

[greglucas]

I realize that "over" is what PROJ uses, but I wonder if there is a more descriptive word we should use for these classes since it doesn't really mean much to me when first looking at this. "extend map beyond 360 degrees" might have some more description about how it ignores the (-180, 180) wrap and where the limits are and what happens when you go beyond them in this case.

[malektronic]

I admit I hesitated on this point. I shall rewrite and clarify the description, but I'm unsure about the parameter name. Possible options could be extend_lon, no_wrap, repeat_lons, cycle_lon, over_lon....

[greglucas]

What happens with +over when supplied with central longitude? Do we need to put the bounds check in the _determine_longitude_bounds() function call?

[malektronic]

Good question. Naively, I left this alone since it ceases to have any meaning when the longitudes are set by the extent of the map with over=True. When over=False everything is as normal, but with over=True it churns and doesn't appear to work, I got tired of waiting and killed it. I wonder if this is because I have hard set the bounds, instead of setting them to the data? A simple fix would be to force central_longitude to be equal to 0 when over=True. I'll dig into it, including looking at _determine_longitude_bounds().

[greglucas]

Could you explain what this is used for? I have to read the code currently because there is no overview docstring.

The name should be CamelCase

Suggested change

	class NaturalEarthFeature_ext(NaturalEarthFeature):
	class NaturalEarthFeatureExtended(NaturalEarthFeature):

[malektronic]

It's to enable the extension of the coastlines beyond 360 degrees. I found that trying to do that within FeatureArtist broke other functionality, so I decided it was best to leave FeatureArtist agnostic with regard to whether +over was set. (I have a branch called extend_in_feature_artist in my cartopy-extlon repository, if you want to see how I tried to implement it.) As a result, FeatureArtist makes no assumptions about the data it is fed, and the unmodified NaturalEarthFeature will work, but will only cover 360 degrees of the map. So the subclass NaturalEarthFeature_ext will repeat the coastlines as many times as necessary to cover the extent defined by the minimum and maximum longitudes set when plotting. (See commit 116bb16. ) Since NaturalEarthFeature is invoked to produce hard-coded feature variables in __init__.py, I could not easily get it to read in the extent beforehand, thus my decision to explicitly call it in coastlines(). It is probably better to implement the changes within the original NaturalEarthFeature, but I wanted to tread lightly for now, so I left it as a separate class, thinking that since the user doesn't see it it should be transparent for them. As for the name, I will amend it, but I'll wait until we have settled on a name for the over variable.

Incidentally, the problem I had with my original approach of extending features in FeatureArtist informed most of my decisions henceforth, and it was why I settled on implementing over on a case-by-case basis for each feature (I have not done them all, since I wasn't sure how this would be received, but I intend to). This speaks somewhat to your comment about playing "whack-a-mole" - I'll say more on that as a direct reply to that comment, but this is part of why it has ended up this way. It has the benefit of not requiring FeatureArtist to care about the feature, so a feature could conceivably only cover a single cycle, in case that was desired, but the flip side is that features which are intended to repeat must be written to be "over-aware", so to speak.

[greglucas]

I feel like we shouldn't need to pass "over" to our features to explicitly get a different feature for the two cases. It seems like we've designed something wrong with so many (-180, 180) assumptions if that is the case. Is there a way we can handle this differently or more generically for any feature to play nicely with the +over map-types?

[malektronic]

Part of it is that my priority was to not break what was already there, and the exceptions allowed me to keep the old code as-is to a large extent. But extending the globe beyond 360 degrees is sort of unnatural, in that there is no such thing as a 361st degree - it is just a place-holder for repeating the 1st degree, so I'm not sure how we could completely eliminate the check. A routine like NightShade, or even coastlines, must first process [0, 360] before then repeating the information, and repetition is a different kind of functionality to the original calculation. I guess this would be better determined by checking if the extents were over 360 degrees, but the code logic wouldn't change.

I had tried modifying FeatureArtist to perform the duplication, but this led to a whole host of problems, not least of which is that since there are different types of feature with different types of data, they do not always behave the same by the time FeatureArtist is invoked. On top of that, some "features" do not pass through FeatureArtist at all! I do think that there is a more elegant way, using a function perhaps, but I haven't yet gone through all the features to see how they work.

I guess this also goes back to your comment about the central_longitude parameter and the automatic resetting of boundaries to [-180, 180].

Regarding the [-180,180] assumption, I did come across this a few times but I was actually surprised how little I needed to modify the code, so far from being "wrong" this suggests the design is actually quite good!

[malektronic]

Mostly general comments about this feeling like we could potentially be playing whack-a-mole with the +over needing to be added / special cased everywhere. I wonder if there is an easier way to handle this somehow and try to remove some of our hard-coded assumptions on (-180, 180) instead.

Thank you for taking the time to look at these modifications in so much detail, it's much appreciated. I hope we can improve it to everyone's satisfaction.

Your comment about central_longitude and being "nice" to users by automatically wrapping are relevant here, so I'll have a think about what can be changed. What surprised me was actually how little modification the code needed to get this to work, so I guess I had a different perspective!

I mention in an earlier reply that the way FeatureArtist works made it exceedingly difficult to use it to duplicate the features beyond a 360-degree Earth. There is also a philosophical choice - do we change FeatureArtist, and how it deals with geometries, so that it can provide this function? Or is it better that it stays true to its intended function, which is to simply draw what it is asked to, and let the user, or the creator of the stock features, provide the data they intend? After wrestling with this throughout this process, I decided on the latter (though I may have been influenced by what was doable). Furthermore, not everything that is plotted passes through FeatureArtist, so those must be explicitly modified anyway - e.g. stock_img(), which calls a modified imshow().

So far as I can tell (and I obviously have not tested everything), the only things that still don't play well with over (appart from central_longitude) are features which I have not yet looked at. So these are the only places over would need to be added - meanwhile, provided I haven't missed anything, the engine seems to be behaving, barring a poor choice of explicit extent limits. As for the individual features, it seems reasonable to provide a generic function (or set of functions) to extend the features, to avoid having to do so explicitly within the feature code. This, I should think, would be best left until I implement (or try to implement) over with all of the features.

[greglucas]

What surprised me was actually how little modification the code needed to get this to work, so I guess I had a different perspective!

I agree that this doesn't seem that bad, so nice work chasing it down and thank you for submitting the PR and work to the main repo too!

The concern I have is for future expansion/maintainability because I think this would mean any new features/methods etc. would need both a "normal" and an "over" handler.

I mention in an earlier reply that the way FeatureArtist works made it exceedingly difficult to use it to duplicate the features beyond a 360-degree Earth. There is also a philosophical choice - do we change FeatureArtist, and how it deals with geometries, so that it can provide this function? Or is it better that it stays true to its intended function, which is to simply draw what it is asked to, and let the user, or the creator of the stock features, provide the data they intend? After wrestling with this throughout this process, I decided on the latter (though I may have been influenced by what was doable). Furthermore, not everything that is plotted passes through FeatureArtist, so those must be explicitly modified anyway - e.g. stock_img(), which calls a modified imshow().

We have some utility functions that add cyclic points. I wonder if this should be a similar helper function that is something like (bad name, but hopefully illustrative) repeat_geometry_to_extended_longitudes() that would repeat the input geometries shifted by the wrap point.
https://github.com/SciTools/cartopy/blob/main/lib/cartopy/util.py

Another thought for this to be more fundamental is that this needs to be handled in the transform_* functions. If a user passes in data that spans from (170 - 190) in longitude do we expect that to show up at (-190, -170) and (530 - 550) or does a user need to explicitly add the duplicate data points into the extended longitude regime before calling the functions?

It might be good to have some other opinions on how this should look because it is a pretty fundamental change somewhere in the function stack. I have not given this enough thought at this point to have preferences one way or the other, I'm simply trying to put some ideas/thoughts down in text in case that sparks something for you or anyone else at this point.

[malektronic]

I agree that this doesn't seem that bad, so nice work chasing it down and thank you for submitting the PR and work to the main repo too!

Thank you, it's kind of you to say so!

The concern I have is for future expansion/maintainability because I think this would mean any new features/methods etc. would need both a "normal" and an "over" handler.

It may be that I was too zealous in keeping as much of the original code as possible. I will see if I can simplify the routines.

We have some utility functions that add cyclic points. I wonder if this should be a similar helper function that is something like (bad name, but hopefully illustrative) repeat_geometry_to_extended_longitudes() that would repeat the input geometries shifted by the wrap point. https://github.com/SciTools/cartopy/blob/main/lib/cartopy/util.py

I think so, I'll look into what that will require.

Another thought for this to be more fundamental is that this needs to be handled in the transform_* functions. If a user passes in data that spans from (170 - 190) in longitude do we expect that to show up at (-190, -170) and (530 - 550) or does a user need to explicitly add the duplicate data points into the extended longitude regime before calling the functions?

This is more subtle than I realised, but I think the present behaviour is actually what we want. The transform functions, if I understand correctly, are in a basic sense dressed-up versions of the Proj Transformer.from_crs() method. Thus they behave just the same. Given that Proj itself is "over-aware", I think we are safe to stick with their choices. They have even included a force_over flag. The behaviour is not, however, always what I expected:

From PlateCarree to PlateCarree
| src over | tgt over | force_over | input_lon | output |
|----------+----------+------------+-----------+--------|
| False    | False    | False      |       500 |    500 |
| False    | True     | False      |       500 |    140 |
| True     | False    | False      |       500 |    140 |
| True     | True     | False      |       500 |    500 |
| False    | False    | True       |       500 |    500 |
| False    | True     | True       |       500 |    500 |
| True     | False    | True       |       500 |    500 |
| True     | True     | True       |       500 |    500 |

From PlateCarree to Mercator
| src over | tgt over | force_over | input_lon |   output |
|----------+----------+------------+-----------+----------|
| False    | False    | False      |       500 | 15584729 |
| False    | True     | False      |       500 | 15584729 |
| True     | False    | False      |       500 | 15584729 |
| True     | True     | False      |       500 | 55659745 |
| False    | False    | True       |       500 | 55659745 |
| False    | True     | True       |       500 | 55659745 |
| True     | False    | True       |       500 | 55659745 |
| True     | True     | True       |       500 | 55659745 |

I've tried this with Mercator to Mercator and Mercator to PlateCarree and the results are equivelent: as long as the source and target CRSs are different, we get what I would expect, which is that if over is false everywhere, then whatever our input longitude we get an output within [0, 360] degrees. It is when we have the same CRS for source and target that we get the unexpected (to me at least) result that our output longitude is the same as the input regardless of whether or not the input lies beyond 360 degrees (first line in first table). I don't think this is a problem, however: after all, if over is false everywhere, then why would an input be greater than 360, and why would you transform to and from the same CRS in the first place? I guess it's an edge case that we should be aware of, but apart from that everything is as I would have designed it. I do not think the transform functions should be duplicating longitudes (e.g. returning a list like (-350, 10, 370)) , they should remain one-to-one. In other words, data duplication is proper to the data, and should remain so. Aside from risking over-complication, it makes certain types of chart, with non-repeated data, impossible. Excuse this long-winded justification, but hopefully it's thorough enough that someone can poke holes in it.

It might be good to have some other opinions on how this should look because it is a pretty fundamental change somewhere in the function stack. I have not given this enough thought at this point to have preferences one way or the other, I'm simply trying to put some ideas/thoughts down in text in case that sparks something for you or anyone else at this point.

That would be welcome. Thanks to these exchanges I am beginning to understand the problem a lot better.

[greglucas]

It is when we have the same CRS for source and target that we get the unexpected (to me at least) result that our output longitude is the same as the input regardless of whether or not the input lies beyond 360 degrees (first line in first table).

I think we might have some logic that says if things are the same don't even call the transform function to speed up those transforms.

One other thought that I had which I don't know if you've looked into is whether you can place axes next to each other with no margin between them and replicate this capability without creating new geometries but rather plotting them on different GeoAxes instances and make them look like they are tiled properly.

As a stretch I'm thinking this could even be interesting for repeating (with a flip) in the south-north direction to make arbitrarily large tiled maps. My thought is that there would be some TiledGeoAxes class that would handle the duplication for a user and stitch things together under the hood and keep track of the 3 separate GeoAxes instances under the hood. If a user called my_map.plot(x, y) it would then pass that to all N tiles it is keeping track of, so you would have N objects with the same lat/lon pairs instead of just one that is a different set of lat/lons.

I think that could be quite complex, but wanted to write that down as a thought too.

[malektronic]

One other thought that I had which I don't know if you've looked into is whether you can place axes next to each other with no margin between them and replicate this capability without creating new geometries but rather plotting them on different GeoAxes instances and make them look like they are tiled properly.

I like this idea.

Back when I first thought about doing this (some years ago now) I initially thought of simply replicating the figure, but I quickly realised that this would most likely require learning the inner workings of matplotlib, which seemed too difficult at the time. Honestly I may have tried going down that route had I not seen #1633 (comment), but I don't think I'd have come up with the Tiled class as you describe it - that seems worth pursuing.

There's a fair bit to recommend this approach: data need not be duplicated, Cartopy need not necessarily be modified (to be determined, I suppose), and we are no longer bound by +over's limitations (+/- 572.95 degrees, only cylindrical projections). I also like the idea of wrapping vertically. Also, such code could be used in other applications, potentially.

The two advantages of exploiting +over that I see are that 1) it is easier to implement, and enables Cartopy to more fully exploit Proj which underlies so much of it and 2) it enables non-repeated data. (1) is mainly my bias really, I've done it so it's there, notwithstanding the work still to do to adequately address the comment's you've made. But there's no reason I can think of that (2) fundamentally can't also be possible in a tiled system, as it can be made to understand both the standard [0, 360] and extended coordinates (with an offset related to the axes geometry) depending on a parameter.

I really would like the functionality to exist one way or another so I'm happy to contribute in either direction - or both!

Over Christmas I hope to find time to get stuck into this a bit more and now that I understand the code - and matplotlib - a bit better I will try to see if I can figure out what will be necessary for a tiling system.