Fixed an issue causing geotransform2mapinfo() return a 'South' mapinfo instead...

Fixed an issue causing geotransform2mapinfo() return a 'South' mapinfo instead of a 'North' one for a GeoTransform
  in the West of the central meridian.
Signed-off-by: Daniel Scheffler <danschef@gfz-potsdam.de>

HISTORY.rst

@@ -2,6 +2,13 @@
 History
 =======
 
+0.16.3 (2020-12-02)
+-------------------
+
+* Fixed an issue causing geotransform2mapinfo() return a 'South' mapinfo instead of a 'North' one for a GeoTransform
+  in the West of the central meridian.
+
+
 0.16.2 (2020-11-18)
 -------------------
 

py_tools_ds/geo/map_info.py

@@ -140,9 +140,10 @@ class Geocoding(object):
                     self.units = proj4['unit'] if 'unit' in proj4 else self.units
 
                 if self.prj_name == 'UTM':
-                    self.utm_zone = srs.GetUTMZone()
+                    utm_zone_srs = srs.GetUTMZone()
+                    self.utm_zone = abs(utm_zone_srs)  # always positive
                     self.utm_north_south = \
-                        'North' if transform_any_prj(prj, 4326, self.ul_x_map, self.ul_y_map)[0] >= 0. else 'South'
+                        'North' if utm_zone_srs >= 0. else 'South'
 
             del srs
 

tests/test_geo/test_map_info.py

@@ -36,13 +36,15 @@ from py_tools_ds.geo.map_info import geotransform2mapinfo, mapinfo2geotransform
 from py_tools_ds.geo.projection import EPSG2WKT
 
 
-geotransform_utm = (331185.0, 30.0, -0.0, 5840115.0, -0.0, -30.0)
-geotransform_utm_rotated = (331185.0, 12.202099292274006, 27.406363729278027,
-                            5840115.0, 27.406363729278027, -12.202099292274006)
+geotransform_utm_NE = (331185.0, 30.0, -0.0, 5840115.0, -0.0, -30.0)
+geotransform_utm_SE = (331185.0, 30.0, -0.0, -5840115.0, -0.0, -30.0)
+geotransform_utm_NE_rotated = (331185.0, 12.202099292274006, 27.406363729278027,
+                               5840115.0, 27.406363729278027, -12.202099292274006)
 geotransform_local = (0, 1, 0, 0, 0, -1)
 geotransform_local_rotated = (0.0, 6.123233995736766e-17, 1.0, 0.0, 1.0, -6.123233995736766e-17)
-map_info_utm = ['UTM', 1, 1, 331185.0, 5840115.0, 30.0, 30.0, 33, 'North', 'WGS-84']
-map_info_utm_rotated = ['UTM', 1, 1, 331185.0, 5840115.0, 30.0, 30.0, 33, 'North', 'WGS-84', 'rotation=66.00000']
+map_info_utm_NE = ['UTM', 1, 1, 331185.0, 5840115.0, 30.0, 30.0, 33, 'North', 'WGS-84']
+map_info_utm_SE = ['UTM', 1, 1, 331185.0, -5840115.0, 30.0, 30.0, 33, 'South', 'WGS-84']
+map_info_utm_NE_rotated = ['UTM', 1, 1, 331185.0, 5840115.0, 30.0, 30.0, 33, 'North', 'WGS-84', 'rotation=66.00000']
 map_info_local = ['Arbitrary', 1, 1, 0, 0, 1, 1, 0, 'North']
 map_info_local_rotated = ['Arbitrary', 1, 1, 0, 0, 1, 1, 0, 'North', 'rotation=90.00000']
 
@@ -78,10 +80,45 @@ class Test_geotransform2mapinfo(unittest.TestCase):
 
     # TODO implement test in case of geographic prj
 
-    def test_UTM_gt_prj(self):
-        map_info = geotransform2mapinfo(gt=geotransform_utm, prj=wkt_utm)
+    def test_UTM_gt_prj_NE(self):
+        map_info = geotransform2mapinfo(gt=geotransform_utm_NE, prj=wkt_utm)
         self.assertTrue(isinstance(map_info, list))
