Create Cartopy projection from pyproj.Proj

[fmaussion]

Maybe I missed it, but since the issues #153 and #710 are still open I thought I would ask this too:

is it possible to create a cartopy projection from a pyproj.Proj instance?

[dopplershift]

No, but I'd love to have it too.

[fmaussion]

I've come up with a quick n dirty solution here

There should be a more generic way to do it (both pyproj and cartopy use the proj4 engine) but I don't know cartopy well enough to help here...

[graziano-giuliani]

I know that this is crude, but...
proposed_pj4.txt

[pelson]

Just to be clear on this, I am against using a proj4 string as a constructor of a "Proj4Projection" - the solution I'd be looking for is to construct a proper cartopy Projection instance from the proj4 string. For example:

>>> projection = ccrs.from_proj('+ellps=WGS84 +proj=robin +lon_0=40 +no_defs')
>>> print(projection)
Robinson(central_longitude=40)

[djhoese]

I'd be willing to work on this if everyone can agree on how it should be done. To do the above behavior where from_proj returns the proper CRS projection class it would probably be easiest if each individual class had a class attribute mapping keyword arguments to PROJ.4 parameters.

I was thinking the from_proj could try to find a matching CRS.Projection class and if it can't it could fall back to a generic PROJ4Projection type class. I noticed that the _EPSGProjection class is already almost a generic PROJ.4 class.

[fmaussion]

I cannot comment on the details here, but I can only emphasize again that a proj4 string -> cartopy converter would be extremely useful...

[pelson]

I'd be willing to work on this if everyone can agree on how it should be done. To do the above behavior where from_proj returns the proper CRS projection class it would probably be easiest if each individual class had a class attribute mapping keyword arguments to PROJ.4 parameters.

In general, I'm completely OK with that (and potentially even automatically adding those kwargs through to the __init__), but kwargs is not 1:1 proj4 arguments. For example https://github.com/SciTools/cartopy/blob/master/lib/cartopy/crs.py#L1277 'lon_0', 180 + pole_longitude. There will have to be special cases, and I'd like those special cases to belong to the respective classes, rather than having if crs == Robinson type conditions in a function.

[djhoese]

Without thinking about it too much maybe a from_proj classmethod on each CRS class? Then a global utility function from_proj that calls the individual from_proj methods.

[fmaussion]

but kwargs is not 1:1 proj4 arguments

Somewhat OT: what was the rationale behind creating cartopy-specific semantics for projection kwargs and names, instead of relying on PROJ4 conventions?

[pelson]

Somewhat OT: what was the rationale behind creating cartopy-specific semantics for projection kwargs and names, instead of relying on PROJ4 conventions?

It isn't something that was done lightly, but given there are a number of cases where the proj4 kwargs aren't the defacto standard within other fields. Rotated pole is a good example, where the pole longitude is more common than the central longitude in environmental science.

Without thinking about it too much maybe a from_proj classmethod on each CRS class?

Yep, sounds sensible. Another complication that I've remembered from having a 30 minute go at doing this several years ago... projection classes aren't 1:1 with proj.4 projections either. Case in point: NorthPolarStereographic. I can completely live with that though - we quickly descend get into a conversation about what is a projection and what is an instance of a projection 😄.

If you're really keen to have a go at this, I'd encourage you to produce something rough and ready, rather than polishing too much - it might take several iterations until we know we are on the right lines, and only then would I be looking to add some polish to the code. 👍

[djhoese]

@pelson Sounds good. I've started on the basics. I'll make a PR hopefully in the next hour or so.

[moccand]

Why not just "extending" the cartopy.crs.epsg(code) in order to offer two diffrent uses :

• EPSG code + call to epsg.io by pyepsg (as it is now)
• PROJ4 string made from whatever source we want

in _epsg.py there is :

class _EPSGProjection(ccrs.Projection):
    def __init__(self, code):
        import pyepsg
        (...)
        projection = pyepsg.get(code)
        (...)
        proj4_str = projection.as_proj4().strip()
        (...)

so the use of proj4's style str is already there.

Am i wrong ?

[djhoese]

@moccand As mentioned in my PR for adding this functionality, see #1023 (comment)

So yes, you're on the right track, but in the current EPSG code the projection object (projection = pyepsg.get(code)) is still used later on (https://github.com/SciTools/cartopy/blob/master/lib/cartopy/_epsg.py#L69).

