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

3
4
# py_tools_ds - A collection of geospatial data analysis tools that simplify standard
# operations when handling geospatial raster and vector data as well as projections.
5
#
6
7
8
9
# Copyright (C) 2016-2021
# - Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de)
# - Helmholtz Centre Potsdam - GFZ German Research Centre for Geosciences Potsdam,
#   Germany (https://www.gfz-potsdam.de/)
10
11
12
13
14
#
# 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).
#
15
16
17
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
18
#
19
#   http://www.apache.org/licenses/LICENSE-2.0
20
#
21
22
23
24
25
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
26

Daniel Scheffler's avatar
Daniel Scheffler committed
27
import re
28
import warnings
29
from pyproj import CRS
Daniel Scheffler's avatar
Daniel Scheffler committed
30
from typing import Union  # noqa F401  # flake8 issue
31
from osgeo import osr
Daniel Scheffler's avatar
Daniel Scheffler committed
32

33
34
35
36
from ..environment import gdal_env

__author__ = "Daniel Scheffler"

37

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

Daniel Scheffler's avatar
Daniel Scheffler committed
41

42
43
def get_proj4info(proj=None):
    # type: (Union[str, int]) -> str
Daniel Scheffler's avatar
Daniel Scheffler committed
44
    """Return PROJ4 formatted projection info for the given projection.
45

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

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


def proj4_to_dict(proj4):
54
    # type: (str) -> dict
Daniel Scheffler's avatar
Daniel Scheffler committed
55
56
    """Convert a PROJ4-like string into a dictionary.

Daniel Scheffler's avatar
Daniel Scheffler committed
57
58
    :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
def dict_to_proj4(proj4dict):
    # type: (dict) -> str
Daniel Scheffler's avatar
Daniel Scheffler committed
64
65
    """Convert a PROJ4-like dictionary into a PROJ4 string.

Daniel Scheffler's avatar
Daniel Scheffler committed
66
67
    :param proj4dict:   <dict> the PROJ4-like dictionary
    """
68
    return CRS.from_dict(proj4dict).to_proj4()
Daniel Scheffler's avatar
Daniel Scheffler committed
69
70


Daniel Scheffler's avatar
Daniel Scheffler committed
71
72
def proj4_to_WKT(proj4str):
    # type: (str) -> str
Daniel Scheffler's avatar
Daniel Scheffler committed
73
74
    """Convert a PROJ4-like string into a WKT string.

Daniel Scheffler's avatar
Daniel Scheffler committed
75
76
    :param proj4str:   <dict> the PROJ4-like string
    """
77
    return CRS.from_proj4(proj4str).to_wkt()
Daniel Scheffler's avatar
Daniel Scheffler committed
78
79


Daniel Scheffler's avatar
Daniel Scheffler committed
80
def prj_equal(prj1, prj2):
Daniel Scheffler's avatar
Daniel Scheffler committed
81
    # type: (Union[None, int, str], Union[None, int, str]) -> bool
Daniel Scheffler's avatar
Daniel Scheffler committed
82
    """Check if the given two projections are equal.
83
84
85
86

    :param prj1: projection 1 (WKT or 'epsg:1234' or <EPSG_int>)
    :param prj2: projection 2 (WKT or 'epsg:1234' or <EPSG_int>)
    """
87
88
89
    if prj1 is None and prj2 is None or prj1 == prj2:
        return True
    else:
90
91
92
        # CRS.equals was added in pyproj 2.5
        crs1 = CRS.from_user_input(prj1)
        crs2 = CRS.from_user_input(prj2)
93

94
        return crs1.equals(crs2)
Daniel Scheffler's avatar
Daniel Scheffler committed
95
96
97


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

101
    :param prj: accepts EPSG, Proj4 and WKT projections
102
    """
103
104
    if prj is None:
        return None
105

106
    crs = CRS.from_user_input(prj)
107

108
    return 'projected' if crs.is_projected else 'geographic' if crs.is_geographic else None
Daniel Scheffler's avatar
Daniel Scheffler committed
109
110


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

    return srs.IsLocal()


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


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


147
148
def WKT2EPSG(wkt):
    # type: (str) -> Union[int, None]
Daniel Scheffler's avatar
Daniel Scheffler committed
149
150
    """Transform a WKT string to an EPSG code.

151
    :param wkt:  WKT definition
Daniel Scheffler's avatar
Daniel Scheffler committed
152
153
    :returns:    EPSG code
    """
154
155
    if not isinstance(wkt, str):
        raise TypeError("'wkt' must be a string. Received %s." % type(wkt))
156
157
    if not wkt:
        return None
158

159
160
161
162
163
164
165
166
167
168
169
    ccrs = CRS.from_wkt(wkt.replace('\n', '').replace('\r', '').replace(' ', ''))\

    if not ccrs.is_bound:
        epsg = ccrs.to_epsg()
    else:
        epsg = ccrs.source_crs.to_epsg()

    if epsg is None:
        warnings.warn('Could not find a suitable EPSG code for the input WKT string.', RuntimeWarning)
    else:
        return epsg
Daniel Scheffler's avatar
Daniel Scheffler committed
170
171


172
173
def get_UTMzone(prj):
    # type: (str) -> Union[int, None]
174
    if isProjectedOrGeographic(prj) == 'projected':
Daniel Scheffler's avatar
Daniel Scheffler committed
175
        srs = osr.SpatialReference()
176
        srs.ImportFromWkt(prj)
Daniel Scheffler's avatar
Daniel Scheffler committed
177
178
        return srs.GetUTMZone()
    else:
179
180
181
182
        return None


def get_prjLonLat(fmt='wkt'):
Daniel Scheffler's avatar
Daniel Scheffler committed
183
    # type: (str) -> Union[str, dict]
Daniel Scheffler's avatar
Daniel Scheffler committed
184
185
    """Return standard geographic projection (EPSG 4326) in the WKT or PROJ4 format.

186
187
    :param fmt:     <str> target format - 'WKT' or 'PROJ4'
    """
188
189
    if not re.search('wkt', fmt, re.I) or re.search('Proj4', fmt, re.I):
        raise ValueError(fmt, 'Unsupported output format.')
190

191
192
193
194
    if re.search('wkt', fmt, re.I):
        return CRS.from_epsg(4326).to_wkt()
    else:
        return CRS.from_epsg(4326).to_proj4()