projection.py 7.64 KB
Newer Older
Daniel Scheffler's avatar
Daniel Scheffler committed
1
2
3
4
5
# -*- coding: utf-8 -*-
__author__ = "Daniel Scheffler"

import re
import os
Daniel Scheffler's avatar
Daniel Scheffler committed
6
import pyproj
7
from typing import Union
Daniel Scheffler's avatar
Daniel Scheffler committed
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

# 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 *


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):
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):
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
109
110
111
    # type: (int) -> str
    if EPSG_code is not None:
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(EPSG_code)
        return srs.ExportToProj4()
    else:
        return ''
Daniel Scheffler's avatar
Daniel Scheffler committed
112
113
114


def EPSG2WKT(EPSG_code):
115
116
    # type: (int) -> str
    if EPSG_code is not None:
117
118
119
120
        # GDAL_DATA_tmp = None
        # if 'GDAL_DATA' in os.environ:
        #     GDAL_DATA_tmp = os.environ['GDAL_DATA']
        #     del os.environ['GDAL_DATA']
121
122
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(EPSG_code)
123
124
125
126
127
128
129
130
131
        wkt = srs.ExportToWkt()
        # if GDAL_DATA_tmp:
        #     os.environ['GDAL_DATA'] = GDAL_DATA_tmp
        if not wkt:
            raise EnvironmentError(gdal.GetLastErrorMsg())

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


134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
def _find_epsgfile():
    """Locate the proj.4 epsg file (defaults to '/usr/local/share/proj/epsg')."""
    try:
        epsgfile = os.environ['GDAL_DATA'].replace('/gdal', '/proj/epsg')
        assert os.path.exists(epsgfile)
    except (KeyError, AssertionError):
        try:
            from pyproj import __file__ as pyprojpath
            epsgfile = os.path.join(os.path.dirname(pyprojpath), 'data/epsg')
            assert os.path.exists(epsgfile)
        except (ImportError, AssertionError):
            epsgfile = '/usr/local/share/proj/epsg'
            if not os.path.exists(epsgfile):
                raise RuntimeError('Could not locate epsg file for converting WKT to EPSG code. '
                                   'Please make sure that your GDAL_DATA environment variable is properly set and the '
                                   'pyproj library is installed.')
    return epsgfile


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

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


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


def get_prjLonLat(fmt='wkt'):
212
    # type: (str) -> Union[str, dict]
213
214
215
    """Returns standard geographic projection (EPSG 4326) in the WKT or PROJ4 format.
    :param fmt:     <str> target format - 'WKT' or 'PROJ4'
    """
216
    assert re.search('wkt', fmt, re.I) or re.search('Proj4', fmt, re.I), 'unsupported output format'
217
218
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)
219
    return srs.ExportToWkt() if re.search('wkt', fmt, re.I) else srs.ExportToProj4()