DAS-2361 - Use np.clip to crop to granule

Author: sudha-murthy

hoss/projection_utilities.py

@@ -203,19 +203,35 @@ def get_projected_x_y_extents(
     grid_lats, grid_lons = get_grid_lat_lons(  # pylint: disable=unpacking-non-sequence
         x_values, y_values, crs
     )
+    max_grid_lon = np.max(grid_lons)
+    min_grid_lon = np.min(grid_lons)
+    max_grid_lat = np.max(grid_lats)
+    min_grid_lat = np.min(grid_lats)
 
     geographic_resolution = get_geographic_resolution(grid_lons, grid_lats)
-
     resolved_points = get_resolved_geojson(
         geographic_resolution, shape_file=shape_file, bounding_box=bounding_box
     )
-    filtered_points = [
-        (
-            max(min(lon, np.nanmax(grid_lons)), np.nanmin(grid_lons)),
-            max(min(lat, np.nanmax(grid_lats)), np.nanmin(grid_lats)),
-        )
-        for lon, lat in resolved_points
-    ]
+    req_lons, req_lats = zip(*resolved_points)
+    min_req_lon = np.min(req_lons)
+    max_req_lon = np.max(req_lons)
+    min_req_lat = np.min(req_lats)
+    max_req_lat = np.max(req_lats)
+    # check if all bbox points are outside the granule.
+    if (
+        min_req_lon > max_grid_lon
+        or max_req_lon < min_grid_lon
+        or min_req_lat > max_grid_lat
+        or max_req_lat < min_grid_lat
+    ):
+        # do not crop. it is outside spatial area
+        filtered_points = resolved_points
+    else:
+        clipped_lons = np.clip(req_lons, min_grid_lon, max_grid_lon)
+        clipped_lats = np.clip(req_lats, min_grid_lat, max_grid_lat)
+        clipped_points = list(zip(clipped_lons, clipped_lats))
+        # print(f'clipped_points={clipped_points}')
+        filtered_points = clipped_points
     return get_x_y_extents_from_geographic_points(filtered_points, crs)
 
 

tests/unit/test_projection_utilities.py

@@ -394,35 +394,60 @@ def test_get_projected_x_y_extents_whole_earth(self):
                 'grid_mapping_name': 'lambert_azimuthal_equal_area',
             }
         )
-        expected_output1 = {
+        expected_output = {
             'x_min': -12702459.818865139,
             'x_max': 12702459.818865139,
             'y_min': -12702440.623773243,
             'y_max': 12702440.710450241,
         }
         with self.subTest('Whole Earth LAEA - Bounding box input'):
-            x_y_extents1 = get_projected_x_y_extents(
-                x_values, y_values, crs, bounding_box=whole_earth_bbox
-            )
-            self.assertDictEqual(x_y_extents1, expected_output1)
             self.assertDictEqual(
-                x_y_extents1,
                 get_projected_x_y_extents(
                     x_values, y_values, crs, bounding_box=whole_earth_bbox
                 ),
-                expected_output1,
+                expected_output,
             )
 
         with self.subTest('Whole Earth LAEA - Shape file input'):
-            x_y_extents2 = get_projected_x_y_extents(
-                x_values, y_values, crs, shape_file=polygon_path
-            )
             self.assertDictEqual(
-                x_y_extents2,
                 get_projected_x_y_extents(
                     x_values, y_values, crs, shape_file=polygon_path
                 ),
-                expected_output1,
+                expected_output,
+            )
+
+    def test_get_projected_x_y_extents_invalid_bbox(self):
+        """Ensure that the an empty extent is returned when the bbox is totally
+        outside the granule extent.
+
+        """
+        bbox_outside_granule = BBox(-180.0, -85.0, 180.0, -75.0)
+        # polygon_path = 'tests/geojson_examples/polygon_whole_earth.geo.json'
+
+        x_values = np.linspace(-8982000, 8982000, 500)
+        y_values = np.linspace(8982000, -8982000, 500)
+
+        crs = CRS.from_cf(
+            {
+                'false_easting': 0.0,
+                'false_northing': 0.0,
+                'longitude_of_central_meridian': 0.0,
+                'latitude_of_projection_origin': 90.0,
+                'grid_mapping_name': 'lambert_azimuthal_equal_area',
+            }
+        )
+        expected_output = {
+            'x_min': -12702459.818865139,
+            'x_max': 12702459.818865139,
+            'y_min': -12702440.623773243,
+            'y_max': 12702440.710450241,
+        }
+        with self.subTest('Spatial area outside granule LAEA - Bounding box input'):
+            self.assertDictEqual(
+                get_projected_x_y_extents(
+                    x_values, y_values, crs, bounding_box=bbox_outside_granule
+                ),
+                expected_output,
             )
 
     def test_get_projected_x_y_variables(self):