Coverage for py_tools_ds/geo/map_info.py: 89%

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 math 

29import warnings 

30import os 

31from typing import Union # noqa F401 # flake8 issue 

32from tempfile import TemporaryDirectory 

33 

34from pyproj import CRS 

35import numpy as np 

36from osgeo import gdal, osr 

37 

38from .projection import isLocal 

39 

40__author__ = "Daniel Scheffler" 

41 

42 

43class Geocoding(object): 

44 def __init__(self, mapinfo=None, gt=None, prj=''): 

45 # type: (Union[list, tuple], Union[list, tuple], str) -> None 

46 """Get an instance of the Geocoding object. 

47 

48 :param mapinfo: ENVI map info, e.g., ['UTM', 1, 1, 192585.0, 5379315.0, 30.0, 30.0, 41, 'North', 'WGS-84'] 

49 :param gt: GDAL GeoTransform, e.g. (249885.0, 30.0, 0.0, 4578615.0, 0.0, -30.0) 

50 :param prj: GDAL Projection - WKT Format 

51 """ 

52 # FIXME this class only supports UTM and Geographic Lat/Lon projections 

53 self.prj_name = 'Arbitrary' 

54 self.ul_x_map = 0. 

55 self.ul_y_map = 0. 

56 self.ul_x_px = 1. 

57 self.ul_y_px = 1. 

58 self.gsd_x = 1. 

59 self.gsd_y = 1. 

60 self.rot_1_deg = 0. 

61 self.rot_1_rad = 0. 

62 self.rot_2_deg = 0. 

63 self.rot_2_rad = 0. 

64 self.utm_zone = 0 

65 self.utm_north_south = 'North' 

66 self.datum = '' 

67 self.units = '' 

68 

69 if mapinfo: 

70 if gt or prj: 

71 warnings.warn("'gt' and 'prj' are not respected if mapinfo is given!") 

72 self.from_mapinfo(mapinfo) 

73 elif gt: 

74 self.from_geotransform_projection(gt, prj) 

75 

76 def from_geotransform_projection(self, gt, prj): 

77 # type: (Union[list, tuple], str) -> 'Geocoding' 

78 """Create Geocoding object from GDAL GeoTransform + WKT projection string. 

79 

80 HOW COMPUTATION OF RADIANTS WORKS: 

81 Building on top of the computation within self.to_geotransform(): 

82 gt[1] = math.cos(rotation_rad) * gsd_x 

83 gt[2] = math.sin(rotation_rad) * gsd_x 

84 

85 -> we have to solve this equation system to get rotation_rad: 

86 gsd_x = gt[2] / math.sin(rotation_rad) 

87 gt[1] = math.cos(rotation_rad) * gt[2] / math.sin(rotation_rad) 

88 gt[1] * math.sin(rotation_rad) = math.cos(rotation_rad) * gt[2] 

89 math.sin(rotation_rad) / math.cos(rotation_rad) = gt[2] / gt[1] 

90 math.tan(rotation_rad) = gt[2] / gt[1] 

91 rotation_rad = math.atan(gt[2] / gt[1]) 

92 

93 :param gt: GDAL GeoTransform, e.g. (249885.0, 30.0, 0.0, 4578615.0, 0.0, -30.0) 

94 :param prj: GDAL Projection - WKT Format 

95 :return: instance of Geocoding 

96 """ 

97 if gt not in [None, [0, 1, 0, 0, 0, -1], (0, 1, 0, 0, 0, -1)]: 

98 # validate input geotransform 

99 if not isinstance(gt, (list, tuple)): 

100 raise TypeError("'gt' must be a list or a tuple. Received type %s." % type(gt)) 

101 if len(gt) != 6: 

102 raise ValueError("'gt' must contain 6 elements.") 

103 

104 self.ul_x_map = float(gt[0]) 

105 self.ul_y_map = float(gt[3]) 

106 

107 # handle rotations 

108 if float(gt[2]) == 0: 

109 # no rotation. use default angles from init 

