Add ability to specify custom bounds to Projections

[djhoese]

Rationale

See #813 for more info. The summary is that I (and others) want to be able to take an arbitrary pyproj CRS object with our own bounds and provide it to cartopy. Currently a Projection instance will try to use the static "area_of_use" provided by the PROJ/pyproj definition of a CRS. This is fine, but there are many cases where datasets only cover a specific region of a CRS. By allowing users to specify the bounds for a Projection you can reuse the same Projection instance in multiple locations without worrying about extents on the plot. This closely follows with the way that Pyresample users resample data and then plot various pieces of data sharing the same geolocation information.

This also adds an optional transform_bounds keyword argument to Projection as a convenience when the user wants to provide bounds in lon/lat degrees and have them transformed to CRS coordinates.

Implications

The bounds and transform_bounds keyword arguments are optional so should not effect existing users at all. I would say "advanced" users seeking this feature will benefit.

There may be cases for certain CRSes where the basic logic provided in Projection isn't good enough (I assume that's why some of the Projection subclasses exist, for example the _RectangularProjection and _CylindricalProjection classes). Having this new flexibility in the Projection class may make users think any projection can be provided to the Projection class and maybe that isn't true. I honestly don't know all of the edge cases in cartopy.

Partial Example

import pyproj
import cartopy.crs as ccrs
pcrs = pyproj.CRS.from_user_input({'proj': 'geos', 'lon_0': 0.0, 'a': 6378169.00, 'b': 6356583.80, 'h': 35785831.00, 'units': 'm'})
crs = ccrs.Projection(pcrs, bounds=[-2000.0, 2000.0, -2000.0, 2000.0])

# use it in normal cartopy usage

CC @snowman2

[djhoese]

Assuming people are comfortable with this, let me know if/where I should add documentation in the sphinx docs.

[snowman2]

If the pin for pyproj is moved to 3.1+ Transformer.transform_bounds could potentially be useful.

[djhoese]

If the pin for pyproj is moved to 3.1+

Obviously up to cartopy devs, but maybe in a future PR?

[snowman2]

maybe in a future PR?

There may be other places in the code that would likely benefit from that as well, so a separate PR would be a better fit for that.

[greglucas]

Currently this only helps with updating rectangular extents? Maybe this would be addressed by the transformer_bounds suggestion... To use transformer_bounds could we add in some hasattr()s with comments to begin using it before the pin would be increased?
See this PR for an example of some desired behavior: #1665
I think the primary issue is when we calculate the boundary (which is in the init of that projection).

Thoughts on if we should add **kwargs and pass them through projection subclasses? That way one could call PlateCarree(bounds=...) right away without having to call Projection(PlateCarree(), bounds=...)

[djhoese]

Currently this only helps with updating rectangular extents? Maybe this would be addressed by the transformer_bounds suggestion... To use transformer_bounds could we add in some hasattr()s with comments to begin using it before the pin would be increased?

I'm fine with adding this. Without knowing all of the edge cases, do you have an example of what you mean beyond that linked PR? I'm used to requesting bounds in X/Y projection coordinates. The issue(s) so are talking about only exist if someone provides lon/lat bounds, right? For example, non-eqc/geographic projections would have different transformed bounds depending on if you transformed (lon_min, lat_min) or (lon_min, lat_max). That's the concern, right?

Thoughts on if we should add **kwargs and pass them through projection subclasses? That way one could call PlateCarree(bounds=...) right away without having to call Projection(PlateCarree(), bounds=...)

Obviously not a cartopy maintainer so feel free to ignore me, but I personally would try to avoid **kwargs as much as possible. Maybe it is OK in this case if the top parent class complains about kwargs that it doesn't understand. Otherwise, you end up with users with typos in their kwargs or using kwargs meant for other interfaces and being confused why the kwarg is having no effect. That concern aside, I agree that this bounds functionality would be great to have on all Projection (or all CRS) classes.

[greglucas]

Taking a step back here as I'm a little confused now. Do you truly want bounds on the Projection, or are you actually looking for bounds/extents on each of your swaths/images? I took a quick look at the pyresample docs and found this in a few places under the Cartopy section plt.imshow(result, transform=crs, extent=crs.bounds, origin='upper', cmap='RdBu_r'), but it looks like that is the only thing that crs.bounds is being used for and it could be obtained from area_def.bounds or similar objects?

My concern in this case is that putting the bounds on Projection objects could lead to issues when adding multiple of these images to the same map because then the bounds of the map/axes would only be defined by one of the image projections. An alternative idea would be to subclass a feature or image source, such as LocatedImage or RasterSource, and then within Cartopy if someone calls ax.imshow(SubclassedImage), we could grab the Projection and extent from that object directly rather than requiring the user to know transform=crs, extent=crs.bounds ....

Just a thought I had, I may be way off base and missing something pretty fundamental too...

[djhoese]

You aren't wrong and it has been a while since I originally wanted this type of interaction between pyresample and cartopy, so it is very possible I'm forgetting some of the use cases. Looking at an old SciPy tutorial I made for Satpy (which uses pyresample underneath), I have code snippets like this:

import matplotlib.pyplot as plt

crs = scn[my_channel].attrs['area'].to_cartopy_crs()
plt.figure()
ax = plt.axes(projection=crs)

my_data = scn[my_channel]
my_data.plot.imshow(transform=crs)
ax.coastlines()
ax.gridlines()

Full notebook: https://github.com/pytroll/tutorial-satpy-half-day/blob/main/notebooks/07_cartopy_geoviews.ipynb

So I think the idea was that the user doesn't have to know extents. At least in use cases like this the extents on the map and the image data are the same thing; we're just viewing the image. Now of course looking at this with fresh eyes and more experience I could see someone wanting to say "I want to display this image but with a margin of 10km on each side" in which case the person would need to access the bounds of the crs.

My concern in this case is that putting the bounds on Projection objects could lead to issues when adding multiple of these images to the same map because then the bounds of the map/axes would only be defined by one of the image projections.

I'm not sure I understand what would be wrong with this. Plus it is fairly typical for our users to have all their data resampled to the same "area" already so the plotted data and the map projection are all the same (makes it so cartopy shouldn't have to resample or transform anything).

The main thing here that I think I need to figure out is what is the cartopy analog for a pyresample AreaDefinition. An AreaDefinition is CRS + bounds/extents + shape. So maybe the new and preferred way of plotting pyresample results should be something like:

crs = area_def.crs  # pyproj CRS
extents = area_def.bounds  # xmin, xmax, ymin, ymax
# use these in matplotlib/cartopy calls

@mraspaud Do you have any opinion on this?

I suppose the more generic use case for this feature is that some Projections won't have bounds so this keyword argument let's you, the knowledgeable user, define them.

[mraspaud]

I would just like to mention the area_of_use.bounds attribute in pyproj CRSs that is used for epsg codes that have bounds. Could that be helpful?

[djhoese]

@mraspaud That is what self.area_of_use is and currently exists in cartopy 0.20.0 based on the work @snowman2 did. I use it here as the default if bounds isn't specified.

[CLAassistant]

Thank you for your submission! We really appreciate it. Like many open source projects, we ask that you sign our Contributor License Agreement before we can accept your contribution.
_{You have signed the CLA already but the status is still pending? Let us recheck it.}