Merge pull request #2804 from ameraner/fix_li_1dacc_latloncomputation

Fix LI L2 accumulated products `'with_area_definition': False` 1-d coordinates computation

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)])

satpy/readers/li_l2_nc.py

@@ -107,7 +107,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."""

satpy/tests/reader_tests/_li_test_utils.py

@@ -522,20 +522,24 @@ 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"
 
     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]),
@@ -549,12 +553,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

satpy/tests/reader_tests/test_li_l2_nc.py

@@ -598,7 +598,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"]
@@ -633,6 +632,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)
@@ -646,6 +646,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 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 (called when with_area_definition==False)
+        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)
+        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"),