Fix geostationary bbox having inf values #493
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -36,10 +36,9 @@ | |
| from pyresample._spatial_mp import Cartesian, Cartesian_MP, Proj_MP | ||
| from pyresample.area_config import create_area_def | ||
| from pyresample.boundary import AreaDefBoundary, Boundary, SimpleBoundary | ||
| from pyresample.utils import check_slice_orientation, load_cf_area | ||
| from pyresample.utils.proj4 import ( | ||
| get_geostationary_height, | ||
| proj4_dict_to_str, | ||
| proj4_radius_parameters, | ||
| ) | ||
|
|
@@ -1556,7 +1555,7 @@ def _get_geo_boundary_sides(self, frequency=None): | |
| # Ensure an even number of vertices for side creation | ||
| if (frequency % 2) != 0: | ||
| frequency = frequency + 1 | ||
| lons, lats = get_geostationary_bounding_box_in_lonlats(self, nb_points=frequency) | ||
| # Retrieve dummy sides for GEO (side1 and side3 always of length 2) | ||
| side02_step = int(frequency / 2) - 1 | ||
| lon_sides = [lons[slice(0, side02_step + 1)], | ||
|
|
@@ -2772,40 +2771,49 @@ def get_geostationary_angle_extent(geos_area): | |
| aeq = 1 - req ** 2 / (h ** 2) | ||
| ap_ = 1 - rp ** 2 / (h ** 2) | ||
|
|
||
| x_angle = np.arccos(np.sqrt(aeq)) | ||
| y_angle = np.arccos(np.sqrt(ap_)) | ||
| return x_angle, y_angle | ||
|
|
||
|
|
||
| def get_geostationary_bounding_box_in_proj_coords(geos_area, nb_points=50): | ||
| """Get the bbox in geos projection coordinates of the valid pixels inside `geos_area`. | ||
|
|
||
| Parameters | ||
| ---------- | ||
| nb_points : Number of points on the polygon. | ||
|
|
||
| """ | ||
| x, y = get_full_geostationary_bounding_box_in_proj_coords(geos_area, nb_points) | ||
| ll_x, ll_y, ur_x, ur_y = geos_area.area_extent | ||
|
|
||
| from shapely.geometry import Polygon | ||
| geo_bbox = Polygon(np.vstack((x, y)).T) | ||
| area_bbox = Polygon(((ll_x, ll_y), (ll_x, ur_y), (ur_x, ur_y), (ur_x, ll_y))) | ||
| intersection = area_bbox.intersection(geo_bbox) | ||
|
Comment on lines
+2790
to
+2793
There was a problem hiding this comment. There is a very good chance I'm missing something, but could this be simplified (logically, maybe not in code) to min/max each of the coordinates in the full geostationary bounding box with the extents of the area? So if the full bbox says a point is at Not sure this is accurate, or makes sense, or is a good idea, but it is the first thing that comes to mind. There was a problem hiding this comment. That's basically what we had before, in some way. The problem is to get the corner points of the intersection between the disc and the truncation, especially if only one corner of the truncated area is in space. There was a problem hiding this comment. Ok so you're saying if the disk polygon was computed and has points partially in space that clipping to the extents of the area wouldn't effect some of those points? How does this current intersection do that? I thought the problem in In your comment above you said:
When you say "cut", do you mean clipping the disk coordinates to the extent of the area or do you mean how the coordinates were calculated for half the disk and then replicated for the other half? There was a problem hiding this comment. ok, so the polygon for the full disk looks like this: Now, if we just clip y at 2e6 like it is in main, we get this: If you look carefully, you see the last point on each side of the y=2e6 line ends up in space, because we 'project' the extreme points of the equator onto a latitude where the sphere is not as wide. So how do we avoid this? There was a problem hiding this comment. Ok, I think I got it. One more alternative solution: mask out (remove) any points not within the bbox. The main point here is that our extent bbox is not a complicated polygon so we only have two limits to check: disk_points = ...
area_points = disk_points[
(disk_points[:, 0] >= min(extent[0], extent[2])) &
(disk_points[:, 0] <= max(extent[0], extent[2])) &
(disk_points[:, 1] >= min(extent[1], extent[3])) &
(disk_points[:, 1] <= min(extent[1], extent[3]))]There was a problem hiding this comment. The only other thing, which is at this point getting ridiculous and we should just use shapely, would be to combine my two suggestions and every outer disk pixel gets clipped to the min/max of the inner disk pixels (the pixels inside the area's bbox). There was a problem hiding this comment. Nevermind again, this last suggestion doesn't fix the issue. You still end up with the points being no closer to the extent of the area. Ok, shapely it is. Maybe at this point we say shapely is required? What is it being used for now? Did @ghiggi want to depend on it more heavily in the future for spherical stuff? There was a problem hiding this comment. Yes, I plan to use it for several tasks. Could make sense to have shapely as a requirement. There was a problem hiding this comment.
Gradient search resampler is using There was a problem hiding this comment. Made shapely a hard requirement |
||
| try: | ||
| x, y = intersection.boundary.xy | ||
| except NotImplementedError: | ||
| return [], [] | ||
| return np.asanyarray(x[:-1]), np.asanyarray(y[:-1]) | ||
|
|
||
|
|
||
| def get_full_geostationary_bounding_box_in_proj_coords(geos_area, nb_points=50): | ||
| """Get the bbox in geos projection coordinates of the full disk in `geos_area` projection. | ||
|
|
||
| Args: | ||
| nb_points: Number of points on the polygon | ||
| """ | ||
| x_max_angle, y_max_angle = get_geostationary_angle_extent(geos_area) | ||
| h = get_geostationary_height(geos_area.crs) | ||
|
|
||
| # generate points around the north hemisphere in satellite projection | ||
| # make it a bit smaller so that we stay inside the valid area | ||
| points_around = np.linspace(-np.pi, np.pi, nb_points, endpoint=False) | ||
| x = np.cos(points_around) * (x_max_angle - 0.0001) | ||
| y = -np.sin(points_around) * (y_max_angle - 0.0001) | ||
| x *= h | ||
| y *= h | ||
|
|
||
| return x, y | ||
|
|
||
|
|
||
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.
Is it a problem for any of the other interfaces (the bounds and polygon stuff) that the end points are not the same?
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.
The tests didn't break, so I think we're good. Also I always thought that was the way we did it in pyresample.