It has been a long time since I've looked at this, but iirc this is essentially what I was doing. Making a PROJ.4 CRS object that the EPSG code could use. It got more complicated when @pelson pointed out that they would rather have "fully qualified" CRS objects returned rather than a single object that can be all projections.

[moccand]

I will take time to read previous posts.

regarding projection.domain_of_validity(), it returns the WGS84 bounds ("extent") as shown in epsg.io website or in

http://epsg.io/?q=2154&format=json&trans=1

for EPSG:2154 WGS84 bounds are :
-9.86 41.15
10.38 51.56

{"status": "ok", "number_result": 1, "results": [{"code": "2154", "kind": "CRS-PROJCRS", "bbox": [51.56, -9.86, 41.15, 10.38], "wkt": "PROJCS["RGF93 / Lambert-93",GEOGCS["RGF93",DATUM["Reseau_Geodesique_Francais_1993",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6171"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4171"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",49],PARAMETER["standard_parallel_2",44],PARAMETER["latitude_of_origin",46.5],PARAMETER["central_meridian",3],PARAMETER["false_easting",700000],PARAMETER["false_northing",6600000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","2154"]]", "unit": "metre", "proj4": "+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs", "name": "RGF93 / Lambert-93", "area": "France - onshore and offshore, mainland and Corsica.", "default_trans": 1671, "trans": [{"proj4": "", "bbox": [51.56, -9.86, 41.15, 10.38], "area": "France - onshore and offshore, mainland and Corsica.", "code_trans": 1671, "accuracy": 1.0, "wkt": "", "unit": "metre", "name": "RGF93 to WGS 84 (1)"}, {"proj4": "", "bbox": [51.56, -9.86, 41.15, 10.38], "area": "France - onshore and offshore, mainland and Corsica.", "code_trans": 8573, "accuracy": "", "wkt": "", "unit": "metre", "name": "RGF93 to WGS 84 (1)"}], "accuracy": 1.0}]}

>>> pyepsg.get(2154)
<ProjectedCRS: 2154, RGF93 / Lambert-93>
>>> pyepsg.get(2154).as_proj4().strip()
'+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs'
>>> pyepsg.get(2154).domain_of_validity()
[-9.86, 10.38, 41.15, 51.56]

setting both proj4 string and extents (in WGS84 or in targeted system with transformation) would not be so far from what openlayers (+ proj4js) needs :

https://openlayers.org/en/latest/examples/wms-image-custom-proj.html

      proj4.defs('EPSG:21781',
        '+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 ' +
          '+x_0=600000 +y_0=200000 +ellps=bessel ' +
          '+towgs84=660.077,13.551,369.344,2.484,1.783,2.939,5.66 +units=m +no_defs');
      register(proj4);

      var projection = new Projection({
        code: 'EPSG:21781',
        extent: [485869.5728, 76443.1884, 837076.5648, 299941.7864]
      });

[rabernat]

👍 to the increased interop btw cartopy and pyproj. Thanks to those who are working on this.

[ecoon]

It looks like the corresponding pull request made a ton of progress on this. It would be a shame to lose that! Thanks all (+1/bump/pretty please) :)

[vedal]

Thank you for the work already done in this direction! Any further activity on this?

[djhoese]

@amundvedal The last I heard on this, at least on the "wish list" level of things, was the idea to try and make pyproj's relatively new CRS objects usable in Cartopy. This isn't implemented as far as I know but is something @snowman2 and @dopplershift and I talked about at the last SciPy conference. I think there is a separate issue talking about it but I can't find it right now.

Edit: Found it #1477

[vedal]

@djhoese Thank you for the answer! I found the following away around my issue, and posting it here to help others in a similar situation:

Issue:
I wanted to plot some geopandasLineString as paths using gv.Path, but their crs EPSG:3006 was not supported by geoviews.

Solution:
I found the function GeoDataFrame.to_crs that could convert my LineString to epsg=3857, which is the same as cartopy.crs.GOOGLE_MERCATOR. Thus, I was able to plot the LineString correctly

Code:

import geopandas
from cartopy.crs import GOOGLE_MERCATOR
paths = gv.Path(gpd.read_file('shapefile.shp').to_crs(epsg=3857), crs=GOOGLE_MERCATOR)

[djhoese]

