-
Notifications
You must be signed in to change notification settings - Fork 365
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
WIP: Add from_proj4
function to create CRS from PROJ.4 definitions
#1023
Conversation
def threshold(self): | ||
x0, x1, y0, y1 = self.bounds | ||
return min(x1 - x0, y1 - y0) / 100. | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
W391 blank line at end of file
lib/cartopy/crs.py
Outdated
elif crs_class is None: | ||
return _PROJ4Projection(proj4_dict, globe=globe, bounds=bounds) | ||
|
||
return crs_class.from_proj4(proj4_dict) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
W292 no newline at end of file
@pelson CLA has been sent. Took me a bit to get used to all the code checks, but when you have time take a look and tell me what you hate. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is a really great start. We need to flesh out the from_proj4
classmethod some more, but really appreciate you showing this early to help steer this in the right direction.
I'm not sure if I should show you this or not, but I went digging an found a commit way back when that attempted to convert cartopy<>proj4 strings. I didn't actually follow through with it but the diff can be seen at pelson/cartopy@b15391d...pelson:repr_from_proj4_init. Would be interested in getting your take on it, and welcome you to lift whatever you like from there into this (or subsequent) pull request.
lib/cartopy/_proj4.py
Outdated
class _PROJ4Projection(ccrs.Projection): | ||
def __init__(self, proj4_terms, globe=None, bounds=None): | ||
terms = get_proj4_dict(proj4_terms) | ||
globe = _globe_from_proj4(terms) if globe is None else globe |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
We need to be careful in case there are globe terms in proj4 and a globe is provided.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Oh good point. I'll fix that.
"""Create a `Globe` object from PROJ.4 parameters.""" | ||
globe_terms = filter(lambda term: term[0] in _GLOBE_PARAMS, | ||
proj4_terms.items()) | ||
globe = ccrs.Globe(**{_GLOBE_PARAMS[name]: value for name, value in |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Style choice... I'd probably get rid of the filter:
globe = ccrs.Globe(**{_GLOBE_PARAMS[name]: value for name, value in proj4_terms.items()
if name in _GLOBE_PARAMS})
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This was copied directly from the EPSG class ;)
other_terms.append(term) | ||
super(_PROJ4Projection, self).__init__(other_terms, globe) | ||
|
||
# FIXME: Can we guess at the bounds if not provided? Maybe transform |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yep, we can do some automatic bounds calculation IFF the target projection doesn't define them itself.
I wouldn't want some projections to be allowed to set its own bounds though (e.g. PlateCarree).
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Isn't this how EPSG works? Isn't this also what a lot of the CRS classes already do? They limit or calculate what the bounds are inside the CRS class?
lib/cartopy/_proj4.py
Outdated
terms = get_proj4_dict(proj4_terms) | ||
globe = _globe_from_proj4(terms) if globe is None else globe | ||
|
||
other_terms = [] |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Call these projection_terms.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Will do. Most of this was copied from the EPSG code.
lib/cartopy/_proj4.py
Outdated
|
||
other_terms = [] | ||
for term in terms.items(): | ||
if term[0] not in _GLOBE_PARAMS: |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Might be nice to allow globe_from_proj4 to manipulate the terms (and it can therefore remove the globe ones).
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Makes sense, especially when you consider the comment above about filtering out global parameters.
lib/cartopy/crs.py
Outdated
@@ -102,6 +102,11 @@ class Projection(six.with_metaclass(ABCMeta, CRS)): | |||
'MultiPolygon': '_project_multipolygon', | |||
} | |||
|
|||
@classmethod | |||
def from_proj4(cls, proj4_dict, **kwargs): | |||
raise NotImplementedError("This projection can not be created from a " |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
We can be friendly by replacing This projection
with the cls.__name__
.
if 'x_0' in proj4_dict: | ||
p_kwargs['false_easting'] = float(proj4_dict['x_0']) | ||
if 'y_0' in proj4_dict: | ||
p_kwargs['false_northing'] = float(proj4_dict['y_0']) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
There is repeated handling of things like false_northing
happening. Can we factorise this?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I thought about that, but for this initial commit I didn't do it because it felt like I would be trying to hard to reduce duplicate code and I wasn't sure if all CRS classes supported the false easting and northing.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Not all projections do support false_easting / false_northing, but those that do all map to y_0
. In some respects, if the proj4_terms
includes this, then whether or not the projection supports it I'd be tempted to try passing the keyword on to the __init__
anyway. At least then the user gets a message about false_northing
not being a valid keyword.
What would be missing from this approach is the user knowing that y_0
was the term that mapped to false_northing
. Perhaps we could make use of inspect
to determine ahead-of-time what the __init__
supports?
lib/cartopy/crs.py
Outdated
crs_class = proj_to_crs.get(proj) | ||
|
||
# special cases | ||
if proj == 'stere' and float(proj4_dict['lat_0']) < 0: |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Can we push that into the Stereographic class?
lib/cartopy/crs.py
Outdated
raise ValueError("'bounds' must be specified for unknown CRS " | ||
"'{}'".format(proj)) | ||
elif crs_class is None: | ||
return _PROJ4Projection(proj4_dict, globe=globe, bounds=bounds) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm 👎 for creating a projection that cartopy hasn't implemented. (i.e. returning a generic _PROJ4Projection
instance)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think we're going to disagree on this, but in the end it is your call. You mentioned in the issue that cartopy doesn't follow PROJ.4 exactly (not a 1:1 mapping) because some scientific fields have existing expectations about how to describe a CRS. In a similar line of thinking I would argue that some fields use PROJ.4 to describe everything so limiting a user to projections that have a custom CRS class seems odd. There will be projections that are available that cartopy doesn't support, if all of the abstract methods are implemented shouldn't a user be able to use whatever projection they want? If not, why support EPSG at all?
I understand wanting to have a nice CRS class name and everything, but if a PROJ.4 generic fallback class is available and the functionality is the same then does it really matter? I guess it would effectively stop development on new CRS classes for cartopy though. </rant>
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I should also mentioned that I'm the maintainer of https://github.com/pytroll/satpy and https://github.com/pytroll/pyresample. The reason I wanted to add this functionality was to make it easier for our users to take our AreaDefinition
(PROJ.4 projection + bounds + number of pixels) and create a cartopy plot from them.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Also if this doesn't get added to cartopy I will likely just add it to pyresample.
|
||
|
||
class _EPSGProjection(ccrs.Projection): | ||
class _EPSGProjection(_PROJ4Projection): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Seems highly likely that we will remove _EPSGProjection and replace it with a function that returns an appropriate proj4 string > a real cartopy.crc.CRS.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
To be clear: nothing for you to do here, just sharing my thoughts.
lib/cartopy/crs.py
Outdated
from cartopy._proj4 import get_proj4_dict, _PROJ4Projection | ||
proj4_dict = get_proj4_dict(proj4_terms) | ||
|
||
proj_to_crs = { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ideally I'd like this to be derived from the projection classes, not a static dictionary, if possible.
My branch does this at pelson/cartopy@b15391d...pelson:repr_from_proj4_init#diff-2af330a440b389503b4c32a5fbd26471R54 with:
+ def all_subclasses(cls):
+ return cls.__subclasses__() + [g for s in cls.__subclasses__()
+ for g in all_subclasses(s)]
+
+ potentials = []
+ for cls in all_subclasses(CRS):
+ if getattr(cls, '_proj4_proj', None) == proj:
+ potentials.append(cls)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ok, I can do this. I wasn't sure how "hacky" you wanted to get. ;)
lib/cartopy/crs.py
Outdated
@@ -2006,3 +2066,28 @@ def epsg(code): | |||
""" | |||
import cartopy._epsg | |||
return cartopy._epsg._EPSGProjection(code) | |||
|
|||
|
|||
def from_proj4(proj4_terms, globe=None, bounds=None): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Let's put the implementation in the _proj4 module and make this a stub function that calls the appropriate function there (crs.py
is already waaaay too long 😄).
from_proj4
function to create CRS from PROJ.4 definitionsfrom_proj4
function to create CRS from PROJ.4 definitions
Thanks a lot! Looking forward to use this ;) Is this going to work with EPSG projection IDs as well? |
@fmaussion Right now, this PR probably won't handle this. As @pelson mentioned the |
@pelson I've done just about every change suggested. Let me know what you think. It seems overall better and somehow worse than before. |
@pelson Do you think I should continue implementing the rest of the projections or is there something that still feels off to you about how I implemented this? |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Do you think I should continue implementing the rest of the projections or is there something that still feels off to you about how I implemented this?
Yes, I think you are heading in the right direction. Wherever possible I'd like to keep the proj4 detail to a minimum in crs.py
and the CRS
classes in general. That means putting much of the mapping from proj4 terms to CRS terms in _proj4.py
. The example we have discussed is false_northing
(but almost all other terms apply too) - ideally each CRS should not be implementing that mapping itself, but is should have the power to override the mapping if it really must.
I'm really excited about moving this forwards upon my return ( 🛫 🌞 🌴 🍸 🛬 ).
Thanks again for moving it forwards!
def __init__(self, proj4_terms, globe=None, bounds=None): | ||
terms = get_proj4_dict(proj4_terms) | ||
projection_terms, globe_terms = _split_globe_parameters(terms) | ||
if globe is None: |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think we need to raise an exception here is globe is not None and globe_terms
is not empty.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I was following the base CRS class which uses Globe()
as its default. I can add an exception, but with this default in mind, does that change your opinion?
return globe | ||
|
||
|
||
class _PROJ4Projection(ccrs.Projection): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think I'm 👎 on this class existing altogether. (much as I am the EPSG class that already does).
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Well right now the EPSG class is based on this class. So either I revert the EPSG class and remove this or I leave this in. I was thinking this could be made public even. It doesn't need to be documented even, but if you want the functionality removed then ok.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
FYI I've added this class to pyresample to convert pyresamples AreaDefinition
s to CRS objects: pytroll/pyresample#102
lib/cartopy/_proj4.py
Outdated
return min(x1 - x0, y1 - y0) / 100. | ||
|
||
|
||
def all_subclasses(cls): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This may as well be private I think.
lib/cartopy/_proj4.py
Outdated
for g in all_subclasses(s)] | ||
|
||
|
||
def from_proj4(proj4_terms, globe=None, bounds=None): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Is there a situation where we need to pass globe through here? It adds redundancy that we need to check against if it remains.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
You're right. Now that the default isn't to use the _Proj4Projection
class as a last resort both of these keywords don't need to be here.
def from_proj4(proj4_terms, globe=None, bounds=None): | ||
proj4_dict = get_proj4_dict(proj4_terms) | ||
|
||
if not PROJ_TO_CRS: |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
PROJ_TO_CRS definitely needs a docstring - it wasn't clear to me that it was a mapping of the proj4 proj
term to coordinate reference system.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
With regards to the case of projection specialisation, do you have an idea in mind for handling the case where multiple CRS classes implement the same projection. This is true of cases such as NorthPolarStereographic (https://github.com/SciTools/cartopy/blob/master/lib/cartopy/crs.py#L1376).
There is potentially a case to be made that such "projections" shouldn't exist at all - they are perhaps better modelled as instances of their base projections...
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I guess I'll find out when I get there. I already have it implemented for Stereographic where the North and South cases are handled inside the class method.
if cls_proj is not None and cls_proj not in PROJ_TO_CRS: | ||
PROJ_TO_CRS[cls_proj] = crs_class | ||
|
||
proj = proj4_dict['proj'] |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
May want to be slightly more defensive here - what if proj
isn't in the dict.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Then you get a KeyError because your PROJ.4 string doesn't make sense...but ok.
@@ -1134,6 +1141,35 @@ def __init__(self, central_longitude=-96.0, central_latitude=39.0, | |||
self._x_limits = bounds[0], bounds[2] | |||
self._y_limits = bounds[1], bounds[3] | |||
|
|||
@classmethod | |||
def from_proj4(cls, proj4_dict, **kwargs): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'd want to assert that the proj
parameter is set correctly (in this case, lcc
).
if 'x_0' in proj4_dict: | ||
p_kwargs['false_easting'] = float(proj4_dict['x_0']) | ||
if 'y_0' in proj4_dict: | ||
p_kwargs['false_northing'] = float(proj4_dict['y_0']) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Not all projections do support false_easting / false_northing, but those that do all map to y_0
. In some respects, if the proj4_terms
includes this, then whether or not the projection supports it I'd be tempted to try passing the keyword on to the __init__
anyway. At least then the user gets a message about false_northing
not being a valid keyword.
What would be missing from this approach is the user knowing that y_0
was the term that mapped to false_northing
. Perhaps we could make use of inspect
to determine ahead-of-time what the __init__
supports?
@pelson I think I see what you are saying now. Maybe the
|
I'd love to move this forward somehow. What are the next steps? Any help needed? |
@knedlsepp I think I need @pelson to look at this again and decide what the best interface is for this and re-review what I have already. After that I can fix up anything that needs fixing (merge conflicts, etc). |
Ugh and I changed my username from |
@pelson Do you need me to do anything special to get the CLA checker fixed? Where would you like this PR to go from here? |
@djhoese I have updated the contributors .json with your new username 👍
@pelson Could you review the code again please? |
I've found it really hard to get the bandwidth to do much on this recently, so apologies for the silence.
|
@pelson Maybe it has been too long for me, but I think most of these things are already handled as you want them. I'm not sure if you are saying you like the way they are done or want them changed.
Are you saying you want the Projection subclasses to handle all of the transform? I have implemented a
This is already handled by special handling in the base
So if any proj parameters are provided that cartopy doesn't know how to handle then an exception is raised? To be clear, exception over warning or something similar?
The class needs to exist to support the
I'm not sure I understand. What type of |
Latest commit is just a copy from the pyresample project where a contributor brought up a specific issue (see PR and related issue): pytroll/pyresample#161 |
This would be a highly welcome functionality. Is there a plan to add this soon? |
I think this is just waiting on feedback from the maintainers. That said, at the last SciPy conference @dopplershift, @snowman2, and I talked about perhaps making cartopy more directly compatible with pyproj's @snowman2 Any further ideas on that? |
I think it would mostly entail using the |
I'm interested in moving to complete reliance on PyProj and eliminating CartoPy's code that links to libproj. It would help simplify many things. I think that needs to happen for 0.19, though, since I'd like to get 0.18 out just as soon as I can clear the cycles. |
I think now that #1808 has been merged this PR can be closed. My understanding is that we can do I'll let one of the maintainers close this. @snowman2 please correct me if I'm wrong. |
That is correct. It is important to note that |
In pyresample we have "AreaDefinition" objects which are essentially a pyproj CRS + extents + size. Is it possible to create a cartopy CRS with a pyproj CRS and the x/y limits? |
The area of use can be set with a spatial reference code, WKT, or PROJ JSON. I think this would be the simplest way to support custom CRS with an area of use: However, pyproj will need to be updated to pass area of use to the PROJ JSON. |
Rationale
See #813 for more details. Many users and other software libraries use PROJ.4 definitions to define their data projections. To convert from these definitions to a cartopy CRS requires knowing the name of the CRS class and mapping PROJ.4 parameters to keyword arguments of the CRS class. This PR attempts to simplify this process by adding a
from_proj4
function similar to theepsg
function. The first commit of this PR only implements the functionality for the Stereographic and LambertConformal projections.Implications
This PR separates out the PROJ.4 parsing code from the
_EPSGProjection
class and adds a_proj4.py
python module to hold a new_PROJ4Projection
class as a generic fallback class for any PROJ.4 strings that couldn't be mapped to a specific CRS class.This PR also requires adding a new
classmethod
to allProjection
subclasses calledfrom_proj4
.Checklist
from_proj4
implemented and tests writtenProjection
subclass)?