Merge branch 'enhancement/revise_projection' into 'master'

Enhancement/revise projection

See merge request !21

.gitlab-ci.yml

@@ -24,6 +24,7 @@ test_py_tools_ds:
     - source activate ci_env
     - export GDAL_DATA=/root/miniconda3/envs/ci_env/share/gdal
     - export PYTHONPATH=$PYTHONPATH:/root  # /root <- directory needed later
+    - conda install -c conda-forge 'pyproj>=2.1'  # FIXME remove as soon as docker container is rebuilt
     - make nosetests
     - make docs
   artifacts:
@@ -66,7 +67,7 @@ test_py_tools_ds_install:
 #    - conda config --set channel_priority strict  # otherwise gdal or libgdal may be installed from defaults channel
 
     # resolve some requirements with conda
-    - conda install --yes -q -c conda-forge numpy gdal pyproj shapely geopandas
+    - conda install --yes -q -c conda-forge numpy gdal 'pyproj>=2.1.0' shapely geopandas
 
     # run installer
     - python setup.py install

README.rst

@@ -62,7 +62,7 @@ Using conda_, the recommended approach is:
     # create virtual environment for py_tools_ds, this is optional
     conda create -y -q -c conda-forge --name py_tools_ds python=3
     source activate py_tools_ds
-    conda install -c conda-forge numpy gdal pyproj shapely scikit-image pandas
+    conda install -c conda-forge numpy gdal 'pyproj>=2.1.0' shapely scikit-image pandas
 
  Then install py_tools_ds using the pip installer:
 

py_tools_ds/geo/coord_trafo.py

@@ -21,8 +21,6 @@
 # You should have received a copy of the GNU Lesser General Public License along
 # with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-from functools import partial
-
 import pyproj
 import numpy as np
 from shapely.ops import transform
@@ -35,8 +33,6 @@ except ImportError:
     import osr
     import gdal
 
-from .projection import get_proj4info
-
 __author__ = "Daniel Scheffler"
 
 
@@ -73,16 +69,9 @@ def transform_any_prj(prj_src, prj_tgt, x, y):
     :param y:
     :return:
     """
-    from pyproj import __version__ as ver
-    if ver.startswith('2') and int(ver.split('.')[1]) >= 1:
-        transformer = pyproj.Transformer.from_crs(get_proj4info(proj=prj_src),
-                                                  get_proj4info(proj=prj_tgt))
-        x, y = transformer.transform(x, y)
-    else:
-        prj_src = pyproj.Proj(get_proj4info(proj=prj_src))
-        prj_tgt = pyproj.Proj(get_proj4info(proj=prj_tgt))
-        x, y = pyproj.transform(prj_src, prj_tgt, x, y)
-    return x, y
+    transformer = pyproj.Transformer.from_crs(pyproj.CRS.from_user_input(prj_src),
+                                              pyproj.CRS.from_user_input(prj_tgt))
+    return transformer.transform(x, y)
 
 
 def transform_coordArray(prj_src, prj_tgt, Xarr, Yarr, Zarr=None):
@@ -313,8 +302,8 @@ def reproject_shapelyGeometry(shapelyGeometry, prj_src, prj_tgt):
     :param prj_src:         GDAL projection as WKT string or EPSG code ('epsg:1234' or <EPSG_int>)
     :param prj_tgt:         GDAL projection as WKT string or EPSG code ('epsg:1234' or <EPSG_int>)
     """
+    project = pyproj.Transformer.from_proj(
+        pyproj.CRS.from_user_input(prj_src),
+        pyproj.CRS.from_user_input(prj_tgt))
 
-    project = partial(pyproj.transform, pyproj.Proj(get_proj4info(proj=prj_src)),
-                      pyproj.Proj(get_proj4info(proj=prj_tgt)))
-
-    return transform(project, shapelyGeometry)  # apply projection
+    return transform(project.transform, shapelyGeometry)  # apply projection

py_tools_ds/geo/map_info.py

@@ -27,6 +27,7 @@ import warnings
 import os
 from typing import Union  # noqa F401  # flake8 issue
 
+from pyproj import CRS
 import numpy as np
 try:
     from osgeo import osr
@@ -37,7 +38,7 @@ except ImportError:
     import gdal
 
 from .coord_trafo import transform_any_prj
-from .projection import get_proj4info, proj4_to_dict, isLocal
+from .projection import isLocal
 
 __author__ = "Daniel Scheffler"
 
@@ -130,14 +131,17 @@ class Geocoding(object):
             if isLocal(prj):
                 self.prj_name = 'Arbitrary'
             else:
