"""Tests :class:`pyrate.common.raster_datasets.geo_datasets.DataSetAccess`.""" # Standard library from math import degrees from math import pi from math import radians # Generic testing from unittest import TestCase # Geometry from rasterio.windows import intersect # Hypothesis testing from hypothesis import given import hypothesis.strategies as st # Numeric testing from numpy import rad2deg from numpy.testing import assert_array_almost_equal # Own Geometry from pyrate.plan.geometry.helpers import meters2rad from pyrate.plan.geometry.helpers import rad2meters from pyrate.plan.geometry import PolarLocation # Test environment helper from ... import _open_test_geo_dataset class TestGeoDataset(TestCase): """Ensure that the :class:`pyrate.plan.graph.generate.geo_datasets.DataSetAccess` works correctly. Uses the *Earth2014* dataset. """ #: The resolution of the dataset in arc-minutes DATASET_RESOLUTION = 30 #: The maximal distance of two data points in the dataset in degrees MAX_POINT_DISTANCE_DEG = DATASET_RESOLUTION / 60 #: The maximal distance of two data points in the dataset in meters (at the equator) MAX_POINT_DISTANCE = rad2meters(radians(MAX_POINT_DISTANCE_DEG)) # Handle the context manager, see https://stackoverflow.com/a/11180583/3753684 def run(self, result=None) -> None: with _open_test_geo_dataset() as dataset: self.dataset = dataset # pylint: disable=attribute-defined-outside-init super().run(result) @given( st.floats(min_value=-pi / 2 * 0.75, max_value=+pi / 2 * 0.75), st.floats(min_value=-pi * 0.75, max_value=+pi * 0.75), st.floats(min_value=0.001, max_value=1000_000.0), ) def test_bounding_window_center_of_dataset( self, latitude: float, longitude: float, radius: float ) -> None: """Tests that the bounding box is correctly calculated if the query point is not at the border.""" # pylint: disable=too-many-locals win_1, win_2 = self.dataset.get_bounding_windows_around(latitude, longitude, radius) self.assertIsNone(win_2, "only one window shall be returned") # Get the geographical extends in degrees, not radians: left, bottom, right, top = self.dataset.dataset.window_bounds(win_1) # Check the position of the window latitude_deg, longitude_deg = degrees(latitude), degrees(longitude) self.assertLessEqual(bottom, latitude_deg + 1e-12) self.assertGreaterEqual(top, latitude_deg - 1e-12) self.assertAlmostEqual( (top + bottom) / 2, latitude_deg, delta=self.MAX_POINT_DISTANCE_DEG, msg="window should be vertically centered around the given center point", ) self.assertLessEqual(left, longitude_deg + 1e-12) self.assertGreaterEqual(right, longitude_deg - 1e-12) self.assertAlmostEqual( (right + left) / 2, longitude_deg, delta=self.MAX_POINT_DISTANCE_DEG, msg="window should be horizontally centered around the given center point", ) # Check the size of the window self.assertLessEqual(left, right) self.assertLessEqual(bottom, top) # Distances are uniform along longitudes radius_deg = degrees(meters2rad(radius)) self.assertGreaterEqual(top - bottom, 2 * radius_deg - 1e-12) self.assertLessEqual(top - bottom, 2 * radius_deg + 3 * self.MAX_POINT_DISTANCE_DEG) # Distances are more complicated along the latitudes left_side_center = PolarLocation(latitude=(top + bottom) / 2, longitude=left) right_side_center = PolarLocation(latitude=(top + bottom) / 2, longitude=right) # Use approximate=True for a spherical model distance_horizontal = left_side_center.distance(right_side_center, approximate=True) # The rough checking of the bounding boxes is very coarse and was determined by fiddling until it # works # This part of the test guarantees that the windows is roughly the right size, but pinning it down to # exact number is hard since we round the discrete window bounds, map them to geographical coordinates # and then have to deal with floating point inaccuracies self.assertGreaterEqual(distance_horizontal, 2 * radius - 6 * self.MAX_POINT_DISTANCE) self.assertLessEqual(distance_horizontal, 2 * radius + 4 * self.MAX_POINT_DISTANCE) @given( st.floats(min_value=-pi / 2 * 0.95, max_value=+pi / 2 * 0.95), st.one_of( [st.floats(min_value=-pi, max_value=-pi * 0.995), st.floats(min_value=+pi * 0.995, max_value=+pi)] ), st.floats(min_value=200_000.0, max_value=1_000_000.0), ) def test_bounding_window_left_and_right_side( self, latitude: float, longitude: float, radius: float ) -> None: """Tests that the bounding box is correctly calculated if the query point is at the border. Very high latitudes (near the poles) are not tested as there might be a single (albeit very wide) window being returned. For the same reason, only moderate radii are tested. """ window_1, window_2 = self.dataset.get_bounding_windows_around(latitude, longitude, radius) # We need a plain assert here for type checking assert window_2 is not None, "two windows should be returned" self.assertFalse(intersect(window_1, window_2), "windows may never overlap") # Also test the same as in :meth:`~test_bounding_window_intersection_empty` self.assertEqual(window_1.height, window_2.height) self.assertGreaterEqual(window_1.height, 1) self.assertGreaterEqual(window_1.width, 1) self.assertGreaterEqual(window_2.width, 1) self.assertLessEqual(window_1.height, self.dataset.dataset.height) self.assertLessEqual(window_1.width + window_2.width, self.dataset.dataset.width) @given( st.floats(min_value=-pi / 2, max_value=+pi / 2), st.floats(min_value=-pi, max_value=+pi), st.floats(min_value=0.001, max_value=100_000_000.0), ) def test_bounding_window_general_properties( self, latitude: float, longitude: float, radius: float ) -> None: """Tests some more general properties that should hold for all windows and window pairs. In particular, it makes sure that even if a window pair is very wide, the intersection is always empty. """ window_1, window_2 = self.dataset.get_bounding_windows_around(latitude, longitude, radius) # Make sure that everything in window_1 is rounded self.assertIsInstance(window_1.col_off, int) self.assertIsInstance(window_1.row_off, int) self.assertIsInstance(window_1.height, int) self.assertIsInstance(window_1.width, int) # Test some general properties of window_1 self.assertTrue(radius == 0 or window_1.height >= 1) self.assertLessEqual(window_1.height, self.dataset.dataset.height) self.assertTrue(radius == 0 or window_1.width >= 1) self.assertLessEqual(window_1.width, self.dataset.dataset.width) if window_2 is not None: self.assertGreater(radius, 0) # Make sure that everything in window_2 is rounded self.assertIsInstance(window_2.col_off, int) self.assertIsInstance(window_2.row_off, int) self.assertIsInstance(window_2.height, int) self.assertIsInstance(window_2.width, int) # Test some general properties of window_2 in relation to window_1 self.assertEqual(window_1.height, window_2.height) self.assertGreaterEqual(window_2.width, 1) self.assertLessEqual(window_1.width + window_2.width, self.dataset.dataset.width) self.assertFalse(intersect(window_1, window_2), "windows may never overlap") @given( st.floats(min_value=-pi / 2 * 0.75, max_value=+pi / 2 * 0.75), st.floats(min_value=-pi * 0.75, max_value=+pi * 0.75), st.floats(min_value=0.001, max_value=1000_000.0), ) def test_meshgrid_generation(self, latitude: float, longitude: float, radius: float) -> None: """Tests that msehgrids are generated correctly no matter whether radians are used or not. Uses the data generation of :meth:`~test_bounding_window_center_of_dataset`. """ window, window_empty = self.dataset.get_bounding_windows_around(latitude, longitude, radius) self.assertIsNone(window_empty, "only one window shall be returned") lats_deg, lons_deg = self.dataset.lat_lon_meshgrid_for(window, window_empty, radians=False) lats_rad, lons_rad = self.dataset.lat_lon_meshgrid_for(window, window_empty, radians=True) assert_array_almost_equal(rad2deg(lats_rad), lats_deg) assert_array_almost_equal(rad2deg(lons_rad), lons_deg)