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

import re
Daniel Scheffler's avatar
Daniel Scheffler committed
4
import pyproj
Daniel Scheffler's avatar
Daniel Scheffler committed
5
from typing import Union  # noqa F401  # flake8 issue
Daniel Scheffler's avatar
Daniel Scheffler committed
6
7
8
9
10
11
12
13
14

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

15
16
17
18
from ..environment import gdal_env

__author__ = "Daniel Scheffler"

19

20
21
22
# try to set GDAL_DATA if unnot set or invalid
gdal_env.try2set_GDAL_DATA()

Daniel Scheffler's avatar
Daniel Scheffler committed
23

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

Daniel Scheffler's avatar
Daniel Scheffler committed
29
    :param ds:      <gdal.Dataset> the gdal dataset to get PROJ4 info for
30
    :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
31
32
33
    """
    assert ds or proj, "Specify at least one of the arguments 'ds' or 'proj'"
    srs = osr.SpatialReference()
34
35
36
37
    proj = ds.GetProjection() if ds else proj

    if isinstance(proj, str) and proj.startswith('epsg:'):
        proj = EPSG2WKT(int(proj.split(':')[1]))
38
    elif isinstance(proj, int):
39
40
41
        proj = EPSG2WKT(proj)

    srs.ImportFromWkt(proj)
Daniel Scheffler's avatar
Daniel Scheffler committed
42
43
44
45
    return srs.ExportToProj4()


def proj4_to_dict(proj4):
46
    # type: (str) -> dict
Daniel Scheffler's avatar
Daniel Scheffler committed
47
48
49
    """Converts a PROJ4-like string into a dictionary.
    :param proj4:   <str> the PROJ4-like string
    """
50
51
    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
52
53


Daniel Scheffler's avatar
Daniel Scheffler committed
54
55
56
57
58
59
60
61
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
62
63
64
65
66
67
68
69
70
71
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
72
def prj_equal(prj1, prj2):
Daniel Scheffler's avatar
Daniel Scheffler committed
73
    # type: (Union[None, int, str], Union[None, int, str]) -> bool
74
75
76
77
78
    """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>)
    """
79
80
81
    if prj1 is None and prj2 is None or prj1 == prj2:
        return True
    else:
82
        return get_proj4info(proj=prj1) == get_proj4info(proj=prj2)
Daniel Scheffler's avatar
Daniel Scheffler committed
83
84
85


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

89
    :param prj: accepts EPSG, Proj4 and WKT projections
90
    """
91
92
    if prj is None:
        return None
93

Daniel Scheffler's avatar
Daniel Scheffler committed
94
95
96
97
98
99
100
101
    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:
102
103
        raise RuntimeError('Unknown input projection.')

Daniel Scheffler's avatar
Daniel Scheffler committed
104
105
106
    return 'projected' if srs.IsProjected() else 'geographic' if srs.IsGeographic() else None


107
108
109
110
111
112
113
114
115
116
117
118
119
120
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)
121
    elif 'GEOGCS' in prj or 'PROJCS' in prj or 'LOCAL_CS' in prj:
122
123
        srs.ImportFromWkt(prj)
    else:
124
        raise RuntimeError('Unknown input projection: \n%s' % prj)
125
126
127
128

    return srs.IsLocal()


Daniel Scheffler's avatar
Daniel Scheffler committed
129
def EPSG2Proj4(EPSG_code):
130
131
132
133
    # type: (int) -> str
    if EPSG_code is not None:
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(EPSG_code)
134
135
136
137
138
139
        proj4 = srs.ExportToProj4()

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

        return proj4
140
141
    else:
        return ''
Daniel Scheffler's avatar
Daniel Scheffler committed
142
143
144


def EPSG2WKT(EPSG_code):
145
146
147
148
    # type: (int) -> str
    if EPSG_code is not None:
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(EPSG_code)
149
        wkt = srs.ExportToWkt()
150

151
152
153
154
155
156
        if not wkt:
            raise EnvironmentError(gdal.GetLastErrorMsg())

        return wkt
    else:
        return ''
