projection.py 6.82 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
Daniel Scheffler's avatar
Daniel Scheffler committed
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

# custom
try:
    from osgeo import ogr
    from osgeo import osr
    from osgeo import gdal
    from osgeo import gdalnumeric
    from osgeo.gdalconst import *
except ImportError:
    import ogr
    import osr
    import gdal
    import gdalnumeric
    from gdalconst import *

21
22
23
24
25
26
27
from ..environment import gdal_env

__author__ = "Daniel Scheffler"

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

Daniel Scheffler's avatar
Daniel Scheffler committed
28

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

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

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


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


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


def isProjectedOrGeographic(prj):
90
    # type: (Union[str, int, dict]) -> Union[str, None]
91
92
    """

93
    :param prj: accepts EPSG, Proj4 and WKT projections
94
    """
95
96
    if prj is None:
        return None
Daniel Scheffler's avatar
Daniel Scheffler committed
97
98
99
100
101
102
103
104
105
106
107
108
109
    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):
110
111
112
113
    # type: (int) -> str
    if EPSG_code is not None:
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(EPSG_code)
114
115
116
117
118
119
        proj4 = srs.ExportToProj4()

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

        return proj4
120
121
    else:
        return ''
Daniel Scheffler's avatar
Daniel Scheffler committed
122
123
124


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

131
132
133
134
135
136
        if not wkt:
            raise EnvironmentError(gdal.GetLastErrorMsg())

        return wkt
    else:
        return ''
Daniel Scheffler's avatar
Daniel Scheffler committed
137
138


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

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


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


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

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

    return out