projection.py 8.2 KB
Newer Older
Daniel Scheffler's avatar
Daniel Scheffler committed
1
2
# -*- coding: utf-8 -*-

3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# 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/>.

Daniel Scheffler's avatar
Daniel Scheffler committed
24
import re
Daniel Scheffler's avatar
Daniel Scheffler committed
25
import pyproj
Daniel Scheffler's avatar
Daniel Scheffler committed
26
from typing import Union  # noqa F401  # flake8 issue
Daniel Scheffler's avatar
Daniel Scheffler committed
27
28
29
30
31
32
33
34
35

# custom
try:
    from osgeo import osr
    from osgeo import gdal
except ImportError:
    import osr
    import gdal

36
37
38
39
from ..environment import gdal_env

__author__ = "Daniel Scheffler"

40

Daniel Scheffler's avatar
Daniel Scheffler committed
41
# try to set GDAL_DATA if not set or invalid
42
43
gdal_env.try2set_GDAL_DATA()

Daniel Scheffler's avatar
Daniel Scheffler committed
44

45
def get_proj4info(ds=None, proj=None):
46
    # type: (gdal.Dataset, Union[str, int]) -> str
47
    """Returns PROJ4 formatted projection info for the given gdal.Dataset or projection respectively,
Daniel Scheffler's avatar
Daniel Scheffler committed
48
    e.g. '+proj=utm +zone=43 +datum=WGS84 +units=m +no_defs '
49

Daniel Scheffler's avatar
Daniel Scheffler committed
50
    :param ds:      <gdal.Dataset> the gdal dataset to get PROJ4 info for
51
    :param proj:    <str,int> the projection to get PROJ4 formatted info for (WKT or 'epsg:1234' or <EPSG_int>)
Daniel Scheffler's avatar
Daniel Scheffler committed
52
53
54
    """
    assert ds or proj, "Specify at least one of the arguments 'ds' or 'proj'"
    srs = osr.SpatialReference()
55
56
57
58
    proj = ds.GetProjection() if ds else proj

    if isinstance(proj, str) and proj.startswith('epsg:'):
        proj = EPSG2WKT(int(proj.split(':')[1]))
59
    elif isinstance(proj, int):
60
61
62
        proj = EPSG2WKT(proj)

    srs.ImportFromWkt(proj)
63
    return srs.ExportToProj4().strip()
Daniel Scheffler's avatar
Daniel Scheffler committed
64
65
66


def proj4_to_dict(proj4):
67
    # type: (str) -> dict
Daniel Scheffler's avatar
Daniel Scheffler committed
68
69
70
    """Converts a PROJ4-like string into a dictionary.
    :param proj4:   <str> the PROJ4-like string
    """
71
72
    items = [item for item in proj4.replace('+', '').split(' ') if '=' in item]
    return {k: v for k, v in [kv.split('=') for kv in items]}
Daniel Scheffler's avatar
Daniel Scheffler committed
73
74


Daniel Scheffler's avatar
Daniel Scheffler committed
75
76
77
78
79
80
81
82
def dict_to_proj4(proj4dict):
    # type: (dict) -> str
    """Converts a PROJ4-like dictionary into a PROJ4 string.
    :param proj4dict:   <dict> the PROJ4-like dictionary
    """
    return pyproj.Proj(proj4dict).srs


Daniel Scheffler's avatar
Daniel Scheffler committed
83
84
85
86
87
88
89
90
91
92
def proj4_to_WKT(proj4str):
    # type: (str) -> str
    """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()


Daniel Scheffler's avatar
Daniel Scheffler committed
93
def prj_equal(prj1, prj2):
Daniel Scheffler's avatar
Daniel Scheffler committed
94
    # type: (Union[None, int, str], Union[None, int, str]) -> bool
95
96
97
98
99
    """Checks if the given two projections are equal.

    :param prj1: projection 1 (WKT or 'epsg:1234' or <EPSG_int>)
    :param prj2: projection 2 (WKT or 'epsg:1234' or <EPSG_int>)
    """
100
101
102
    if prj1 is None and prj2 is None or prj1 == prj2:
        return True
    else:
103
        return get_proj4info(proj=prj1) == get_proj4info(proj=prj2)
Daniel Scheffler's avatar
Daniel Scheffler committed
104
105
106


def isProjectedOrGeographic(prj):
Daniel Scheffler's avatar
Daniel Scheffler committed
107
    # type: (Union[str, int, dict]) -> Union[str, None]
108
109
    """

110
    :param prj: accepts EPSG, Proj4 and WKT projections
