projection.py 5.54 KB
Newer Older
Daniel Scheffler's avatar
Daniel Scheffler committed
1
2
# -*- coding: utf-8 -*-

3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# py_tools_ds
#
# Copyright (C) 2019  Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de)
#
# This software was developed within the context of the GeoMultiSens project funded
# by the German Federal Ministry of Education and Research
# (project grant code: 01 IS 14 010 A-C).
#
# This program is free software: you can redistribute it and/or modify it under
# the terms of the GNU Lesser General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
# details.
#
# You should have received a copy of the GNU Lesser General Public License along
# with this program.  If not, see <http://www.gnu.org/licenses/>.

Daniel Scheffler's avatar
Daniel Scheffler committed
24
import re
25
from pyproj import CRS
Daniel Scheffler's avatar
Daniel Scheffler committed
26
from typing import Union  # noqa F401  # flake8 issue
Daniel Scheffler's avatar
Daniel Scheffler committed
27
28
29
30
31
32
33
34
35

# custom
try:
    from osgeo import osr
    from osgeo import gdal
except ImportError:
    import osr
    import gdal

36
37
38
39
from ..environment import gdal_env

__author__ = "Daniel Scheffler"

40

Daniel Scheffler's avatar
Daniel Scheffler committed
41
# try to set GDAL_DATA if not set or invalid
42
43
gdal_env.try2set_GDAL_DATA()

Daniel Scheffler's avatar
Daniel Scheffler committed
44

45
def get_proj4info(ds=None, proj=None):
46
    # type: (gdal.Dataset, Union[str, int]) -> str
47
    """Returns PROJ4 formatted projection info for the given gdal.Dataset or projection respectively,
Daniel Scheffler's avatar
Daniel Scheffler committed
48
    e.g. '+proj=utm +zone=43 +datum=WGS84 +units=m +no_defs '
49

Daniel Scheffler's avatar
Daniel Scheffler committed
50
    :param ds:      <gdal.Dataset> the gdal dataset to get PROJ4 info for
51
    :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
52
53
54
    """
    assert ds or proj, "Specify at least one of the arguments 'ds' or 'proj'"
    srs = osr.SpatialReference()
55
56
57
58
    proj = ds.GetProjection() if ds else proj

    if isinstance(proj, str) and proj.startswith('epsg:'):
        proj = EPSG2WKT(int(proj.split(':')[1]))
59
    elif isinstance(proj, int):
60
61
62
        proj = EPSG2WKT(proj)

    srs.ImportFromWkt(proj)
63
    return srs.ExportToProj4().strip()
Daniel Scheffler's avatar
Daniel Scheffler committed
64
65
66


def proj4_to_dict(proj4):
67
    # type: (str) -> dict
Daniel Scheffler's avatar
Daniel Scheffler committed
68
69
70
    """Converts a PROJ4-like string into a dictionary.
    :param proj4:   <str> the PROJ4-like string
    """
71
    return CRS.from_proj4(proj4).to_dict()
Daniel Scheffler's avatar
Daniel Scheffler committed
72
73


Daniel Scheffler's avatar
Daniel Scheffler committed
74
75
76
77
78
def dict_to_proj4(proj4dict):
    # type: (dict) -> str
    """Converts a PROJ4-like dictionary into a PROJ4 string.
    :param proj4dict:   <dict> the PROJ4-like dictionary
    """
79
    return CRS.from_dict(proj4dict).to_proj4()
Daniel Scheffler's avatar
Daniel Scheffler committed
80
81


Daniel Scheffler's avatar
Daniel Scheffler committed
82
83
84
85
86
def proj4_to_WKT(proj4str):
    # type: (str) -> str
    """Converts a PROJ4-like string into a WKT string.
    :param proj4str:   <dict> the PROJ4-like string
    """
87
    return CRS.from_proj4(proj4str).to_wkt()
Daniel Scheffler's avatar
Daniel Scheffler committed
88
89


Daniel Scheffler's avatar
Daniel Scheffler committed
90
def prj_equal(prj1, prj2):
Daniel Scheffler's avatar
Daniel Scheffler committed
91
    # type: (Union[None, int, str], Union[None, int, str]) -> bool
92
93
94
95
96
    """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>)
    """