110 self.gsd_x = float(gt[1]) 

111 else: 

112 self.rot_1_rad = math.atan(gt[2] / gt[1]) 

113 self.rot_1_deg = math.degrees(self.rot_1_rad) 

114 self.gsd_x = gt[2] / math.sin(self.rot_1_rad) 

115 

116 if float(gt[4]) == 0: 

117 # no rotation. use default angles from init 

118 self.gsd_y = float(abs(gt[5])) 

119 else: 

120 self.rot_2_rad = math.atan(gt[4] / gt[5]) 

121 self.rot_2_deg = math.degrees(self.rot_2_rad) 

122 self.gsd_y = gt[4] / math.sin(self.rot_2_rad) 

123 

124 # handle projection 

125 srs = osr.SpatialReference() 

126 srs.ImportFromWkt(prj) 

127 

128 if isLocal(prj): 

129 self.prj_name = 'Arbitrary' 

130 else: 

131 with warnings.catch_warnings(): 

132 warnings.filterwarnings("ignore", message="You will likely lose important projection information") 

133 

134 # get prj_name and datum 

135 proj4 = CRS(prj).to_dict() # FIXME avoid to convert to proj4 

136 self.prj_name = \ 

137 'Geographic Lat/Lon' if proj4['proj'] == 'longlat' else \ 

138 'UTM' if proj4['proj'] == 'utm' else proj4['proj'] 

139 # FIXME there is no 'datum' key in case of, e.g., LAEA projection # still true? 

140 # -> use CRS.datum.name instead? 

141 self.datum = 'WGS-84' if proj4['datum'] == 'WGS84' else proj4['datum'] # proj4['ellps']? 

142 self.units = proj4['unit'] if 'unit' in proj4 else self.units 

143 

144 if self.prj_name == 'UTM': 

145 utm_zone_srs = srs.GetUTMZone() 

146 self.utm_zone = abs(utm_zone_srs) # always positive 

147 self.utm_north_south = \ 

148 'North' if utm_zone_srs >= 0. else 'South' 

149 

150 del srs 

151 

152 return self 

153 

154 def from_mapinfo(self, mapinfo): 

155 # type: (Union[list, tuple]) -> 'Geocoding' 

156 """Create Geocoding object from ENVI map info. 

157 

158 :param mapinfo: ENVI map info, e.g., ['UTM', 1, 1, 192585.0, 5379315.0, 30.0, 30.0, 41, 'North', 'WGS-84'] 

159 :return: instance of Geocoding 

160 """ 

161 if mapinfo: 

162 # validate input map info 

163 if not isinstance(mapinfo, (list, tuple)): 

164 raise TypeError("'mapinfo' must be a list or a tuple. Received type %s." % type(mapinfo)) 

165 

166 def assert_mapinfo_length(min_len): 

167 if len(mapinfo) < min_len: 

168 raise ValueError("A map info of type '%s' must contain at least %s elements. Received %s." 

169 % (mapinfo[0], min_len, len(mapinfo))) 

170 

171 assert_mapinfo_length(10 if mapinfo[0] == 'UTM' else 9 if mapinfo[0] == 'Arbitrary' else 8) 

172 

173 # parse mapinfo 

174 self.prj_name = mapinfo[0] 

175 self.ul_x_px, self.ul_y_px, self.ul_x_map, self.ul_y_map, self.gsd_x = (float(i) for i in mapinfo[1:6]) 

176 self.gsd_y = float(abs(mapinfo[6])) 

177 

178 if self.prj_name == 'UTM': 

179 self.utm_zone = mapinfo[7] 

180 self.utm_north_south = mapinfo[8] 

181 self.datum = mapinfo[9] 

182 else: 

183 self.datum = mapinfo[7] 

184 

185 # handle rotation 

186 for i in mapinfo: 

187 if isinstance(i, str) and re.search('rotation', i, re.I): 

188 self.rot_1_deg = float(i.split('=')[1].strip()) 

189 self.rot_2_deg = self.rot_1_deg 

