Geos area reduction #115
Geos area reduction #115
Conversation
pyresample/geometry.py
Outdated
| @@ -1369,6 +1370,127 @@ def proj4_string(self): | |||
| items = self.proj_dict.items() | |||
| return '+' + ' +'.join([t[0] + '=' + str(t[1]) for t in items]) | |||
|
|
|||
| def get_area_slices(self, area_to_cover): | |||
| """Compute the slice to read from an *area* based on an *area_to_cover*.""" | |||
There was a problem hiding this comment.
E501 line too long (83 > 79 characters)
pyresample/geometry.py
Outdated
| # Intersection only required for two different projections | ||
| if area_to_cover.proj_dict['proj'] == self.proj_dict['proj']: | ||
| logger.debug('Projections for data and slice areas are' | ||
| ' identical: {}'.format(area_to_cover.proj_dict['proj'])) |
There was a problem hiding this comment.
E501 line too long (82 > 79 characters)
|
Couldn't we use pyproj to transform X/Y coordinates between the source area and the crop area? Why do the 'geos' angles matter? Shouldn't we be able to use lon/lat as intermediate coordinates (or do the transform directly using pyproj's functions for doing that)? Lastly, |
|
I wrote this a long time ago, but iirc it was because I wanted to go directly from the instrument scanning angle to the lon/lat coordinates. But you're right, we could use proj and just divide the x and y projection coordinates by the altitude of the satellite to get the angle. About |
|
I was thinking of not using satellite height at all. We should be able to use transform (https://jswhit.github.io/pyproj/pyproj-module.html#transform). This would only require the source PROJ.4, the target PROJ.4, and the extents of the two areas. We could convert from the crop area extents to the source area extents and use those to figure out the indexes of the area to crop. If any of the crop area extents can't be transformed (1e30 iirc) then we use the source area extent value for that coordinate. I thought maybe |
|
Summary of conversation with @mraspaud on slack: The boundary objects already handle everything that would need to be handled if we used the X/Y PROJ.4 transforms. They also use spherical coordinates (I think) so it should work in all cases. |
|
|
||
| from trollsched.boundary import AreaDefBoundary, Boundary | ||
|
|
||
| data_boundary = Boundary(*get_geostationary_bounding_box(self)) |
There was a problem hiding this comment.
How hard would it be to make this work for other projections? Couldn't you just use AreaDefBoundary if not geos?
There was a problem hiding this comment.
I'm afraid that every projection has it's own boundary to undefined/space pixels, eg in mollweide, polar stereographic. And to make this efficient, this boundary should be computed, vs eg search iteratively.
There was a problem hiding this comment.
Then the use of AreaDefBoundary a couple lines below this shouldn't be allowed by that same reasoning, right?
There was a problem hiding this comment.
Additionally, if the ROI is geos then shouldn't that also use the above logic Boundary(*get_geostation_bounding_box(...))?
There was a problem hiding this comment.
Yes, this is true theoretically. However, I believe that users will mostly provide a well define ROI
Implement area-slicing for geos-based areas.
git diff origin/develop **/*py | flake8 --diff