geotransform2mapinfo, mapinfo2geotransform: added compatibility to local...

geotransform2mapinfo, mapinfo2geotransform: added compatibility to local coordinate systems. Added test_map_info module.

py_tools_ds/geo/map_info.py

@@ -13,6 +13,7 @@ __author__ = "Daniel Scheffler"
 
 
 def geotransform2mapinfo(gt, prj):
+    # type: (Union[list, tuple, None], Union[str, None]) -> list
     """Builds an ENVI geo info from given GDAL GeoTransform and Projection (compatible with UTM and LonLat projections).
     :param gt:  GDAL GeoTransform, e.g. (249885.0, 30.0, 0.0, 4578615.0, 0.0, -30.0)
     :param prj: GDAL Projection - WKT Format
@@ -20,39 +21,46 @@ def geotransform2mapinfo(gt, prj):
     :rtype:     list
     """
 
-    if gt in [None, [0, 1, 0, 0, 0, -1]]:
-        return None
+    if gt in [None, [0, 1, 0, 0, 0, -1], (0, 1, 0, 0, 0, -1)]:
+        return ['Arbitrary', 1, 1, 0, 0, 1, 1, 0, 'North']
 
     if gt[2] != 0 or gt[4] != 0:  # TODO
         raise NotImplementedError('Currently rotated datasets are not supported.')
 
+    UL_X, UL_Y, GSD_X, GSD_Y = gt[0], gt[3], gt[1], gt[5]
+
     srs = osr.SpatialReference()
     srs.ImportFromWkt(prj)
-    Proj4 = [i[1:] for i in srs.ExportToProj4().split()]
-    Proj4_proj = [v.split('=')[1] for i, v in enumerate(Proj4) if '=' in v and v.split('=')[0] == 'proj'][0]
-    Proj4_ellps = \
-        [v.split('=')[1] for i, v in enumerate(Proj4) if '=' in v and v.split('=')[0] in ['ellps', 'datum']][0]
-    proj = 'Geographic Lat/Lon' if Proj4_proj == 'longlat' else 'UTM' if Proj4_proj == 'utm' else Proj4_proj
-    ellps = 'WGS-84' if Proj4_ellps == 'WGS84' else Proj4_ellps
-    UL_X, UL_Y, GSD_X, GSD_Y = gt[0], gt[3], gt[1], gt[5]
 
-    def utm2wgs84(utmX, utmY): return transform_utm_to_wgs84(utmX, utmY, srs.GetUTMZone())
+    if prj and not srs.IsLocal():
+        Proj4 = [i[1:] for i in srs.ExportToProj4().split()]
+        Proj4_proj = [v.split('=')[1] for i, v in enumerate(Proj4) if '=' in v and v.split('=')[0] == 'proj'][0]
+        Proj4_ellps = \
+            [v.split('=')[1] for i, v in enumerate(Proj4) if '=' in v and v.split('=')[0] in ['ellps', 'datum']][0]
+        proj = 'Geographic Lat/Lon' if Proj4_proj == 'longlat' else 'UTM' if Proj4_proj == 'utm' else Proj4_proj
+        ellps = 'WGS-84' if Proj4_ellps == 'WGS84' else Proj4_ellps
 
-    def is_UTM_North_South(LonLat): return 'North' if LonLat[1] >= 0. else 'South'
+        def utm2wgs84(utmX, utmY): return transform_utm_to_wgs84(utmX, utmY, srs.GetUTMZone())
+
+        def is_UTM_North_South(LonLat): return 'North' if LonLat[1] >= 0. else 'South'
+
+        north_south = is_UTM_North_South(utm2wgs84(UL_X, UL_Y))
+        map_info = [proj, 1, 1, UL_X, UL_Y, GSD_X, abs(GSD_Y), ellps] if proj != 'UTM' else \
+            ['UTM', 1, 1, UL_X, UL_Y, GSD_X, abs(GSD_Y), srs.GetUTMZone(), north_south, ellps]
+    else:
+        map_info = ['Arbitrary', 1, 1, UL_X, UL_Y, GSD_X, abs(GSD_Y), 0, 'North']
 
-    map_info = [proj, 1, 1, UL_X, UL_Y, GSD_X, abs(GSD_Y), ellps] if proj != 'UTM' else \
-        ['UTM', 1, 1, UL_X, UL_Y, GSD_X, abs(GSD_Y), srs.GetUTMZone(), is_UTM_North_South(utm2wgs84(UL_X, UL_Y)), ellps]
     return map_info
 
 
 def mapinfo2geotransform(map_info):
