From 71afcda5608938aa49439497b76ca2eea7343a64 Mon Sep 17 00:00:00 2001
From: BENR0 <mail2benr@gmail.com>
Date: Wed, 16 Nov 2022 12:19:52 +0100
Subject: [PATCH] feat_swathdefinition_with_gcp

---
 satpy/etc/readers/modis_l1b.yaml | 20 ++++++++++++++++-
 satpy/readers/hdfeos_base.py     |  4 ++++
 satpy/readers/yaml_reader.py     | 37 ++++++++++++++++++++++++++++++++
 3 files changed, 60 insertions(+), 1 deletion(-)

diff --git a/satpy/etc/readers/modis_l1b.yaml b/satpy/etc/readers/modis_l1b.yaml
index d39ffbb99c..1347c137d9 100644
--- a/satpy/etc/readers/modis_l1b.yaml
+++ b/satpy/etc/readers/modis_l1b.yaml
@@ -403,12 +403,30 @@ datasets:
     name: '36'
     resolution: 1000
     calibration: [brightness_temperature, radiance]
-    coordinates: [longitude, latitude]
+    coordinates: [longitude, latitude, longitude_tiepoints, latitude_tiepoints]
     wavelength:
     - 14.085
     - 14.235
     - 14.385
 
+  longitude_tiepoints:
+    name: longitude_tiepoints
+    resolution:
+     1000:
+        file_type: hdf_eos_data_1000m
+    file_key: Longitude
+    standard_name: longitude_tiepoints
+    units: degree
+
+  latitude_tiepoints:
+    name: latitude_tiepoints
+    resolution:
+      1000:
+        file_type: hdf_eos_data_1000m
+    file_key: Latitude
+    standard_name: latitude_tiepoints
+    units: degree
+
   longitude:
     name: longitude
     resolution:
diff --git a/satpy/readers/hdfeos_base.py b/satpy/readers/hdfeos_base.py
index 3f9e62bf9b..16f43fee4a 100644
--- a/satpy/readers/hdfeos_base.py
+++ b/satpy/readers/hdfeos_base.py
@@ -296,6 +296,8 @@ class HDFEOSGeoReader(HDFEOSBaseFileReader):
     DATASET_NAMES = {
         'longitude': 'Longitude',
         'latitude': 'Latitude',
+        'longitude_tiepoints': 'Longitude',
+        'latitude_tiepoints': 'Latitude',
         'satellite_azimuth_angle': ('SensorAzimuth', 'Sensor_Azimuth'),
         'satellite_zenith_angle': ('SensorZenith', 'Sensor_Zenith'),
         'solar_azimuth_angle': ('SolarAzimuth', 'SolarAzimuth'),
@@ -393,6 +395,8 @@ def get_dataset(self, dataset_id: DataID, dataset_info: dict) -> xr.DataArray:
             # otherwise use the default name for this variable
             data = self._load_ds_by_name(dataset_name)
         if resolution != self.geo_resolution:
+            if dataset_name in ["longitude_tiepoints", "latitude_tiepoints"]:
+                return data
             if in_file_dataset_name is not None:
                 # they specified a custom variable name but
                 # we don't know how to interpolate this yet
diff --git a/satpy/readers/yaml_reader.py b/satpy/readers/yaml_reader.py
index ba1aeadbea..e56b787231 100644
--- a/satpy/readers/yaml_reader.py
+++ b/satpy/readers/yaml_reader.py
@@ -774,6 +774,12 @@ def _make_area_from_coords(self, coords):
             lons, lats = self._get_lons_lats_from_coords(coords)
             sdef = self._make_swath_definition_from_lons_lats(lons, lats)
             return sdef
+        if len(coords) == 4:
+            lons_tiepoints, lats_tiepoints = coords[2::]
+            lons, lats = self._get_lons_lats_from_coords(coords[0:2])
+            sdef = self._make_tiepoints_definition_from_lons_lats_and_gcps(lons, lats, lons_tiepoints,
+                                                                           lats_tiepoints)
+            return sdef
         if len(coords) != 0:
             raise NameError("Don't know what to do with coordinates " + str(
                 coords))
@@ -809,6 +815,37 @@ def _make_swath_definition_from_lons_lats(self, lons, lats):
                 FileYAMLReader._coords_cache[key] = sdef
         return sdef
 
+    def _make_tiepoints_definition_from_lons_lats_and_gcps(self, lons, lats, lons_tiepoints, lats_tiepoints):
+        """Make a swath definition instance from lons and lats.
+
+        Args:
+            lons (xarray.DataArray): Interpolated Longitudes.
+            lats (xarray.DataArray): Interpolated Latitudes.
+            lons_tiepoints (xarray.DataArray): Longitude tiepoints.
+            lats_tiepoints (xarray.DataArray): Latitude tiepoints.
+
+        Returns:
+            pyresample.geometry.GCPDefinition
+
+        """
+        from pyresample.geometry import GCPDefinition
+        key = None
+        try:
+            key = (lons.data.name, lats.data.name)
+            sdef = FileYAMLReader._coords_cache.get(key)
+        except AttributeError:
+            sdef = None
+        if sdef is None:
+            sdef = GCPDefinition(lons, lats, {"longitude": lons_tiepoints, "latitude": lats_tiepoints})
+            sensor_str = '_'.join(self.info['sensors'])
+            shape_str = '_'.join(map(str, lons.shape))
+            sdef.name = "{}_{}_{}_{}".format(sensor_str, shape_str,
+                                             lons.attrs.get('name', lons.name),
+                                             lats.attrs.get('name', lats.name))
+            if key is not None:
+                FileYAMLReader._coords_cache[key] = sdef
+        return sdef
+
     def _load_dataset_with_area(self, dsid, coords, **kwargs):
         """Load *dsid* and its area if available."""
         file_handlers = self._get_file_handlers(dsid)