With #1808 merged I think the original issue here can be closed, but it also doesn't cover all cases. I'm wondering if I can get @snowman2 and the rest of the community here to help me add a little more functionality to completely handle my use case (and I think @fmaussion's).

The Projection class is now a subclass of the pyproj CRS class so doing Projection(some_pyproj_crs) works. In my pyresample library we want to take any pyproj CRS but use our own custom bounds. The use case is loading some geographic data from a data file that has a limited set of extents and use those in the Projection definition. So my questions for @snowman2 and the rest of you are:

1. Is it a bad idea to limit a Projection/CRS instance in cartopy to a subset of the projection space? Would this cause issues with drawing certain cartopy features?
2. Previously, .bounds/.boundary/.threshold information had to be in radians for geographic projections. Should these now be in degrees since pyproj defaults to keeping things in degrees (I think).
3. Would an optional bounds keyword argument to Projection be good enough?

class Projection(CRS, metaclass=ABCMeta):
    ...

    def __init__(self, *args, bounds=None, transform_bounds=True, **kwargs):
        super().__init__(*args, **kwargs)
        self.bounds = None
        bounds = bounds if bounds is not None else self.area_of_use
        if bounds and transform_bounds:
            # existing transformation logic

This may kind of break the idea of a "Projection" though, but also fits decently well in my opinion.

Side note: The existing self.bounds = and self.threshold = will have issues if bounds are flipped which is true for some satellite data (SEVIRI).

[snowman2]

I like the idea of 3 (bounds kwarg for Projection). Adding in 3 and tinkering a bit may help answer 1 & 2.

[greglucas]

I don't think that Cartopy actually restricts the transforms based on the bounds, so setting the bounds would only affect the drawn boundary. For example, if you limit the extent of your projection but then start zooming/panning around and your dataset was outside of the initial extent it will come into view later. If you do want to restrict the bounds for the transforms maybe you'd need to update the area_of_use that pyproj uses?

It looks like the bounds currently have to be specified in degrees since that is what the area_of_use is in?

cartopy/lib/cartopy/crs.py

Line 651 in 019f013

# Convert lat/lon bounds to projected bounds.

3 seems pretty reasonable and wouldn't be very invasive.

[djhoese]

I'm testing something to be included in pyresample since it currently doesn't work with cartopy 0.20.0. We previously had a custom subclass of the cartopy CRS which sparked my work in #1023. Here's the replacement that I'll try to submit as a PR tomorrow if I can find the time:

class _Projection(ccrs.Projection):
    def __init__(self,
                 crs: pyproj.CRS,
                 bounds: list[float, float, float, float] = None,
                 transform_bounds: bool = False):
        # XXX: Skip the base cartopy Projection class __init__
        super(ccrs.Projection, self).__init__(crs)
        if bounds is None and self.area_of_use is not None:
            bounds = (
                self.area_of_use.west,
                self.area_of_use.east,
                self.area_of_use.south,
                self.area_of_use.north,
            )
            transform_bounds = True

        self.bounds = bounds
        if bounds is not None and transform_bounds:
            # Convert lat/lon bounds to projected bounds.
            # Geographic area of the entire dataset referenced to WGS 84
            # NB. We can't use a polygon transform at this stage because
            # that relies on the existence of the map boundary... the very
            # thing we're trying to work out! ;-)
            x0, x1, y0, y1 = bounds
            lons = np.array([x0, x0, x1, x1])
            lats = np.array([y0, y1, y1, y0])
            points = self.transform_points(self.as_geodetic(), lons, lats)
            x = points[:, 0]
            y = points[:, 1]
            self.bounds = (x.min(), x.max(), y.min(), y.max())
        if self.bounds is not None:
            x0, x1, y0, y1 = self.bounds
            self.threshold = min(abs(x1 - x0), abs(y1 - y0)) / 100.

Seems to work for pyresample's use cases so far. The transform_bounds keyword argument is a convenience so that users who know their bounds in lon/lat space don't have to use pyproj to transform the points to the CRS space, let cartopy do it for you.

[djhoese]

Also as to answer number 2, it seems that for a CRS for EPSG:4326 that the bounds are stored as degrees so I think radians being used is no longer needed (at least I had to use them in pyresample's code) for geographic CRS objects or at least the ones that don't explicitly define themselves as being in radians.

[djhoese]

See #1888