190 self.rot_1_rad = math.radians(self.rot_1_deg) 

191 self.rot_2_rad = math.radians(self.rot_2_deg) 

192 

193 return self 

194 

195 def to_geotransform(self): 

196 # type: () -> tuple 

197 """Return GDAL GeoTransform list using the attributes of the Geocoding instance. 

198 

199 For equations, see: 

200 https://gis.stackexchange.com/questions/229952/rotate-envi-hyperspectral-imagery-with-gdal/229962 

201 

202 :return: GDAL GeoTransform, e.g. [249885.0, 30.0, 0.0, 4578615.0, 0.0, -30.0] 

203 """ 

204 # handle pixel coordinates of UL unequal to (1/1) 

205 ul_map_x = self.ul_x_map if self.ul_x_px == 1 else (self.ul_x_map - (self.ul_x_px * self.gsd_x - self.gsd_x)) 

206 ul_map_y = self.ul_y_map if self.ul_y_px == 1 else (self.ul_y_map - (self.ul_y_px * self.gsd_y - self.gsd_y)) 

207 

208 # handle rotation and pixel sizes 

209 gsd_x, rot_1 = (self.gsd_x, 0) if self.rot_1_deg == 0 else \ 

210 (math.cos(self.rot_1_rad) * self.gsd_x, math.sin(self.rot_1_rad) * self.gsd_x) 

211 gsd_y, rot_2 = (self.gsd_y, 0) if self.rot_2_deg == 0 else \ 

212 (math.cos(self.rot_2_rad) * self.gsd_y, math.sin(self.rot_2_rad) * self.gsd_y) 

213 

214 return ul_map_x, gsd_x, rot_1, ul_map_y, rot_2, -gsd_y 

215 

216 def to_mapinfo(self): 

217 """Return ENVI map info list using the attributes of the Geocoding instance. 

218 

219 :return: ENVI map info, e.g. [ UTM , 1 , 1 , 256785.0 , 4572015.0 , 30.0 , 30.0 , 43 , North , WGS-84 ] 

220 """ 

221 mapinfo = [self.prj_name, self.ul_x_px, self.ul_y_px, self.ul_x_map, self.ul_y_map, self.gsd_x, abs(self.gsd_y)] 

222 

223 # add UTM infos 

224 if self.prj_name in ['UTM', 'Arbitrary']: 

225 mapinfo.extend([self.utm_zone, self.utm_north_south]) 

226 

227 # add datum 

228 if self.prj_name != 'Arbitrary': 

229 mapinfo.append(self.datum) 

230 

231 # add rotation 

232 if self.rot_1_deg != 0.: 

233 mapinfo.append('rotation=%.5f' % self.rot_1_deg) 

234 

235 return mapinfo 

236 

237 

238def geotransform2mapinfo(gt, prj): 

239 # type: (Union[list, tuple, None], Union[str, None]) -> list 

240 """Build an ENVI geo info from given GDAL GeoTransform and Projection (compatible with UTM and LonLat projections). 

241 

242 :param gt: GDAL GeoTransform, e.g. (249885.0, 30.0, 0.0, 4578615.0, 0.0, -30.0) 

243 :param prj: GDAL Projection - WKT Format 

244 :returns: ENVI geo info, e.g. [ UTM , 1 , 1 , 256785.0 , 4572015.0 , 30.0 , 30.0 , 43 , North , WGS-84 ] 

245 :rtype: list 

246 """ 

247 try: 

248 return Geocoding(gt=gt, prj=prj).to_mapinfo() 

249 

250 except KeyError: # KeyError: 'datum' - in case of, e.g., ETRS/LAEA projection 

251 if int(gdal.__version__[0]) < 3: 

252 # noinspection PyTypeChecker 

253 prj = CRS(prj).to_wkt(version="WKT1_GDAL") 

254 

255 with TemporaryDirectory() as fdir: 

256 fn = "py_tools_ds__geotransform2mapinfo_temp" 

