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


import re
import os
Daniel Scheffler's avatar
Daniel Scheffler committed
7
import pyproj
Daniel Scheffler's avatar
Daniel Scheffler committed
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

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



def get_proj4info(ds=None,proj=None):
    # 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
30
    :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
31
32
33
    """
    assert ds or proj, "Specify at least one of the arguments 'ds' or 'proj'"
    srs = osr.SpatialReference()
34
35
36
37
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]))
    elif isinstance(proj,int):
        proj = EPSG2WKT(proj)

    srs.ImportFromWkt(proj)
Daniel Scheffler's avatar
Daniel Scheffler committed
42
43
44
45
46
47
48
49
50
51
52
53
    return srs.ExportToProj4()


def proj4_to_dict(proj4):
    #type: (str) -> dict
    """Converts a PROJ4-like string into a dictionary.
    :param proj4:   <str> the PROJ4-like string
    """
    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
56
57
58
59
60
61
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
62
63
64
65
66
67
68
69
70
71
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
72
def prj_equal(prj1, prj2):
73
74
75
76
77
78
    #type: (any,any) -> bool
    """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>)
    """
Daniel Scheffler's avatar
Daniel Scheffler committed
79
80
81
82
83
    return get_proj4info(proj=prj1)==get_proj4info(proj=prj2)


def isProjectedOrGeographic(prj):
    """:param prj: accepts EPSG, Proj4 and WKT projections"""
84
    if prj is None: return None
Daniel Scheffler's avatar
Daniel Scheffler committed
85
86
87
88
89
90
91
92
93
94
95
96
97
    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):
98
99
100
101
102
103
104
    # 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
105
106
107


def EPSG2WKT(EPSG_code):
108
109
110
111
112
    # type: (int) -> str
    if EPSG_code is not None:
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(EPSG_code)
        return srs.ExportToWkt()
Daniel Scheffler's avatar
Daniel Scheffler committed
113
114


115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
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=''):
Daniel Scheffler's avatar
Daniel Scheffler committed
135
    """ Transform a WKT string to an EPSG code
136
    :param wkt:  WKT definition
137
    :param epsgfile: the proj.4 epsg file (automatically located if no path is provided)
Daniel Scheffler's avatar
Daniel Scheffler committed
138
139
140
    :returns:    EPSG code
    http://gis.stackexchange.com/questions/20298/is-it-possible-to-get-the-epsg-value-from-an-osr-spatialreference-class-using-th
    """
141
    # FIXME this function returns None if datum=NAD27 but works with datum=WGS84, e.g.:
142
143
144
145
146
    # 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
147
148
149

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


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


def get_prjLonLat(fmt='wkt'):
    # type: (str) -> Any
    """Returns standard geographic projection (EPSG 4326) in the WKT or PROJ4 format.
    :param fmt:     <str> target format - 'WKT' or 'PROJ4'
    """
    assert re.search('wkt',fmt,re.I) or re.search('Proj4',fmt,re.I), 'unsupported output format'
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)
199
    return srs.ExportToWkt() if re.search('wkt',fmt,re.I) else srs.ExportToProj4()