Coverage for py_tools_ds/geo/projection.py: 62%

1# -*- coding: utf-8 -*- 

2 

3# py_tools_ds - A collection of geospatial data analysis tools that simplify standard 

4# operations when handling geospatial raster and vector data as well as projections. 

5# 

6# Copyright (C) 2016-2021 

7# - Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de) 

8# - Helmholtz Centre Potsdam - GFZ German Research Centre for Geosciences Potsdam, 

9# Germany (https://www.gfz-potsdam.de/) 

10# 

11# This software was developed within the context of the GeoMultiSens project funded 

12# by the German Federal Ministry of Education and Research 

13# (project grant code: 01 IS 14 010 A-C). 

14# 

15# Licensed under the Apache License, Version 2.0 (the "License"); 

16# you may not use this file except in compliance with the License. 

17# You may obtain a copy of the License at 

18# 

19# http://www.apache.org/licenses/LICENSE-2.0 

20# 

21# Unless required by applicable law or agreed to in writing, software 

22# distributed under the License is distributed on an "AS IS" BASIS, 

23# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 

24# See the License for the specific language governing permissions and 

25# limitations under the License. 

26 

27import re 

28import warnings 

29from pyproj import CRS 

30from typing import Union # noqa F401 # flake8 issue 

31from osgeo import osr 

32 

33from ..environment import gdal_env 

34 

35__author__ = "Daniel Scheffler" 

36 

37 

38# try to set GDAL_DATA if not set or invalid 

39gdal_env.try2set_GDAL_DATA() 

40 

41 

42def get_proj4info(proj=None): 

43 # type: (Union[str, int]) -> str 

44 """Return PROJ4 formatted projection info for the given projection. 

45 

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>) 

49 """ 

50 return CRS.from_user_input(proj).to_proj4().strip() 

51 

52 

53def proj4_to_dict(proj4): 

54 # type: (str) -> dict 

55 """Convert a PROJ4-like string into a dictionary. 

56 

57 :param proj4: <str> the PROJ4-like string 

58 """ 

59 return CRS.from_proj4(proj4).to_dict() 

60 

61 

62def dict_to_proj4(proj4dict): 

63 # type: (dict) -> str 

64 """Convert a PROJ4-like dictionary into a PROJ4 string. 

65 

66 :param proj4dict: <dict> the PROJ4-like dictionary 

67 """ 

68 return CRS.from_dict(proj4dict).to_proj4() 

69 

70 

71def proj4_to_WKT(proj4str): 

72 # type: (str) -> str 

73 """Convert a PROJ4-like string into a WKT string. 

74 

75 :param proj4str: <dict> the PROJ4-like string 

76 """ 

77 return CRS.from_proj4(proj4str).to_wkt() 

78 

79 

80def prj_equal(prj1, prj2): 

81 # type: (Union[None, int, str], Union[None, int, str]) -> bool 

82 """Check if the given two projections are equal. 

83 

84 :param prj1: projection 1 (WKT or 'epsg:1234' or <EPSG_int>) 

85 :param prj2: projection 2 (WKT or 'epsg:1234' or <EPSG_int>) 

86 """ 

87 if prj1 is None and prj2 is None or prj1 == prj2: 

88 return True 

89 else: 

90 # CRS.equals was added in pyproj 2.5 

91 crs1 = CRS.from_user_input(prj1) 

92 crs2 = CRS.from_user_input(prj2) 

93 

94 return crs1.equals(crs2) 

95 

96 

97def isProjectedOrGeographic(prj): 

98 # type: (Union[str, int, dict]) -> Union[str, None] 

99 """ 

100 

101 :param prj: accepts EPSG, Proj4 and WKT projections 

102 """ 

103 if prj is None: 

104 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 

109 

110 

111def isLocal(prj): 

112 # type: (Union[str, int, dict]) -> Union[bool, None] 

113 """ 

114 

115 :param prj: accepts EPSG, Proj4 and WKT projections 

116 """ 

117 if not prj: 

118 return True 

119 

120 srs = osr.SpatialReference() 

121 if prj.startswith('EPSG:'): 

122 srs.ImportFromEPSG(int(prj.split(':')[1])) 

123 elif prj.startswith('+proj='): 

124 srs.ImportFromProj4(prj) 

125 elif 'GEOGCS' in prj or \ 

126 'GEOGCRS' in prj or \ 

127 'PROJCS' in prj or \ 

128 'PROJCS' in prj or \ 

129 'LOCAL_CS' in prj: 

130 srs.ImportFromWkt(prj) 

131 else: 

132 raise RuntimeError('Unknown input projection: \n%s' % prj) 

133 

134 return srs.IsLocal() 

135 

136 

137def EPSG2Proj4(EPSG_code): 

138 # type: (int) -> str 

139 return CRS.from_epsg(EPSG_code).to_proj4() if EPSG_code is not None else '' 

140 

141 

142def 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 '' 

145 

146 

147def WKT2EPSG(wkt): 

148 # type: (str) -> Union[int, None] 

149 """Transform a WKT string to an EPSG code. 

150 

151 :param wkt: WKT definition 

152 :returns: EPSG code 

153 """ 

154 if not isinstance(wkt, str): 

155 raise TypeError("'wkt' must be a string. Received %s." % type(wkt)) 

156 if not wkt: 

157 return None 

158 

159 wkt = ''.join([line.strip().replace('\n', '').replace('\r', '') 

160 for line in wkt.splitlines()]) 

161 ccrs = CRS.from_wkt(wkt) 

162 

163 if not ccrs.is_bound: 

164 epsg = ccrs.to_epsg() 

165 else: 

166 epsg = ccrs.source_crs.to_epsg() 

167 

168 if epsg is None: 

169 warnings.warn('Could not find a suitable EPSG code for the input WKT string.', RuntimeWarning) 

170 else: 

171 return epsg 

172 

173 

174def get_UTMzone(prj): 

175 # type: (str) -> Union[int, None] 

176 if isProjectedOrGeographic(prj) == 'projected': 

177 srs = osr.SpatialReference() 

178 srs.ImportFromWkt(prj) 

179 return srs.GetUTMZone() 

180 else: 

181 return None 

182 

183 

184def get_prjLonLat(fmt='wkt'): 

185 # type: (str) -> Union[str, dict] 

186 """Return standard geographic projection (EPSG 4326) in the WKT or PROJ4 format. 

187 

188 :param fmt: <str> target format - 'WKT' or 'PROJ4' 

189 """ 

190 if not re.search('wkt', fmt, re.I) or re.search('Proj4', fmt, re.I): 

191 raise ValueError(fmt, 'Unsupported output format.') 

192 

193 if re.search('wkt', fmt, re.I): 

194 return CRS.from_epsg(4326).to_wkt() 

195 else: 

196 return CRS.from_epsg(4326).to_proj4()