-
Notifications
You must be signed in to change notification settings - Fork 301
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
Switch to pyproj for projection to CF NetCDF grid mapping #1038
Conversation
I'd like to use the opportunity to add an oblique mercator test. I'll try to provide something during the day. |
I learned while working on the short term solution last night that I'm pretty sure the projection tests aren't actually comparing any of the |
I'll have a look now |
you're right, the attrs weren't checked... |
Codecov Report
@@ Coverage Diff @@
## master #1038 +/- ##
==========================================
+ Coverage 89.47% 89.53% +0.05%
==========================================
Files 198 200 +2
Lines 29210 29356 +146
==========================================
+ Hits 26136 26284 +148
+ Misses 3074 3072 -2
Continue to review full report at Codecov.
|
02986d1
to
10337c1
Compare
So I've just rebased this on to master and I have a few things I think we need to decide on before merging this:
Edit: Also see pyproj4/pyproj#515 |
|
Here is the result of a before/after comparison with omerc: Before
After
It's basically fine, but there are some parameters missing (both in the before and after), since CF 1.8 says:
|
It's better with pyproj master, but geographic_crs_name is still missing
|
|
http://cfconventions.org/Data/cf-conventions/cf-conventions-1.8/cf-conventions.html it is there, in table F.1 |
Hm, but for your projection it isn't defined by the OGC WKT spec, so maybe that's ok? The CF docs don't say anything about "unknown" when it isn't defined. I assume that means you just don't include it if it isn't defined which is what pyproj is doing here. What if you do have a OGC WKT GEOGCS defined? Does pyproj include it then? |
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.
Very nice development! And thanks for getting in touch with the pyproj developers.
- Agree with Martin.
- Never used that projection.
- See my comment.
satpy/writers/cf_writer.py
Outdated
grid_mapping['name'] = grid_mapping['grid_mapping_name'] | ||
except AttributeError: | ||
grid_mapping = mappings[area.proj_dict['proj']](area) | ||
grid_mapping['name'] = area.proj_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.
Does that mean that the name
attribute would depend on the pyproj version? If so, should we align those names?
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.
My suggestion for 3 is not implemented. Right now this does not depend on the pyproj version because proj
is statically defined and unchanging (merc
, geos
, etc). However, the name
attribute is not required by CF for the grid mapping variable so we can remove this portion and instead use area.area_id
as the name of the variable.
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 see. Using area.area_id
sounds good!
@djhoese I don't think there is much more to do here, I'm glad with the result of this actually. If there are improvements to make, I suppose they should go in pyproj, no here, right ? |
There still need to be some changes for what |
I'm good with depending in a new pyproj. But maybe we skip this PR for v0.21 then. |
I'm OK with that. Are you saying that due to the time it would take me to finish this or because you don't want this writer to depend on pyproj 2+? |
I was thinking time constraints, I'd like a release within a week. |
# Conflicts: # .travis.yml # appveyor.yml # satpy/tests/writer_tests/test_cf.py
DeepCode's analysis on #11ae11 found:1 minor issue. ✔️ 1 issue were fixed. 👉 View analysis in DeepCode’s Dashboard |
@mraspaud This is ready for re-review. The big change is that the CRS conversion now completely depends on pyproj. The biggest change here on the user side is that there is no longer a warning if a CRS can't be converted "properly". This is because pyproj follows what the current CF standard is (not sure it is fully documented but I saw it on the mailing list) which is to provide a @sfinkens This also implements the grid mapping variable naming talked about above. Basically, it now names the grid mapping variable to |
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 looks good to me as afaics, but we need to check how this performs against the cf-checker.
Moreover, @sfinkens are you going to be affected by these changes ?
satpy/tests/writer_tests/test_cf.py
Outdated
with mock.patch('satpy.writers.cf_writer.warnings.warn') as warn: | ||
res, grid_mapping = area2gridmapping(ds) | ||
warn.assert_called() | ||
proj_dict = pyresample.geometry.proj4_str_to_dict(res.attrs['grid_proj4']) | ||
self.assertEqual(proj_dict['lon_0'], 4.535) | ||
self.assertEqual(proj_dict['lat_0'], 46.0) | ||
self.assertEqual(proj_dict['o_lon_p'], -5.465) | ||
self.assertEqual(proj_dict['o_lat_p'], 90.0) | ||
self.assertEqual(proj_dict['proj'], 'ob_tran') | ||
self.assertEqual(proj_dict['o_proj'], 'stere') | ||
self.assertEqual(proj_dict['ellps'], 'WGS84') | ||
self.assertEqual(grid_mapping.attrs['name'], 'proj4') | ||
# with mock.patch('satpy.writers.cf_writer.warnings.warn') as warn: | ||
res, grid_mapping = area2gridmapping(ds) | ||
# warn.assert_called() | ||
# proj_dict = pyresample.geometry.proj4_str_to_dict(res.attrs['grid_proj4']) | ||
# self.assertEqual(proj_dict['lon_0'], 4.535) | ||
# self.assertEqual(proj_dict['lat_0'], 46.0) | ||
# self.assertEqual(proj_dict['o_lon_p'], -5.465) | ||
# self.assertEqual(proj_dict['o_lat_p'], 90.0) | ||
# self.assertEqual(proj_dict['proj'], 'ob_tran') | ||
# self.assertEqual(proj_dict['o_proj'], 'stere') | ||
# self.assertEqual(proj_dict['ellps'], 'WGS84') | ||
self.assertIn('crs_wkt', grid_mapping.attrs) | ||
self.assertEqual(res.attrs['grid_mapping'], 'cosmo7') |
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.
@loreclem I believe you implemented this, would you mind checking if you are ok with this ?
The tests were failing previously because |
Ok I just pushed a commit with some more big changes after @peters77 pointed out that this PR wasn't actually working in a real world case. Turns out the tests weren't actually testing very much. The lesson from us (the developers) should be "just because you can mock something in a test, doesn't mean you should". There are ways with mocks to declare them with The other big change in this last commit: the Lastly, |
Regarding my comment on geographic_crs_name, an issue has been opened in pyproj: |
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, this looks quite good now. I still have a few comments, but nothing serious.
@@ -759,14 +754,14 @@ def _gm_matches(gmapping, expected): | |||
'grid_mapping_name': 'geostationary', | |||
'semi_major_axis': a, | |||
'semi_minor_axis': b, | |||
'sweep_axis': None, | |||
'name': 'geos'}) | |||
# 'sweep_angle_axis': 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.
In my offline test, I get 'sweep_angle_axis': 'y'
. Should this be added ?
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 sweep angle axis is only guaranteed with pyproj 2.5+ since that version now bases the CF output on the WKT version of the CRS instead of the PROJ string. Ok with requiring pyproj 2.5+ for this writer?
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 fine with it, but it might be too recent. Could we just check that sweep_angle_axis
is y
if it is there ?
satpy/tests/writer_tests/test_cf.py
Outdated
# warn.assert_called() | ||
# proj_dict = pyresample.geometry.proj4_str_to_dict(res.attrs['grid_proj4']) | ||
# self.assertEqual(proj_dict['lon_0'], 4.535) | ||
# self.assertEqual(proj_dict['lat_0'], 46.0) | ||
# self.assertEqual(proj_dict['o_lon_p'], -5.465) | ||
# self.assertEqual(proj_dict['o_lat_p'], 90.0) | ||
# self.assertEqual(proj_dict['proj'], 'ob_tran') | ||
# self.assertEqual(proj_dict['o_proj'], 'stere') | ||
# self.assertEqual(proj_dict['ellps'], 'WGS84') |
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 isn't useful anymore and should be removed
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.
Yes, I just wanted to make sure it was obvious that this is no longer being done. I can remove it.
# self.assertEqual(proj_dict['proj'], 'ob_tran') | ||
# self.assertEqual(proj_dict['o_proj'], 'stere') | ||
# self.assertEqual(proj_dict['ellps'], 'WGS84') | ||
self.assertIn('crs_wkt', grid_mapping.attrs) |
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.
Should we copy the generated wkt and test against it ?
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.
Maybe not the whole string but I'd be comfortable with checking if certain parts are there. Just in case the WKT changes a bit I don't want to have to update these tests every time.
satpy/writers/cf_writer.py
Outdated
# modifies dataarray in-place | ||
area2lonlat(dataarray) |
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.
For style, I'd rather avoid in-place modifications. Can we make this return a modified copy of the array and do:
# modifies dataarray in-place | |
area2lonlat(dataarray) | |
dataarray = area2lonlat(dataarray) |
?
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 know we had this before, but since we're changing this anyway...
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.
Sure. The way this is set up is that area2lonlat
and area2gridmapping
are given a copy by area2cf
. So these two functions are in-place but are only called with a copy. Still want the change? Not a big deal.
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.
yes, I'd rather have each function make a copy.
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 the lonlat case we add coordinates, does that still only need a shallow copy?
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 so, but don't take my word for it.
satpy/writers/cf_writer.py
Outdated
res.extend(area2gridmapping(dataarray)) | ||
res.append(area2gridmapping(dataarray)) |
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.
Same here, can we have area2gridmapping
also return a modified copy of the array ? In this case we probably just need a shallow copy as we just add an attribute.
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.
LGTM
This is the long term solution for fixing #1029 (FYI @RickKohrs).
This switches the CF writer to depend on pyproj 2.0+ for converting the PROJ information to a CF compatible set of attributes. This is a no brainer in the long run as it lets the projection-specific library handling the projection logic (and people who are much more up to date on changes in CF and PROJ).
However, this has some important changes and this PR is kind of 50/50 on deciding how it handles them. I'm not sure whether we want to completely switch to depending on pyproj and remove all previous methods for "unknown" projections or not. Pyproj also doesn't provide any default way of "naming" the grid_mapping variable. I made it default to the same name as the grid_mapping_name attribute.
I think I might make a short term solution for a 0.19.1 release of Satpy to fix #1029 but wanted to make this for now so people could see it and discuss.
flake8 satpy