Added first algorithm to warp sensor to map geometry using GDAL, in-memory.

Signed-off-by: Daniel Scheffler <danschef@gfz-potsdam.de>

sensormapgeo/transformer_2d.py

@@ -391,6 +391,115 @@ class SensorMapGeometryTransformer(object):
 
         return data_mapgeo, out_gt, out_prj
 
+    def to_map_geometry(self, data: np.ndarray,
+                        tgt_prj: Union[str, int] = None,
+                        tgt_extent: Tuple[float, float, float, float] = None,
+                        tgt_res: Tuple = None,
+                        tgt_coordgrid: Tuple[Tuple, Tuple] = None,
+                        area_definition: AreaDefinition = None) -> Tuple[np.ndarray, tuple, str]:
+        """Transform the input sensor geometry array into map geometry.
+
+        :param data:            numpy array (representing sensor geometry) to be warped to map geometry
+        :param tgt_prj:         target projection (WKT or 'epsg:1234' or <EPSG_int>)
+        :param tgt_extent:      extent coordinates of output map geometry array (LL_x, LL_y, UR_x, UR_y) in the tgt_prj
+        :param tgt_res:         target X/Y resolution (e.g., (30, 30))
+        :param tgt_coordgrid:   target coordinate grid ((x, x), (y, y)):
+                                if given, the output extent is moved to this coordinate grid
+        :param area_definition: an instance of pyresample.geometry.AreaDefinition;
+                                OVERRIDES tgt_prj, tgt_extent and tgt_res; saves computation time
+        """
+        from osgeo import osr
+
+        if self.lons.ndim > 2 >= data.ndim:
+            raise ValueError(data.ndim, "'data' must at least have %d dimensions because of %d longiture array "
+                                        "dimensions." % (self.lons.ndim, self.lons.ndim))
+
+        if data.shape[:2] != self.lons.shape[:2]:
+            raise ValueError(data.shape, 'Expected a sensor geometry data array with %d rows and %d columns.'
+                             % self.lons.shape[:2])
+
+        tgt_wkt = CRS(tgt_prj).to_wkt(version='WKT1_GDAL')
+        # tgt_wkt = osr.GetUserInputAsWKT(tgt_prj)
+        wgs84_wkt = osr.GetUserInputAsWKT('WGS84')
+
+        gdal.UseExceptions()
+
+        p_lons_tmp = '/vsimem/lons.tif'
+        p_lats_tmp = '/vsimem/lats.tif'
+        p_data_tmp = '/vsimem/data.tif'
+        p_data_mapgeo_vrt = '/vsimem/data_mapgeo.vrt'
+
+        # save numpy arrays to temporary tif files
+        drv = gdal.GetDriverByName('GTiff')
+        rows, cols = self.lons.shape
+
+        ds_lons = drv.Create(p_lons_tmp, cols, rows, 1, gdal.GDT_Float64)
+        ds_lons.GetRasterBand(1).WriteArray(self.lons)
+        del ds_lons
+
+        ds_lats = drv.Create(p_lats_tmp, cols, rows, 1, gdal.GDT_Float64)
+        ds_lats.GetRasterBand(1).WriteArray(self.lats)
+        del ds_lats
+
+        ds_swath_data = drv.Create(p_data_tmp, cols, rows, 1, gdal.GDT_Float64)
+        ds_swath_data.GetRasterBand(1).WriteArray(data)
+
+        # add geolocation information to input data
+        ds_swath_data.SetMetadata(
+            dict(
+                SRS=wgs84_wkt,
+                X_DATASET=p_lons_tmp,
+                Y_DATASET=p_lats_tmp,
+                X_BAND='1',
+                Y_BAND='1',
+                PIXEL_OFFSET='0',
+                LINE_OFFSET='0',
+                PIXEL_STEP='1',
+                LINE_STEP='1'
+            ),
+            'GEOLOCATION'
+        )
+        del ds_swath_data
+
+        # TODO try to avoid GDAL cuts off the corners of the output image
+        # warp from geolocation arrays and read the result
+        tgt_epsg = CRS(tgt_prj).to_epsg()
+        tgt_extent = tgt_extent or self._get_target_extent(tgt_epsg)
+
+        if tgt_coordgrid:
+            tgt_extent = _move_extent_to_coordgrid(tgt_extent, *tgt_coordgrid)
+
+        xmin, ymin, xmax, ymax = tgt_extent
+
+        self.area_definition = self.compute_areadefinition_sensor2map(
+            data, tgt_prj=tgt_prj, tgt_extent=tgt_extent, tgt_res=tgt_res)
+
+        warpOptions = gdal.WarpOptions(
+            format='VRT',
+            # width=self.area_definition.width,
+            # height=self.area_definition.height,
+            outputBounds=tgt_extent,
+            # outputBounds=(xmin -1000, ymin -1000, xmax + 1000, ymax + 1000),
+            outputBoundsSRS=tgt_wkt,  # FIXME not needed?
+            xRes=30, yRes=30,  # FIXME hardcoded
+            targetAlignedPixels=True,
+            errorThreshold=0,
+            srcSRS=wgs84_wkt,
+            dstSRS=tgt_wkt,
+            geoloc=True,
+            resampleAlg='cubic'  # FIXME hardcoded
+        )
+
+        gdal.Warp(p_data_mapgeo_vrt, p_data_tmp, options=warpOptions)
+        ds_mapgeo = gdal.Open(p_data_mapgeo_vrt)
+
+        data_mapgeo = ds_mapgeo.ReadAsArray()
+        out_gt = ds_mapgeo.GetGeoTransform()
+        out_prj = ds_mapgeo.GetProjection()
+        del ds_mapgeo
+
+        return data_mapgeo, out_gt, out_prj
+
     def to_sensor_geometry(self, data: np.ndarray,
                            src_prj: Union[str, int],
                            src_extent: Tuple[float, float, float, float]) -> np.ndarray:

tests/test_sensormapgeo.py

@@ -139,6 +139,11 @@ class Test_SensorMapGeometryTransformer:
                 # GeoArray(dem_map_geo, dem_gt, dem_prj)\
                 #     .save(os.path.join(tests_path, 'test_output', 'resampled_pyresample_ll.bsq'))
 
+                # TODO: remove this
+                from geoarray import GeoArray
+                gA = GeoArray(dem_map_geo, dem_gt, dem_prj)
+                gA.show_map()
+
                 self.check_result_mapgeo(
                     dem_map_geo,
                     (SMGT.area_definition.height,
@@ -183,6 +188,11 @@ class Test_SensorMapGeometryTransformer:
                 # GeoArray(dem_map_geo, dem_gt, dem_prj)\
                 #     .save(os.path.join(tests_path, 'test_output', 'resampled_pyresample_bilinear_16n.bsq'))
 
+                # TODO: remove this
+                from geoarray import GeoArray
+                gA = GeoArray(dem_map_geo, dem_gt, dem_prj)
+                gA.show_map()
+
                 self.check_result_mapgeo(
                     dem_map_geo,
                     (366, 976),