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

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

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


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


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


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

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

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

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


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

    return srs.IsLocal()


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

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

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


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

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

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


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

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


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


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

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

    return out