Daniel Scheffler's avatar
Daniel Scheffler committed
157
158


159
def WKT2EPSG(wkt, epsgfile=''):
Daniel Scheffler's avatar
Daniel Scheffler committed
160
    # type: (str) -> Union[int, None]
Daniel Scheffler's avatar
Daniel Scheffler committed
161
    """ Transform a WKT string to an EPSG code
162
    :param wkt:  WKT definition
163
    :param epsgfile: the proj.4 epsg file (automatically located if no path is provided)
Daniel Scheffler's avatar
Daniel Scheffler committed
164
165
166
    :returns:    EPSG code
    http://gis.stackexchange.com/questions/20298/is-it-possible-to-get-the-epsg-value-from-an-osr-spatialreference-class-using-th
    """
167
    # FIXME this function returns None if datum=NAD27 but works with datum=WGS84, e.g.:
168
169
170
171
172
    # 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
173

174
175
176
    if not isinstance(wkt, str):
        raise TypeError("'wkt' must be a string. Received %s." % type(wkt))
    code = None  # default
177
178
    if not wkt:
        return None
Daniel Scheffler's avatar
Daniel Scheffler committed
179
    p_in = osr.SpatialReference()
180
    s = p_in.ImportFromWkt(wkt)
Daniel Scheffler's avatar
Daniel Scheffler committed
181
    if s == 5:  # invalid WKT
182
        raise Exception('Received an invalid WKT string: %s' % wkt)
Daniel Scheffler's avatar
Daniel Scheffler committed
183
184
185
    if p_in.IsLocal():
        raise Exception('The given WKT is a local coordinate system.')
    cstype = 'GEOGCS' if p_in.IsGeographic() else 'PROJCS'
186
    p_in.AutoIdentifyEPSG()
Daniel Scheffler's avatar
Daniel Scheffler committed
187
    an = p_in.GetAuthorityName(cstype)
188
    assert an in [None, 'EPSG'], "No EPSG code found. Found %s instead." % an
Daniel Scheffler's avatar
Daniel Scheffler committed
189
190
191
192
    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:
193
            epsgfile = epsgfile or gdal_env.find_epsgfile()
194
            with open(epsgfile) as f:
Daniel Scheffler's avatar
Daniel Scheffler committed
195
196
197
198
199
200
                for line in f:
                    if line.find(p_out) != -1:
                        m = re.search('<(\\d+)>', line)
                        if m:
                            code = m.group(1)
                            break
201
                code = int(code) if code else None  # match or no match
Daniel Scheffler's avatar
Daniel Scheffler committed
202
    else:
Daniel Scheffler's avatar
Daniel Scheffler committed
203
        code = int(ac)
204
    return code
Daniel Scheffler's avatar
Daniel Scheffler committed
205
206


207
def get_UTMzone(ds=None, prj=None):
Daniel Scheffler's avatar
Daniel Scheffler committed
208
    assert ds or prj, "Specify at least one of the arguments 'ds' or 'prj'"
209
    if isProjectedOrGeographic(prj) == 'projected':
Daniel Scheffler's avatar
Daniel Scheffler committed
210
211
212
213
        srs = osr.SpatialReference()
        srs.ImportFromWkt(ds.GetProjection() if ds else prj)
        return srs.GetUTMZone()
    else:
214
215
216
217
        return None


def get_prjLonLat(fmt='wkt'):
Daniel Scheffler's avatar
Daniel Scheffler committed
218
    # type: (str) -> Union[str, dict]
219
220
221
    """Returns standard geographic projection (EPSG 4326) in the WKT or PROJ4 format.
    :param fmt:     <str> target format - 'WKT' or 'PROJ4'
    """
222
    assert re.search('wkt', fmt, re.I) or re.search('Proj4', fmt, re.I), 'unsupported output format'
223
224
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)
225
226
227
228
229
230
    out = srs.ExportToWkt() if re.search('wkt', fmt, re.I) else srs.ExportToProj4()

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

    return out