projection.py 6.77 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
5
from typing import Union, TypeVar, TYPE_CHECKING
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
if TYPE_CHECKING:
    _T_prj1 = TypeVar(Union[None, int, str])
    _T_prj2 = TypeVar(Union[str, int, dict])

23
24
25
# try to set GDAL_DATA if unnot set or invalid
gdal_env.try2set_GDAL_DATA()

Daniel Scheffler's avatar
Daniel Scheffler committed
26

27
def get_proj4info(ds=None, proj=None):
Daniel Scheffler's avatar
Daniel Scheffler committed
28
29
30
31
    # 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
32
    :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
33
34
35
    """
    assert ds or proj, "Specify at least one of the arguments 'ds' or 'proj'"
    srs = osr.SpatialReference()
36
37
38
39
    proj = ds.GetProjection() if ds else proj

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

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


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


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


def isProjectedOrGeographic(prj):
88
    # type: (_T_prj2) -> TypeVar(Union[str, None])
89
90
    """

91
    :param prj: accepts EPSG, Proj4 and WKT projections
92
    """
93
94
    if prj is None:
        return None
Daniel Scheffler's avatar
Daniel Scheffler committed
95
96
97
98
99
100
101
102
103
104
105
106
107
    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 Exception('Unknown input projection.')
    return 'projected' if srs.IsProjected() else 'geographic' if srs.IsGeographic() else None


def EPSG2Proj4(EPSG_code):
108
109
110
111
    # type: (int) -> str
    if EPSG_code is not None:
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(EPSG_code)
112
113
114
115
116
117
        proj4 = srs.ExportToProj4()

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

        return proj4
118
119
    else:
        return ''
Daniel Scheffler's avatar
Daniel Scheffler committed
120
121
122


def EPSG2WKT(EPSG_code):
123
124
125
126
    # type: (int) -> str
    if EPSG_code is not None:
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(EPSG_code)
127
        wkt = srs.ExportToWkt()
128

129
130
131
132
133
134
        if not wkt:
            raise EnvironmentError(gdal.GetLastErrorMsg())

        return wkt
    else:
        return ''
Daniel Scheffler's avatar
Daniel Scheffler committed
135
136


137
def WKT2EPSG(wkt, epsgfile=''):
138
    # type: (str) -> TypeVar(Union[int, None])
Daniel Scheffler's avatar
Daniel Scheffler committed
139
    """ Transform a WKT string to an EPSG code
140
    :param wkt:  WKT definition
141
    :param epsgfile: the proj.4 epsg file (automatically located if no path is provided)
Daniel Scheffler's avatar
Daniel Scheffler committed
142
143
144
    :returns:    EPSG code
    http://gis.stackexchange.com/questions/20298/is-it-possible-to-get-the-epsg-value-from-an-osr-spatialreference-class-using-th
    """
145
    # FIXME this function returns None if datum=NAD27 but works with datum=WGS84, e.g.:
146
147
148
149
150
    # 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
151

152
153
154
    if not isinstance(wkt, str):
        raise TypeError("'wkt' must be a string. Received %s." % type(wkt))
    code = None  # default
155
156
    if not wkt:
        return None
Daniel Scheffler's avatar
Daniel Scheffler committed
157
    p_in = osr.SpatialReference()
158
    s = p_in.ImportFromWkt(wkt)
Daniel Scheffler's avatar
Daniel Scheffler committed
159
    if s == 5:  # invalid WKT
160
        raise Exception('Received an invalid WKT string: %s' % wkt)
Daniel Scheffler's avatar
Daniel Scheffler committed
161
162
163
    if p_in.IsLocal():
        raise Exception('The given WKT is a local coordinate system.')
    cstype = 'GEOGCS' if p_in.IsGeographic() else 'PROJCS'
164
    p_in.AutoIdentifyEPSG()
Daniel Scheffler's avatar
Daniel Scheffler committed
165
    an = p_in.GetAuthorityName(cstype)
166
    assert an in [None, 'EPSG'], "No EPSG code found. Found %s instead." % an
Daniel Scheffler's avatar
Daniel Scheffler committed
167
168
169
170
    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:
171
            epsgfile = epsgfile or gdal_env.find_epsgfile()
172
            with open(epsgfile) as f:
Daniel Scheffler's avatar
Daniel Scheffler committed
173
174
175
176
177
178
                for line in f:
                    if line.find(p_out) != -1:
                        m = re.search('<(\\d+)>', line)
                        if m:
                            code = m.group(1)
                            break
179
                code = int(code) if code else None  # match or no match
Daniel Scheffler's avatar
Daniel Scheffler committed
180
    else:
Daniel Scheffler's avatar
Daniel Scheffler committed
181
        code = int(ac)
182
    return code
Daniel Scheffler's avatar
Daniel Scheffler committed
183
184


185
def get_UTMzone(ds=None, prj=None):
Daniel Scheffler's avatar
Daniel Scheffler committed
186
    assert ds or prj, "Specify at least one of the arguments 'ds' or 'prj'"
187
    if isProjectedOrGeographic(prj) == 'projected':
Daniel Scheffler's avatar
Daniel Scheffler committed
188
189
190
191
        srs = osr.SpatialReference()
        srs.ImportFromWkt(ds.GetProjection() if ds else prj)
        return srs.GetUTMZone()
    else:
192
193
194
195
        return None


def get_prjLonLat(fmt='wkt'):
196
    # type: (str) -> TypeVar(Union[str, dict])
197
198
199
    """Returns standard geographic projection (EPSG 4326) in the WKT or PROJ4 format.
    :param fmt:     <str> target format - 'WKT' or 'PROJ4'
    """
200
    assert re.search('wkt', fmt, re.I) or re.search('Proj4', fmt, re.I), 'unsupported output format'
201
202
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)
203
204
205
206
207
208
    out = srs.ExportToWkt() if re.search('wkt', fmt, re.I) else srs.ExportToProj4()

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

    return out