97
98
99
    if prj1 is None and prj2 is None or prj1 == prj2:
        return True
    else:
100
        return get_proj4info(proj=prj1) == get_proj4info(proj=prj2)
Daniel Scheffler's avatar
Daniel Scheffler committed
101
102
103


def isProjectedOrGeographic(prj):
Daniel Scheffler's avatar
Daniel Scheffler committed
104
    # type: (Union[str, int, dict]) -> Union[str, None]
105
106
    """

107
    :param prj: accepts EPSG, Proj4 and WKT projections
108
    """
109
110
    if prj is None:
        return None
111

112
    crs = CRS.from_user_input(prj)
113

114
    return 'projected' if crs.is_projected else 'geographic' if crs.is_geographic else None
Daniel Scheffler's avatar
Daniel Scheffler committed
115
116


117
118
119
120
121
122
123
124
125
126
127
128
129
130
def isLocal(prj):
    # type: (Union[str, int, dict]) -> Union[bool, None]
    """

    :param prj: accepts EPSG, Proj4 and WKT projections
    """
    if not prj:
        return True

    srs = osr.SpatialReference()
    if prj.startswith('EPSG:'):
        srs.ImportFromEPSG(int(prj.split(':')[1]))
    elif prj.startswith('+proj='):
        srs.ImportFromProj4(prj)
131
    elif 'GEOGCS' in prj or 'PROJCS' in prj or 'LOCAL_CS' in prj:
132
133
        srs.ImportFromWkt(prj)
    else:
134
        raise RuntimeError('Unknown input projection: \n%s' % prj)
135
136
137
138

    return srs.IsLocal()


Daniel Scheffler's avatar
Daniel Scheffler committed
139
def EPSG2Proj4(EPSG_code):
140
    # type: (int) -> str
141
    return CRS.from_epsg(EPSG_code).to_proj4() if EPSG_code is not None else ''
Daniel Scheffler's avatar
Daniel Scheffler committed
142
143
144


def EPSG2WKT(EPSG_code):
145
    # type: (int) -> str
146
    return CRS.from_epsg(EPSG_code).to_wkt(pretty=False) if EPSG_code is not None else ''
Daniel Scheffler's avatar
Daniel Scheffler committed
147
148


149
150
def WKT2EPSG(wkt):
    # type: (str) -> Union[int, None]
Daniel Scheffler's avatar
Daniel Scheffler committed
151
    """ Transform a WKT string to an EPSG code
152
    :param wkt:  WKT definition
Daniel Scheffler's avatar
Daniel Scheffler committed
153
154
    :returns:    EPSG code
    """
155
156
    if not isinstance(wkt, str):
        raise TypeError("'wkt' must be a string. Received %s." % type(wkt))
157
158
    if not wkt:
        return None
159
160

    return CRS.from_wkt(wkt.replace('\n', '').replace('\r', '').replace(' ', '')).to_epsg()
Daniel Scheffler's avatar
Daniel Scheffler committed
161
162


163
def get_UTMzone(ds=None, prj=None):
Daniel Scheffler's avatar
Daniel Scheffler committed
164
    assert ds or prj, "Specify at least one of the arguments 'ds' or 'prj'"
165
    if isProjectedOrGeographic(prj) == 'projected':
Daniel Scheffler's avatar
Daniel Scheffler committed
166
167
168
169
        srs = osr.SpatialReference()
        srs.ImportFromWkt(ds.GetProjection() if ds else prj)
        return srs.GetUTMZone()
    else:
170
171
172
173
        return None


def get_prjLonLat(fmt='wkt'):
Daniel Scheffler's avatar
Daniel Scheffler committed
174
    # type: (str) -> Union[str, dict]
175
176
177
    """Returns standard geographic projection (EPSG 4326) in the WKT or PROJ4 format.
    :param fmt:     <str> target format - 'WKT' or 'PROJ4'
    """
178
179
    if not re.search('wkt', fmt, re.I) or re.search('Proj4', fmt, re.I):
        raise ValueError(fmt, 'Unsupported output format.')
180

181
182
183
184
    if re.search('wkt', fmt, re.I):
        return CRS.from_epsg(4326).to_wkt()
    else:
        return CRS.from_epsg(4326).to_proj4()