-        self.assertEqual(map_info, map_info_utm)
+        self.assertEqual(map_info, map_info_utm_NE)
+
+    def test_UTM_gt_prj_NW(self):
+        gt = (394185.0, 30.0, -0.0, 4572315.0, -0.0, -30.0)
+        prj = EPSG2WKT(32610)  # UTM zone 10 North
+        exp = ['UTM', 1.0, 1.0, 394185.0, 4572315.0, 30.0, 30.0, 10, 'North', 'WGS-84']
+
+        map_info = geotransform2mapinfo(gt=gt, prj=prj)
+        self.assertTrue(isinstance(map_info, list))
+        self.assertEqual(map_info, exp)
+
+    def test_UTM_gt_prj_SW(self):
+        # option 1: projection says 'South', gt has a positive y origin (e.g., Sentinel-2 data)
+        gt = (394185.0, 30.0, -0.0, 4572315.0, -0.0, -30.0)
+        prj = EPSG2WKT(32710)  # UTM zone 10 South
+        exp = ['UTM', 1.0, 1.0, 394185.0, 4572315.0, 30.0, 30.0, 10, 'South', 'WGS-84']
+
+        map_info = geotransform2mapinfo(gt=gt, prj=prj)
+        self.assertTrue(isinstance(map_info, list))
+        self.assertEqual(map_info, exp)
+
+        # option 2: projection says 'North', gt has a negative y origin (e.g., Landsat-5 data)
+        gt = (394185.0, 30.0, -0.0, -4572315.0, -0.0, -30.0)
+        prj = EPSG2WKT(32610)  # UTM zone 10 North
+        exp = ['UTM', 1.0, 1.0, 394185.0, -4572315.0, 30.0, 30.0, 10, 'North', 'WGS-84']
+
+        map_info = geotransform2mapinfo(gt=gt, prj=prj)
+        self.assertTrue(isinstance(map_info, list))
+        self.assertEqual(map_info, exp)
+
+    def test_UTM_gt_prj_SE(self):
+        prj = EPSG2WKT(32733)  # UTM zone 33 South
+
+        map_info = geotransform2mapinfo(gt=geotransform_utm_SE, prj=prj)
+        self.assertTrue(isinstance(map_info, list))
+        self.assertEqual(map_info, map_info_utm_SE)
 
     def test_gt_is_none(self):
         # test gt=None
@@ -97,7 +134,7 @@ class Test_geotransform2mapinfo(unittest.TestCase):
 
     def test_prj_is_empty(self):
         exp_map_info = ['Arbitrary', 1, 1, 331185.0, 5840115.0, 30.0, 30.0, 0, 'North']
-        map_info = geotransform2mapinfo(gt=geotransform_utm, prj='')
+        map_info = geotransform2mapinfo(gt=geotransform_utm_NE, prj='')
         self.assertTrue(isinstance(map_info, list))
         self.assertEqual(map_info, exp_map_info)
 
@@ -109,14 +146,14 @@ class Test_geotransform2mapinfo(unittest.TestCase):
 
     def test_prj_is_not_wgs84_rotation(self):
         exp_map_info = ['Lambert Azimuthal Equal Area', 1.0, 1.0, 331185.0, 5840115.0, 30.0, 30.0, 'rotation=66']
-        map_info = geotransform2mapinfo(gt=geotransform_utm_rotated, prj=EPSG2WKT(3035))
+        map_info = geotransform2mapinfo(gt=geotransform_utm_NE_rotated, prj=EPSG2WKT(3035))
         self.assertTrue(isinstance(map_info, list))
         self.assertEqual(map_info, exp_map_info)
 
     def test_gt_contains_rotation(self):
-        map_info = geotransform2mapinfo(gt=geotransform_utm_rotated, prj=wkt_utm)
+        map_info = geotransform2mapinfo(gt=geotransform_utm_NE_rotated, prj=wkt_utm)
         self.assertTrue(isinstance(map_info, list))
-        self.assertEqual(map_info, map_info_utm_rotated)
+        self.assertEqual(map_info, map_info_utm_NE_rotated)
 
     def test_gt_contains_rotation_prj_is_local(self):
         map_info = geotransform2mapinfo(gt=geotransform_local_rotated, prj='')
@@ -127,11 +164,15 @@ class Test_geotransform2mapinfo(unittest.TestCase):
 class Test_mapinfo2geotransform(unittest.TestCase):
 
     def test_map_info_is_valid(self):
-        gt = mapinfo2geotransform(map_info_utm)
+        gt = mapinfo2geotransform(map_info_utm_NE)
+        self.assertTrue(isinstance(gt, (tuple, list)))
+        self.assertEqual(gt, geotransform_utm_NE)
+
+        gt = mapinfo2geotransform(map_info_utm_SE)
         self.assertTrue(isinstance(gt, (tuple, list)))
-        self.assertEqual(gt, geotransform_utm)
+        self.assertEqual(gt, geotransform_utm_SE)
 
-        # test
+        # test arbitrary map info
         gt = mapinfo2geotransform(['Arbitrary', 1, 1, 5, -7, 1, 1, 0, 'North'])
         self.assertTrue(isinstance(gt, (tuple, list)))
         self.assertEqual(gt, (5, 1, 0, -7, 0, -1))
@@ -147,9 +188,9 @@ class Test_mapinfo2geotransform(unittest.TestCase):
         self.assertEqual(gt, (4526026.0, 10.0, 0.0, 3284919.5, 0.0, -10.0))
 
     def test_map_info_contains_rotation(self):
-        gt = mapinfo2geotransform(map_info_utm_rotated)
+        gt = mapinfo2geotransform(map_info_utm_NE_rotated)
         self.assertTrue(isinstance(gt, (tuple, list)))
-        self.assertEqual(gt, geotransform_utm_rotated)
+        self.assertEqual(gt, geotransform_utm_NE_rotated)
 
     def test_map_info_is_local_contains_rotation(self):
         gt = mapinfo2geotransform(map_info_local_rotated)