projection.py 6.69 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
Daniel Scheffler's avatar
Daniel Scheffler committed
92
93
94
95
96
97
98
99
100
101
102
103
104
    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):
105
106
107
108
    # type: (int) -> str
    if EPSG_code is not None:
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(EPSG_code)
109
110
111
112
113
114
        proj4 = srs.ExportToProj4()

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

        return proj4
115
116
    else:
        return ''
Daniel Scheffler's avatar
Daniel Scheffler committed
117
118
119


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

126
127
128
129
130
131
        if not wkt:
            raise EnvironmentError(gdal.GetLastErrorMsg())

        return wkt
    else:
        return ''
Daniel Scheffler's avatar
Daniel Scheffler committed
132
133


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

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


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


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

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

    return out