Skip to content
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

Fix spill over of ocean/and features in cartopy plots in case of geostationary full disc plot. #601

Merged
merged 6 commits into from
Jul 30, 2024

Conversation

BENR0
Copy link
Contributor

@BENR0 BENR0 commented Jun 14, 2024

In case of geostationary full disc plots the boundary of the cartopy CRS needs to be the actual full disc otherwise ocean/land features spill over.

Since the to_cartopy_crs method in pyresample uses a custom Projection class this is a hot fix using code from the Cartopy Geostationary class.

  • Tests added
  • Tests passed
  • Fully documented

Copy link

codecov bot commented Jun 14, 2024

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 93.97%. Comparing base (56f4506) to head (9bf142c).
Report is 116 commits behind head on main.

Additional details and impacted files
@@            Coverage Diff             @@
##             main     #601      +/-   ##
==========================================
- Coverage   94.01%   93.97%   -0.04%     
==========================================
  Files          92       86       -6     
  Lines       13836    13811      -25     
==========================================
- Hits        13008    12979      -29     
- Misses        828      832       +4     
Flag Coverage Δ
unittests 93.97% <100.00%> (-0.04%) ⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

@coveralls
Copy link

Coverage Status

coverage: 93.685% (+0.009%) from 93.676%
when pulling e16203a on BENR0:feat_html_repr
into db94a87 on pytroll:main.


# For geostationary full disc the boundary needs to be the actuall boundary of the earth disc otherwise
# when ocean or land features are added to the cartopy plot these spill over.
if "Geostationary Satellite" in crs.to_wkt():
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would be great if this if block for geos was put into a separate function/method. Or at least invert the condition and put the much smaller else: block in the main if block and return immediately.

@@ -69,3 +70,70 @@ def __init__(self,
if self.bounds is not None:
x0, x1, y0, y1 = self.bounds
self.threshold = min(abs(x1 - x0), abs(y1 - y0)) / 100.

# For geostationary full disc the boundary needs to be the actuall boundary of the earth disc otherwise
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What if it's only part of the disc, but including some boundary?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this gets into the difference between how I thought and wanted the cartopy projection classes to work (as an AreaDefinition clone) versus how they are intended to work. In all my usage and examples of this class I use it to hold the projection information but also the literal bounds of the data. Cartopy uses bounds to mean the actual physical (?) bounds of the projection space, not of the data. Cartopy/matplotlib use xlims/ylims for data bounds.

So I wonder if we could take the produced geostationary bound here and clip it to the x and y limits passed by the user. I'm not sure that is always accurate though.

Or...we completely rewrite these interfaces and/or examples so that they use upstream cartopy projection objects and includes an extra step to set the xlim/ylims of the plot to match the data.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is exactly the issue. And while I think it makes sense to rethink this for the display of the html representation of the AreaDefinition this currently seems to work.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok @BENR0, but doesn't the current method (_geos_boundary) get the full disk geostationary boundary? So if I try to do the HTML representations of a subset of the full disk (ex. ABI CONUS or mesoscale sector) then my representation would be much bigger than intended/expected?

@coveralls
Copy link

coveralls commented Jul 17, 2024

Coverage Status

coverage: 93.677% (+0.001%) from 93.676%
when pulling 9bf142c on BENR0:feat_html_repr
into db94a87 on pytroll:main.

Comment on lines +143 to +146
@property
def boundary(self):
"""Return boundary."""
return self._boundary
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is this needed?

Copy link
Member

@djhoese djhoese left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm OK merging this if it improves what you're seeing and trying to fix. I still have a question on non-full disk geostationary images and what they look like with (or without) this PR. I made one other comment about if the boundary property is strictly needed, but that's not a big deal.

@mraspaud if you wake up and think this looks fine as is then feel free to merge and make a release. If this is not ready then I'd be fine with a release of main as is and save this for a bug fix release shortly after.

@BENR0
Copy link
Contributor Author

BENR0 commented Jul 30, 2024

Yes this should fix the problem with the geostationary plots. I am on vacation this week and can not make plots and post them here. I would also be fine with merging this later if you want to make a release this week.

@djhoese djhoese added the bug label Jul 30, 2024
@djhoese
Copy link
Member

djhoese commented Jul 30, 2024

Ok I've made a couple changes and tested things in a jupyter notebook. It looks like ABI CONUS shows as expected which was my main concern.

First I added a context manager to remove the PROJ warnings when getting the PROJ dict. The context manager was created by me for cases exactly like this where we know that PROJ.4 is inexact, but it is the representation we want to use/display.

The other fix was converting extent floats to float as the ones I was using from the Satpy abi_l1b reader were np.float64 so they would show up as np.float64(123456.78) which was just ugly and unnecessary especially since you were already rounding them.

@djhoese djhoese merged commit 70b37ab into pytroll:main Jul 30, 2024
24 of 26 checks passed
@djhoese
Copy link
Member

djhoese commented Jul 30, 2024

I'll try to make a release including this PR later tonight. Have to go for now.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants