Open ghiggi opened 2 years ago
@ghiggi thanks a lot for putting all the effort! Some preliminary remarks:
As such, I strongly suggest you set this PR aside, and start making instead small incremental PRs, adding at most a couple of features/bug fixes at the time, with accompanying tests and documentation. The changes you propose here look very promising, and it would be a pity they don't get through because of having too much to chew at once...
No need to ask for forgiveness for what looks like great work toward solving some long-standing problems, but I second @mraspaud that smaller bites help to avoid choking. Ten small PRs would be easier to handle than one big one. Would it be possible to split it?
If the change to a more shapely-compatible set of names is where we want to go (I'm ok with that) then I'd suggest either:
This way we can put deprecation warnings on all the old stuff and eventually remove it. Heck you could even put it in the future
subpackage of pyresample where all the pyresample 2.0 stuff is supposed to go.
Without reading the actual changes, I worry about testing for this. Some of the tests that were added for the existing interfaces were made because of edge cases with the implementation. Verifying that a new interface doesn't suffer from those or similar issues might be difficult.
In addition to agreeing with the many smaller PRs idea, I'd also suggest thinking about making multiple commits. This entire PR is one commit. That's 2000 lines of code where you never had a breaking/pause point. Separate commits make it easier to review your changes and issues you've run into as you develop something so we can work on it. You could even consider Test Driven Development (TDD) where you right tests first and then write the user code that makes that test pass. It may make it easier to make smaller commits and pull requests this way too since you'd be dividing work into user use cases.
I definitely understand your concerns and I could start to do small incremental PRs instead of this huge one. But it would be almost impossible to perform these PRs to spherical.py
. I modified spherical.py
in this PR with the intention of providing you with an overview of where the main changes happen and discussing your opinions.
I am aware that it exists also the old spherical_geometry.py
file, but if you agree with the need for the refactoring I propose, I think it would be nice if we add a new file called spherical_future.py
(or let me know the name) where I start adding incremental PRs that you can review progressively. I could proceed in the following order:
SPoint
SArc
SLine
SExtent
SPolygon
and SMultiPolygon
Regarding the integration of such classes defined in spherical_future.py
into the current workflow, I would foresee the following:
geometry.py
the SphPolygon
is currently obtained by creating a Boundary
/AreaBoundary
/AreaDefBoundary
object and then calling the contour_poly
method. SwathDefinition
crop/get_area_slices
methods, I essentially come up adding a method boundary
to AreaDefinition
and SwathDefinition
classes, which returns an AreaBoundary
object. I then added the method polygon
to AreaBoundary
which enables the retrieval of the boundary SPolygon
. It will look like area_def.boundary().polygon()
instead of the current AreaDefBoundary(area_def).contour_poly
.spherical_future.py
is ready, we can then progressively think to replace the usage of contours_poly
in the code with polygon
. spherical.py
will remain available for back-compatibility reasons, while the new computations will be performed with the S<geom>
class objects defined in spherical_future.py
.When the above will be in place, I will be able to add the PR to enhance get_area_slices
also for SwathDefinition. I come up with a quite efficient implementation that can return in about 30 seconds the intersection slices of 1 day of LEO orbits with a GEO disk, with a custom specification of the minimum_overlap_fraction
required to define the intersection (see for value 0.9 and 0 in the below pictures) .
Let me know what you think and the path forward !!
Others may disagree with this, but here is my gut reaction in order to:
I think you should:
pyresample/future/spherical/
(https://github.com/pytroll/pyresample/tree/main/pyresample/future). You already have these shapely utils, but are they used by any other module? I'd say put them in pyresample/future/spherical/shapely_utils.py
or something like that. You could go Java-style and put each spherical object in its own spherical/
submodule, but maybe there is a more consolidated structure that would work (ex. primitives.py
for lines, points, and arcs; polygons.py
for extent and polygons).pyresample/test/test_futures/test_spherical.py
(or test_spherical
sub-directory).pyresample/future
module and be in line with existing Pyresample 2.0 design concepts (#181, https://github.com/pytroll/pyresample/milestone/3). Although as I say that I realize none of those issues include the decisions/discussions we had made about geometry objects...oops.Questions:
future
, should we make more drastic changes to the spherical polygon API?We appreciate all you've done and I think these will all be really good improvements for pyresample. Thanks for helping and hopefully we can get these things figured out sooner rather than later.
I agree with all that you suggest !
spherical.py
into futures/spherical.py
. After this first merge, I can start to do incremental PR to this file, so that you can easily track the changes, and you can easily see the point at which it breaks the compatibility. Answering already your question:
area_def.boundary().line()
and area_def.boundary().polygon()
But I think people might be interested in accessing some of the high classes with something like from pyresample.spherical import SPolygon, SLine, SPoint
to perform operations that currently shapely does not enable. As an example, with less than 10 lines of code is possible to simulate/generate SwathDefinition objects using pyorbital
satellite orbits. Let me know @mraspaud @djhoese @gerritholl. Have a nice weekend.
- So, I think that for starting, we can do a copy of current
spherical.py
intofutures/spherical.py
. After this first merge, I can start to do incremental PR to this file, so that you can easily track the changes, and you can easily see the point at which it breaks the compatibility.
No, please do not duplicate code! Make calls to existing functionality from functions with new names, but please do not copy...
If there are bug fixes that can be applied to the existing interfaces :+1:. Make those PRs. I see you already made one for the inplace operations in the intersection code.
Otherwise, for these new interfaces, start a PR where you take what you've done here (which seems almost entirely new code, right?) and put it in pyresample/future/spherical/polygons.py
. Add imports to pyresample/future/spherical/__init__.py
to support the imports you are hoping for (I would have suggested this anyway). Make a commit, push, and start the PR from here. You could even just modify this existing one if you really want. Then add tests with inspiration from the existing tests and any other information you know since you've developed this (what are the edge cases, what new functionality is possible by the user). Make multiple commits during this step and push. Then start refactoring into multiple modules if it makes sense. If you or anyone else doesn't see the benefit of this then we can move it all to pyresample/future/spherical.py
.
Given that this seems like a new and separate interface from the existing functionality (correct me if I'm wrong) then starting from the point of the existing interfaces doesn't have a huge benefit in my eyes. So I say just make the new modules and I'll review it when it is ready.
@ghiggi I talked with Martin today in a video call about this (we had the video call about something else, but then talked about this PR). I thought this PR was more a complete rewrite than it was so I have some suggestions. Let me know what you think.
I think you could make a lot of little PRs out of this PR. I'm not entirely sure as I write this where everything should go. Possible PRs (with tests!):
__init__
could be updated with a DeprecationWarning saying to use one of the new classes.__contains__
/within
/intersection
, one for plotting, etc. If any of this can be done from a base class for all spherical classes that would be great. Like a central def plot
method. Or make plot
a separate utility function outside of the classes (Martin wanted this, I'm less sure).Some quick considerations/thoughts:
SCoordinate
and CCoordinate
are essentially unvaried except for the plot method and some improved docs.Arc
, I can add the methods to spherical.py
without breaking back-compatibility.SPolygon
is almost a complete rewrite and I will put it in /future/spherical/polygon.py
. Really no sense to add PRs to SphPolygon
. I will do PR "per topic". SLine
and SExtent
are completely new and I will place them into /future/spherical/line.py
and /future/spherical/extent.py
.SArc
and SPoint
will be done in /future/__init__
by importing spherical.Arc
and SCoordinate
.SPoint
on Cartopy is necessary to use ax.scatter
, because shapely Point/MultiPoint
can not be plotted AFAIK. Plotting LineString
requires the specification of facecolor='none'
otherwise it plots polygons ... we will see ... In my mind I plan to commit 1/2 day per week on these PRs, to give you time to review.
I will also start parallelly to add incremental small PR for the Boundary
retrieval and SwathDefinition methods. I am not sure I see how to add new methods in /future/geometry/swath
(monkey patching?). Maybe I will need an example ;)
At a given point, we should discuss also how to implement a function to test if a point is contained in a polygon. This will be required to enable the API to work on all edge cases. I couldn't figure out which maths are required ...
That sounds like it should work. For the future SwathDefinition, we'd have to make a subclass but I'm not sure that is even the right approach yet. As long as it doesn't break backwards compatibility, could you define new methods or change existing methods in the current SwathDefinition that use the future spherical stuff?
@mraspaud @djhoese. It's a long time since the last iterations on this refactoring.
Unfortunately, I never found the time to advance and implement the tests ...
I was thinking ... if I come to the PCW meeting in November in Athens, do you think it could make sense to address (part of) this this refactoring during the week? I think that in-person discussion will facilitate the design and implementation of the required PRs. Let me know ;-)
It will definitely help, yes! Live discussions, but also dedicated time will allow you to work on this efficiently.
In such a case, would you agree to have a meeting already in mid of October to briefly discuss again the structure of the PRs? In this way, I could already do some preparatory PRs to be efficient as possible during the week in Athens.
Works for me!
Hey @mraspaud. Do you have time in the week between the 17 and 21 October? Do @djhoese want to participate too? Have a nice week :)
I'd be OK joining any conversations, but given the time zone differences don't let my schedule stop you from meeting.
Edit: It may be best to have this conversation on slack.
This PR provides a refactoring of
spherical.py
, introduces new methods and solves long-standing issues. This PR is required for a future PR addressing the missing methods/option for cropping/slicing ofSwathDefinition
andAreaDefinition
. It is also likely thatpyorbital
andpytrollschedule
can benefit of this PR.The refactoring of
spherical.py
tries to stick to the naming/conventions employed by theshapely
library, which is commonly used to perform planar geometrical computation. The PR will introduce backward incompatibilities![ ] Closes
https://github.com/pytroll/satpy/issues/383
https://github.com/pytroll/pyresample/issues/329
[ ] Step forward toward solving
[ ] Tests added
[ ] Tests passed
[ ] Passes
git diff origin/main **/*py | flake8 --diff
[ ] Fully documented
Here below I summarize briefly the content of this very long (sorry!) PR. It implement 6 main classes:
SExtent
[NEW]SPoint
[previously called SCoordinate]SArc
[previously called SArc]SLine
[NEW]SPolygon
[previously called SphPolygon]SMultiPolygon
[NEW]SPoint
plot
to plot with Cartopy.SArc
midpoint
,segmentize
,to_line
,to_shapely
andplot
.intersections
has been renamed_great_circle_intersections
.intersection
has been renamedintersection_point
._hash__
has been added to possibly LRU cacheintersection_point
._great_circle_intersections
/intersections
which was causing https://github.com/pytroll/satpy/issues/383 has been resolved.SLine
intersects
,intersection_points
,segmentize
,buffer
,get_parallel_lines
,to_shapely
andplot
.SPolygon/SMultiPolygon
within
,contains
,equals
,equals_exact
,intersects
,disjoint
,overlap_fraction
,overlap_rate
,intersection_points
,buffer
,segmentize
,to_line
,to_shapely
andplot
.SphPolygon
, it does not accept aradius
argument and performs all computations on a unit sphere.area
can be calculated more accurately withpyproj
specifying a customellips
oid. A method for computing theperimeter
is also provided.is_valid_geometry
property.intersection
, it now does not fail when the intersection with itself, providesall
the intersection also in cases of multiple intersections (which returns a SMultiPolygon), and deal correctly with the antimeridian crossing cases.union
, if theSPolygon
(s)/SMultiPolygon
(s) are not intersecting it now returns the original geometries in a singleSMultiPolygon
(instead ofNone
). Cascaded union is performed in case of unioning multiple SPolygons.SMultiPolygons
contain multipleSPolygon
(s) with the constraint that they must not overlap each other!SExtent
SPolygon
taking into account if it crosses the antimeridian. As an example, the GOES-17 full disk would be represented by a list of 2 extent:[[-180.0, -57.76, -79.14, 79.14], [143.76, 180.0, -74.79, 74.79]]
so that the extents does not cross the antimeridian and shapely can be employed for the computations.SExtent
can be composed of an unbounded number of extents, under the conditions that the extents do not overlap each other. All the extents can be plotted usingplot
.intersects
,disjoint
,within
,contains
,touches
andequals
employs shapely in the background for fast planar rectangle geometry computation. These methods enable to speed up some of theSPolygon
spherical computations and the upcoming PR enablingSwathDefinition.crop
,SwathDefinition.get_area_slices
andSwathDefinition.get_extent_slices
.I didn't implement yet carefully tests as I wait for your opinion before proceeding. But I used the code to run intersection/union and get MODIS and GPM swath slices with GEO AreaDefinition over 1 year of data and everything looks fine ;) As a note of caution, in case of polygons encloses the poles, the
SExtent
of the SPolygon is currently wrongly calculated as I didn't manage to implement a function that checks if a point (i.e. the pole) is contained in the polygon. I AM LOOKING FOR YOUR HELP WITH THIS.I am sorry for this long PR landing on your head, I hope you can forgive me @mraspaud @djhoese @gerritholl.