projection.py 5.22 KB
Newer Older
Daniel Scheffler's avatar
Daniel Scheffler committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
# -*- coding: utf-8 -*-
__author__ = "Daniel Scheffler"


import re
import os

# 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
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
37
38
39
40
    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
41
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]}


def prj_equal(prj1, prj2):
54
55
56
57
58
59
    #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
60
61
62
63
64
    return get_proj4info(proj=prj1)==get_proj4info(proj=prj2)


def isProjectedOrGeographic(prj):
    """:param prj: accepts EPSG, Proj4 and WKT projections"""
65
    if prj is None: return None
Daniel Scheffler's avatar
Daniel Scheffler committed
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
    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):
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(EPSG_code)
    return srs.ExportToProj4()


def EPSG2WKT(EPSG_code):
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(EPSG_code)
    return srs.ExportToWkt()


90
def WKT2EPSG(wkt, epsg=os.environ['GDAL_DATA'].replace('/gdal', '/proj/epsg')):
Daniel Scheffler's avatar
Daniel Scheffler committed
91
    """ Transform a WKT string to an EPSG code
92
    :param wkt:  WKT definition
Daniel Scheffler's avatar
Daniel Scheffler committed
93
94
95
96
    :param epsg: the proj.4 epsg file (defaults to '/usr/local/share/proj/epsg')
    :returns:    EPSG code
    http://gis.stackexchange.com/questions/20298/is-it-possible-to-get-the-epsg-value-from-an-osr-spatialreference-class-using-th
    """
97
98
    # FIXME this function returns None if datum=NAD27 but works with datum=WGS84, e.g.:
    # FIXME {PROJCS["UTM_Zone_33N",GEOGCS["GCS_North_American_1927",DATUM["D_North_American_1927",SPHEROID["Clarke_1866",6378206.4,294.9786982]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",15.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]}
Daniel Scheffler's avatar
Daniel Scheffler committed
99
100
    code = None #default
    p_in = osr.SpatialReference()
101
    s    = p_in.ImportFromWkt(wkt)
Daniel Scheffler's avatar
Daniel Scheffler committed
102
103
104
105
106
    if s == 5:  # invalid WKT
        raise Exception('Invalid WKT.')
    if p_in.IsLocal():
        raise Exception('The given WKT is a local coordinate system.')
    cstype = 'GEOGCS' if p_in.IsGeographic() else 'PROJCS'
107
    p_in.AutoIdentifyEPSG()
Daniel Scheffler's avatar
Daniel Scheffler committed
108
109
110
111
112
113
114
115
116
117
118
119
120
    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:
            with open(epsg) as f:
                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
121
                code = int(code) if code else None # match or no match
Daniel Scheffler's avatar
Daniel Scheffler committed
122
    else:
Daniel Scheffler's avatar
Daniel Scheffler committed
123
        code = int(ac)
124
    return code
Daniel Scheffler's avatar
Daniel Scheffler committed
125
126
127
128
129
130
131
132
133


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:
134
135
136
137
138
139
140
141
142
143
144
145
        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)
    return srs.ExportToWkt() if re.search('wkt',fmt,re.I) else srs.ExportToProj4()