-    # type: (list) -> Union[list, None]
+    # type: (Union[list, None]) -> list
     """Builds GDAL GeoTransform tuple from an ENVI geo info.
 
     :param map_info: ENVI geo info (list), e.g., ['UTM', 1, 1, 192585.0, 5379315.0, 30.0, 30.0, 41, 'North', 'WGS-84']
     :returns:        GDAL GeoTransform, e.g. [249885.0, 30.0, 0.0, 4578615.0, 0.0, -30.0]
     """
-    if map_info:
+    if map_info and map_info[0] != 'Arbitrary':
         # FIXME rotated datasets are currently not supported -> function must return rotation at gt[2] and gt[4]
         # validate input map info
         exp_len = 10 if map_info[0] == 'UTM' else 8
@@ -67,6 +75,9 @@ def mapinfo2geotransform(map_info):
 
         return [ULmapX, float(map_info[5]), 0., ULmapY, 0., -float(map_info[6])]
 
+    else:
+        return [0, 1, 0, 0, 0, -1]
+
 
 def get_corner_coordinates(gdal_ds=None, gt=None, cols=None, rows=None):
     """Returns (ULxy, LLxy, LRxy, URxy) in the same coordinate units like the given geotransform."""

tests/test_map_info.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+"""
+test_map_info
+-------------
+
+Tests for `py_tools_ds.geo.map_info` module.
+"""
+
+import unittest
+
+from py_tools_ds.geo.map_info import geotransform2mapinfo, mapinfo2geotransform
+
+
+geotransform_utm = [331185.0, 30.0, -0.0, 5840115.0, -0.0, -30.0]
+map_info_utm = ['UTM', 1, 1, 331185.0, 5840115.0, 30.0, 30.0, 33, 'North', 'WGS-84']
+geotransform_local = [0, 1, 0, 0, 0, -1]
+map_info_local = ['Arbitrary', 1, 1, 0, 0, 1, 1, 0, 'North']
+
+wkt_utm = \
+    """
+    PROJCS["WGS 84 / UTM zone 33N",
+           GEOGCS["WGS 84",
+                  DATUM["WGS_1984",
+                        SPHEROID["WGS 84", 6378137, 298.257223563,
+                                 AUTHORITY["EPSG", "7030"]],
+                        AUTHORITY["EPSG", "6326"]],
+                  PRIMEM["Greenwich", 0,
+                         AUTHORITY["EPSG", "8901"]],
+                  UNIT["degree", 0.0174532925199433,
+                       AUTHORITY["EPSG", "9122"]],
+                  AUTHORITY["EPSG", "4326"]],
+           PROJECTION["Transverse_Mercator"],
+           PARAMETER["latitude_of_origin", 0],
+           PARAMETER["central_meridian", 15],
+           PARAMETER["scale_factor", 0.9996],
+           PARAMETER["false_easting", 500000],
+           PARAMETER["false_northing", 0],
+           UNIT["metre", 1,
+                AUTHORITY["EPSG", "9001"]],
+           AXIS["Easting", EAST],
+           AXIS["Northing", NORTH],
+           AUTHORITY["EPSG", "32633"]]
+    """
+
+
+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)
+        self.assertTrue(isinstance(map_info, list))
+        self.assertEqual(map_info, map_info_utm)
+
+    def test_gt_is_arbitrary_or_none(self):
+        # test gt=None
+        map_info = geotransform2mapinfo(gt=None, prj=wkt_utm)
+        self.assertTrue(isinstance(map_info, list))
+        self.assertEqual(map_info, map_info_local)
+
+        # test gt=[0, 1, 0, 0, 0, -1]
+        map_info = geotransform2mapinfo(gt=geotransform_local, prj=wkt_utm)
+        self.assertTrue(isinstance(map_info, list))
+        self.assertEqual(map_info, map_info_local)
+
+    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='')
+        self.assertTrue(isinstance(map_info, list))
+        self.assertEqual(map_info, exp_map_info)
+
+
+class Test_mapinfo2geotransform(unittest.TestCase):
+
+    def test_UTM(self):
+        gt = mapinfo2geotransform(map_info_utm)
+        self.assertTrue(isinstance(gt, list))
+        self.assertEqual(gt, geotransform_utm)
+
+    def test_map_info_is_empty(self):
+        gt = mapinfo2geotransform(None)
+        self.assertTrue(isinstance(gt, list))
+        self.assertEqual(gt, geotransform_local)