Replaced old code with simple pyproj calls.

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

py_tools_ds/geo/projection.py

@@ -22,7 +22,7 @@
 # with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 import re
-import pyproj
+from pyproj import CRS
 from typing import Union  # noqa F401  # flake8 issue
 
 # custom
@@ -68,8 +68,7 @@ def proj4_to_dict(proj4):
     """Converts a PROJ4-like string into a dictionary.
     :param proj4:   <str> the PROJ4-like string
     """
-    items = [item for item in proj4.replace('+', '').split(' ') if '=' in item]
-    return {k: v for k, v in [kv.split('=') for kv in items]}
+    return CRS.from_proj4(proj4).to_dict()
 
 
 def dict_to_proj4(proj4dict):
@@ -77,7 +76,7 @@ def dict_to_proj4(proj4dict):
     """Converts a PROJ4-like dictionary into a PROJ4 string.
     :param proj4dict:   <dict> the PROJ4-like dictionary
     """
-    return pyproj.Proj(proj4dict).srs
+    return CRS.from_dict(proj4dict).to_proj4()
 
 
 def proj4_to_WKT(proj4str):
@@ -85,9 +84,7 @@ def proj4_to_WKT(proj4str):
     """Converts a PROJ4-like string into a WKT string.
     :param proj4str:   <dict> the PROJ4-like string
     """
-    srs = osr.SpatialReference()
-    srs.ImportFromProj4(proj4str)
-    return srs.ExportToWkt()
+    return CRS.from_proj4(proj4str).to_wkt()
 
 
 def prj_equal(prj1, prj2):
@@ -112,17 +109,9 @@ def isProjectedOrGeographic(prj):
     if prj is None:
         return None
 
-    srs = osr.SpatialReference()
-    if prj.startswith('EPSG:'):
-        srs.ImportFromEPSG(int(prj.split(':')[1]))
-    elif prj.startswith('+proj='):
-        srs.ImportFromProj4(prj)
-    elif prj.startswith('GEOGCS') or prj.startswith('PROJCS'):
-        srs.ImportFromWkt(prj)
-    else:
-        raise RuntimeError('Unknown input projection.')
+    crs = CRS.from_user_input(prj)
 
-    return 'projected' if srs.IsProjected() else 'geographic' if srs.IsGeographic() else None
+    return 'projected' if crs.is_projected else 'geographic' if crs.is_geographic else None
 
 
 def isLocal(prj):
@@ -149,80 +138,26 @@ def isLocal(prj):
 
 def EPSG2Proj4(EPSG_code):
     # type: (int) -> str
-    if EPSG_code is not None:
-        srs = osr.SpatialReference()
-        srs.ImportFromEPSG(EPSG_code)
-        proj4 = srs.ExportToProj4()
-
-        if not proj4:
-            raise EnvironmentError(gdal.GetLastErrorMsg())
-
-        return proj4
-    else:
-        return ''
+    return CRS.from_epsg(EPSG_code).to_proj4() if EPSG_code is not None else ''
 
 
 def EPSG2WKT(EPSG_code):
     # type: (int) -> str
-    if EPSG_code is not None:
-        srs = osr.SpatialReference()
-        srs.ImportFromEPSG(EPSG_code)
-        wkt = srs.ExportToWkt()
-
-        if not wkt:
-            raise EnvironmentError(gdal.GetLastErrorMsg())
-
-        return wkt
-    else:
-        return ''
+    return CRS.from_epsg(EPSG_code).to_wkt(pretty=False) if EPSG_code is not None else ''
 
 
-def WKT2EPSG(wkt, epsgfile=''):
-    # type: (str, str) -> Union[int, None]
+def WKT2EPSG(wkt):
+    # type: (str) -> Union[int, None]
     """ Transform a WKT string to an EPSG code
     :param wkt:  WKT definition
-    :param epsgfile: the proj.4 epsg file (automatically located if no path is provided)
     :returns:    EPSG code
-    http://gis.stackexchange.com/questions/20298/is-it-possible-to-get-the-epsg-value-from-an-osr-spatialreference-class-using-th
     """