257 fP_bsq = os.path.join(fdir, fn + '.bsq') 

258 fP_hdr = os.path.join(fdir, fn + '.hdr') 

259 

260 ds_out = gdal.GetDriverByName('ENVI').Create(fP_bsq, 2, 2, 1, gdal.GDT_Int32) 

261 ds_out.SetGeoTransform(gt) 

262 ds_out.SetProjection(prj) 

263 ds_out.FlushCache() 

264 del ds_out 

265 

266 with open(fP_hdr, 'r') as inF: 

267 content = inF.read() 

268 if 'map info' in content: 

269 res = re.search("map info = {(.*?)}", content, re.I).group(1) 

270 map_info = [i.strip() for i in res.split(',')] 

271 

272 for i, ele in enumerate(map_info): 

273 try: 

274 map_info[i] = float(ele) 

275 except ValueError: 

276 pass 

277 else: 

278 map_info = ['Arbitrary', 1.0, 1.0, 0.0, 0.0, 1.0, 1.0] 

279 

280 return map_info 

281 

282 

283def mapinfo2geotransform(map_info): 

284 # type: (Union[list, None]) -> tuple 

285 """Build GDAL GeoTransform tuple from an ENVI geo info. 

286 

287 :param map_info: ENVI geo info (list), e.g., ['UTM', 1, 1, 192585.0, 5379315.0, 30.0, 30.0, 41, 'North', 'WGS-84'] 

288 :returns: GDAL GeoTransform, e.g. [249885.0, 30.0, 0.0, 4578615.0, 0.0, -30.0] 

289 """ 

290 try: 

291 return Geocoding(mapinfo=map_info).to_geotransform() 

292 

293 except (KeyError, ValueError): # KeyError: 'datum' - in case of, e.g., ETRS/LAEA projection 

294 with TemporaryDirectory() as fdir: 

295 fn = "py_tools_ds__mapinfo2geotransform_temp" 

296 fP_bsq = os.path.join(fdir, fn + '.bsq') 

297 fP_hdr = os.path.join(fdir, fn + '.hdr') 

298 

299 ds_out = gdal.GetDriverByName('ENVI').Create(fP_bsq, 2, 2, 1) 

300 ds_out.GetRasterBand(1).WriteArray(np.array([[1, 2], [2, 3]])) 

301 ds_out.FlushCache() 

302 del ds_out 

303 

304 with open(fP_hdr, 'r') as InHdr: 

305 lines = InHdr.readlines() 

306 

307 lines.append('map info = { %s }\n' % ', '.join([str(i) for i in map_info])) 

308 

309 with open(fP_hdr, 'w') as OutHdr: 

310 OutHdr.writelines(lines) 

311 

312 ds = gdal.Open(fP_bsq) 

313 gt = ds.GetGeoTransform() 

314 del ds 

315 

316 return gt 

317 

318 

319def get_corner_coordinates(gdal_ds=None, gt=None, cols=None, rows=None): 

320 """Return (ULxy, LLxy, LRxy, URxy) in the same coordinate units like the given geotransform.""" 

321 assert gdal_ds or (gt and cols and rows), \ 

322 "GEOP.get_corner_coordinates: Missing argument! Please provide either 'gdal_ds' or 'gt', 'cols' AND 'rows'." 

323 

324 gdal_ds_GT = gdal_ds.GetGeoTransform() if gdal_ds else gt 

325 ext = [] 

326 xarr = [0, gdal_ds.RasterXSize if gdal_ds else cols] 

327 yarr = [0, gdal_ds.RasterYSize if gdal_ds else rows] 

328 

329 for px in xarr: 

330 for py in yarr: 

331 x = gdal_ds_GT[0] + (px * gdal_ds_GT[1]) + (py * gdal_ds_GT[2]) 

332 y = gdal_ds_GT[3] + (px * gdal_ds_GT[4]) + (py * gdal_ds_GT[5]) 

333 ext.append([x, y]) 

334 yarr.reverse() 

335 del gdal_ds_GT 

336 

337 return ext