projection.py 5.17 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

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

34
35
36
37
from ..environment import gdal_env

__author__ = "Daniel Scheffler"

38

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

Daniel Scheffler's avatar
Daniel Scheffler committed
42

43
44
45
46
def get_proj4info(proj=None):
    # type: (Union[str, int]) -> str
    """Returns PROJ4 formatted projection info for the given projection.

Daniel Scheffler's avatar
Daniel Scheffler committed
47
    e.g. '+proj=utm +zone=43 +datum=WGS84 +units=m +no_defs '
48

49
    :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
50
    """
51
    return CRS.from_user_input(proj).to_proj4().strip()
Daniel Scheffler's avatar
Daniel Scheffler committed
52
53
54


def proj4_to_dict(proj4):
55
    # type: (str) -> dict
Daniel Scheffler's avatar
Daniel Scheffler committed
56
57
58
    """Converts a PROJ4-like string into a dictionary.
    :param proj4:   <str> the PROJ4-like string
    """
59
    return CRS.from_proj4(proj4).to_dict()
Daniel Scheffler's avatar
Daniel Scheffler committed
60
61


Daniel Scheffler's avatar
Daniel Scheffler committed
62
63
64
65
66
def dict_to_proj4(proj4dict):
    # type: (dict) -> str
    """Converts a PROJ4-like dictionary into a PROJ4 string.
    :param proj4dict:   <dict> the PROJ4-like dictionary
    """
67
    return CRS.from_dict(proj4dict).to_proj4()
Daniel Scheffler's avatar
Daniel Scheffler committed
68
69


Daniel Scheffler's avatar
Daniel Scheffler committed
70
71
72
73
74
def proj4_to_WKT(proj4str):
    # type: (str) -> str
    """Converts a PROJ4-like string into a WKT string.
    :param proj4str:   <dict> the PROJ4-like string
    """
75
    return CRS.from_proj4(proj4str).to_wkt()
Daniel Scheffler's avatar
Daniel Scheffler committed
76
77


Daniel Scheffler's avatar
Daniel Scheffler committed
78
def prj_equal(prj1, prj2):
Daniel Scheffler's avatar
Daniel Scheffler committed
79
    # type: (Union[None, int, str], Union[None, int, str]) -> bool
80
81
82
83
84
    """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>)
    """
85
86
87
    if prj1 is None and prj2 is None or prj1 == prj2:
        return True
    else:
88
89
90
91
        crs1 = CRS.from_user_input(prj1)
        crs2 = CRS.from_user_input(prj2)

        return crs1.equals(crs2)
Daniel Scheffler's avatar
Daniel Scheffler committed
92
93
94


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

98
    :param prj: accepts EPSG, Proj4 and WKT projections
99
    """
100
101
    if prj is None:
        return None
102

103
    crs = CRS.from_user_input(prj)
104

105
    return 'projected' if crs.is_projected else 'geographic' if crs.is_geographic else None
Daniel Scheffler's avatar
Daniel Scheffler committed
106
107


108
109
110
111
112
113
114
115
116
117
118
119
120
121
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)
Daniel Scheffler's avatar
Daniel Scheffler committed
122
123
124
125
126
    elif 'GEOGCS' in prj or \
         'GEOGCRS' in prj or \
         'PROJCS' in prj or \
         'PROJCS' in prj or \
         'LOCAL_CS' in prj:
127
128
        srs.ImportFromWkt(prj)
    else:
129
        raise RuntimeError('Unknown input projection: \n%s' % prj)
130
131
132
133

    return srs.IsLocal()


Daniel Scheffler's avatar
Daniel Scheffler committed
134
def EPSG2Proj4(EPSG_code):
135
    # type: (int) -> str
136
    return CRS.from_epsg(EPSG_code).to_proj4() if EPSG_code is not None else ''
Daniel Scheffler's avatar
Daniel Scheffler committed
137
138
139


def EPSG2WKT(EPSG_code):
140
    # type: (int) -> str
141
    return CRS.from_epsg(EPSG_code).to_wkt(pretty=False) if EPSG_code is not None else ''
Daniel Scheffler's avatar
Daniel Scheffler committed
142
143


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

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


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


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

176
177
178
179
    if re.search('wkt', fmt, re.I):
        return CRS.from_epsg(4326).to_wkt()
    else:
        return CRS.from_epsg(4326).to_proj4()