-    # FIXME this function returns None if datum=NAD27 but works with datum=WGS84, e.g.:
-    # FIXME {PROJCS["UTM_Zone_33N",GEOGCS["GCS_North_American_1927",DATUM["D_North_American_1927",SPHEROID
-    # FIXME ["Clarke_1866",6378206.4,294.9786982]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],
-    # FIXME PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],
-    # FIXME PARAMETER["Central_Meridian",15.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],
-    # FIXME UNIT["Meter",1.0]]}
-
     if not isinstance(wkt, str):
         raise TypeError("'wkt' must be a string. Received %s." % type(wkt))
-    code = None  # default
     if not wkt:
         return None
-    p_in = osr.SpatialReference()
-    s = p_in.ImportFromWkt(wkt)
-    if s == 5:  # invalid WKT
-        raise Exception('Received an invalid WKT string: %s' % wkt)
-    if p_in.IsLocal():
-        raise Exception('The given WKT is a local coordinate system.')
-    cstype = 'GEOGCS' if p_in.IsGeographic() else 'PROJCS'
-    p_in.AutoIdentifyEPSG()
-    an = p_in.GetAuthorityName(cstype)
-    assert an in [None, 'EPSG'], "No EPSG code found. Found %s instead." % an
-    ac = p_in.GetAuthorityCode(cstype)
-    if ac is None:  # try brute force approach by grokking proj epsg definition file
-        p_out = p_in.ExportToProj4()
-        if p_out:
-            epsgfile = epsgfile or gdal_env.find_epsgfile()
-            with open(epsgfile) as f:
-                for line in f:
-                    if line.find(p_out) != -1:
-                        m = re.search('<(\\d+)>', line)
-                        if m:
-                            code = m.group(1)
-                            break
-                code = int(code) if code else None  # match or no match
-    else:
-        code = int(ac)
-    return code
+
+    return CRS.from_wkt(wkt.replace('\n', '').replace('\r', '').replace(' ', '')).to_epsg()
 
 
 def get_UTMzone(ds=None, prj=None):
@@ -240,12 +175,10 @@ def get_prjLonLat(fmt='wkt'):
     """Returns standard geographic projection (EPSG 4326) in the WKT or PROJ4 format.
     :param fmt:     <str> target format - 'WKT' or 'PROJ4'
     """
-    assert re.search('wkt', fmt, re.I) or re.search('Proj4', fmt, re.I), 'unsupported output format'
-    srs = osr.SpatialReference()
-    srs.ImportFromEPSG(4326)
-    out = srs.ExportToWkt() if re.search('wkt', fmt, re.I) else srs.ExportToProj4()
+    if not re.search('wkt', fmt, re.I) or re.search('Proj4', fmt, re.I):
+        raise ValueError(fmt, 'Unsupported output format.')
 
-    if not out:
-        raise EnvironmentError(gdal.GetLastErrorMsg())
-
-    return out
+    if re.search('wkt', fmt, re.I):
+        return CRS.from_epsg(4326).to_wkt()
+    else:
+        return CRS.from_epsg(4326).to_proj4()

tests/test_geo/test_projection.py

@@ -68,7 +68,7 @@ class Test_WKT2EPSG(unittest.TestCase):
         self.wkt_utm = wkt_utm
 
     def test_UTM_wkt(self):
-        epsg = WKT2EPSG(self.wkt_utm, epsgfile='')
+        epsg = WKT2EPSG(self.wkt_utm)
         self.assertTrue(isinstance(epsg, int))
 
 
@@ -80,4 +80,4 @@ class Test_EPSG2WKT(unittest.TestCase):
     def test_UTM_epsg(self):
         wkt = EPSG2WKT(self.epsg_utm)
         self.assertTrue(isinstance(wkt, str), "EPSG2WKT returned a %s object instead of a string!" % type(wkt))
-        self.assertNotEquals(wkt, "", msg="EPSG2WKT returned an empty WKT string!")
+        self.assertNotEqual(wkt, "", msg="EPSG2WKT returned an empty WKT string!")