-                # get prj_name and datum
-                proj4 = proj4_to_dict(get_proj4info(proj=prj))  # type: dict
-                self.prj_name = \
-                    'Geographic Lat/Lon' if proj4['proj'] == 'longlat' else \
-                    'UTM' if proj4['proj'] == 'utm' else proj4['proj']
-                # FIXME there is no 'datum' key in case of, e.g., LAEA projection
-                self.datum = 'WGS-84' if proj4['datum'] == 'WGS84' else proj4['datum']  # proj4['ellps']?
-                self.units = proj4['unit'] if 'unit' in proj4 else self.units
+                with warnings.catch_warnings():
+                    warnings.filterwarnings("ignore", message="You will likely lose important projection information")
+
+                    # get prj_name and datum
+                    proj4 = CRS(prj).to_dict()  # FIXME avoid to convert to proj4
+                    self.prj_name = \
+                        'Geographic Lat/Lon' if proj4['proj'] == 'longlat' else \
+                        'UTM' if proj4['proj'] == 'utm' else proj4['proj']
+                    # FIXME there is no 'datum' key in case of, e.g., LAEA projection  # still true?
+                    self.datum = 'WGS-84' if proj4['datum'] == 'WGS84' else proj4['datum']  # proj4['ellps']?
+                    self.units = proj4['unit'] if 'unit' in proj4 else self.units
 
                 if self.prj_name == 'UTM':
                     self.utm_zone = srs.GetUTMZone()
@@ -249,11 +253,16 @@ def geotransform2mapinfo(gt, prj):
         fdir = os.path.join(os.path.abspath(os.curdir))
 
         try:
+            if int(gdal.__version__[0]) < 3:
+                prj = CRS(prj).to_wkt(version="WKT1_GDAL")  # noqa
+
             ds_out = gdal.GetDriverByName('ENVI').Create(fn_bsq, 2, 2, 1, gdal.GDT_Int32)
             ds_out.SetGeoTransform(gt)
             ds_out.SetProjection(prj)
             ds_out.FlushCache()
-            del ds_out
+
+            # noinspection PyUnusedLocal
+            ds_out = None
 
             with open(fn_hdr, 'r') as inF:
                 content = inF.read()
@@ -296,7 +305,9 @@ def mapinfo2geotransform(map_info):
             ds_out = gdal.GetDriverByName('ENVI').Create(fn_bsq, 2, 2, 1)
             ds_out.GetRasterBand(1).WriteArray(np.array([[1, 2], [2, 3]]))
             ds_out.FlushCache()
-            del ds_out
+
+            # noinspection PyUnusedLocal
+            ds_out = None
 
             with open(fn_hdr, 'r') as InHdr:
                 lines = InHdr.readlines()

py_tools_ds/geo/projection.py

@@ -22,16 +22,14 @@
 # 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
 try:
     from osgeo import osr
-    from osgeo import gdal
 except ImportError:
     import osr
-    import gdal
 
 from ..environment import gdal_env
 
@@ -42,25 +40,15 @@ __author__ = "Daniel Scheffler"
 gdal_env.try2set_GDAL_DATA()
 
 
-def get_proj4info(ds=None, proj=None):
-    # type: (gdal.Dataset, Union[str, int]) -> str
-    """Returns PROJ4 formatted projection info for the given gdal.Dataset or projection respectively,
+def get_proj4info(proj=None):
+    # type: (Union[str, int]) -> str
+    """Returns PROJ4 formatted projection info for the given projection.
+
     e.g. '+proj=utm +zone=43 +datum=WGS84 +units=m +no_defs '
 
