-
Notifications
You must be signed in to change notification settings - Fork 286
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
Add lambert azimuthal equal area coord_system #1738
Conversation
X offset from planar origin in metres. | ||
|
||
* false_northing | ||
Y offset from planar origin in metres. |
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.
Please could you add "Defaults to 0." for these arguments.
👍 Now it just needs a new Cartopy release... |
It probably also needs an entry to the what's new. |
Yes, good point. 👍 |
951b4a3
to
5441678
Compare
@bjlittle , this is the pull request I was talking about that is now needed. |
(iris.coord_systems.TransverseMercator, | ||
iris.coord_systems.LambertAzimuthalEqualArea)): | ||
units = 'm' | ||
|
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 block was removed in PR #1943 but I don't agree that it should have been removed.
Section 4.1 and 4.2 of the cf conventions says:
Coordinates of latitude with respect to a rotated pole should be given units of degrees
also 'm' is a perfectly fine unit (see the example in this section).
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.
@lbdreyer Interesting that swapping this back doesn't cause any test failures ...
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.
In fact, if you look at Example 5.6, that's a rotated_lat_lon grid with units of "degrees". Also Example 5.10 is a tranverse_mercator grid with units of "m".
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.
Right, I think I must have gotten confused trying to rebase a PR that was over a year old!
On closer inspection I now see that, in the case of tranverse_mercator we're just passing the units from the coord. Although this has it's disadvantages (i.e. it means that a user could produce a netcdf file with a badly defined coordinate:
e.g. DimCoord(array([30, 60]), standard_name='latitude', units=Unit('m'), coord_system=RotatedGeogCS(10.0, 20.0))
But I would prefer for users to have this flexibility, and cf says that there is no default value for lat or lon units.
As the CF conventions recommend a specific unit for GeogCS and RotatedGeogCS I think you could argue that we should 'enforce' these:
The recommended unit of latitude is degrees_north
and
Coordinates of latitude with respect to a rotated pole should be given units of degrees
Failing because of this: |
That failure was fixed in #2144, targeted at the v1.10.x branch. This change needs to make its way onto master, and then you can rebase onto that to get passing tests here (or we can trust that this feature doesn't cause the failure). There is another test failure too (an error actually), not sure what that one is caused by. |
I am going to add a test for the coordinate units discussed in the above comment (I'm hoping that's not a controversial change), and then rebase. |
5441678
to
0af35ac
Compare
def check(self, coord_system, units, expected_units): | ||
coord = iris.coords.DimCoord([30, 45], 'latitude', units=units, | ||
coord_system=coord_system) | ||
saver = Saver('filename', 'NETCDF4') |
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.
Argh! This is wrong. I was going to mock it out.
I had made _cf_coord_identity a static method in the hopes that I wouldn't have to instantiate the Saver class to test it, but I was getting an error
>>> coord = iris.coords.DimCoord([45], 'latitude', units='m')
>>> r = Saver._cf_coord_identity(coord)
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
TypeError: unbound method _cf_coord_identity() must be called with Saver instance as first argument (got DimCoord instance instead)
35883f7
to
9e3a3b8
Compare
@bjlittle The tests are now passing and the static method is working correctly. |
Hi @lbdreyer i'm afraid something has occurred on master such that your branch now has conflicts please rebase and fix so we can see a plausible change with test results are you still waiting on a cartopy release for this to be able to work properly? mark |
9e3a3b8
to
cde0ac4
Compare
I have rebased and the tests are passing again.
No, the change that this required was in the Cartopy 0.14 release |
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.
Some minor changes, and some feedback from @djkirkham and @marqh would be welcome regarding rotated pole and (CF spec) units of degrees
(iris.coord_systems.TransverseMercator, | ||
iris.coord_systems.LambertAzimuthalEqualArea)): | ||
units = 'm' | ||
|
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.
@lbdreyer Interesting that swapping this back doesn't cause any test failures ...
@@ -1327,7 +1328,10 @@ def _cf_coord_identity(self, coord): | |||
elif coord.standard_name == "longitude": | |||
units = 'degrees_east' | |||
|
|||
return (coord.standard_name, coord.long_name, units) | |||
elif isinstance(coord.coord_system, iris.coord_systems.RotatedGeogCS): | |||
units = 'degrees' |
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.
@lbdreyer As per the CF spec sections 4.1 and 4.2.
@djkirkham and @marqh This reverts your changes from #1943. Can we all agree this is what we want for rotated pole units, as per CF.
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.
If I remember correctly, we took this out because the units may be radians, say, so changing to degrees would be incorrect.
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.
Thanks @djkirkham, that makes sense. I will remove this.
self.semi_major_axis = 6377563.396 | ||
self.semi_minor_axis = 6356256.909 | ||
self.ellipsoid = GeogCS(self.semi_major_axis, self.semi_minor_axis) | ||
self.laea_cs = LambertAzimuthalEqualArea( |
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.
@lbdreyer Add in the false_easting
and false_northing
... thanks!
def check_call(self, coord_system, units, expected_units): | ||
coord = iris.coords.DimCoord([30, 45], 'latitude', units=units, | ||
coord_system=coord_system) | ||
result = Saver._cf_coord_identity(coord) |
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.
@lbdreyer Hence the need for @staticmethod
, okay ... 😄
class Test__cf_coord_identity(tests.IrisTest): | ||
def check_call(self, coord_system, units, expected_units): | ||
coord = iris.coords.DimCoord([30, 45], 'latitude', units=units, | ||
coord_system=coord_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.
@lbdreyer Ideally, to give full coverage, could you also test longitude
also ... thanks!
self.check_call(coord_system=crs, units='degrees', | ||
expected_units='degrees') | ||
|
||
def test_crs_with_no_default_units(self): |
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.
@lbdreyer This is really a pass-thru
rather than default
behaviour ... what do you think?
latitude_of_projection_origin=52, | ||
longitude_of_projection_origin=10, | ||
false_easting=100, | ||
false_northing=200) |
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.
@lbdreyer Did you want to pass-thru a choice of ellipsoid
here ...
5bbe26e
to
fa72ee7
Compare
@lbdreyer Great work! Thanks 👍 |
Depends on SciTools/iris-test-data#37 and SciTools/cartopy#619