From 370467917a2fc373721614b56ebb463433d84a51 Mon Sep 17 00:00:00 2001 From: andream Date: Tue, 21 May 2024 18:00:58 +0200 Subject: [PATCH 1/6] add *= -1 in test --- satpy/tests/reader_tests/test_li_l2_nc.py | 1 + 1 file changed, 1 insertion(+) diff --git a/satpy/tests/reader_tests/test_li_l2_nc.py b/satpy/tests/reader_tests/test_li_l2_nc.py index 5e9d0ff563..c5e02f93d0 100644 --- a/satpy/tests/reader_tests/test_li_l2_nc.py +++ b/satpy/tests/reader_tests/test_li_l2_nc.py @@ -627,6 +627,7 @@ def test_coords_generation(self, filetype_infos): projection = Proj(proj_dict) azimuth_vals = azimuth.values * point_height elevation_vals = elevation.values * point_height + azimuth_vals *= -1 lon_ref, lat_ref = projection(azimuth_vals, elevation_vals, inverse=True) # Convert to float32: lon_ref = lon_ref.astype(np.float32) From dd478273ab63abc56dc6d397c32183df94d27f5f Mon Sep 17 00:00:00 2001 From: andream Date: Tue, 21 May 2024 18:01:35 +0200 Subject: [PATCH 2/6] add *= -1 in reader code --- satpy/readers/li_base_nc.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/satpy/readers/li_base_nc.py b/satpy/readers/li_base_nc.py index eba9548985..cefbcc7e55 100644 --- a/satpy/readers/li_base_nc.py +++ b/satpy/readers/li_base_nc.py @@ -371,6 +371,9 @@ def inverse_projection(self, azimuth, elevation, proj_dict): azimuth = azimuth.values * point_height elevation = elevation.values * point_height + # In the MTG world, azimuth is defined as positive towards west, while proj expects it positive towards east + azimuth *= -1 + lon, lat = projection(azimuth, elevation, inverse=True) return np.stack([lon.astype(azimuth.dtype), lat.astype(elevation.dtype)]) From 766c4ece53fb8d7e4cb19ae5c00e2d7c023ad0bf Mon Sep 17 00:00:00 2001 From: andream Date: Tue, 21 May 2024 19:51:37 +0200 Subject: [PATCH 3/6] add a test to check 1-d and 2-d consistency --- satpy/tests/reader_tests/_li_test_utils.py | 16 ++++++---- satpy/tests/reader_tests/test_li_l2_nc.py | 34 ++++++++++++++++++---- 2 files changed, 38 insertions(+), 12 deletions(-) diff --git a/satpy/tests/reader_tests/_li_test_utils.py b/satpy/tests/reader_tests/_li_test_utils.py index 32107006fc..051743b238 100644 --- a/satpy/tests/reader_tests/_li_test_utils.py +++ b/satpy/tests/reader_tests/_li_test_utils.py @@ -521,9 +521,13 @@ def accumulation_dimensions(nacc, nobs): def fci_grid_definition(axis, nobs): """FCI grid definition on X or Y axis.""" + scale_factor = 5.58871526031607e-5 + add_offset = -0.15561777642350116 if axis == "X": long_name = "azimuth angle encoded as column" standard_name = "projection_x_coordinate" + scale_factor *= -1 + add_offset *= -1 else: long_name = "zenith angle encoded as row" standard_name = "projection_y_coordinate" @@ -531,10 +535,10 @@ def fci_grid_definition(axis, nobs): return { "format": "i2", "shape": ("pixels",), - "add_offset": -0.155619516, + "add_offset": add_offset, "axis": axis, "long_name": long_name, - "scale_factor": 5.58878e-5, + "scale_factor": scale_factor, "standard_name": standard_name, "units": "radian", "valid_range": np.asarray([1, 5568]), @@ -548,12 +552,12 @@ def mtg_geos_projection(): "format": "i4", "shape": ("accumulations",), "grid_mapping_name": "geostationary", - "inverse_flattening": 298.2572221, + "inverse_flattening": 298.257223563, "latitude_of_projection_origin": 0, "longitude_of_projection_origin": 0, - "perspective_point_height": 42164000, - "semi_major_axis": 6378169, - "semi_minor_axis": 6356583.8, + "perspective_point_height": 3.57864e7, + "semi_major_axis": 6378137.0, + "semi_minor_axis": 6356752.31424518, "sweep_angle_axis": "y", "long_name": "MTG geostationary projection", "default_data": lambda: -2147483647 diff --git a/satpy/tests/reader_tests/test_li_l2_nc.py b/satpy/tests/reader_tests/test_li_l2_nc.py index c5e02f93d0..bba36a2155 100644 --- a/satpy/tests/reader_tests/test_li_l2_nc.py +++ b/satpy/tests/reader_tests/test_li_l2_nc.py @@ -592,7 +592,6 @@ def test_generate_coords_called_once(Self, filetype_infos): def test_coords_generation(self, filetype_infos): """Compare daskified coords generation results with non-daskified.""" - # Prepare dummy (but somewhat realistic) arrays of azimuth/elevation values. products = ["li_l2_af_nc", "li_l2_afr_nc", "li_l2_afa_nc"] @@ -601,11 +600,10 @@ def test_coords_generation(self, filetype_infos): handler = LIL2NCFileHandler("filename", {}, extract_filetype_info(filetype_infos, prod)) # Get azimuth/elevation arrays from handler - azimuth = handler.get_measured_variable(handler.swath_coordinates["azimuth"]) - azimuth = handler.apply_use_rescaling(azimuth) - - elevation = handler.get_measured_variable(handler.swath_coordinates["elevation"]) - elevation = handler.apply_use_rescaling(elevation) + x = handler.get_measured_variable(handler.swath_coordinates["azimuth"]) + azimuth = handler.apply_use_rescaling(x) + y = handler.get_measured_variable(handler.swath_coordinates["elevation"]) + elevation = handler.apply_use_rescaling(y) # Initialize proj_dict proj_var = handler.swath_coordinates["projection"] @@ -641,6 +639,30 @@ def test_coords_generation(self, filetype_infos): np.testing.assert_equal(lon, lon_ref) np.testing.assert_equal(lat, lat_ref) + + def test_coords_and_grid_consistency(self, filetype_infos): + """Compare computed latlon coords for 1-d version with latlon from areadef as for the gridded version.""" + handler = LIL2NCFileHandler("filename", {}, extract_filetype_info(filetype_infos, "li_l2_af_nc"), + with_area_definition=True) + + # Get azimuth/elevation arrays from handler + x = handler.get_measured_variable(handler.swath_coordinates["azimuth"]) + y = handler.get_measured_variable(handler.swath_coordinates["elevation"]) + + handler.generate_coords_from_scan_angles() + lon = handler.internal_variables["longitude"].values + lat = handler.internal_variables["latitude"].values + + dsid = make_dataid(name="flash_accumulation") + area_def = handler.get_area_def(dsid) + rows = (LI_GRID_SHAPE[0] - y.astype(int)) + cols = x.astype(int) - 1 + lon_areadef, lat_areadef = area_def.get_lonlat_from_array_coordinates(cols, rows) + + np.testing.assert_allclose(lon, lon_areadef, rtol=1e-3) + np.testing.assert_allclose(lat, lat_areadef, rtol=1e-3) + + def test_get_area_def_acc_products(self, filetype_infos): """Test retrieval of area def for accumulated products.""" handler = LIL2NCFileHandler("filename", {}, extract_filetype_info(filetype_infos, "li_l2_af_nc"), From 3c0d46fc6e01a23feb06819ce9245294f5d1c894 Mon Sep 17 00:00:00 2001 From: andream Date: Wed, 22 May 2024 10:20:02 +0200 Subject: [PATCH 4/6] rearrange and add comments to test --- satpy/tests/reader_tests/test_li_l2_nc.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/satpy/tests/reader_tests/test_li_l2_nc.py b/satpy/tests/reader_tests/test_li_l2_nc.py index bba36a2155..76b9e74e04 100644 --- a/satpy/tests/reader_tests/test_li_l2_nc.py +++ b/satpy/tests/reader_tests/test_li_l2_nc.py @@ -639,30 +639,30 @@ def test_coords_generation(self, filetype_infos): np.testing.assert_equal(lon, lon_ref) np.testing.assert_equal(lat, lat_ref) - def test_coords_and_grid_consistency(self, filetype_infos): """Compare computed latlon coords for 1-d version with latlon from areadef as for the gridded version.""" handler = LIL2NCFileHandler("filename", {}, extract_filetype_info(filetype_infos, "li_l2_af_nc"), - with_area_definition=True) + with_area_definition=True) - # Get azimuth/elevation arrays from handler + # Get cols/rows arrays from handler x = handler.get_measured_variable(handler.swath_coordinates["azimuth"]) y = handler.get_measured_variable(handler.swath_coordinates["elevation"]) + cols = x.astype(int) - 1 + rows = (LI_GRID_SHAPE[0] - y.astype(int)) + # compute lonlat from 1-d coords generation handler.generate_coords_from_scan_angles() lon = handler.internal_variables["longitude"].values lat = handler.internal_variables["latitude"].values + # compute lonlat from 2-d areadef dsid = make_dataid(name="flash_accumulation") area_def = handler.get_area_def(dsid) - rows = (LI_GRID_SHAPE[0] - y.astype(int)) - cols = x.astype(int) - 1 lon_areadef, lat_areadef = area_def.get_lonlat_from_array_coordinates(cols, rows) np.testing.assert_allclose(lon, lon_areadef, rtol=1e-3) np.testing.assert_allclose(lat, lat_areadef, rtol=1e-3) - def test_get_area_def_acc_products(self, filetype_infos): """Test retrieval of area def for accumulated products.""" handler = LIL2NCFileHandler("filename", {}, extract_filetype_info(filetype_infos, "li_l2_af_nc"), From 2799a3225e784cddd48f2f2ed989be4fec965ded Mon Sep 17 00:00:00 2001 From: andream Date: Wed, 22 May 2024 11:56:18 +0200 Subject: [PATCH 5/6] restore coordinates variable naming in coords test --- satpy/tests/reader_tests/test_li_l2_nc.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/satpy/tests/reader_tests/test_li_l2_nc.py b/satpy/tests/reader_tests/test_li_l2_nc.py index 76b9e74e04..24280539a1 100644 --- a/satpy/tests/reader_tests/test_li_l2_nc.py +++ b/satpy/tests/reader_tests/test_li_l2_nc.py @@ -600,10 +600,11 @@ def test_coords_generation(self, filetype_infos): handler = LIL2NCFileHandler("filename", {}, extract_filetype_info(filetype_infos, prod)) # Get azimuth/elevation arrays from handler - x = handler.get_measured_variable(handler.swath_coordinates["azimuth"]) - azimuth = handler.apply_use_rescaling(x) - y = handler.get_measured_variable(handler.swath_coordinates["elevation"]) - elevation = handler.apply_use_rescaling(y) + azimuth = handler.get_measured_variable(handler.swath_coordinates["azimuth"]) + azimuth = handler.apply_use_rescaling(azimuth) + + elevation = handler.get_measured_variable(handler.swath_coordinates["elevation"]) + elevation = handler.apply_use_rescaling(elevation) # Initialize proj_dict proj_var = handler.swath_coordinates["projection"] From 1b52311852fecbb1cd1daef1233fa755315f298e Mon Sep 17 00:00:00 2001 From: andream Date: Wed, 22 May 2024 12:01:45 +0200 Subject: [PATCH 6/6] update comments and error message --- satpy/readers/li_l2_nc.py | 2 +- satpy/tests/reader_tests/test_li_l2_nc.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/satpy/readers/li_l2_nc.py b/satpy/readers/li_l2_nc.py index 4fe0826380..891372d596 100644 --- a/satpy/readers/li_l2_nc.py +++ b/satpy/readers/li_l2_nc.py @@ -73,7 +73,7 @@ def get_area_def(self, dsid): if var_with_swath_coord and self.with_area_def: return get_area_def("mtg_fci_fdss_2km") - raise NotImplementedError("Area definition is not supported for accumulated products.") + raise NotImplementedError("Area definition is not supported for non-accumulated products.") def is_var_with_swath_coord(self, dsid): """Check if the variable corresponding to this dataset is listed as variable with swath coordinates.""" diff --git a/satpy/tests/reader_tests/test_li_l2_nc.py b/satpy/tests/reader_tests/test_li_l2_nc.py index 24280539a1..1a198a7831 100644 --- a/satpy/tests/reader_tests/test_li_l2_nc.py +++ b/satpy/tests/reader_tests/test_li_l2_nc.py @@ -651,7 +651,7 @@ def test_coords_and_grid_consistency(self, filetype_infos): cols = x.astype(int) - 1 rows = (LI_GRID_SHAPE[0] - y.astype(int)) - # compute lonlat from 1-d coords generation + # compute lonlat from 1-d coords generation (called when with_area_definition==False) handler.generate_coords_from_scan_angles() lon = handler.internal_variables["longitude"].values lat = handler.internal_variables["latitude"].values