-    :param ds:      <gdal.Dataset> the gdal dataset to get PROJ4 info for
     :param proj:    <str,int> the projection to get PROJ4 formatted info for (WKT or 'epsg:1234' or <EPSG_int>)
     """
-    assert ds or proj, "Specify at least one of the arguments 'ds' or 'proj'"
-    srs = osr.SpatialReference()
-    proj = ds.GetProjection() if ds else proj
-
-    if isinstance(proj, str) and proj.startswith('epsg:'):
-        proj = EPSG2WKT(int(proj.split(':')[1]))
-    elif isinstance(proj, int):
-        proj = EPSG2WKT(proj)
-
-    srs.ImportFromWkt(proj)
-    return srs.ExportToProj4().strip()
+    return CRS.from_user_input(proj).to_proj4().strip()
 
 
 def proj4_to_dict(proj4):
@@ -68,8 +56,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 +64,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 +72,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):
@@ -100,7 +85,16 @@ def prj_equal(prj1, prj2):
     if prj1 is None and prj2 is None or prj1 == prj2:
         return True
     else:
-        return get_proj4info(proj=prj1) == get_proj4info(proj=prj2)
+        from pyproj import __version__ as ver
+        if ver.startswith('2') and int(ver.split('.')[1]) >= 5:
+            # CRS.equals was added in pyproj 2.5 which does not exist for Python 2.7 in conda-forge channel
+            crs1 = CRS.from_user_input(prj1)
+            crs2 = CRS.from_user_input(prj2)
+
+            return crs1.equals(crs2)
+
+        else:
+            return get_proj4info(proj=prj1) == get_proj4info(proj=prj2)
 
 
 def isProjectedOrGeographic(prj):
@@ -112,17 +106,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):
@@ -139,7 +125,11 @@ def isLocal(prj):
         srs.ImportFromEPSG(int(prj.split(':')[1]))
     elif prj.startswith('+proj='):
         srs.ImportFromProj4(prj)
-    elif 'GEOGCS' in prj or 'PROJCS' in prj or 'LOCAL_CS' in prj:
+    elif 'GEOGCS' in prj or \
+         'GEOGCRS' in prj or \
+         'PROJCS' in prj or \
+         'PROJCS' in prj or \
+         'LOCAL_CS' in prj:
         srs.ImportFromWkt(prj)
     else:
         raise RuntimeError('Unknown input projection: \n%s' % prj)
@@ -149,80 +139,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 +176,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 out:
-        raise EnvironmentError(gdal.GetLastErrorMsg())
+    if not re.search('wkt', fmt, re.I) or re.search('Proj4', fmt, re.I):
+        raise ValueError(fmt, 'Unsupported output format.')
 
-    return out
+    if re.search('wkt', fmt, re.I):
+        return CRS.from_epsg(4326).to_wkt()
+    else:
+        return CRS.from_epsg(4326).to_proj4()

requirements.txt

@@ -6,5 +6,5 @@ six
 scikit-image>=0.16.0
 pandas
 geopandas
-pyproj
+pyproj>=2.1.0
 spectral

setup.py

@@ -36,8 +36,8 @@ version = {}
 with open("py_tools_ds/version.py") as version_file:
     exec(version_file.read(), version)
 
-requirements = ['gdal>=2.1.0', 'numpy', 'shapely', 'six', 'pandas', 'scikit-image>=0.16.0', 'geopandas', 'pyproj',
-                'spectral']
+requirements = ['gdal>=2.1.0', 'numpy', 'shapely', 'six', 'pandas', 'scikit-image>=0.16.0', 'geopandas',
+                'pyproj>=2.1.0', 'spectral']
 setup_requirements = []  # TODO(danschef): put setup requirements (distutils extensions, etc.) here
 test_requirements = requirements + ["coverage", "nose", "nose2", "nose-htmloutput", "rednose"]
 

tests/CI_docker/context/environment_py_tools_ds.yml

@@ -10,7 +10,7 @@ dependencies:
   - gdal>=2.1.0
   - conda-forge::libgdal  # force to use conda-forge for libgdal to avoid package version incompatiblies due to mixed channels (libkea.so.1.4.7: cannot open shared object file: No such file or directory)
   # - kealib    # fix for libkea.so.1.4.7: cannot open shared object file: No such file or directory (not needed as long as libgdal and gdal are conda-forge packages)
-  - pyproj
+  - pyproj>=2.1.0
   - scikit-image>=0.16.0
   - geopandas
   - pip:

tests/test_geo/test_coord_trafo.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+# py_tools_ds
+#
+# Copyright (C) 2019  Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de)
+#
+# This software was developed within the context of the GeoMultiSens project funded
+# by the German Federal Ministry of Education and Research
+# (project grant code: 01 IS 14 010 A-C).
+#
+# This program is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Lesser General Public License as published by the Free
+# Software Foundation, either version 3 of the License, or (at your option) any
+# later version.
+#
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+# details.
+#
+# You should have received a copy of the GNU Lesser General Public License along
+# with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+
+"""
+test_coord_trafo
+----------------
+
+Tests for `py_tools_ds.geo.coord_trafo` module.
+"""
+
+import unittest
+from shapely.geometry import Polygon
+
+from py_tools_ds.geo.coord_trafo import reproject_shapelyGeometry
+
+
+poly_local = Polygon([(5708.2, -3006), (5708, -3262), (5452, -3262), (5452, -3006), (5708, -3006)])
+
+
+class Test_reproject_shapelyGeometry(unittest.TestCase):
+
+    def test_reproject_shapelyGeometry(self):
+        poly_lonlat = reproject_shapelyGeometry(poly_local, 32636, 4326)
+        self.assertTrue(isinstance(poly_lonlat, Polygon))

tests/test_geo/test_map_info.py

@@ -71,6 +71,7 @@ wkt_utm = \
            AXIS["Northing", NORTH],
            AUTHORITY["EPSG", "32633"]]
     """
+wkt_utm = ' '.join(wkt_utm.split())
 
 
 class Test_geotransform2mapinfo(unittest.TestCase):

tests/test_geo/test_projection.py

@@ -65,10 +65,10 @@ wkt_utm = \
 class Test_WKT2EPSG(unittest.TestCase):
 
     def setUp(self):
-        self.wkt_utm = wkt_utm