111
    """
112
113
    if prj is None:
        return None
114

Daniel Scheffler's avatar
Daniel Scheffler committed
115
116
117
118
119
120
121
122
    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:
123
124
        raise RuntimeError('Unknown input projection.')

Daniel Scheffler's avatar
Daniel Scheffler committed
125
126
127
    return 'projected' if srs.IsProjected() else 'geographic' if srs.IsGeographic() else None


128
129
130
131
132
133
134
135
136
137
138
139
140
141
def isLocal(prj):
    # type: (Union[str, int, dict]) -> Union[bool, None]
    """

    :param prj: accepts EPSG, Proj4 and WKT projections
    """
    if not prj:
        return True

    srs = osr.SpatialReference()
    if prj.startswith('EPSG:'):
        srs.ImportFromEPSG(int(prj.split(':')[1]))
    elif prj.startswith('+proj='):
        srs.ImportFromProj4(prj)
142
    elif 'GEOGCS' in prj or 'PROJCS' in prj or 'LOCAL_CS' in prj:
143
144
        srs.ImportFromWkt(prj)
    else:
145
        raise RuntimeError('Unknown input projection: \n%s' % prj)
146
147
148
149

    return srs.IsLocal()


Daniel Scheffler's avatar
Daniel Scheffler committed
150
def EPSG2Proj4(EPSG_code):
151
152
153
154
    # type: (int) -> str
    if EPSG_code is not None:
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(EPSG_code)
155
156
157
158
159
160
        proj4 = srs.ExportToProj4()

        if not proj4:
            raise EnvironmentError(gdal.GetLastErrorMsg())

        return proj4
161
162
    else:
        return ''
Daniel Scheffler's avatar
Daniel Scheffler committed
163
164
165


def EPSG2WKT(EPSG_code):
166
167
168
169
    # type: (int) -> str
    if EPSG_code is not None:
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(EPSG_code)
170
        wkt = srs.ExportToWkt()
171

172
173
174
175
176
177
        if not wkt:
            raise EnvironmentError(gdal.GetLastErrorMsg())

        return wkt
    else:
        return ''
Daniel Scheffler's avatar
Daniel Scheffler committed
178
179


180
def WKT2EPSG(wkt, epsgfile=''):
181
    # type: (str, str) -> Union[int, None]
Daniel Scheffler's avatar
Daniel Scheffler committed
182
    """ Transform a WKT string to an EPSG code
183
    :param wkt:  WKT definition
184
    :param epsgfile: the proj.4 epsg file (automatically located if no path is provided)
Daniel Scheffler's avatar
Daniel Scheffler committed
185
186
187
    :returns:    EPSG code
    http://gis.stackexchange.com/questions/20298/is-it-possible-to-get-the-epsg-value-from-an-osr-spatialreference-class-using-th
    """
188
    # FIXME this function returns None if datum=NAD27 but works with datum=WGS84, e.g.:
189
190
191
192
193
    # 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]]}
Daniel Scheffler's avatar
Daniel Scheffler committed
194

195
196
197
    if not isinstance(wkt, str):
        raise TypeError("'wkt' must be a string. Received %s." % type(wkt))
    code = None  # default
198
199
    if not wkt:
        return None
Daniel Scheffler's avatar
Daniel Scheffler committed
200
    p_in = osr.SpatialReference()
201
    s = p_in.ImportFromWkt(wkt)
Daniel Scheffler's avatar
Daniel Scheffler committed
202
    if s == 5:  # invalid WKT
203
        raise Exception('Received an invalid WKT string: %s' % wkt)
Daniel Scheffler's avatar
Daniel Scheffler committed
204
205
206
    if p_in.IsLocal():
        raise Exception('The given WKT is a local coordinate system.')
    cstype = 'GEOGCS' if p_in.IsGeographic() else 'PROJCS'
207
    p_in.AutoIdentifyEPSG()
Daniel Scheffler's avatar
Daniel Scheffler committed
208
    an = p_in.GetAuthorityName(cstype)
209
    assert an in [None, 'EPSG'], "No EPSG code found. Found %s instead." % an
Daniel Scheffler's avatar
Daniel Scheffler committed
210
211
212
213
    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:
214
            epsgfile = epsgfile or gdal_env.find_epsgfile()
215
            with open(epsgfile) as f:
Daniel Scheffler's avatar
Daniel Scheffler committed
216
217
218
219
220
221
                for line in f:
                    if line.find(p_out) != -1:
                        m = re.search('<(\\d+)>', line)
                        if m:
                            code = m.group(1)
                            break
222
                code = int(code) if code else None  # match or no match
Daniel Scheffler's avatar
Daniel Scheffler committed
223
    else:
Daniel Scheffler's avatar
Daniel Scheffler committed
224
        code = int(ac)
225
    return code
Daniel Scheffler's avatar
Daniel Scheffler committed
226
227


228
def get_UTMzone(ds=None, prj=None):
Daniel Scheffler's avatar
Daniel Scheffler committed
229
    assert ds or prj, "Specify at least one of the arguments 'ds' or 'prj'"
230
    if isProjectedOrGeographic(prj) == 'projected':
Daniel Scheffler's avatar
Daniel Scheffler committed
231
232
233
234
        srs = osr.SpatialReference()
        srs.ImportFromWkt(ds.GetProjection() if ds else prj)
        return srs.GetUTMZone()
    else:
235
236
237
238
        return None


def get_prjLonLat(fmt='wkt'):
Daniel Scheffler's avatar
Daniel Scheffler committed
239
    # type: (str) -> Union[str, dict]
240
241
242
    """Returns standard geographic projection (EPSG 4326) in the WKT or PROJ4 format.
    :param fmt:     <str> target format - 'WKT' or 'PROJ4'
    """
243
    assert re.search('wkt', fmt, re.I) or re.search('Proj4', fmt, re.I), 'unsupported output format'
244
245
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)
246
247
248
249
250
251
    out = srs.ExportToWkt() if re.search('wkt', fmt, re.I) else srs.ExportToProj4()

    if not out:
        raise EnvironmentError(gdal.GetLastErrorMsg())

    return out