diff --git a/ioos_qc/qartod.py b/ioos_qc/qartod.py index bf6b455..8278be9 100644 --- a/ioos_qc/qartod.py +++ b/ioos_qc/qartod.py @@ -15,10 +15,10 @@ isfixedlength, add_flag_metadata, great_circle_distance, - mapdates + mapdates, + distance_from_target ) - L = logging.getLogger(__name__) # noqa @@ -85,10 +85,13 @@ def qartod_compare(vectors : Sequence[Sequence[N]] @add_flag_metadata(standard_name='location_test_quality_flag', long_name='Location Test Quality Flag') -def location_test(lon : Sequence[N], - lat : Sequence[N], - bbox : Tuple[N, N, N, N] = (-180, -90, 180, 90), - range_max : N = None +def location_test(lon: Sequence[N], + lat: Sequence[N], + bbox: Tuple[N, N, N, N] = (-180, -90, 180, 90), + range_max: N = None, + target_lat: [N] = None, + target_lon: [N] = None, + target_range: N = None ) -> np.ma.core.MaskedArray: """Checks that a location is within reasonable bounds. @@ -103,6 +106,11 @@ def location_test(lon : Sequence[N], lat: Latitudes as a numeric numpy array or a list of numbers. bbox: A length 4 tuple expressed in (minx, miny, maxx, maxy) [optional]. range_max: Maximum allowed range expressed in geodesic curve distance (meters). + target_lat: Target Latitude as numeric numpy array or a list of numbers, + it can either same size as lat/lon or a unique values + target_lon: Target Longitude as numeric numpy array or a list of numbers, + it can either same size as lat/lon or a unique values + target_range: Maximum allowed range in geodesic curve distance (meters) away from target position. Returns: A masked array of flag values equal in size to that of the input. @@ -130,6 +138,34 @@ def location_test(lon : Sequence[N], lon = lon.flatten() lat = lat.flatten() + # Handle target inputs + # If any target inputs are provided + if target_lon is not None or target_lat is not None or target_range is not None: + # All target inputs should be there + if target_lon is not None and target_lat is not None and target_range is not None: + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + target_lat = np.ma.masked_invalid(np.array(target_lat).astype(np.float64)) + target_lon = np.ma.masked_invalid(np.array(target_lon).astype(np.float64)) + if type(target_range) not in [int, float]: + raise ValueError('Bad target_range input. target_range should either float or int.') + + elif target_lon is not None and target_lat is not None and target_range is None: + raise ValueError('Missing target_range input if target_lat and target_lon are provided') + else: + raise ValueError('Missing some target inputs') + + if target_lon.shape != target_lat.shape: + raise ValueError( + 'Target_lon ({0.shape}) and target_lat ({1.shape}) are different shapes'.format( + target_lon, target_lat + ) + ) + + # Flatten target_lon and target_lat + target_lon = target_lon.flatten() + target_lat = target_lat.flatten() + # Start with everything as passing (1) flag_arr = np.ma.ones(lon.size, dtype='uint8') @@ -147,7 +183,24 @@ def location_test(lon : Sequence[N], d = great_circle_distance(lat, lon) flag_arr[d > range_max] = QartodFlags.SUSPECT - # Ignore warnings when comparing NaN values even though they are masked + # Distance From Target Test + if target_lat is not None and target_lon is not None and \ + target_range is not None: + # If only one value is given assume to be constant for all positions + if target_lon.size == 1 and target_lat.size == 1: + (target_lon, target_lat) = (target_lon * np.ones(lat.size), target_lat * np.ones(lat.size)) + + # Compute the range from the target location + d_from_target = distance_from_target(lat, lon, target_lat, target_lon) + + # Flag as suspect distances greater than target_range + with np.errstate(invalid='ignore'): + flag_arr[d_from_target > target_range] = QartodFlags.SUSPECT + + # Flag as missing target_location distance + flag_arr[(target_lat.mask | target_lon.mask)] = QartodFlags.MISSING + + # Ignore warnings when comparing NaN values even though they are masked # https://github.com/numpy/numpy/blob/master/doc/release/1.8.0-notes.rst#runtime-warnings-when-comparing-nan-numbers with np.errstate(invalid='ignore'): flag_arr[(lon < bbox.minx) | (lat < bbox.miny) | diff --git a/ioos_qc/utils.py b/ioos_qc/utils.py index 417ef9e..b8d0ee3 100644 --- a/ioos_qc/utils.py +++ b/ioos_qc/utils.py @@ -263,3 +263,10 @@ def great_circle_distance(lat_arr, lon_arr): g = Geod(ellps='WGS84') _, _, dist[1:] = g.inv(lon_arr[:-1], lat_arr[:-1], lon_arr[1:], lat_arr[1:]) return dist + + +def distance_from_target(lat, lon, target_lat, target_lon): + g = Geod(ellps='WGS84') + _, _, dist_to_target = g.inv(lon, lat, target_lon, target_lat) + dist_to_target = np.ma.masked_invalid(dist_to_target.astype(np.float64)) + return dist_to_target diff --git a/tests/test_qartod.py b/tests/test_qartod.py index 8c89b22..3273c68 100644 --- a/tests/test_qartod.py +++ b/tests/test_qartod.py @@ -117,6 +117,18 @@ def test_location_bad_input(self): with self.assertRaises(ValueError): qartod.location_test(lon=70, lat=70, bbox=(1, 2)) + # Wrong target lon + with self.assertRaises(ValueError): + qartod.location_test(lon=70, lat=70, bbox=(1, 2, 3, 4), target_lon='foo', target_lat=70, target_range=3000) + + # Wrong target lat + with self.assertRaises(ValueError): + qartod.location_test(lon=70, lat=70, bbox=(1, 2, 3, 4), target_lon=70, target_lat='bad', target_range=3000) + + # Wrong target range + with self.assertRaises(ValueError): + qartod.location_test(lon=70, lat=70, bbox=(1, 2, 3, 4), target_lon=70, target_lat=70, target_range='300') + def test_location_bbox(self): lon = [80, -78, -71, -79, 500] lat = [None, 50, 59, 10, -60] @@ -155,6 +167,74 @@ def test_location_distance_threshold(self): np.ma.array([1, 1, 3]) ) + def test_location_single_target_threshold(self): + lon = np.array([-71.05, -71.06, -80.0]) + lat = np.array([41.0, 41.02, 45.05]) + + npt.assert_array_equal( + qartod.location_test(lon, lat, target_range=3000.0, target_lon=-71.06, target_lat=41), + np.ma.array([1, 1, 3]) + ) + + def test_location_multiple_target_threshold(self): + lon = np.array([-71.05, -71.06, -80.0]) + lat = np.array([41.0, 41.02, 45.05]) + target_lon = np.array([-71.05, -75.06, -80.0]) + target_lat = np.array([41.0, 41.02, 45.05]) + + npt.assert_array_equal( + qartod.location_test(lon, lat, target_range=3000.0, target_lon=target_lon, target_lat=target_lat), + np.ma.array([1, 3, 1]) + ) + + def test_location_multiple_target_missing_target(self): + lon = np.array([-71.05, -71.06, -80.0, -80.0]) + lat = np.array([41.0, 41.02, 45.05, 45.05]) + target_lon = np.array([-71.05, None, -80.0, None]) + target_lat = np.array([41.0, 41.02, None, None]) + + npt.assert_array_equal( + qartod.location_test(lon, lat, target_range=3000.0, target_lon=target_lon, target_lat=target_lat), + np.ma.array([1, 9, 9, 9]) + ) + + def test_location_target_missing_threshold(self): + lon = np.array([-71.05, -71.06, -80.0]) + lat = np.array([41.0, 41.02, 45.05]) + target_lon = np.array([-71.05, -75.06, -80.0]) + target_lat = np.array([41.0, 41.02, 45.05]) + + with self.assertRaises(ValueError): + qartod.location_test(lon, lat, target_range=None, target_lon=target_lon, target_lat=target_lat) + with self.assertRaises(ValueError): + qartod.location_test(lon, lat, target_range=None, target_lon=target_lon[0], target_lat=target_lat[0]) + + def test_location_target_missing_threshold(self): + lon = np.array([-71.05, -71.06, -80.0]) + lat = np.array([41.0, 41.02, 45.05]) + target_lon = np.array([-71.05, -75.06, -80.0]) + target_lat = np.array([41.0, 41.02, 45.05]) + + with self.assertRaises(ValueError): + qartod.location_test(lon, lat, target_range=3000, target_lat=target_lat) + with self.assertRaises(ValueError): + qartod.location_test(lon, lat, target_range=None, target_lon=target_lon) + with self.assertRaises(ValueError): + qartod.location_test(lon, lat, target_range=3000, target_lat=45) + with self.assertRaises(ValueError): + qartod.location_test(lon, lat, target_range=3000, target_lon=45) + + def test_location_dual_threshold_test(self): + lon = np.array([-71.05, -71.06, -80.0, -80.0]) + lat = np.array([41.0, 41.02, 45.05, 45.05]) + target_lon = np.array([-71.05, -71.06, -80.0, -80.0]) + target_lat = np.array([41.0, 41.02, 45.05, 46.05]) + + npt.assert_array_equal( + qartod.location_test(lon, lat, range_max=3000, + target_range=3000.0, target_lon=target_lon, target_lat=target_lat), + np.ma.array([1, 1, 3, 3]) + ) class QartodGrossRangeTest(unittest.TestCase):