Cartopy: Create Cartopy projection from pyproj.Proj

[fmaussion]

Maybe I missed it, but since the issues https://github.com/SciTools/cartopy/issues/153 and https://github.com/SciTools/cartopy/issues/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 https://github.com/SciTools/cartopy/pull/1023#discussion_r168395702

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 https://github.com/SciTools/cartopy/issues/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)