Add kwargs config option to turn off mitiff corner correction

satpy/tests/writer_tests/test_mitiff.py

@@ -1,4 +1,4 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 # Copyright (c) 2018-2020 Satpy developers
 #
@@ -879,6 +879,40 @@
             proj4_string = w._add_proj4_string(ds1, ds1)
             assert proj4_string == check["proj4"]
 
+    def test_correction_proj4_string(self):
+        """Test correction of proj4 lower left coordinate."""
+        import dask.array as da
+        import xarray as xr
+        from pyresample.geometry import AreaDefinition
+
+        from satpy.writers.mitiff import MITIFFWriter
+        area_def = AreaDefinition(
+            "test",
+            "test",
+            "test",
+            "+proj=merc",
+            100,
+            200,
+            (-1000., -1500., 1000., 1500.),
+        )
+
+        ds1 = xr.DataArray(
+            da.zeros((10, 20), chunks=20),
+            dims=("y", "x"),
+            attrs={"area": area_def}
+        )
+        default_expected_correction = (20.0, 15.0)
+        w = MITIFFWriter(filename="dummy.tif", base_dir=self.base_dir)
+        mitiff_pixel_adjustment = True
+        correction = w._set_correction_size(ds1, mitiff_pixel_adjustment)
+        assert correction == default_expected_correction
+
+        mitiff_pixel_adjustment = False
+        new_expected_correction = (0, 0)
+        w = MITIFFWriter(filename="dummy.tif", base_dir=self.base_dir)
+        correction = w._set_correction_size(ds1, mitiff_pixel_adjustment)
+        assert correction == new_expected_correction
+
     def test_save_dataset_palette(self):
         """Test writer operation as palette."""
         from satpy.writers.mitiff import MITIFFWriter

satpy/writers/mitiff.py

@@ -1,4 +1,4 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 # Copyright (c) 2018, 2019 Satpy developers
 #
@@ -220,72 +220,86 @@
 
         return _image_description
 
-    def _add_proj4_string(self, datasets, first_dataset):
+    def _add_proj4_string(self, datasets, first_dataset, **kwargs):
         import warnings
 
         proj4_string = " Proj string: "
 
         if isinstance(datasets, list):
             area = first_dataset.attrs["area"]
         else:
             area = datasets.attrs["area"]
         # Use pyproj's CRS object to get a valid EPSG code if possible
         # only in newer pyresample versions with pyproj 2.0+ installed
         if hasattr(area, "crs") and area.crs.to_epsg() is not None:
             proj4_string += "+init=EPSG:{}".format(area.crs.to_epsg())
         else:
             # Filter out the PROJ warning of losing projection information
             with warnings.catch_warnings():
                 warnings.filterwarnings("ignore", category=UserWarning,
                                         message=r"You will likely lose important projection information")
                 proj4_string += area.proj_str
 
         x_0 = 0
         y_0 = 0
         # FUTURE: Use pyproj 2.0+ to convert EPSG to PROJ4 if possible
         proj4_string, x_0 = self._convert_epsg_to_proj(proj4_string, x_0)
 
+        proj4_string = self._special_correction_of_proj_string(proj4_string)
+
+        if isinstance(datasets, list):
+            _dataset = first_dataset
+        else:
+            _dataset = datasets
+        mitiff_pixel_adjustment = kwargs.get("mitiff_pixel_adjustment", True)
+        proj4_string = self._append_projection_center(proj4_string, _dataset, x_0, y_0, mitiff_pixel_adjustment)
+        LOG.debug("proj4_string: %s", proj4_string)
+        proj4_string += "\n"
+
+        return proj4_string
+
+    def _special_correction_of_proj_string(self, proj4_string):
         if "geos" in proj4_string:
             proj4_string = proj4_string.replace("+sweep=x ", "")
             if "+a=6378137.0 +b=6356752.31414" in proj4_string:
                 proj4_string = proj4_string.replace("+a=6378137.0 +b=6356752.31414",
                                                     "+ellps=WGS84")
         if "+units=m" in proj4_string:
             proj4_string = proj4_string.replace("+units=m", "+units=km")
 
         if not any(datum in proj4_string for datum in ["datum", "towgs84"]):
             proj4_string += " +towgs84=0,0,0"
 
         if "units" not in proj4_string:
             proj4_string += " +units=km"
-
-        proj4_string = self._append_projection_center(proj4_string, datasets, first_dataset, x_0, y_0)
-        LOG.debug("proj4_string: %s", proj4_string)
-        proj4_string += "\n"
-
         return proj4_string
 
-    def _append_projection_center(self, proj4_string, datasets, first_dataset, x_0, y_0):
-        if isinstance(datasets, list):
-            dataset = first_dataset
-        else:
-            dataset = datasets
+    def _append_projection_center(self, proj4_string, dataset, x_0, y_0, mitiff_pixel_adjustment):
+        corner_correction_x, corner_correction_y = self._set_correction_size(dataset, mitiff_pixel_adjustment)
         if "x_0" not in proj4_string:
             proj4_string += " +x_0=%.6f" % (
                     (-dataset.attrs["area"].area_extent[0] +
-                     dataset.attrs["area"].pixel_size_x) + x_0)
+                     corner_correction_x) + x_0)
             proj4_string += " +y_0=%.6f" % (
                     (-dataset.attrs["area"].area_extent[1] +
-                     dataset.attrs["area"].pixel_size_y) + y_0)
+                     corner_correction_y) + y_0)
         elif "+x_0=0" in proj4_string and "+y_0=0" in proj4_string:
             proj4_string = proj4_string.replace("+x_0=0", "+x_0=%.6f" % (
                     (-dataset.attrs["area"].area_extent[0] +
-                     dataset.attrs["area"].pixel_size_x) + x_0))
+                     corner_correction_x) + x_0))
             proj4_string = proj4_string.replace("+y_0=0", "+y_0=%.6f" % (
                     (-dataset.attrs["area"].area_extent[1] +
-                     dataset.attrs["area"].pixel_size_y) + y_0))
+                     corner_correction_y) + y_0))
         return proj4_string
 
+    def _set_correction_size(self, dataset, mitiff_pixel_adjustment):
+        corner_correction_x = dataset.attrs["area"].pixel_size_x
+        corner_correction_y = dataset.attrs["area"].pixel_size_y
+        if not mitiff_pixel_adjustment:
+            corner_correction_x = 0
+            corner_correction_y = 0
+        return corner_correction_x,corner_correction_y
+
     def _convert_epsg_to_proj(self, proj4_string, x_0):
         if "EPSG:32631" in proj4_string:
             proj4_string = proj4_string.replace("+init=EPSG:32631",
@@ -563,7 +577,7 @@
 
         _image_description += " Map projection: Stereographic\n"
 
-        _image_description += self._add_proj4_string(datasets, first_dataset)
+        _image_description += self._add_proj4_string(datasets, first_dataset, **kwargs)
 
         _image_description += " TrueLat: 60N\n"
         _image_description += " GridRot: 0\n"