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

2 

3# gms_preprocessing, spatial and spectral homogenization of satellite remote sensing data 

4# 

5# Copyright (C) 2020 Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de) 

6# 

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

8# by the German Federal Ministry of Education and Research 

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

10# 

11# This program is free software: you can redistribute it and/or modify it under 

12# the terms of the GNU General Public License as published by the Free Software 

13# Foundation, either version 3 of the License, or (at your option) any later version. 

14# Please note the following exception: `gms_preprocessing` depends on tqdm, which 

15# is distributed under the Mozilla Public Licence (MPL) v2.0 except for the files 

16# "tqdm/_tqdm.py", "setup.py", "README.rst", "MANIFEST.in" and ".gitignore". 

17# Details can be found here: https://github.com/tqdm/tqdm/blob/master/LICENCE. 

18# 

19# This program is distributed in the hope that it will be useful, but WITHOUT 

20# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 

21# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more 

22# details. 

23# 

24# You should have received a copy of the GNU Lesser General Public License along 

25# with this program. If not, see <http://www.gnu.org/licenses/>. 

26 

27"""Module 'metadata' for handling any type of metadata of GeoMultiSens compatible sensor systems.""" 

28 

29from __future__ import (division, print_function, unicode_literals, absolute_import) 

30 

31import collections 

32import datetime 

33import glob 

34import math 

35import os 

36import re 

37import sys 

38import warnings 

39import iso8601 

40import xml.etree.ElementTree as ET 

41from typing import List, TYPE_CHECKING # noqa F401 # flake8 issue 

42 

43import numpy as np 

44import pyproj 

45from matplotlib import dates as mdates 

46from pyorbital import astronomy 

47from natsort import natsorted 

48 

49from py_tools_ds.geo.map_info import geotransform2mapinfo 

50from py_tools_ds.geo.projection import WKT2EPSG 

51from pyrsr import RSR 

52from sicor.options import get_options as get_ac_options 

53 

54from ..options.config import GMS_config as CFG 

55from ..io.input_reader import open_specific_file_within_archive, Solar_Irradiance_reader, SRF_Reader 

56from ..io.output_writer import enviHdr_keyOrder 

57from ..algorithms import geoprocessing as GEOP 

58from ..misc import helper_functions as HLP_F 

59from ..misc import database_tools as DB_T 

60from ..misc.path_generator import path_generator, get_path_ac_options 

61from ..misc.definition_dicts import get_GMS_sensorcode, is_dataset_provided_as_fullScene, datasetid_to_sat_sen 

62from ..misc.exceptions import ACNotSupportedError 

63 

64if TYPE_CHECKING: 

65 from ..model.gms_object import GMS_identifier # noqa F401 # flake8 issue 

66 

67__author__ = 'Daniel Scheffler', 'Robert Behling' 

68 

69 

70class METADATA(object): 

71 def __init__(self, GMS_id): 

72 # private attributes 

73 self._AcqDateTime = None 

74 

75 # unpack GMS_identifier 

76 self.GMS_identifier = GMS_id 

77 self.image_type = GMS_id.image_type 

78 self.Satellite = GMS_id.satellite 

79 self.Sensor = GMS_id.sensor 

80 self.Subsystem = GMS_id.subsystem 

81 self.proc_level = GMS_id.proc_level 

82 self.logger = GMS_id.logger 

83 

84 self.Dataname = '' 

85 self.FolderOrArchive = '' 

86 self.Metafile = "" # File containing image metadata (automatically found) 

87 self.EntityID = "" # ID to identify the original scene 

88 self.SceneID = '' # postgreSQL-database identifier 

89 self.Sensormode = "" 

90 self.gResolution = -99. # resolution [m] 

91 self.AcqDate = "" # YYYY-MM-DD 

92 self.AcqTime = "" # HH:MM:SS 

93 self.Offsets = {} # Dict with offset for each band (radiance) 

94 self.Gains = {} # Dict with gain for each band (radiance) 

95 self.OffsetsRef = {} # Dict with offset for each band for conversion to Reflectance (Landsat8) 

96 self.GainsRef = {} # Dict with gain for each band for conversion to Reflectance (Landsat8) 

97 self.CWL = {} 

98 self.FWHM = {} 

99 self.ProcLCode = "" # processing Level: Sensor specific Code 

100 self.ProcLStr = "" # processing Level: Sensor independent String (raw,rad cal, rad+geom cal, ortho) 

101 self.SunElevation = -99.0 # Sun elevation angle at center of product [Deg] 

102 # Sun azimuth angle at center of product, in degrees from North (clockwise) at the time of the first image line 

103 self.SunAzimuth = -99.0 

104 self.SolIrradiance = [] 

105 self.ThermalConstK1 = {} 

106 self.ThermalConstK2 = {} 

107 self.EarthSunDist = 1.0 

108 # viewing zenith angle of the Sensor (offNadir) [Deg] (usually) in Case of RapidEye "+" 

109 # being East and “-” being West 

110 self.ViewingAngle = -99.0 

111 self.ViewingAngle_arrProv = {} 

112 # viewing azimuth angle. The angle between the view direction of the satellite and a line perpendicular to 

113 # the image or tile center.[Deg] 

114 self.IncidenceAngle = -9999.99 

115 self.IncidenceAngle_arrProv = {} 

116 self.FOV = 9999.99 # field of view of the sensor [Deg] 

117 # Sensor specific quality code: See ro2/behling/Satelliten/doc_allg/ReadOutMetadata/SatMetadata.xls 

118 self.Quality = [] 

119 self.PhysUnit = "DN" 

120 # Factor to get reflectance values in [0-1]. Sentinel2A provides scaling factors for the Level1C 

121 # TOA-reflectance products 

122 self.ScaleFactor = -99 

123 self.CS_EPSG = -99. # EPSG-Code of coordinate system 

124 self.CS_TYPE = "" 

125 self.CS_DATUM = "" 

126 self.CS_UTM_ZONE = -99. 

127 self.WRS_path = -99. 

128 self.WRS_row = -99. 

129 # List of Corner Coordinates in order of Lon/Lat/DATA_X/Data_Y for all given image corners 

130 self.CornerTieP_LonLat = [] 

131 self.CornerTieP_UTM = [] 

132 self.LayerBandsAssignment = [] # List of spectral bands in the same order as they are stored in the rasterfile. 

133 self.additional = [] 

134 self.image_type = 'RSD' 

135 self.orbitParams = {} 

136 self.overpassDurationSec = -99. 

137 self.scene_length = -99. 

138 self.rows = -99. 

139 self.cols = -99. 

140 self.bands = -99. 

141 self.nBands = -99. 

142 self.map_info = [] 

143 self.projection = "" 

144 self.wvlUnit = "" 

145 self.spec_vals = {'fill': None, 'zero': None, 'saturated': None} 

146 

147 self.version_gms_preprocessing = CFG.version 

148 self.versionalias_gms_preprocessing = CFG.versionalias 

149 

150 def read_meta(self, scene_ID, stacked_image, data_folderOrArchive, LayerBandsAssignment=None): 

151 """ 

152 Read metadata. 

153 """ 

154 

155 self.SceneID = scene_ID 

156 self.Dataname = stacked_image 

157 self.FolderOrArchive = data_folderOrArchive 

158 self.LayerBandsAssignment = LayerBandsAssignment if LayerBandsAssignment else [] 

159 

160 if re.search(r"SPOT", self.Satellite, re.I): 

161 self.Read_SPOT_dimap2() 

162 elif re.search(r"Terra", self.Satellite, re.I): 

163 self.Read_ASTER_hdffile(self.Subsystem) 

164 elif re.search(r"Sentinel-2A", self.Satellite, re.I) or re.search(r"Sentinel-2B", self.Satellite, re.I): 

165 self.Read_Sentinel2_xmls() 

166 elif re.search(r"LANDSAT", self.Satellite, re.I): 

167 self.Read_LANDSAT_mtltxt(self.LayerBandsAssignment) 

168 elif re.search(r"RapidEye", self.Satellite, re.I): 

169 self.Read_RE_metaxml() 

170 elif re.search(r"ALOS", self.Satellite, re.I): 

171 self.Read_ALOS_summary() 

172 self.Read_ALOS_LEADER() # for gains and offsets 

173 else: 

174 raise NotImplementedError("%s is not a supported sensor or sensorname is misspelled." % self.Satellite) 

175 

176 return self 

177 

178 def __getstate__(self): 

179 """Defines how the attributes of MetaObj instances are pickled.""" 

180 if self.logger: 

181 self.logger.close() 

182 

183 return self.__dict__ 

184 

185 @property 

186 def AcqDateTime(self): 

187 """Returns a datetime.datetime object containing date, time and timezone (UTC time).""" 

188 if not self._AcqDateTime and self.AcqDate and self.AcqTime: 

189 self._AcqDateTime = datetime.datetime.strptime('%s %s%s' % (self.AcqDate, self.AcqTime, '.000000+0000'), 

190 '%Y-%m-%d %H:%M:%S.%f%z') 

191 

192 return self._AcqDateTime 

193 

194 @AcqDateTime.setter 

195 def AcqDateTime(self, DateTime): 

196 # type: (datetime.datetime) -> None 

197 if isinstance(DateTime, str): 

198 self._AcqDateTime = datetime.datetime.strptime(DateTime, '%Y-%m-%d %H:%M:%S.%f%z') 

199 elif isinstance(DateTime, datetime.datetime): 

200 self._AcqDateTime = DateTime 

201 

202 self.AcqDate = DateTime.strftime('%Y-%m-%d') 

203 self.AcqTime = DateTime.strftime('%H:%M:%S') 

204 

205 @property 

206 def overview(self): 

207 attr2include = \ 

208 ['Dataname', 'Metafile', 'EntityID', 'Satellite', 'Sensor', 'Subsystem', 'gResolution', 'AcqDate', 

209 'AcqTime', 'CWL', 'FWHM', 'Offsets', 'Gains', 'ProcLCode', 'ProcLStr', 'SunElevation', 'SunAzimuth', 

210 'ViewingAngle', 'IncidenceAngle', 'FOV', 'SolIrradiance', 'ThermalConstK1', 'ThermalConstK2', 

211 'EarthSunDist', 'Quality', 'PhysUnit', 'additional', 'GainsRef', 'OffsetsRef', 'CornerTieP_LonLat', 

212 'CS_EPSG', 'CS_TYPE', 'CS_DATUM', 'CS_UTM_ZONE', 'LayerBandsAssignment'] 

213 return collections.OrderedDict((attr, getattr(self, attr)) for attr in attr2include) 

214 

215 @property 

216 def LayerBandsAssignment_full(self): 

217 # type: () -> list 

218 """Return complete LayerBandsAssignment without excluding thermal or panchromatic bands. 

219 

220 NOTE: CFG.sort_bands_by_cwl is respected, so returned list may be sorted by central wavelength 

221 """ 

222 return get_LayerBandsAssignment(self.GMS_identifier, no_thermal=False, no_pan=False, return_fullLBA=True) 

223 

224 @property 

225 def bandnames(self): 

226 return [('Band %s' % i) for i in self.LayerBandsAssignment] 

227 

228 def Read_SPOT_dimap2(self): 

229 """----METHOD_2------------------------------------------------------------ 

230 # read metadata from spot dimap file 

231 """ 

232 

233 # self.default_attr() 

234 if os.path.isdir(self.FolderOrArchive): 

235 glob_res = glob.glob(os.path.join(self.FolderOrArchive, '*/scene01/metadata.dim')) 

236 assert len(glob_res) > 0, 'No metadata.dim file can be found in %s!' % self.FolderOrArchive 

237 self.Metafile = glob_res[0] 

238 dim_ = open(self.Metafile, "r").read() 

239 else: # archive 

240 dim_, self.Metafile = open_specific_file_within_archive(self.FolderOrArchive, '*/scene01/metadata.dim') 

241 

242 # special values 

243 h1 = re.findall(r"<SPECIAL_VALUE_INDEX>([a-zA-Z0-9]*)</SPECIAL_VALUE_INDEX>", dim_, re.I) 

244 h2 = re.findall(r"<SPECIAL_VALUE_TEXT>([a-zA-Z0-9]*)</SPECIAL_VALUE_TEXT>", dim_, re.I) 

245 self.additional.append(["SpecialValues:"]) 

246 for ii, ind in enumerate(h1): 

247 self.additional[0].append(["%s:%s" % (ind, h2[ii])]) 

248 

249 # EntityID 

250 h3 = re.search(r"<SOURCE_ID>([a-zA-Z0-9]*)</SOURCE_ID>", dim_, re.I) 

251 self.EntityID = h3.group(1) 

252 

253 # AcqDate 

254 h4 = re.search(r"<IMAGING_DATE>([0-9-]*)</IMAGING_DATE>", dim_, re.I) 

255 AcqDate = h4.group(1) 

256 

257 # Earth sun distance 

258 self.EarthSunDist = self.get_EarthSunDistance(AcqDate) 

259 

260 # AcqTime 

261 h5 = re.search(r"<IMAGING_TIME>([0-9:]*)</IMAGING_TIME>", dim_, re.I) 

262 AcqTime = h5.group(1) 

263 

264 # set self.AcqDateTime as well as self.AcqDate and self.AcqTime 

265 self.AcqDateTime = datetime.datetime.strptime( 

266 '%s %s%s' % (AcqDate, AcqTime, '.000000+0000'), '%Y-%m-%d %H:%M:%S.%f%z') 

267 

268 # Satellite + Sensor 

269 h6 = re.search(r"<MISSION>([a-zA-Z]*)</MISSION>[a-zA-Z0-9\s]*" 

270 r"<MISSION_INDEX>([0-9]*)</MISSION_INDEX>[a-zA-Z0-9\s]*" 

271 r"<INSTRUMENT>([a-zA-Z]*)</INSTRUMENT>[a-zA-Z0-9\s]*" 

272 r"<INSTRUMENT_INDEX>([0-9]*)</INSTRUMENT_INDEX>[a-zA-Z0-9\s]*" 

273 r"<SENSOR_CODE>([a-zA-Z0-9]*)</SENSOR_CODE>", 

274 dim_, re.I) 

275 self.Satellite = "-".join([h6.group(1), h6.group(2)]) 

276 self.Sensor = "".join([h6.group(3), h6.group(4), h6.group(5)]) 

277 

278 # Angles: incidence angle, sunAzimuth, sunElevation 

279 h7 = re.search(r"<INCIDENCE_ANGLE>(.*)</INCIDENCE_ANGLE>[\s\S]*" 

280 r"<SUN_AZIMUTH>(.*)</SUN_AZIMUTH>[\s\S]" 

281 r"*<SUN_ELEVATION>(.*)</SUN_ELEVATION>", dim_, re.I) 

282 self.IncidenceAngle = float(h7.group(1)) 

283 self.SunAzimuth = float(h7.group(2)) 

284 self.SunElevation = float(h7.group(3)) 

285 

286 # Viewing Angle: not always included in the metadata.dim file 

287 h8 = re.search(r"<VIEWING_ANGLE>(.*)</VIEWING_ANGLE>", dim_, re.I) 

288 if type(h8).__name__ == 'NoneType': 

289 self.ViewingAngle = 0 

290 else: 

291 self.ViewingAngle = float(h8.group(1)) 

292 

293 # Field of View 

294 self.FOV = get_FieldOfView(self.GMS_identifier) 

295 

296 # Units 

297 self.ScaleFactor = 1 

298 self.PhysUnit = "DN" 

299 

300 # ProcLevel 

301 h11 = re.search(r"<PROCESSING_LEVEL>([a-zA-Z0-9]*)</PROCESSING_LEVEL>", dim_, re.I) 

302 self.ProcLCode = h11.group(1) 

303 

304 # Quality 

305 h12 = re.findall(r"<Bad_Pixel>[\s]*" 

306 r"<PIXEL_INDEX>([0-9]*)</PIXEL_INDEX>[\s]*" 

307 r"<BAD_PIXEL_STATUS>([^<]*)</BAD_PIXEL_STATUS>" 

308 r"</Bad_Pixel>", dim_, 

309 re.I) 

310 self.Quality = h12 

311 

312 # Coordinate Reference System 

313 re_CS_EPSG = re.search(r'<HORIZONTAL_CS_CODE>epsg:([0-9]*)</HORIZONTAL_CS_CODE>', dim_, re.I) 

314 re_CS_TYPE = re.search(r'<HORIZONTAL_CS_TYPE>([a-zA-Z0-9]*)</HORIZONTAL_CS_TYPE>', dim_, re.I) 

315 re_CS_DATUM = re.search(r'<HORIZONTAL_CS_NAME>([\w+\s]*)</HORIZONTAL_CS_NAME>', dim_, re.I) 

316 self.CS_EPSG = int(re_CS_EPSG.group(1)) if re_CS_EPSG is not None else self.CS_EPSG 

317 self.CS_TYPE = 'LonLat' if re_CS_TYPE is not None and re_CS_TYPE.group(1) == 'GEOGRAPHIC' else 'UTM' \ 

318 if re_CS_TYPE is not None and re_CS_TYPE.group(1) == 'UTM' else self.CS_TYPE 

319 self.CS_DATUM = 'WGS84' if re_CS_DATUM is not None and re_CS_DATUM.group(1) == 'WGS 84' else self.CS_DATUM 

320 

321 # Corner Coordinates 

322 h121 = re.findall(r"<TIE_POINT_CRS_X>(.*)</TIE_POINT_CRS_X>", dim_, re.I) 

323 h122 = re.findall(r"<TIE_POINT_CRS_Y>(.*)</TIE_POINT_CRS_Y>", dim_, re.I) 

324 if len(h121) == 4 and self.CS_TYPE == 'LonLat': 

325 # Set Corner Tie Point Coordinates (order = UL,UR,LL,LR) 

326 self.CornerTieP_LonLat.append(tuple([float(h121[0]), float(h122[0])])) # UL 

327 self.CornerTieP_LonLat.append(tuple([float(h121[1]), float(h122[1])])) # UR 

328 self.CornerTieP_LonLat.append(tuple([float(h121[3]), float(h122[3])])) # LL 

329 self.CornerTieP_LonLat.append(tuple([float(h121[2]), float(h122[2])])) # LR 

330 # ul_lon,ul_lat = self.inDs.GetGCPs()[0].GCPX,self.inDs.GetGCPs()[0].GCPY # funktioniert nur bei SPOT 

331 # lr_lon,lr_lat = self.inDs.GetGCPs()[2].GCPX,self.inDs.GetGCPs()[2].GCPY 

332 

333 ########################## 

334 # band specific metadata # 

335 ########################## 

336 

337 LBA_full_sorted = natsorted(self.LayerBandsAssignment_full) 

338 

339 # Gains and Offsets 

340 h9 = re.search(r"<Image_Interpretation>[\s\S]*</Image_Interpretation>", dim_, re.I) 

341 h9_ = h9.group(0) 

342 h91 = re.findall(r"<PHYSICAL_UNIT>([^<]*)</PHYSICAL_UNIT>", h9_, re.I) 

343 h92 = re.findall(r"<PHYSICAL_BIAS>([^<]*)</PHYSICAL_BIAS>", h9_, re.I) 

344 h93 = re.findall(r"<PHYSICAL_GAIN>([^<]*)</PHYSICAL_GAIN>", h9_, re.I) 

345 self.additional.append(["Physical Units per band:"]) 

346 for ii, ind in enumerate(h91): # FIXME does not respect sort_bands_by_cwl 

347 self.additional[-1].append(ind) 

348 # Offsets 

349 for ii, (ind, band) in enumerate(zip(h92, LBA_full_sorted)): 

350 self.Offsets[band] = float(ind) 

351 # Gains 

352 for ii, (ind, band) in enumerate(zip(h93, LBA_full_sorted)): 

353 # gains in dimap file represent reciprocal 1/gain 

354 self.Gains[band] = 1 / float(ind) 

355 

356 # Solar irradiance, central wavelengths , full width half maximum per band 

357 self.wvlUnit = 'Nanometers' 

358 # derive number of bands from number of given gains/offsets in header file or from raster data 

359 # noinspection PyBroadException 

360 try: 

361 self.nBands = (np.mean([len(self.Gains), len(self.Offsets)]) 

362 if np.std([len(self.Gains), len(self.Offsets)]) == 0 and len(self.Gains) != 0 else 

363 GEOP.GEOPROCESSING(self.Dataname, self.logger).bands) 

364 except Exception: 

365 self.logger.warning('Unable to get number of bands for the dataset %s Provider values are used for ' 

366 'solar irradiation, CWL and FWHM!.' % self.Dataname) 

367 self.LayerBandsAssignment = get_LayerBandsAssignment(self.GMS_identifier, self.nBands) 

368 

369 # Exact values calculated based in SRFs 

370 self.SolIrradiance, self.CWL, self.FWHM = self.calc_solar_irradiance_CWL_FWHM_per_band() 

371 # Provider values 

372 if not self.SolIrradiance: 

373 h10 = re.search(r"<Solar_Irradiance>[\s\S]*" 

374 r"</Solar_Irradiance>", dim_, re.I) 

375 h10_ = h10.group(0) 

376 h101 = re.findall(r"<SOLAR_IRRADIANCE_VALUE>([^<]*)" 

377 r"</SOLAR_IRRADIANCE_VALUE>", h10_, re.I) 

378 if h101: 

379 self.SolIrradiance = dict(zip(LBA_full_sorted, h101)) 

380 # self.additional.append(["Solar Irradiance per band:"]) 

381 # for ii,ind in enumerate(h101): 

382 # self.additional[-1].append(ind) 

383 else: # Preconfigured Irradiation values 

384 spot_irr_dic = { 

385 'SPOT1a': dict(zip(LBA_full_sorted, [1855., 1615., 1090., 1680.])), 

386 'SPOT1b': dict(zip(LBA_full_sorted, [1845., 1575., 1040., 1690.])), 

387 'SPOT2a': dict(zip(LBA_full_sorted, [1865., 1620., 1085., 1705.])), 

388 'SPOT2b': dict(zip(LBA_full_sorted, [1865., 1615., 1090., 1670.])), 

389 'SPOT3a': dict(zip(LBA_full_sorted, [1854., 1580., 1065., 1668.])), 

390 'SPOT3b': dict(zip(LBA_full_sorted, [1855., 1597., 1067., 1667.])), 

391 'SPOT4a': dict(zip(LBA_full_sorted, [1843., 1568., 1052., 233., 1568.])), 

392 'SPOT4b': dict(zip(LBA_full_sorted, [1851., 1586., 1054., 240., 1586.])), 

393 'SPOT5a': dict(zip(LBA_full_sorted, [1858., 1573., 1043., 236., 1762.])), 

394 'SPOT5b': dict(zip(LBA_full_sorted, [1858., 1575., 1047., 234., 1773.]))} 

395 self.SolIrradiance = spot_irr_dic[get_GMS_sensorcode(self.GMS_identifier)] 

396 # Preconfigured CWLs -- # ref: Guyot & Gu (1994): 'Effect of Radiometric Corrections on NDVI-Determined 

397 # from SPOT-HRV and Landsat-TM Data'; Hill 1990 Comparative Analysis of Landsat-5 TM and SPOT HRV-1 Data for 

398 # Use in Multiple Sensor Approaches ; # resource: SPOT techical report: Resolutions and spectral modes 

399 sensorcode = get_GMS_sensorcode(self.GMS_identifier)[:2] 

400 if sensorcode in ['SPOT1a', 'SPOT1b', 'SPOT2a', 'SPOT2b', 'SPOT3a', 'SPOT3b']: 

401 self.CWL = dict(zip(LBA_full_sorted, [545., 638., 819., 615.])) 

402 elif sensorcode in ['SPOT4a', 'SPOT4b']: 

403 self.CWL = dict(zip(LBA_full_sorted, [540., 650., 835., 1630., 645.])) 

404 elif sensorcode == ['SPOT5a', 'SPOT5b']: 

405 self.CWL = dict(zip(LBA_full_sorted, [540., 650., 835., 1630., 595.])) 

406 

407 self.orbitParams = get_orbit_params(self.GMS_identifier) 

408 self.filter_layerdependent_metadata() 

409 self.spec_vals = get_special_values(self.GMS_identifier) 

410 

411 def Read_LANDSAT_mtltxt(self, LayerBandsAssignment): 

412 """----METHOD_3------------------------------------------------------------ 

413 read metadata from LANDSAT metafile: <dataname>.MTL.txt. Metadatafile of LPGS processing chain 

414 :param LayerBandsAssignment: 

415 """ 

416 

417 # self.default_attr() 

418 self.LayerBandsAssignment = LayerBandsAssignment 

419 self.nBands = len(LayerBandsAssignment) 

420 if os.path.isdir(self.FolderOrArchive): 

421 glob_res = glob.glob(os.path.join(self.FolderOrArchive, '*MTL.txt')) 

422 assert len(glob_res) > 0, 'No *.MTL metadata file can be found in %s!' % self.FolderOrArchive 

423 self.Metafile = glob_res[0] 

424 mtl_ = open(self.Metafile, "r").read() 

425 else: # archive 

426 mtl_, self.Metafile = open_specific_file_within_archive(self.FolderOrArchive, '*MTL.txt') 

427 

428 # Processing Level 

429 h1 = re.search(r'PRODUCT_TYPE = "([a-zA-Z0-9]*)"', mtl_, re.I) 

430 if h1 is None: 

431 h1 = re.search(r'DATA_TYPE = "([a-zA-Z0-9]*)"', mtl_, re.I) 

432 self.ProcLCode = h1.group(1) 

433 

434 # Satellite + Sensor + Sensor Mode 

435 h2 = re.search(r'SPACECRAFT_ID = "([a-zA-Z0-9_]*)"[\s]*' 

436 r'SENSOR_ID = "([a-zA-Z0-9+]*)"[\s]*' 

437 r'SENSOR_MODE = "([\S]*)"', 

438 mtl_, re.I) 

439 if h2: 

440 self.Satellite = 'Landsat-%s' % re.search(r'LANDSAT[\D]*([0-9])', h2.group(1), re.I).group(1) 

441 self.Sensor = h2.group(2) 

442 self.Sensormode = h2.group(3) 

443 else: 

444 h2a = re.search(r'SPACECRAFT_ID = "([a-zA-Z0-9_]*)"', mtl_, re.I) 

445 h2b = re.search(r'SENSOR_ID = "([a-zA-Z0-9_+]*)"', mtl_, re.I) 

446 h2c = re.search(r'SENSOR_MODE = "([a-zA-Z0-9_+]*)"', mtl_, re.I) 

447 self.Satellite = 'Landsat-%s' % re.search(r'LANDSAT[\D]*([0-9])', h2a.group(1), re.I).group(1) 

448 self.Sensor = h2b.group(1) 

449 self.Sensormode = h2c.group(1) if h2c is not None else self.Sensormode # Landsat-8 

450 self.Sensor = 'ETM+' if self.Sensor == 'ETM' else self.Sensor 

451 

452 # EntityID 

453 h2c = re.search(r'LANDSAT_SCENE_ID = "([A-Z0-9]*)"', mtl_, re.I) 

454 if h2c: 

455 self.EntityID = h2c.group(1) 

456 

457 # Acquisition Date + Time 

458 h3 = re.search(r'ACQUISITION_DATE = ([0-9-]*)[\s]*' 

459 r'SCENE_CENTER_SCAN_TIME = "?([0-9:]*)"?', 

460 mtl_, re.I) 

461 if h3 is None: 

462 h3 = re.search(r'DATE_ACQUIRED = ([0-9-]*)[\s]*' 

463 r'SCENE_CENTER_TIME = "?([0-9:]*)"?', 

464 mtl_, re.I) 

465 AcqDate = h3.group(1) 

466 AcqTime = h3.group(2) 

467 

468 # set self.AcqDateTime as well as self.AcqDate and self.AcqTime 

469 self.AcqDateTime = datetime.datetime.strptime( 

470 '%s %s%s' % (AcqDate, AcqTime, '.000000+0000'), '%Y-%m-%d %H:%M:%S.%f%z') 

471 

472 # Earth sun distance 

473 self.EarthSunDist = self.get_EarthSunDistance(self.AcqDate) 

474 

475 # Units 

476 self.ScaleFactor = 1 

477 self.PhysUnit = "DN" 

478 

479 # Angles: incidence angle, sunAzimuth, sunElevation, field of view 

480 h5 = re.search(r"SUN_AZIMUTH = ([\S]*)[\s]*" 

481 r"SUN_ELEVATION = ([\S]*)", 

482 mtl_, re.I) 

483 self.SunAzimuth = float(h5.group(1)) 

484 self.SunElevation = float(h5.group(2)) 

485 self.FOV = get_FieldOfView(self.GMS_identifier) 

486 

487 # Quality 

488 h6 = re.search(r"GROUP = CORRECTIONS_APPLIED[\s\S]*" 

489 r"END_GROUP = CORRECTIONS_APPLIED", 

490 mtl_, re.I) 

491 if h6 is None: 

492 h6 = re.search(r"GROUP = IMAGE_ATTRIBUTES[\s\S]*" 

493 r"END_GROUP = IMAGE_ATTRIBUTES", 

494 mtl_, re.I) 

495 

496 h6_ = h6.group(0) 

497 h61 = (h6_.split("\n")) 

498 x = 0 

499 for i in h61: 

500 if x == 0 or x + 1 == len(h61): 

501 pass 

502 else: 

503 i_ = i.strip().replace('"', "") 

504 self.Quality.append(i_.split(" = ")) 

505 x += 1 

506 

507 # Additional: coordinate system, geom. Resolution 

508 h7 = re.search(r"GROUP = PROJECTION_PARAMETERS[\s\S]*END_GROUP = L1_METADATA_FILE", mtl_, re.I) 

509 h7_ = h7.group(0) 

510 h71 = (h7_.split("\n")) 

511 for x, i in enumerate(h71): 

512 if re.search(r"Group", i, re.I): 

513 pass 

514 else: 

515 i_ = i.strip().replace('"', "") 

516 self.additional.append(i_.split(" = ")) 

517 re_CS_TYPE = re.search(r'MAP_PROJECTION = "([\w+\s]*)"', h7_, re.I) 

518 re_CS_DATUM = re.search(r'DATUM = "([\w+\s]*)"', h7_, re.I) 

519 re_CS_UTM_ZONE = re.search(r'ZONE_NUMBER = ([0-9]*)\n', mtl_, re.I) 

520 re_CS_UTM_ZONE = re.search(r'UTM_ZONE = ([0-9]*)\n', h7_, 

521 re.I) if re_CS_UTM_ZONE is None else re_CS_UTM_ZONE # Landsat8 

522 self.CS_TYPE = 'LonLat' if re_CS_TYPE is not None and re_CS_TYPE.group(1) == 'GEOGRAPHIC' else 'UTM' \ 

523 if re_CS_TYPE is not None and re_CS_TYPE.group(1) == 'UTM' else self.CS_TYPE 

524 self.CS_DATUM = 'WGS84' if re_CS_DATUM is not None and re_CS_DATUM.group(1) == 'WGS84' else self.CS_DATUM 

525 self.CS_UTM_ZONE = int(re_CS_UTM_ZONE.group(1)) if re_CS_UTM_ZONE is not None else self.CS_UTM_ZONE 

526 # viewing Angle 

527 self.ViewingAngle = 0 

528 # Landsat8 

529 h8 = re.search(r"ROLL_ANGLE = ([\S]*)", mtl_, re.I) 

530 if h8: 

531 self.ViewingAngle = float(h8.group(1)) 

532 

533 # Corner Coordinates ## Lon/Lat for all given image corners UL,UR,LL,LR (tuples) 

534 h10_UL = re.findall(r"PRODUCT_UL_CORNER_[A-Z]+ = (.*)[\S]*", mtl_, re.I) 

535 h10_UR = re.findall(r"PRODUCT_UR_CORNER_[A-Z]+ = (.*)[\S]*", mtl_, re.I) 

536 h10_LL = re.findall(r"PRODUCT_LL_CORNER_[A-Z]+ = (.*)[\S]*", mtl_, re.I) 

537 h10_LR = re.findall(r"PRODUCT_LR_CORNER_[A-Z]+ = (.*)[\S]*", mtl_, re.I) 

538 if not h10_UL: # Landsat8 

539 h10_UL = re.findall(r"CORNER_UL_[\S]+ = (.*)[\S]*", mtl_, re.I) 

540 h10_UR = re.findall(r"CORNER_UR_[\S]+ = (.*)[\S]*", mtl_, re.I) 

541 h10_LL = re.findall(r"CORNER_LL_[\S]+ = (.*)[\S]*", mtl_, re.I) 

542 h10_LR = re.findall(r"CORNER_LR_[\S]+ = (.*)[\S]*", mtl_, re.I) 

543 if h10_UL: # Set Corner Tie Point Coordinates (order = UL,UR,LL,LR) 

544 self.CornerTieP_LonLat.append(tuple([float(h10_UL[i]) for i in [1, 0]])) 

545 self.CornerTieP_LonLat.append(tuple([float(h10_UR[i]) for i in [1, 0]])) 

546 self.CornerTieP_LonLat.append(tuple([float(h10_LL[i]) for i in [1, 0]])) 

547 self.CornerTieP_LonLat.append(tuple([float(h10_LR[i]) for i in [1, 0]])) 

548 self.CornerTieP_UTM.append(tuple([float(h10_UL[i]) for i in [2, 3]])) 

549 self.CornerTieP_UTM.append(tuple([float(h10_UR[i]) for i in [2, 3]])) 

550 self.CornerTieP_UTM.append(tuple([float(h10_LL[i]) for i in [2, 3]])) 

551 self.CornerTieP_UTM.append(tuple([float(h10_LR[i]) for i in [2, 3]])) 

552 

553 # WRS path/row 

554 h11_p = re.search(r'WRS_PATH = ([0-9]*)', mtl_, re.I) 

555 if h11_p: 

556 self.WRS_path = h11_p.group(1) 

557 h11_r1 = re.search(r'STARTING_ROW = ([0-9]*)', mtl_, re.I) 

558 h11_r2 = re.search(r'ENDING_ROW = ([0-9]*)', mtl_, re.I) 

559 if h11_r1 is None: # Landsat-8 

560 h11_r = re.search(r'WRS_ROW = ([0-9]*)', mtl_, re.I) 

561 self.WRS_row = int(h11_r.group(1)) 

562 else: 

563 self.WRS_row = int(h11_r1.group(1)) if h11_r1.group(1) == h11_r2.group(1) else self.WRS_row 

564 

565 # Fill missing values 

566 if self.EntityID == '': 

567 self.logger.info('Scene-ID could not be extracted and has to be retrieved from %s metadata database...' 

568 % self.Satellite) 

569 result = DB_T.get_info_from_postgreSQLdb(CFG.conn_database, 'scenes', ['entityID'], 

570 {'id': self.SceneID}) 

571 

572 if len(result) == 1: # e.g. [('LE71950282003121EDC00',)] 

573 self.EntityID = result[0][0] 

574 elif len(result) == 0: 

575 self.logger.warning("Scene ID could not be assigned because no dataset matching to the query " 

576 "parameters could be found in metadata database.") 

577 else: # e.g. [('LE71950282003121EDC00',), ('LE71950282003105ASN00',)] 

578 self.logger.warning("Scene ID could not be assigned because %s datasets matching to the query " 

579 "parameters were found in metadata database." % len(result)) 

580 # if self.EntityID=='': 

581 # dataset = 'LANDSAT_TM' if self.Satellite=='Landsat-5' else \ 

582 # 'LANDSAT_ETM' if self.Satellite=='Landsat-7' else 'LANDSAT_8' if self.Satellite=='Landsat-8' else '' 

583 # if dataset != '': 

584 # webmeta = list(usgsapi.search(dataset,'EE',distance=0,ll={'longitude':self.CornerTieP_LonLat[2][0], \ 

585 # 'latitude':self.CornerTieP_LonLat[2][1]},ur={'longitude':self.CornerTieP_LonLat[1][0], \ 

586 # 'latitude':self.CornerTieP_LonLat[1][1]},start_date=self.AcqDate,end_date=self.AcqDate)) 

587 # if len(webmeta)==1: 

588 # self.EntityID=webmeta[0]['entityId'] 

589 # else: 

590 # sen = {'MSS':'M','TM':'T','ETM+':'E','OLI_TIRS':'C','OLI':'O'}[self.Sensor] 

591 # sat = self.Satellite.split('-')[1] 

592 # path_res = re.search(r'WRS_PATH = ([0-9]+)',mtl_, re.I) 

593 # path_str = "%03d"%int(path_res.group(1)) if path_res!=None else '000' 

594 # row_res = re.search(r'STARTING_ROW = ([0-9]+)',mtl_, re.I) 

595 # if row_res is None: row_res = re.search(r'WRS_ROW = ([0-9]+)',mtl_, re.I) 

596 # row_str = "%03d"%int(row_res.group(1)) if row_res!=None else '000' 

597 # tt = datetime.datetime.strptime(self.AcqDate, '%Y-%m-%d').timetuple() 

598 # julianD = '%d%03d'%(tt.tm_year,tt.tm_yday) 

599 # station_res = re.search(r'GROUND_STATION = "([\S]+)"',mtl_, re.I) 

600 # if station_res is None: station_res = re.search(r'STATION_ID = "([\S]+)"',mtl_, re.I) 

601 # station_str = station_res.group(1) if station_res is not None else 'XXX' 

602 # idWithoutVersion = 'L%s%s%s%s%s%s'%(sen,sat,path_str,row_str,julianD,station_str) 

603 # filt_webmeta = [i for i in webmeta if i['entityId'].startswith(idWithoutVersion)] 

604 # if len(filt_webmeta) == 1: 

605 # self.EntityID = filt_webmeta[0]['entityId'] 

606 

607 ########################## 

608 # band specific metadata # 

609 ########################## 

610 LBA_full_sorted = natsorted(self.LayerBandsAssignment_full) 

611 

612 # Gains and Offsets 

613 h4 = re.search(r"GROUP = MIN_MAX_RADIANCE[\s\S]*END_GROUP = MIN_MAX_PIXEL_VALUE", mtl_, re.I) 

614 h4_ = h4.group(0) 

615 h4max_rad = re.findall(r" LMAX_BAND[0-9]+ = ([\S]*)", h4_, re.I) # space in front IS needed 

616 h4min_rad = re.findall(r" LMIN_BAND[0-9]+ = ([\S]*)", h4_, re.I) # space in front IS needed 

617 h4max_pix = re.findall(r"QCALMAX_BAND[0-9]+ = ([\S]*)", h4_, re.I) 

618 h4min_pix = re.findall(r"QCALMIN_BAND[0-9]+ = ([\S]*)", h4_, re.I) 

619 if not h4max_rad: 

620 h4max_rad = re.findall(r" RADIANCE_MAXIMUM_BAND_[0-9_VCID]+ = ([\S]*)", h4_, 

621 re.I) # space in front IS needed 

622 h4min_rad = re.findall(r" RADIANCE_MINIMUM_BAND_[0-9_VCID]+ = ([\S]*)", h4_, 

623 re.I) # space in front IS needed 

624 h4max_pix = re.findall(r"QUANTIZE_CAL_MAX_BAND_[0-9_VCID]+ = ([\S]*)", h4_, re.I) 

625 h4min_pix = re.findall(r"QUANTIZE_CAL_MIN_BAND_[0-9_VCID]+ = ([\S]*)", h4_, re.I) 

626 # additional: LMAX, LMIN, QCALMAX, QCALMIN 

627 self.additional.append([["LMAX"], ["LMIN"], ["QCALMAX"], ["QCALMIN"]]) 

628 # Offsets + Gains 

629 Gains, Offsets = [], [] 

630 for ii, ind in enumerate(h4min_rad): 

631 Gains.append( 

632 (float(h4max_rad[ii]) - float(h4min_rad[ii])) / (float(h4max_pix[ii]) - float(h4min_pix[ii]))) 

633 Offsets.append(float(ind) - float(h4min_pix[ii]) * Gains[-1]) 

634 self.additional[-1][-4].append(h4max_rad[ii]) 

635 self.additional[-1][-3].append(h4min_rad[ii]) 

636 self.additional[-1][-2].append(h4max_pix[ii]) 

637 self.additional[-1][-1].append(h4min_pix[ii]) 

638 self.Gains = {bN: Gains[i] for i, bN in enumerate(LBA_full_sorted)} 

639 self.Offsets = {bN: Offsets[i] for i, bN in enumerate(LBA_full_sorted)} 

640 

641 # Solar irradiance, central wavelengths , full width half maximum per band 

642 self.wvlUnit = 'Nanometers' 

643 # Exact values calculated based in SRFs 

644 self.SolIrradiance, self.CWL, self.FWHM = self.calc_solar_irradiance_CWL_FWHM_per_band() # 3x dict 

645 # Provider values 

646 if not self.SolIrradiance: # Preconfigured Irradiation and CWL values 

647 if re.search(r"[\D]*5", self.Satellite): 

648 # Landsat 5; resource for center wavelength (6 = thermal) 

649 self.SolIrradiance = {'1': 1957., 

650 '2': 1826., 

651 '3': 1554., 

652 '4': 1036., 

653 '5': 215., 

654 '6': 0.0, 

655 '7': 80.67} 

656 self.CWL = {'1': 485., 

657 '2': 560., 

658 '3': 660., 

659 '4': 830., 

660 '5': 1650., 

661 '6': 11450., 

662 '7': 2215.} 

663 if re.search(r"[\D]*7", self.Satellite): 

664 # Landsat 7; resource for center wavelength: 

665 # https://opticks.org/display/opticksDev/Sensor+Wavelength+Definitions 

666 # 6L(thermal),6H(thermal),B8(pan) 

667 self.SolIrradiance = {'1': 1969., 

668 '2': 1840., 

669 '3': 1551., 

670 '4': 1044., 

671 '5': 225.7, 

672 '6L': 0.0, 

673 '6H': 0.0, 

674 '7': 82.07, 

675 '8': 1368.} 

676 self.CWL = {'1': 482.5, 

677 '2': 565., 

678 '3': 660., 

679 '4': 825., 

680 '5': 1650., 

681 '6L': 11450., 

682 '6H': 11450., 

683 '7': 2215., 

684 '8': 710.} 

685 if re.search(r"[\D]*8", self.Satellite): # Landsat 8 

686 # no irradiation values available 

687 self.CWL = {'1': 443., 

688 '2': 482.6, 

689 '3': 561.3, 

690 '4': 654.6, 

691 '5': 864.6, 

692 '6': 1609.1, 

693 '7': 2201.2, 

694 '8': 591.7, 

695 '9': 1373.5, 

696 '10': 10903.6, 

697 '11': 12003.} 

698 # if None in SolIrradiance: 

699 

700 # Thermal Constants (K1 = [W*m-2*um-1]; K1 = [K]) 

701 if re.search(r"[\D]*5", self.Satellite): 

702 # resource: http://geography.middlebury.edu/data/gg1002/Handouts/Landsat_DN_Temp.pdf 

703 ThermalConstK1 = [0., 0., 0., 0., 0., 607.76, 0.] 

704 ThermalConstK2 = [0., 0., 0., 0., 0., 1260.56, 0.] 

705 if re.search(r"[\D]*7", self.Satellite): 

706 # resource: http://geography.middlebury.edu/data/gg1002/Handouts/Landsat_DN_Temp.pdf 

707 ThermalConstK1 = [0., 0., 0., 0., 0., 666.09, 666.09, 0., 0.] 

708 ThermalConstK2 = [0., 0., 0., 0., 0., 1282.71, 1282.71, 0., 0.] 

709 if re.search(r"[\D]*8", self.Satellite): # Landsat 8 

710 K1_res = re.findall(r"K1_CONSTANT_BAND_[0-9]+ = ([\S]*)", mtl_, re.I) 

711 K2_res = re.findall(r"K2_CONSTANT_BAND_[0-9]+ = ([\S]*)", mtl_, re.I) 

712 if len(K1_res) == 2: 

713 ThermalConstK1 = [0., 0., 0., 0., 0., 0., 0., 0., 0., float(K1_res[0]), float(K1_res[1])] 

714 ThermalConstK2 = [0., 0., 0., 0., 0., 0., 0., 0., 0., float(K2_res[0]), float(K2_res[1])] 

715 else: 

716 self.logger.error('Unable to set thermal constants. Expected to derive 2 values for K1 and K2. ' 

717 'Found %s' % len(K1_res)) 

718 self.ThermalConstK1 = {bN: ThermalConstK1[i] for i, bN in enumerate(LBA_full_sorted)} 

719 self.ThermalConstK2 = {bN: ThermalConstK2[i] for i, bN in enumerate(LBA_full_sorted)} 

720 

721 # reflectance coefficients (Landsat8) 

722 h9 = re.search(r"GROUP = RADIOMETRIC_RESCALING[\s\S]*END_GROUP = RADIOMETRIC_RESCALING", mtl_, re.I) 

723 if h9: 

724 h9_ = h9.group(0) 

725 h9gain_ref = re.findall(r" REFLECTANCE_MULT_BAND_[0-9]+ = ([\S]*)", h9_, re.I) 

726 h9gain_ref_bandNr = re.findall(r" REFLECTANCE_MULT_BAND_([0-9]+) = [\S]*", h9_, re.I) 

727 h9offs_ref = re.findall(r" REFLECTANCE_ADD_BAND_[0-9]+ = ([\S]*)", h9_, re.I) 

728 h9offs_ref_bandNr = re.findall(r" REFLECTANCE_ADD_BAND_([0-9]+) = [\S]*", h9_, re.I) 

729 if h9gain_ref: 

730 GainsRef = [None] * len(self.Gains) 

731 OffsetsRef = [None] * len(self.Offsets) 

732 

733 for ii, ind in zip(h9gain_ref_bandNr, h9gain_ref): 

734 GainsRef[LBA_full_sorted.index(ii)] = ind 

735 for ii, ind in zip(h9offs_ref_bandNr, h9offs_ref): 

736 OffsetsRef[LBA_full_sorted.index(ii)] = ind 

737 

738 self.GainsRef = {bN: GainsRef[i] for i, bN in enumerate(LBA_full_sorted)} 

739 self.OffsetsRef = {bN: OffsetsRef[i] for i, bN in enumerate(LBA_full_sorted)} 

740 

741 self.orbitParams = get_orbit_params(self.GMS_identifier) 

742 self.filter_layerdependent_metadata() 

743 self.spec_vals = get_special_values(self.GMS_identifier) 

744 

745 # mtl.close() 

746 

747 def Read_RE_metaxml(self): 

748 """----METHOD_4------------------------------------------------------------ 

749 read metadata from RapidEye metafile: <dataname>metadata.xml 

750 """ 

751 

752 # self.default_attr() 

753 if os.path.isdir(self.FolderOrArchive): 

754 glob_res = glob.glob(os.path.join(self.FolderOrArchive, '*/*_metadata.xml')) 

755 assert len(glob_res) > 0, 'No *metadata.xml file can be found in %s/*/!' % self.FolderOrArchive 

756 self.Metafile = glob_res[0] 

757 mxml_ = open(self.Metafile, "r").read() 

758 else: # archive 

759 mxml_, self.Metafile = open_specific_file_within_archive(self.FolderOrArchive, '*/*_metadata.xml') 

760 

761 # EntityID 

762 h1 = re.search(r"<[a-z]*:identifier>([\S]*)</[a-z]*:identifier>", mxml_, re.I) 

763 self.EntityID = h1.group(1) if h1 else "Not found" 

764 

765 # Processing Level 

766 h2 = re.search(r"<[a-z]*:productType>([a-zA-Z0-9]*)</[a-z]*:productType>", mxml_, re.I) 

767 self.ProcLCode = h2.group(1) if h2 else "Not found" 

768 

769 # Satellite 

770 h3 = re.search(r"<[a-z]*:serialIdentifier>([a-zA-Z0-9-]*)</[a-z]*:serialIdentifier>", mxml_, re.I) 

771 self.Satellite = 'RapidEye-%s' % re.search(r'[\D]*([0-9])', h3.group(1), re.I).group(1) if h3 else "Not found" 

772 

773 # Sensor (Instrument shortname) 

774 h4 = re.search(r"<[a-z]*:Instrument>[\s]*<eop:shortName>([\S]*)</[a-z]*:shortName>", mxml_, re.I) 

775 self.Sensor = h4.group(1) if h4 else "Not found" 

776 

777 # Acquisition Parameters: Incidence Angle, SunAzimuth, SunElevation, ViewingAngle, FieldOfView, AcqDate, AcqTime 

778 try: 

779 h5 = re.search(r'<[a-z]*:incidenceAngle uom="deg">([\S]*)' 

780 r'</[a-z]*:incidenceAngle>[\s]*' 

781 r'<opt:illuminationAzimuthAngle uom="deg">([\S]*)' 

782 r'</opt:illuminationAzimuthAngle>[\s]*' 

783 r'<opt:illuminationElevationAngle uom="deg">([\S]*)' 

784 r'</opt:illuminationElevationAngle>[\s]*' 

785 r'<re:azimuthAngle uom="deg">([\S]*)' 

786 r'</re:azimuthAngle>[\s]*' 

787 r'<re:spaceCraftViewAngle uom="deg">([\S]*)' 

788 r'</re:spaceCraftViewAngle>[\s]*' 

789 r'<re:acquisitionDateTime>([0-9-]*)T([0-9:]*)[\S]*' 

790 r'</re:acquisitionDateTime>', 

791 mxml_, re.I) 

792 self.IncidenceAngle = float(h5.group(1)) 

793 self.SunAzimuth = float(h5.group(2)) 

794 self.SunElevation = float(h5.group(3)) 

795 # angle from true north at the image or tile center to the scan (line) direction at image center, 

796 # in clockwise positive degrees. 

797 self.additional.append(["spaceCraftAzimuthAngle:", str(round(float(h5.group(4)), 1))]) 

798 self.ViewingAngle = float(h5.group(5)) 

799 self.FOV = get_FieldOfView(self.GMS_identifier) 

800 AcqDate = h5.group(6) 

801 AcqTime = h5.group(7) 

802 

803 except AttributeError: 

804 h5 = re.search(r'<hma:acrossTrackIncidenceAngle uom="deg">([\S]*)' 

805 r'</hma:acrossTrackIncidenceAngle>[\s]*' 

806 r'<ohr:illuminationAzimuthAngle uom="deg">([\S]*)' 

807 r'</ohr:illuminationAzimuthAngle>[\s]*' 

808 r'<ohr:illuminationElevationAngle uom="deg">([\S]*)' 

809 r'</ohr:illuminationElevationAngle>[\s]*' 

810 r'<re:azimuthAngle uom="deg">([\S]*)' 

811 r'</re:azimuthAngle>[\s]*' 

812 r'<re:acquisitionDateTime>([0-9-]*)' 

813 r'T([0-9:]*)[\S]*' 

814 r'</re:acquisitionDateTime>', 

815 mxml_, re.I) # 

816 self.IncidenceAngle = float(h5.group(1)) 

817 self.SunAzimuth = float(h5.group(2)) 

818 self.SunElevation = float(h5.group(3)) 

819 self.additional.append(["spaceCraftAzimuthAngle:", str(round(float(h5.group(4)), 1))]) 

820 self.ViewingAngle = 9999.99 

821 self.FOV = get_FieldOfView(self.GMS_identifier) 

822 AcqDate = h5.group(5) 

823 AcqTime = h5.group(6) 

824 # set self.AcqDateTime as well as self.AcqDate and self.AcqTime 

825 self.AcqDateTime = datetime.datetime.strptime( 

826 '%s %s%s' % (AcqDate, AcqTime, '.000000+0000'), '%Y-%m-%d %H:%M:%S.%f%z') 

827 

828 # Earth sun distance 

829 self.EarthSunDist = self.get_EarthSunDistance(self.AcqDate) 

830 

831 # Additional: Projection 

832 h6 = re.search(r"<re:geodeticDatum>([\S]*)" 

833 r"</re:geodeticDatum>[\s]*" 

834 r"<re:projection>([\S^<\s]*)" 

835 r"</re:projection>", 

836 mxml_, re.I) 

837 try: 

838 self.additional.append(["SpatialReference", 

839 ["geodeticDatum", h6.group(1)], 

840 ["projection", h6.group(2)] 

841 ]) 

842 except AttributeError: # Level1-B data has no geodaticDatum 

843 pass 

844 

845 # Corrections applied + Quality 

846 h7 = re.search(r"<re:radiometricCorrectionApplied>([\w]*)" 

847 r"</re:radiometricCorrectionApplied>[\s]*" 

848 r"<re:radiometricCalibrationVersion>([\S]*)" 

849 r"</re:radiometricCalibrationVersion>[\s]*" 

850 r"<re:geoCorrectionLevel>([\S\s^<]*)" 

851 r"</re:geoCorrectionLevel>[\s]*" 

852 r"<re:elevationCorrectionApplied>([\S]*)" 

853 r"</re:elevationCorrectionApplied>[\s]*" 

854 r"<re:atmosphericCorrectionApplied>([\w]*)" 

855 r"</re:atmosphericCorrectionApplied>[\s]*" 

856 r"<re:productAccuracy>([\S\s^<]*)" 

857 r"</re:productAccuracy>", mxml_, re.I) 

858 # fuer IRIS ihre Daten 

859 if h7 is None: 

860 h7 = re.search(r"<re:radiometricCorrectionApplied>([\w]*)" 

861 r"</re:radiometricCorrectionApplied>[\s]*" 

862 r"<re:geoCorrectionLevel>([\S\s^<]*)" 

863 r"</re:geoCorrectionLevel>[\s]*" 

864 r"<re:elevationCorrectionApplied>([\S]*)" 

865 r"</re:elevationCorrectionApplied>[\s]*" 

866 r"<re:atmosphericCorrectionApplied>([\w]*)" 

867 r"</re:atmosphericCorrectionApplied>", 

868 mxml_, re.I) 

869 self.additional.append( 

870 ["Corrections", 

871 ["radiometricCorrectionApplied", h7.group(1)], 

872 ["geoCorrectionLevel", h7.group(2)], 

873 ["elevationCorrectionApplied", h7.group(3)], 

874 ["atmosphericCorrectionApplied", h7.group(4)] 

875 ]) 

876 else: 

877 self.additional.append( 

878 ["Corrections", 

879 ["radiometricCorrectionApplied", h7.group(1)], 

880 ["radiometricCalibrationVersion", h7.group(2)], 

881 ["geoCorrectionLevel", h7.group(3)], 

882 ["elevationCorrectionApplied", h7.group(4)], 

883 ["atmosphericCorrectionApplied", h7.group(5)] 

884 ]) 

885 self.Quality.append(["geomProductAccuracy[m]:", str( 

886 round(float(h7.group(6)), 1))]) # Estimated product horizontal CE90 uncertainty [m] 

887 

888 # additional. missing lines, suspectlines, binning, shifting, masking 

889 h81 = re.findall(r"<re:percentMissingLines>([^<]*)</re:percentMissingLines>", mxml_, re.I) 

890 h82 = re.findall(r"<re:percentSuspectLines>([^<]*)</re:percentSuspectLines>", mxml_, re.I) 

891 h83 = re.findall(r"<re:binning>([^<]*)</re:binning>", mxml_, re.I) 

892 h84 = re.findall(r"<re:shifting>([^<]*)</re:shifting>", mxml_, re.I) 

893 h85 = re.findall(r"<re:masking>([^<]*)</re:masking>", mxml_, re.I) 

894 

895 self.Quality.append( 

896 ["MissingLines[%]perBand", h81]) # Percentage of missing lines in the source data of this band 

897 # Percentage of suspect lines (lines that contained downlink errors) in the source data for the band 

898 self.Quality.append(["SuspectLines[%]perBand", h82]) 

899 self.additional.append( 

900 ["binning_perBand", h83]) # Indicates the binning used (across track x along track) [1x1,2x2,3x3,1x2,2x1] 

901 self.additional.append( 

902 ["shifting_perBand", h84]) # Indicates the sensor applied right shifting [none, 1bit, 2bits, 3bits, 4bits] 

903 self.additional.append(["masking_perBand", h85]) # Indicates the sensor applied masking [111, 110, 100, 000] 

904 

905 # Units 

906 self.ScaleFactor = 1 

907 self.PhysUnit = "DN" 

908 

909 # Coordinate Reference System 

910 re_CS_EPSG = re.search( 

911 r'<re:ProductInformation>[\s\S]*' 

912 r'<re:epsgCode>([0-9]*)</re:epsgCode>[\s\S]*' 

913 r'</re:ProductInformation>', 

914 mxml_, re.I) 

915 re_CS_TYPE = re.search( 

916 r'<re:ProductInformation>[\s\S]*' 

917 r'<re:projection>([\s\S]*)</re:projection>[\s\S]*' 

918 r'</re:ProductInformation>', 

919 mxml_, re.I) 

920 re_CS_DATUM = re.search( 

921 r'<re:ProductInformation>[\s\S]*' 

922 r'<re:geodeticDatum>([\w+\s]*)</re:geodeticDatum>[\s\S]*' 

923 r'</re:ProductInformation>', 

924 mxml_, re.I) 

925 re_CS_UTM_ZONE = re.search( 

926 r'<re:ProductInformation>[\s\S]*' 

927 r'<re:projectionZone>([0-9]*)</re:projectionZone>[\s\S]*' 

928 r'</re:ProductInformation>', 

929 mxml_, re.I) 

930 self.CS_EPSG = int(re_CS_EPSG.group(1)) if re_CS_EPSG else self.CS_EPSG 

931 self.CS_TYPE = 'LonLat' if re_CS_TYPE and not re.search(r'UTM', re_CS_TYPE.group(1)) else 'UTM' \ 

932 if re_CS_TYPE and re.search(r'UTM', re_CS_TYPE.group(1)) else self.CS_TYPE 

933 self.CS_DATUM = 'WGS84' if re_CS_DATUM is not None and re_CS_DATUM.group(1) == 'WGS_1984' else self.CS_DATUM 

934 self.CS_UTM_ZONE = int(re_CS_UTM_ZONE.group(1)) if re_CS_UTM_ZONE is not None else self.CS_UTM_ZONE 

935 

936 # Corner Coordinates ## Lon/Lat for all given image corners UL,UR,LL,LR (tuples) 

937 h10_UL = re.findall( 

938 r'<re:topLeft>' 

939 r'<re:latitude>(.*)</re:latitude>' 

940 r'<re:longitude>(.*)</re:longitude>' 

941 r'</re:topLeft>', 

942 mxml_, re.I) 

943 h10_UR = re.findall( 

944 r'<re:topRight>' 

945 r'<re:latitude>(.*)</re:latitude>' 

946 r'<re:longitude>(.*)</re:longitude>' 

947 r'</re:topRight>', 

948 mxml_, re.I) 

949 h10_LL = re.findall( 

950 r'<re:bottomLeft>' 

951 r'<re:latitude>(.*)</re:latitude>' 

952 r'<re:longitude>(.*)</re:longitude>' 

953 r'</re:bottomLeft>', 

954 mxml_, re.I) 

955 h10_LR = re.findall( 

956 r'<re:bottomRight>' 

957 r'<re:latitude>(.*)</re:latitude>' 

958 r'<re:longitude>(.*)</re:longitude>' 

959 r'</re:bottomRight>', 

960 mxml_, re.I) 

961 if h10_UL: # Set Corner Tie Point Coordinates (order = UL,UR,LL,LR) 

962 self.CornerTieP_LonLat.append(tuple([float(h10_UL[0][1]), float(h10_UL[0][0])])) 

963 self.CornerTieP_LonLat.append(tuple([float(h10_UR[0][1]), float(h10_UR[0][0])])) 

964 self.CornerTieP_LonLat.append(tuple([float(h10_LL[0][1]), float(h10_LL[0][0])])) 

965 self.CornerTieP_LonLat.append(tuple([float(h10_LR[0][1]), float(h10_LR[0][0])])) 

966 

967 ########################## 

968 # band specific metadata # 

969 ########################## 

970 

971 LBA_full_sorted = natsorted(self.LayerBandsAssignment_full) 

972 

973 # Gains + Offsets 

974 h9 = re.findall(r"<re:radiometricScaleFactor>([^<]*)</re:radiometricScaleFactor>", 

975 mxml_, re.I) 

976 self.Gains = dict(zip(LBA_full_sorted, [float(gain) for gain in h9])) 

977 self.Offsets = dict(zip(LBA_full_sorted, [0, 0, 0, 0, 0])) 

978 

979 # Solar irradiance, central wavelengths , full width half maximum per band 

980 self.wvlUnit = 'Nanometers' 

981 # derive number of bands from number of given gains/offsets in header file or from raster data 

982 # noinspection PyBroadException 

983 try: 

984 self.nBands = (np.mean([len(self.Gains), len(self.Offsets)]) 

985 if np.std([len(self.Gains), len(self.Offsets)]) == 0 and len(self.Gains) != 0 else 

986 GEOP.GEOPROCESSING(self.Dataname, self.logger).bands) 

987 except Exception: 

988 self.logger.warning('Unable to get number of bands for the dataset %s Provider values are used for ' 

989 'solar irradiation, CWL and FWHM!.' % self.Dataname) 

990 self.LayerBandsAssignment = get_LayerBandsAssignment(self.GMS_identifier, self.nBands) 

991 

992 # Exact values calculated based in SRFs 

993 self.SolIrradiance, self.CWL, self.FWHM = self.calc_solar_irradiance_CWL_FWHM_per_band() 

994 # Provider values 

995 if not self.SolIrradiance: 

996 # Preconfigured Irradiation values 

997 self.SolIrradiance = dict(zip(LBA_full_sorted, [1997.8, 1863.5, 1560.4, 1395.0, 1124.4])) 

998 # Preconfigured CWLs 

999 self.CWL = dict(zip(LBA_full_sorted, [475., 555., 657., 710., 805.])) 

1000 

1001 self.orbitParams = get_orbit_params(self.GMS_identifier) 

1002 self.filter_layerdependent_metadata() 

1003 self.spec_vals = get_special_values(self.GMS_identifier) 

1004 

1005 def Read_ASTER_hdffile(self, subsystem): 

1006 """#----METHOD_5------------------------------------------------------------ 

1007 read metadata from ASTER hdf 

1008 input: 

1009 hdffile: 

1010 subsystem: 

1011 output: 

1012 :param subsystem: 

1013 """ 

1014 

1015 # self.default_attr() 

1016 dat_ = open(self.FolderOrArchive, "r").read() if sys.version_info[0] < 3 else \ 

1017 open(self.FolderOrArchive, "rb").read().decode('latin-1') 

1018 

1019 # Split meta from raster data 

1020 meta = re.search(r"GROUP[\s]*=[\s]" 

1021 r"ASTERGENERICMETADATA[\s\S]*?" 

1022 r"END_GROUP[\s]*=[\s]" 

1023 r"INVENTORYMETADATA", 

1024 dat_, re.I) 

1025 meta_ = meta.group(0) if meta else '' 

1026 genericmeta = re.search(r"GROUP[\s]*=[\s]" 

1027 r"ASTERGENERICMETADATA[\s\S]*?" 

1028 r"END_GROUP[\s]*=[\s]" 

1029 r"ASTERGENERICMETADATA", 

1030 meta_, re.I) 

1031 inventorymeta = re.search(r"GROUP[\s]*=[\s]" 

1032 r"INVENTORYMETADATA[\s\S]*?" 

1033 r"END_GROUP[\s]*=[\s]" 

1034 r"INVENTORYMETADATA", 

1035 meta_, re.I) 

1036 gcsgenericmeta = re.search(r"GROUP[\s]*=[\s]" 

1037 r"GDSGENERICMETADATA[\s\S]*?" 

1038 r"END_GROUP[\s]*=[\s]" 

1039 r"GDSGENERICMETADATA", 

1040 meta_, re.I) 

1041 vnirmeta = re.search(r"GROUP[\s]*=[\s]" 

1042 r"PRODUCTSPECIFICMETADATAVNIR[\s\S]*?" 

1043 r"END_GROUP[\s]*=[\s]" 

1044 r"PRODUCTSPECIFICMETADATAVNIR", 

1045 meta_, re.I) 

1046 swirmeta = re.search(r"GROUP[\s]*=[\s]" 

1047 r"PRODUCTSPECIFICMETADATASWIR[\s\S]*?" 

1048 r"END_GROUP[\s]*=[\s]" 

1049 r"PRODUCTSPECIFICMETADATASWIR", 

1050 meta_, re.I) 

1051 tirmeta = re.search(r"GROUP[\s]*=[\s]" 

1052 r"PRODUCTSPECIFICMETADATATIR[\s\S]*?" 

1053 r"END_GROUP[\s]*=[\s]" 

1054 r"PRODUCTSPECIFICMETADATATIR", 

1055 meta_, re.I) 

1056 genericmeta_ = genericmeta.group(0) if genericmeta else '' 

1057 inventorymeta_ = inventorymeta.group(0) if inventorymeta else '' 

1058 gcsgenericmeta_ = gcsgenericmeta.group(0) if gcsgenericmeta else '' 

1059 vnirmeta_ = vnirmeta.group(0) if vnirmeta else '' 

1060 swirmeta_ = swirmeta.group(0) if swirmeta else '' 

1061 tirmeta_ = tirmeta.group(0) if tirmeta else '' 

1062 h_ = '\n\n\n'.join([genericmeta_, inventorymeta_, gcsgenericmeta_, vnirmeta_, swirmeta_, tirmeta_]) 

1063 

1064 # with open("./testing/out/ASTER_HDFmeta__h_.txt", "w") as allMetaOut: 

1065 # allMetaOut.write(h_) 

1066 

1067 # EntityID 

1068 h1 = re.search(r"OBJECT[\s]*=[\s]*" 

1069 r"IDOFASTERGDSDATAGRANULE[\s\S]*" 

1070 r"END_OBJECT[\s]*=[\s]*" 

1071 r"IDOFASTERGDSDATAGRANULE", 

1072 h_, re.I) 

1073 h11 = re.search(r'VALUE[\s]*=[\s]*"([\s\S]*)"', h1.group(0), re.I) 

1074 self.EntityID = [h11.group(1), re.search(r"\"(ASTL1A[\s0-9]*)\"", h_).group(1)] 

1075 

1076 # Viewing Angle 

1077 h2 = re.search(r"GROUP[\s]*=[\s]*" 

1078 r"POINTINGANGLES[\s\S]*" 

1079 r"END_GROUP[\s]*=[\s]*" 

1080 r"POINTINGANGLES", 

1081 h_, re.I) 

1082 h21 = re.findall(r"VALUE[\s]*=[\s]*([-]*[0-9.0-9]+)", h2.group(0), re.I) 

1083 self.additional.append(["ViewingAngles", "VNIR", float(h21[0]), "SWIR", float(h21[1]), "TIR", float(h21[2])]) 

1084 if max(float(h21[0]), float(h21[1]), float(h21[2])) - min(float(h21[0]), float(h21[1]), float(h21[2])) < 1: 

1085 self.ViewingAngle = float(h21[0]) 

1086 else: 

1087 self.ViewingAngle = -99.0 

1088 

1089 # additional GainMode 

1090 h3 = re.search(r"GROUP[\s]*=[\s]*" 

1091 r"GAININFORMATION[\s\S]*" 

1092 r"END_GROUP[\s]*=[\s]*" 

1093 r"GAININFORMATION", 

1094 genericmeta_, re.I) 

1095 h31 = re.findall(r'VALUE[\s]*=[\s]*[\S]?\"[0-9A-Z]*\", \"([A-Z]*)\"', h3.group(0), re.I) 

1096 gains = {'HGH': [170.8, 179.0, 106.8, 27.5, 8.8, 7.9, 7.55, 5.27, 4.02, 'N/A', 'N/A', 'N/A', 'N/A', 'N/A'], 

1097 'NOR': [427.0, 358.0, 218.0, 55.0, 17.6, 15.8, 15.1, 10.55, 8.04, 28.17, 27.75, 26.97, 23.30, 21.38], 

1098 'LOW': [569.0, 477.0, 290.0, 73.3, 23.4, 21.0, 20.1, 14.06, 10.72, 'N/A', 'N/A', 'N/A', 'N/A', 'N/A'], 

1099 'LO1': [569.0, 477.0, 290.0, 73.3, 23.4, 21.0, 20.1, 14.06, 10.72, 'N/A', 'N/A', 'N/A', 'N/A', 'N/A'], 

1100 'LO2': ['N/A', 'N/A', 'N/A', 73.3, 103.5, 98.7, 83.8, 62.0, 67.0, 'N/A', 'N/A', 'N/A', 'N/A', 'N/A'], 

1101 'OFF': "OFF"} 

1102 self.additional.append([["GainMode:"], ["Max_radiance:"]]) 

1103 for x, i in enumerate(h31[:15]): 

1104 self.additional[-1][-2].append(i) 

1105 self.additional[-1][-1].append(gains[i][x]) 

1106 

1107 # Units 

1108 self.ScaleFactor = 1 

1109 self.PhysUnit = "DN" 

1110 

1111 # Satellite 

1112 self.Satellite = 'Terra' 

1113 

1114 # Sensor 

1115 h10 = re.search(r"OBJECT[\s]*=[\s]*" 

1116 r"INSTRUMENTSHORTNAME[\s\S]*" 

1117 r"END_OBJECT[\s]*=[\s]*" 

1118 r"INSTRUMENTSHORTNAME", h_, re.I) 

1119 self.Sensor = re.search(r"VALUE[\s]*=[\s]*[\"]([A-Za-z]*)\"", h10.group(0), re.I).group(1) 

1120 

1121 # Subsystem 

1122 h5 = re.search(r"GROUP[\s]*=[\s]*" 

1123 r"OBSERVATIONMODE[\s\S]*" 

1124 r"END_GROUP[\s]*=[\s]*" 

1125 r"OBSERVATIONMODE", h_, re.I) 

1126 avail_subsystems = re.findall(r'VALUE[\s]*=[\s]*[(]\"([a-zA-Z0-9]*)\", \"([ONF]*)\"', h5.group(0), re.I) 

1127 if subsystem not in ['VNIR1', 'VNIR2', 'SWIR', 'TIR']: 

1128 raise ValueError('Unexpected subsystem >%s<. Unable to specify the correct ASTER subsystem to be ' 

1129 'processed.' % subsystem) 

1130 else: 

1131 if subsystem == 'VNIR1' and avail_subsystems[0][1] == 'ON': 

1132 self.nBands = 3 

1133 if subsystem == 'VNIR2' and avail_subsystems[1][1] == 'ON': 

1134 self.nBands = 1 

1135 if subsystem == 'SWIR' and avail_subsystems[2][1] == 'ON': 

1136 self.nBands = 6 

1137 if subsystem == 'TIR' and avail_subsystems[3][1] == 'ON': 

1138 self.nBands = 5 

1139 self.Subsystem = subsystem 

1140 

1141 # Field of view (requires Satellite, Sensor, Subsystem) 

1142 self.FOV = get_FieldOfView(self.GMS_identifier) 

1143 

1144 # Ground resolution 

1145 re_res_GSD = re.findall(r'OBJECT[\s]*=[\s]*' 

1146 r'SPATIALRESOLUTION[\s\S]*' 

1147 r'VALUE[\s]*=[\s]*[\S]{1}([1-9][0-9]), ([1-9][0-9]), ([1-9][0-9])[\s\S]*' 

1148 r'END_OBJECT[\s]*=[\s]*' 

1149 r'SPATIALRESOLUTION', h_, re.I) 

1150 idx_subsystem = \ 

1151 0 if subsystem[:4] == 'VNIR' else \ 

1152 1 if subsystem[:4] == 'SWIR' else \ 

1153 2 if subsystem[:4] == 'TIR' else None 

1154 self.gResolution = float(re_res_GSD[0][idx_subsystem]) \ 

1155 if re_res_GSD and idx_subsystem is not None else self.gResolution 

1156 

1157 # Flight direction 

1158 h6 = re.search(r"OBJECT[\s]*=[\s]*" 

1159 r"FLYINGDIRECTION[\s\S]*\"([ASDE]*?)\"", 

1160 h_, re.I) 

1161 Fdir = {'AS': "Ascending", 'DE': "Descending"} 

1162 self.additional.append(["Flight Direction", Fdir[h6.group(1)]]) 

1163 

1164 # SunAzimuth SunElevation 

1165 h7 = re.search(r"OBJECT[\s]*=[\s]*" 

1166 r"SOLARDIRECTION[\s\S]*" 

1167 r"END_OBJECT[\s]*=[\s]*" 

1168 r"SOLARDIRECTION", 

1169 h_, re.I) 

1170 h71 = re.search(r"VALUE[\s]*=[\s]*[\S]{1}([0-9.]*), ([0-9.]*)", h7.group(0), re.I) 

1171 self.SunAzimuth = float(h71.group(1)) 

1172 self.SunElevation = float(h71.group(2)) 

1173 

1174 # data Quality 

1175 h8 = re.findall(r"GROUP[\s]*=[\s]*" 

1176 r"DATAQUALITY[1-9][0-4BN]?[^.]*" 

1177 r"END_GROUP[\s]*=[\s]*" 

1178 r"DATAQUALITY[1-9][0-4BN]?", 

1179 h_, re.I) 

1180 h81 = re.findall(r'VALUE[\s]*=[\s]*([^\n]*)', "\n".join(h8), re.I) 

1181 self.Quality.append(["NoOfBadPixel:"]) 

1182 for i in h81: 

1183 self.Quality[-1].append(i) 

1184 

1185 # ProdLevel 

1186 h9 = re.search(r"OBJECT[\s]*=[\s]*" 

1187 r"SHORTNAME[\s\S]*" 

1188 r"END_OBJECT[\s]*=[\s]*" 

1189 r"SHORTNAME", h_, re.I) 

1190 self.ProcLCode = re.search(r"VALUE[\s]*=[\s]*[\"]([0-9A-Za-z]*)\"", h9.group(0), re.I).group(1) 

1191 

1192 # AcqTime 

1193 h11 = re.search(r"OBJECT[\s]*=[\s]*" 

1194 r"TIMEOFDAY[\s\S]*" 

1195 r"END_OBJECT[\s]*=[\s]*" 

1196 r"TIMEOFDAY", h_, re.I) 

1197 h111 = re.search(r"VALUE[\s]*=[\s]*[\"]([\d][\d][\d][\d][\d][\d])[\S]*\"", h11.group(0), re.I) 

1198 if type(h111).__name__ != 'NoneType': 

1199 AcqTime = "%s:%s:%s" % (h111.group(1)[:2], h111.group(1)[2:4], h111.group(1)[4:6]) 

1200 else: 

1201 h112 = re.search(r"VALUE[\s]*=[\s]*[\"]([\d:]*)[\S]*\"", h11.group(0), re.I) 

1202 AcqTime = h112.group(1) 

1203 

1204 # AcqDate 

1205 h12 = re.search(r"OBJECT[\s]*=[\s]*" 

1206 r"CALENDARDATE[\s\S]*" 

1207 r"END_OBJECT[\s]*=[\s]*" 

1208 r"CALENDARDATE", h_, re.I) 

1209 h121 = re.search(r"VALUE[\s]*=[\s]*[\"]([\d]*)\"", h12.group(0), re.I) 

1210 if type(h121).__name__ != 'NoneType': 

1211 AcqDate = "%s-%s-%s" % (h121.group(1)[:4], h121.group(1)[4:6], h121.group(1)[6:8]) 

1212 else: 

1213 h121 = re.search(r"VALUE[\s]*=[\s]*[\"]([\d-]*)\"", h12.group(0), re.I) 

1214 AcqDate = h121.group(1) 

1215 

1216 # set self.AcqDateTime as well as self.AcqDate and self.AcqTime 

1217 self.AcqDateTime = datetime.datetime.strptime( 

1218 '%s %s%s' % (AcqDate, AcqTime, '.000000+0000'), '%Y-%m-%d %H:%M:%S.%f%z') 

1219 

1220 # Earth sun distance 

1221 self.EarthSunDist = self.get_EarthSunDistance(self.AcqDate) 

1222 

1223 # data Quality 

1224 h13 = re.search(r"GROUP[\s]*=[\s]*" 

1225 r"QASTATS[\s\S]*" 

1226 r"END_GROUP[\s]*=[\s]*" 

1227 r"QASTATS", 

1228 h_, re.I) 

1229 h131 = re.findall(r'VALUE[\s]*=[\s]*([0-9.0-9]*)', h13.group(0), re.I) 

1230 self.Quality.append(["PercMissingData:", h131[0]]) # The percentage of missing data of the scene [%] 

1231 self.Quality.append(["PercOutOfBoundsData:", h131[1]]) # The percentage of out of bounds data of the scene [%] 

1232 self.Quality.append(["PercInterpolatedData:", h131[2]]) # The percentage of interpolated data of the scene [%] 

1233 

1234 # Coordinate Reference System (L1B+) 

1235 topgroupname = 'PRODUCTSPECIFICMETADATA' + subsystem[:4] 

1236 # read only the information from first band in current ASTER subsystem 

1237 # (the others should have the same parameters) 

1238 re_res_procParm = re.search(r'GROUP[\s]*=[\s]*%s[\s\S]*' 

1239 r'GROUP[\s]*=[\s]*(PROCESSINGPARAMETERS%s[\s\S]*' 

1240 r'END_GROUP[\s]*=[\s]*PROCESSINGPARAMETERS%s)[\s\S]*' 

1241 r'END_GROUP[\s]*=[\s]*%s' 

1242 % (topgroupname, self.LayerBandsAssignment[0], 

1243 self.LayerBandsAssignment[0], topgroupname), 

1244 h_, re.I) 

1245 if type(re_res_procParm).__name__ != 'NoneType': 

1246 re_res_procParm = re_res_procParm.group(1) 

1247 re_CS_TYPE = re.search(r'OBJECT[\s]*=[\s]*' 

1248 r'MPMETHOD[\s\S]*' 

1249 r'VALUE[\s]*=[\s]*"([\s\S]*)"[\s\S]*' 

1250 r'END_OBJECT[\s]*=[\s]*' 

1251 r'MPMETHOD', 

1252 re_res_procParm, re.I) 

1253 re_CS_UTM_ZONE = re.search(r'OBJECT[\s]*=[\s]*' 

1254 r'UTMZONECODE[\s\S]*' 

1255 r'VALUE[\s]*=[\s]*([0-9][0-9]?)[\s\S]*' 

1256 r'END_OBJECT[\s]*=[\s]*' 

1257 r'UTMZONECODE', 

1258 re_res_procParm, re.I) 

1259 self.CS_TYPE = 'LonLat' if re_CS_TYPE and not re.search(r'UTM', re_CS_TYPE.group(1)) else 'UTM' \ 

1260 if re_CS_TYPE and re.search(r'UTM', re_CS_TYPE.group(1)) else self.CS_TYPE 

1261 self.CS_DATUM = 'WGS84' if self.CS_TYPE == 'UTM' else self.CS_DATUM 

1262 self.CS_UTM_ZONE = int(re_CS_UTM_ZONE.group(1)) if re_CS_UTM_ZONE else self.CS_UTM_ZONE 

1263 

1264 # Corner Coordinates ## Lon/Lat for all given image corners UL,UR,LL,LR (tuples) 

1265 re_res_sceneInfo = re.search(r'GROUP[\s]*=[\s]*' 

1266 r'SCENEINFORMATION[\s\S]*' 

1267 r'END_GROUP[\s]*=[\s]*' 

1268 r'SCENEINFORMATION', h_, 

1269 re.I).group(0) 

1270 re_res_UL = re.findall(r'OBJECT[\s]*=[\s]*' 

1271 r'UPPERLEFT[\s\S]*' 

1272 r'VALUE[\s]*=[\s]*[\S]{1}([0-9.]*), ([0-9.]*)[\s\S]*' 

1273 r'END_OBJECT[\s]*=[\s]*' 

1274 r'UPPERLEFT', 

1275 re_res_sceneInfo, re.I) 

1276 re_res_UR = re.findall(r'OBJECT[\s]*=[\s]*' 

1277 r'UPPERRIGHT[\s\S]*' 

1278 r'VALUE[\s]*=[\s]*[\S]{1}([0-9.]*), ([0-9.]*)[\s\S]*' 

1279 r'END_OBJECT[\s]*=[\s]*' 

1280 r'UPPERRIGHT', 

1281 re_res_sceneInfo, re.I) 

1282 re_res_LL = re.findall(r'OBJECT[\s]*=[\s]*' 

1283 r'LOWERLEFT[\s\S]*' 

1284 r'VALUE[\s]*=[\s]*[\S]{1}([0-9.]*), ([0-9.]*)[\s\S]*' 

1285 r'END_OBJECT[\s]*=[\s]*' 

1286 r'LOWERLEFT', 

1287 re_res_sceneInfo, re.I) 

1288 re_res_LR = re.findall(r'OBJECT[\s]*=[\s]*' 

1289 r'LOWERRIGHT[\s\S]*' 

1290 r'VALUE[\s]*=[\s]*[\S]{1}([0-9.]*), ([0-9.]*)[\s\S]*' 

1291 r'END_OBJECT[\s]*=[\s]*' 

1292 r'LOWERRIGHT', 

1293 re_res_sceneInfo, re.I) 

1294 if re_res_UL: # Set Corner Tie Point Coordinates (order = UL,UR,LL,LR) 

1295 self.CornerTieP_LonLat.append(tuple([float(re_res_UL[0][1]), float(re_res_UL[0][0])])) 

1296 self.CornerTieP_LonLat.append(tuple([float(re_res_UR[0][1]), float(re_res_UR[0][0])])) 

1297 self.CornerTieP_LonLat.append(tuple([float(re_res_LL[0][1]), float(re_res_LL[0][0])])) 

1298 self.CornerTieP_LonLat.append(tuple([float(re_res_LR[0][1]), float(re_res_LR[0][0])])) 

1299 

1300 ########################## 

1301 # band specific metadata # 

1302 ########################## 

1303 

1304 LBA_full_sorted = natsorted(self.LayerBandsAssignment_full) 

1305 

1306 # Gains/Offsets 

1307 h4 = re.findall(r"GROUP[\s]*=[\s]*" 

1308 r"UNITCONVERSIONCOEFF[1-9][0-4BN]?[^(]*" 

1309 r"END_GROUP[\s]*=[\s]*" 

1310 r"UNITCONVERSIONCOEFF[1-9][0-4BN]?", 

1311 h_, re.I) 

1312 h41 = re.findall(r'VALUE[\s]*=[\s]*([0-9.0-9]+)', "\n".join(h4), re.I) 

1313 h42 = re.findall(r'VALUE[\s]*=[\s]*([-]+[0-9.0-9]+)', "\n".join(h4), re.I) 

1314 self.Gains = dict(zip(LBA_full_sorted, [float(i) for i in h41])) 

1315 self.Offsets = dict(zip(LBA_full_sorted, [float(i) for i in h42])) 

1316 

1317 # Solar irradiance, central wavelengths , full width half maximum per band 

1318 self.wvlUnit = 'Nanometers' 

1319 self.LayerBandsAssignment = get_LayerBandsAssignment(self.GMS_identifier, self.nBands) 

1320 

1321 # Exact values calculated based in SRFs 

1322 self.SolIrradiance, self.CWL, self.FWHM = self.calc_solar_irradiance_CWL_FWHM_per_band() 

1323 # Provider values 

1324 if not self.SolIrradiance: 

1325 # Preconfigured Irradiation values 

1326 self.SolIrradiance = dict(zip(LBA_full_sorted, [1845.99, 1555.74, 1119.47] + [None] * 12)) 

1327 # Preconfigured CWLs 

1328 self.CWL = dict(zip(LBA_full_sorted, [556., 659., 807., 804., 1656., 2167., 2208., 2266., 2336., 2400., 

1329 8291., 8634., 9075., 10657., 11318])) 

1330 

1331 # Thermal Constants (K1 = [W*m-2*um-1]; K1 = [K]) 

1332 # Preconfigured values; resource: http://www.pancroma.com/ASTER%20Surface%20Temperature%20Analysis.html 

1333 self.ThermalConstK1 = dict(zip(LBA_full_sorted, 

1334 [0.] * 10 + [3040.136402, 2482.375199, 1935.060183, 866.468575, 641.326517])) 

1335 self.ThermalConstK2 = dict(zip(LBA_full_sorted, 

1336 [0.] * 10 + [1735.337945, 1666.398761, 1585.420044, 1350.069147, 1271.221673])) 

1337 

1338 self.orbitParams = get_orbit_params(self.GMS_identifier) 

1339 self.filter_layerdependent_metadata() 

1340 self.spec_vals = get_special_values(self.GMS_identifier) 

1341 

1342 def Read_ALOS_summary(self): 

1343 """----METHOD_6------------------------------------------------------------ 

1344 read metadata from ALOS summary.txt 

1345 """ 

1346 

1347 # self.default_attr() 

1348 if os.path.isdir(self.FolderOrArchive): 

1349 glob_res = glob.glob(os.path.join(self.FolderOrArchive, '*/*data*/summary.txt')) 

1350 assert len(glob_res) > 0, 'No summary.txt file can be found in %s/*/*data*/!' % self.FolderOrArchive 

1351 self.Metafile = glob_res[0] 

1352 sum_ = open(self.Metafile, "r").read() 

1353 else: # archive 

1354 sum_, self.Metafile = open_specific_file_within_archive(self.FolderOrArchive, '*/*data*/summary.txt') 

1355 

1356 # EntityID 

1357 h1 = re.search(r'Scs_SceneID="([a-zA-Z0-9]*)"', sum_, re.I) 

1358 self.EntityID = h1.group(1) 

1359 

1360 # Satellite, Sensor, ProcLevel, 

1361 h5 = re.search(r'Lbi_Satellite="([\S]*)"[\s]*' 

1362 r'Lbi_Sensor="([\S-]*)"[\s]*' 

1363 r'Lbi_ProcessLevel="([\S-]*)"', 

1364 sum_, re.I) 

1365 self.Satellite = h5.group(1) 

1366 self.Sensor = h5.group(2) 

1367 self.ProcLCode = h5.group(3) 

1368 

1369 # Ground resolution 

1370 re_res_GSD = re.search(r'Pds_PixelSpacing="([\d]*)"', sum_, re.I) 

1371 self.gResolution = float(re_res_GSD.group(1)) if re_res_GSD else self.gResolution 

1372 

1373 # Additional map info + resampling 

1374 try: # no mapinfo or resampling for Level1A 

1375 hx = re.search(r'Ps_(ResamplingMethod="[\S]*")[\s]*' 

1376 r'Pds_(UTM_ZoneNo="[\S]*")[\s]*' 

1377 r'Pds_(MapDirection="[\S]*")[\s]*' 

1378 r'Pds_(PixelSpacing="[\d]*")', 

1379 sum_, re.I) 

1380 self.additional.append(["Map_Info: " + ', '.join(hx.groups())]) 

1381 except AttributeError: 

1382 pass 

1383 

1384 # SunElevation, SunAzimuth, ViewingAngle =(Img_Pointing_Angle), 

1385 # IncidenceAngle =(Img_SceneCenterAngle), Field of View 

1386 h2 = re.search(r'Img_SunAngleElevation="([\d.]*)"[\s]*' 

1387 r'Img_SunAngleAzimuth="([\d.]*)"[\s]*' 

1388 r'Img_PointingAngle="([\d.]*)"[\s]*' 

1389 r'Img_SceneCenterAngle="([\S]*)"', 

1390 sum_, re.I) 

1391 self.SunElevation = float(h2.group(1)) 

1392 self.SunAzimuth = float(h2.group(2)) 

1393 self.ViewingAngle = float(h2.group(3)) 

1394 incAngle = h2.group(4) 

1395 if incAngle.lower().startswith("l"): # L means incidence angle to the left. Set to negativ values 

1396 self.IncidenceAngle = float("-" + incAngle[1:]) 

1397 elif incAngle.lower().startswith("r"): # R means incidence angle to the right. Set to positiv values 

1398 self.IncidenceAngle = float(incAngle[1:]) 

1399 else: # Sign ("L" or "R") will not be added in case of zero. 

1400 self.IncidenceAngle = float(incAngle) 

1401 self.FOV = get_FieldOfView(self.GMS_identifier) 

1402 

1403 self.ScaleFactor = 1 

1404 self.PhysUnit = "DN" 

1405 

1406 # Additional: Exposure coefficients, Saturation coefficients of each band 

1407 h4 = re.search(r'Img_ExposureOfBand1="([\d.-]*)"[\s]*' 

1408 r'Img_ExposureOfBand2="([\d.-]*)"[\s]*' 

1409 r'Img_ExposureOfBand3="([\d.-]*)"[\s]*' 

1410 r'Img_ExposureOfBand4="([\d.-]*)"[\s]*' 

1411 r'Img_SaturationLevelOfBand1="([\d.-]*)"[\s]*' 

1412 r'Img_SaturationLevelOfBand2="([\d.-]*)"[\s]*' 

1413 r'Img_SaturationLevelOfBand3="([\d.-]*)"[\s]*' 

1414 r'Img_SaturationLevelOfBand4="([\d.-]*)"', 

1415 sum_, re.I) 

1416 

1417 self.additional.append(["Exposure Coeffients: B1:'%s',B2:'%s',B3:'%s',B4:'%s'" 

1418 % (h4.group(1), h4.group(2), h4.group(3), h4.group(4))]) 

1419 self.additional.append(["Saturation coeffients: B1:'%s',B2:'%s',B3:'%s',B4:'%s'" 

1420 % (h4.group(5), h4.group(6), h4.group(7), h4.group(8))]) 

1421 

1422 # AcqDate + AcqTime 

1423 h6 = re.search(r'Lbi_ObservationDate="([\d]*)"', sum_, re.I) 

1424 AcqDate = "%s-%s-%s" % (h6.group(1)[:4], h6.group(1)[4:6], h6.group(1)[6:8]) 

1425 re_AcqTime = re.search(r'Img_SceneCenterDateTime="[0-9]* ([0-9][0-9]:[0-9][0-9]:.*)"', sum_, re.I) 

1426 AcqTime_dec = re_AcqTime.group(1) if re_AcqTime else self.AcqTime 

1427 AcqTime = AcqTime_dec.split('.')[0] 

1428 

1429 # set self.AcqDateTime as well as self.AcqDate and self.AcqTime 

1430 if AcqTime_dec: 

1431 self.AcqDateTime = datetime.datetime.strptime( 

1432 '%s %s%s' % (AcqDate, AcqTime_dec, '+0000'), '%Y-%m-%d %H:%M:%S.%f%z') 

1433 else: 

1434 self.AcqDateTime = datetime.datetime.strptime('%s%s' % (AcqDate, '+0000'), '%Y-%m-%d%z') 

1435 self.AcqTime = AcqTime 

1436 

1437 # Earth sun distance 

1438 self.EarthSunDist = self.get_EarthSunDistance(self.AcqDate) 

1439 

1440 # Quality (clouds) 

1441 # cloudsperc = {'0': "0 to 2%", '1': "3 to 10%", '1': "11 to 20%", '1': "21 to 30%", '1': "31 to 40%", 

1442 # '1': "41 to 50%", 

1443 # '1': "51 to 60%", '1': "61 to 70%", '1': "71 to 80%", '1': "81 to 90%", '1': "91 to 100%", 

1444 # '99': "No assessment"} FIXME F601 dictionary key '1' repeated with different values 

1445 # h7 = re.search(r'Img_CloudQuantityOfAllImage="([\d]*)"', sum_, re.I) 

1446 # self.Quality.append(["CloudCoverage: " + cloudsperc[h7.group(1)]]) 

1447 

1448 # Corner Coordinates ## Lon/Lat for all given image corners UL,UR,LL,LR (tuples) 

1449 h10_UL = re.findall(r'Img_[\s\S]*SceneLeftTopLatitude="(.*)"[\s\S]' 

1450 r'Img_[\s\S]*SceneLeftTopLongitude="(.*)"', 

1451 sum_, re.I) 

1452 h10_UR = re.findall(r'Img_[\s\S]*SceneRightTopLatitude="(.*)"[\s\S]' 

1453 r'Img_[\s\S]*SceneRightTopLongitude="(.*)"', 

1454 sum_, re.I) 

1455 h10_LL = re.findall(r'Img_[\s\S]*SceneLeftBottomLatitude="(.*)"[\s\S]' 

1456 r'Img_[\s\S]*SceneLeftBottomLongitude="(.*)"', 

1457 sum_, re.I) 

1458 h10_LR = re.findall( 

1459 r'Img_[\s\S]*SceneRightBottomLatitude="(.*)"[\s\S]' 

1460 r'Img_[\s\S]*SceneRightBottomLongitude="(.*)"', 

1461 sum_, re.I) 

1462 if h10_UL: # Set Corner Tie Point Coordinates (order = UL,UR,LL,LR) 

1463 self.CornerTieP_LonLat.append(tuple([float(h10_UL[0][1]), float(h10_UL[0][0])])) 

1464 self.CornerTieP_LonLat.append(tuple([float(h10_UR[0][1]), float(h10_UR[0][0])])) 

1465 self.CornerTieP_LonLat.append(tuple([float(h10_LL[0][1]), float(h10_LL[0][0])])) 

1466 self.CornerTieP_LonLat.append(tuple([float(h10_LR[0][1]), float(h10_LR[0][0])])) 

1467 

1468 # Coordinate Reference System 

1469 re_res_CS_UTM_ZONE = re.search(r'Pds_UTM_ZoneNo="([0-9][0-9]?)"', sum_, re.I) 

1470 self.CS_UTM_ZONE = int(re_res_CS_UTM_ZONE.group(1)) if re_res_CS_UTM_ZONE else self.CS_UTM_ZONE 

1471 if self.CS_UTM_ZONE != -99.: 

1472 self.CS_TYPE = 'UTM' 

1473 self.CS_EPSG = self.CS_EPSG if self.CornerTieP_LonLat == [] else int('326' + str(self.CS_UTM_ZONE)) \ 

1474 if self.CornerTieP_LonLat[0][1] > 0.0 else int('327' + str(self.CS_UTM_ZONE)) 

1475 self.CS_DATUM = 'WGS84' 

1476 else: 

1477 self.CS_TYPE = 'LonLat' if self.CornerTieP_LonLat != [] and -180. <= self.CornerTieP_LonLat[0][0] <= 180. \ 

1478 and -90. <= self.CornerTieP_LonLat[0][1] <= 90. else self.CS_TYPE 

1479 self.CS_EPSG = 4326 if self.CS_TYPE == 'LonLat' else self.CS_TYPE 

1480 self.CS_DATUM = 'WGS84' if self.CS_TYPE == 'LonLat' else self.CS_DATUM 

1481 

1482 ########################## 

1483 # band specific metadata # 

1484 ########################## 

1485 

1486 LBA_full_sorted = natsorted(self.LayerBandsAssignment_full) 

1487 

1488 # GainMode with corresponding coefficients + Offsets 

1489 gains_AVNIR = {'1': ['N/A', 'N/A', 'N/A', 'N/A'], '2': [0.5880, 0.5730, 0.5020, 0.5570], 

1490 '3': [0.9410, 0.9140, 1.2040, 0.835], '4': [1.412, 1.373, 'N/A', 0.8350]} 

1491 # gains_PRISM = {'1': ['N/A', 'N/A', 'N/A', 'N/A'], '2': [0.501, 'N/A', 'N/A', 'N/A'], 

1492 # '3': ['N/A', 'N/A', 'N/A', 'N/A'], '4': ['N/A', 'N/A', 'N/A', 'N/A']} 

1493 

1494 h3 = re.search(r'Img_GainModeBand1="([1-4])"[\s]*' 

1495 r'Img_GainModeBand2="([1-4])"[\s]*' 

1496 r'Img_GainModeBand3="([1-4])"' 

1497 r'[\s]*Img_GainModeBand4="([1-4])"', 

1498 sum_, re.I) 

1499 self.additional.append( 

1500 ["GainMode: B1:'%s',B2:'%s',B3:'%s',B4:'%s'" % (h3.group(1), h3.group(2), h3.group(3), h3.group(4))]) 

1501 self.Gains = dict(zip(LBA_full_sorted, [gains_AVNIR[h3.group(1)][0], gains_AVNIR[h3.group(2)][1], 

1502 gains_AVNIR[h3.group(3)][2], gains_AVNIR[h3.group(4)][3]])) 

1503 self.Offsets = dict(zip(LBA_full_sorted, [0, 0, 0, 0])) # only gains are required for DN2radiance calculation 

1504 

1505 # Solar irradiance, central wavelengths, full width half maximum per band 

1506 self.wvlUnit = 'Nanometers' 

1507 # derive number of bands from number of given gains/offsets in header file or from raster data 

1508 # noinspection PyBroadException 

1509 try: 

1510 self.nBands = (np.mean([len(self.Gains), len(self.Offsets)]) 

1511 if np.std([len(self.Gains), len(self.Offsets)]) == 0 and len(self.Gains) != 0 else 

1512 GEOP.GEOPROCESSING(self.Dataname, self.logger).bands) 

1513 except Exception: 

1514 self.logger.warning('Unable to get number of bands for the dataset %s Provider values are used for ' 

1515 'solar irradiation, CWL and FWHM!.' % self.Dataname) 

1516 self.LayerBandsAssignment = get_LayerBandsAssignment(self.GMS_identifier, self.nBands) 

1517 

1518 # Exact values calculated based in SRFs 

1519 self.SolIrradiance, self.CWL, self.FWHM = self.calc_solar_irradiance_CWL_FWHM_per_band() 

1520 # Provider values 

1521 if not self.SolIrradiance: 

1522 # Preconfigured Irradiation values 

1523 # = Thullier. source: https://earth.esa.int/c/document_library/get_file?folderId=19584&name=DLFE-262.pdf 

1524 self.SolIrradiance = dict(zip(LBA_full_sorted, [1943.3, 1813.7, 1562.3, 1076.5])) 

1525 # Preconfigured CWLs 

1526 self.CWL = dict(zip(LBA_full_sorted, [460., 560., 650., 825.])) 

1527 # http://www.isprs-ann-photogramm-remote-sens-spatial-inf-sci.net/I-7/291/2012/isprsannals-I-7-291-2012.pdf 

1528 self.FWHM = dict(zip(LBA_full_sorted, [80., 80., 80., 130.])) 

1529 

1530 self.filter_layerdependent_metadata() 

1531 self.orbitParams = get_orbit_params(self.GMS_identifier) 

1532 self.spec_vals = get_special_values(self.GMS_identifier) 

1533 

1534 def Read_ALOS_LEADER(self): 

1535 """Read metadata from ALOS leader file. binary. 

1536 

1537 For exact information content see: 

1538 file:///misc/ro2/behling/Satelliten/ALOS/doc/ALOS Product Format description.pdf 

1539 """ 

1540 

1541 # self.default_attr() 

1542 if os.path.isdir(self.FolderOrArchive): 

1543 glob_res = glob.glob(os.path.join(self.FolderOrArchive, '*/*data*/LED-*')) 

1544 assert len(glob_res) > 0, 'No LED* file can be found in %s/*/*data*/!' % self.FolderOrArchive 

1545 self.Metafile = glob_res[0] 

1546 lead_ = open(self.Metafile, "rb").read() 

1547 else: # archive 

1548 lead_, self.Metafile = \ 

1549 open_specific_file_within_archive(self.FolderOrArchive, '*/*data*/LED-*', read_mode='rb') 

1550 

1551 # Gains & offsets 

1552 LBA_full_sorted = natsorted(self.LayerBandsAssignment_full) 

1553 self.Gains = dict(zip(LBA_full_sorted, [float(lead_[4680 * 3 + 2703:4680 * 3 + 2710]), 

1554 float(lead_[4680 * 3 + 2719:4680 * 3 + 2726]), 

1555 float(lead_[4680 * 3 + 2735:4680 * 3 + 2742]), 

1556 float(lead_[4680 * 3 + 2751:4680 * 3 + 2758])])) 

1557 self.Offsets = dict(zip(LBA_full_sorted, [float(lead_[4680 * 3 + 2711:4680 * 3 + 2718]), 

1558 float(lead_[4680 * 3 + 2727:4680 * 3 + 2734]), 

1559 float(lead_[4680 * 3 + 2743:4680 * 3 + 2750]), 

1560 float(lead_[4680 * 3 + 2759:4680 * 3 + 2766])])) 

1561 

1562 def Read_Sentinel2_xmls(self): 

1563 """Read metadata from Sentinel-2 generic xml and granule xml""" 

1564 

1565 # query database to get entityid 

1566 assert self.SceneID and self.SceneID != -9999, "Read_Sentinel2_xmls(): Missing scene ID. " 

1567 res = DB_T.get_info_from_postgreSQLdb(CFG.conn_database, 'scenes', ['entityid'], {'id': self.SceneID}) 

1568 assert len(res) != 0, \ 

1569 "Invalid SceneID given - no corresponding scene with the ID=%s found in database.\n" % self.SceneID 

1570 assert len(res) == 1, "Error in database. The sceneid %s exists more than once. \n" % self.SceneID 

1571 

1572 S2AB = 'S2A' if re.search(r"Sentinel-2A", self.Satellite, re.I) else 'S2B' 

1573 S2ABgranuleID = res[0][0] 

1574 

1575 ################# 

1576 # get XML roots # 

1577 ################# 

1578 

1579 if os.path.isdir(self.FolderOrArchive): 

1580 # metadata has to be read from folder 

1581 ##################################### 

1582 

1583 # get xml_Scene_root (contains scene metadata) 

1584 glob_res = glob.glob(os.path.join(self.FolderOrArchive, '%s*.xml' % S2AB)) 

1585 if not glob_res: 

1586 # new style packaging 

1587 glob_res = glob.glob(os.path.join(self.FolderOrArchive, 'MTD*.xml')) 

1588 assert len(glob_res) > 0, 'No %s*.xml or MTD*.xml file can be found in %s/!' % (S2AB, self.FolderOrArchive) 

1589 self.Metafile = glob_res[0] 

1590 xml_Scene_root = ET.parse(glob_res[0]).getroot() # xml_parser from file 

1591 

1592 # get xml_GR_root (contains granule metadata) 

1593 glob_res = glob.glob(os.path.join(self.FolderOrArchive, 'GRANULE/' + S2ABgranuleID + '/%s*.xml' % S2AB)) 

1594 if not glob_res: 

1595 # new style packaging 

1596 glob_res = glob.glob(os.path.join(self.FolderOrArchive, 'GRANULE/' + S2ABgranuleID + '/MTD*.xml')) 

1597 assert len(glob_res) > 0, \ 

1598 'No /GRANULE/<%sgranuleID>/%s*.xml or MTD*.xml file can be found in %s/!' \ 

1599 % (S2AB, S2AB, self.FolderOrArchive) 

1600 self.Metafile = self.Metafile + ", " + glob_res[0] 

1601 xml_GR_root = ET.parse(glob_res).getroot() # xml_parser from file 

1602 

1603 else: 

1604 # metadata has to be read from within archive 

1605 ############################################# 

1606 

1607 # get xml_Scene_root and xml_GR_root (contain scene and granule metadata) 

1608 try: 

1609 # old style packaging 

1610 xml_SC_str_, self.Metafile = \ 

1611 open_specific_file_within_archive(self.FolderOrArchive, '*.SAFE/%s*.xml' % S2AB) 

1612 xml_GR_str_, Metafile_ = \ 

1613 open_specific_file_within_archive(self.FolderOrArchive, '*.SAFE/GRANULE/' + 

1614 S2ABgranuleID + '/%s*.xml' % S2AB) 

1615 except RuntimeError: # wrong matching expression 

1616 # new style packaging 

1617 xml_SC_str_, self.Metafile = open_specific_file_within_archive(self.FolderOrArchive, '*.SAFE/MTD*.xml') 

1618 xml_GR_str_, Metafile_ = open_specific_file_within_archive(self.FolderOrArchive, '*.SAFE/GRANULE/' + 

1619 S2ABgranuleID + '/MTD*.xml') 

1620 

1621 xml_Scene_root = ET.fromstring(xml_SC_str_) 

1622 xml_GR_root = ET.fromstring(xml_GR_str_) 

1623 self.Metafile = self.Metafile + ", " + Metafile_ 

1624 

1625 ################################### 

1626 # EXTRACT METADATA FROM XML ROOTS # 

1627 ################################### 

1628 

1629 # define Sentinel 2 metadata (hard coded) 

1630 self.Sensor = "MSI" 

1631 

1632 # extract metadata from xml_Scene_root # 

1633 ######################################## 

1634 

1635 # get the namespace within the xml_Scene_root 

1636 m = re.match(r'{(.*)\}', xml_Scene_root.tag) # extract namespace from "{https://....}Level-1C_Tile_ID" 

1637 assert m, 'XML Namespace could not be identified within Sentinel-2 metadata XML file.' 

1638 namespace = m.group(1) 

1639 

1640 self.EntityID = \ 

1641 xml_Scene_root.find(".//Datatake").attrib['datatakeIdentifier'] # FIXME tileID (Granule) oder scene ID??? 

1642 self.Satellite = xml_Scene_root.find(".//SPACECRAFT_NAME").text 

1643 

1644 # AcqDate, AcqTime 

1645 # NOTE ths corresponds to the while data take, not to the granule, BUT has th same value as the one from granule 

1646 # DateTime = xml_Scene_root.find(".//PRODUCT_START_TIME").text 

1647 # ''' PRODUCT_START_TIME = Sensing Time of the first line of the first scene in the product. 

1648 # Alternative: 'DATATAKE_SENSING_START': Sensing start time of the Datatake''' 

1649 # self.AcqDate = DateTime[:10] 

1650 # self.AcqTime = DateTime[11:19] 

1651 

1652 self.ScaleFactor = int(xml_Scene_root.find(".//QUANTIFICATION_VALUE").text) 

1653 self.PhysUnit = "TOA_Reflectance in [0-%d]" % self.ScaleFactor 

1654 

1655 # Flight direction 

1656 Fdir = {'ASCENDING': "Ascending", 'DESCENDING': "Descending"} 

1657 self.additional.append(["Flight Direction", Fdir[xml_Scene_root.find(".//SENSING_ORBIT_DIRECTION").text]]) 

1658 

1659 self.ProcLCode = xml_Scene_root.find(".//PROCESSING_LEVEL").text 

1660 

1661 # extract metadata from xml_GR_root # 

1662 ##################################### 

1663 

1664 # get the namespace within the xml_GR_root 

1665 m = re.match(r'{(.*)\}', xml_GR_root.tag) # extract namespace from "{https://....}Level-1C_Tile_ID" 

1666 assert m, 'XML Namespace could not be identified within metadata XML file.' 

1667 namespace = m.group(1) # e.g., "https://psd-12.sentinel2.eo.esa.int/PSD/S2_PDI_Level-1C_Tile_Metadata.xsd" 

1668 

1669 # set self.AcqDateTime as well as self.AcqDate and self.AcqTime 

1670 self.AcqDateTime = iso8601.parse_date(xml_GR_root.find(".//SENSING_TIME").text) 

1671 

1672 # SunAngles 

1673 self.SunElevation = 90 - float(xml_GR_root.find(".//Mean_Sun_Angle/ZENITH_ANGLE").text) # mean angle of granule 

1674 self.SunAzimuth = float(xml_GR_root.find(".//Mean_Sun_Angle/AZIMUTH_ANGLE").text) # mean angle of granule 

1675 

1676 # coordinate system 

1677 geo_codings = HLP_F.find_in_xml_root(namespace, xml_GR_root, 'Geometric_Info', "Tile_Geocoding") 

1678 self.CS_EPSG = int(geo_codings.find(".//HORIZONTAL_CS_CODE").text.split(":")[1]) 

1679 CooSys = geo_codings.find(".//HORIZONTAL_CS_NAME").text 

1680 self.CS_DATUM = "WGS84" if re.search(r"wgs84", CooSys, re.I).group(0) else self.CS_DATUM 

1681 self.CS_TYPE = "UTM" if re.search(r"utm", CooSys, re.I).group(0) else self.CS_TYPE 

1682 if self.CS_TYPE == "UTM": 

1683 tmp = CooSys.split(" ")[-1] 

1684 self.CS_UTM_ZONE = \ 

1685 int(tmp[:-1]) if tmp[-1] == 'N' else \ 

1686 -int(tmp[:-1]) if tmp[-1] == 'S' else self.CS_UTM_ZONE 

1687 

1688 # corner coords 

1689 subsytem_Res_dic = {"%s10" % S2AB: 10, "%s20" % S2AB: 20, "%s60" % S2AB: 60} 

1690 if self.CS_TYPE == 'UTM': 

1691 spatial_samplings = {float(size.get("resolution")): {key: int(size.find(key).text) 

1692 for key in ["NROWS", "NCOLS"]} for size in 

1693 geo_codings.findall("Size")} 

1694 for geo in geo_codings.findall("Geoposition"): 

1695 spatial_samplings[float(geo.get("resolution"))].update( 

1696 {key: float(geo.find(key).text) for key in ["ULX", "ULY", "XDIM", "YDIM"]}) 

1697 

1698 ss_sub = spatial_samplings[subsytem_Res_dic[self.Subsystem]] 

1699 LRX = ss_sub['ULX'] + ss_sub['NCOLS'] * ss_sub['XDIM'] 

1700 LRY = ss_sub['ULY'] + ss_sub['NROWS'] * ss_sub['YDIM'] 

1701 self.CornerTieP_UTM = [(ss_sub['ULX'], ss_sub['ULY']), (LRX, ss_sub['ULY']), 

1702 (ss_sub['ULX'], LRY), (LRX, LRY)] # (x,y) for UL,UR,LL,LR 

1703 

1704 # geometricResolution 

1705 self.gResolution = subsytem_Res_dic[self.Subsystem] 

1706 

1707 # determine metadata from extracted metadata values 

1708 self.EarthSunDist = self.get_EarthSunDistance(self.AcqDate) 

1709 

1710 # Quality flags # FIXME does not work (at least with old data) 

1711 Quality_temp = (xml_Scene_root.find(".//Technical_Quality_Assessment")) 

1712 self.Quality.append(["DEGRADED_ANC_DATA_PERCENTAGE", Quality_temp.find("./DEGRADED_ANC_DATA_PERCENTAGE").text]) 

1713 self.Quality.append(["DEGRADED_MSI_DATA_PERCENTAGE", Quality_temp.find("./DEGRADED_MSI_DATA_PERCENTAGE").text]) 

1714 Quality_temp2 = xml_Scene_root.find(".//Quality_Inspections") 

1715 quality_flags = ["SENSOR_QUALITY_FLAG", "GEOMETRIC_QUALITY_FLAG", "GENERAL_QUALITY_FLAG", 

1716 "FORMAT_CORRECTNESS_FLAG", "RADIOMETRIC_QUALITY_FLAG"] 

1717 

1718 try: 

1719 # layout example: <SENSOR_QUALITY_FLAG>PASSED</SENSOR_QUALITY_FLAG> 

1720 for ql in quality_flags: 

1721 self.Quality.append([ql, Quality_temp2.find("./" + ql).text]) 

1722 except AttributeError: 

1723 # since ~11/2017 the quality checks layout in the XML has changed: 

1724 # layout example: <quality_check checkType="SENSOR_QUALITY">PASSED</quality_check> 

1725 elements = Quality_temp2.findall('quality_check') 

1726 checkTypeValDict = {ele.attrib['checkType']: ele.text for ele in elements} 

1727 for ql in quality_flags: 

1728 self.Quality.append([ql, checkTypeValDict[ql.split('_FLAG')[0]]]) 

1729 

1730 ########################## 

1731 # band specific metadata # 

1732 ########################## 

1733 

1734 LBA_full_sorted = natsorted(self.LayerBandsAssignment_full) 

1735 

1736 # ATTENTION Gains are only provided for 12 bands! I don't know why? 

1737 Gains = [float(ele.text) for ele in xml_Scene_root.findall(".//PHYSICAL_GAINS")] 

1738 Gains = Gains if len(Gains) == 13 else [1] + Gains 

1739 self.Gains = dict(zip(LBA_full_sorted, Gains)) 

1740 # FIXME assuming that the first band at 443nm has been left out here IS POSSIBLY WRONG 

1741 # FIXME (could also be band 8A oder band 9 (water vapour)) 

1742 

1743 # Solar irradiance, central wavelengths, full width half maximum per band 

1744 self.wvlUnit = 'Nanometers' 

1745 self.LayerBandsAssignment = get_LayerBandsAssignment(self.GMS_identifier) 

1746 

1747 # Exact values calculated based in SRFs 

1748 self.SolIrradiance, self.CWL, self.FWHM = self.calc_solar_irradiance_CWL_FWHM_per_band() 

1749 # Provider values 

1750 if not self.SolIrradiance: 

1751 # get from xml file 

1752 self.SolIrradiance = \ 

1753 dict(zip(LBA_full_sorted, 

1754 [float(ele.text) for ele in 

1755 xml_Scene_root.find(".//Solar_Irradiance_List").findall("SOLAR_IRRADIANCE")])) 

1756 # Preconfigured CWLs 

1757 self.CWL = dict(zip(LBA_full_sorted, 

1758 [float(ele.text) for ele in 

1759 xml_Scene_root.find(".//Spectral_Information_List").findall(".//CENTRAL")])) 

1760 

1761 # SensorAngles 

1762 meta_temp = {} 

1763 branch = HLP_F.find_in_xml_root(namespace, xml_GR_root, *("Geometric_Info", "Tile_Angles")) 

1764 meta_temp["bandId2bandName"] = {int(ele.get("bandId")): ele.text.split("_")[-2] for ele in 

1765 xml_GR_root.findall(".//MASK_FILENAME") if ele.get("bandId") is not None} 

1766 meta_temp["bandName2bandId"] = {bandName: bandId for bandId, bandName in 

1767 meta_temp["bandId2bandName"].items()} 

1768 meta_temp["bandIds"] = sorted(list(meta_temp["bandId2bandName"].keys())) 

1769 meta_temp["viewing_zenith_detectors"] = \ 

1770 {bandId: {bf.get("detectorId"): HLP_F.get_values_from_xml(HLP_F.find_in_xml(bf, *("Zenith", "Values_List"))) 

1771 for bf in branch.findall("Viewing_Incidence_Angles_Grids[@bandId='%i']" % bandId)} 

1772 for bandId in meta_temp["bandIds"]} 

1773 meta_temp["viewing_zenith"] = HLP_F.stack_detectors(meta_temp["viewing_zenith_detectors"]) 

1774 

1775 meta_temp["viewing_azimuth_detectors"] = {bandId: {bf.get("detectorId"): HLP_F.get_values_from_xml( 

1776 HLP_F.find_in_xml(bf, *("Azimuth", "Values_List"))) for bf in branch.findall( 

1777 "Viewing_Incidence_Angles_Grids[@bandId='%i']" % bandId)} for bandId in meta_temp["bandIds"]} 

1778 meta_temp["viewing_azimuth"] = HLP_F.stack_detectors(meta_temp["viewing_azimuth_detectors"]) 

1779 # mean values of all mean band values # FIXME 

1780 self.ViewingAngle = np.mean([np.nanmean(i) for i in meta_temp['viewing_zenith'].values()]) 

1781 self.IncidenceAngle = np.mean([np.nanmean(i) for i in meta_temp['viewing_azimuth'].values()]) 

1782 # 5000m step arrays per band. dict: key = bandindex, value = stacked_array for all detectors 

1783 

1784 self.ViewingAngle_arrProv = meta_temp['viewing_zenith'] 

1785 self.IncidenceAngle_arrProv = meta_temp['viewing_azimuth'] 

1786 LBA2Id_dic = {'1': 0, '2': 1, '3': 2, '4': 3, '5': 4, '6': 5, '7': 6, '8': 7, '8A': 8, '9': 9, '10': 10, 

1787 '11': 11, '12': 12} 

1788 

1789 def filter_dic(AngleArr): return {LBAn: AngleArr[LBA2Id_dic[LBAn]] for LBAn in self.LayerBandsAssignment} 

1790 

1791 self.ViewingAngle_arrProv = filter_dic(self.ViewingAngle_arrProv) 

1792 self.IncidenceAngle_arrProv = filter_dic(self.IncidenceAngle_arrProv) 

1793 

1794 self.FOV = get_FieldOfView(self.GMS_identifier) 

1795 self.orbitParams = get_orbit_params(self.GMS_identifier) 

1796 self.filter_layerdependent_metadata() 

1797 self.spec_vals = get_special_values(self.GMS_identifier) 

1798 

1799 def add_rasObj_dims_projection_physUnit(self, rasObj, dict_LayerOptTherm, temp_logger=None): 

1800 self.rows = rasObj.rows 

1801 self.cols = rasObj.cols 

1802 self.bands = rasObj.bands 

1803 if rasObj.projection != '': 

1804 self.map_info = geotransform2mapinfo(rasObj.geotransform, rasObj.projection) 

1805 self.projection = rasObj.projection 

1806 self.gResolution = abs(rasObj.geotransform[1]) 

1807 self.CS_EPSG = WKT2EPSG(rasObj.projection) 

1808 

1809 dict_conv_physUnit = {'Rad': "W * m-2 * sr-1 * micrometer-1", 

1810 'TOA_Ref': 'TOA_Reflectance in [0-%d]' % CFG.scale_factor_TOARef, 

1811 'BOA_Ref': 'BOA_Reflectance in [0-%d]' % CFG.scale_factor_BOARef, 

1812 'Temp': 'Degrees Celsius with scale factor = 100'} 

1813 if list(set(dict_LayerOptTherm.values())) == ['optical']: 

1814 self.PhysUnit = dict_conv_physUnit[CFG.target_radunit_optical] 

1815 elif list(set(dict_LayerOptTherm.values())) == ['thermal']: 

1816 self.PhysUnit = dict_conv_physUnit[CFG.target_radunit_thermal] 

1817 elif sorted(list(set(dict_LayerOptTherm.values()))) == ['optical', 'thermal']: 

1818 self.PhysUnit = ['Optical bands: %s' % dict_conv_physUnit[CFG.target_radunit_optical], 

1819 'Thermal bands: %s' % dict_conv_physUnit[CFG.target_radunit_thermal]] 

1820 else: 

1821 logger = self.logger if hasattr(self, 'logger') else temp_logger 

1822 assert logger, "ERROR: Physical unit could not be determined due to unexpected 'dict_LayerOptTherm'. " \ 

1823 "Got %s." % dict_LayerOptTherm 

1824 logger.error("Physical unit could not be determined due to unexpected 'dict_LayerOptTherm'. Got %s." 

1825 % dict_LayerOptTherm) 

1826 

1827 def to_odict(self): 

1828 # type: () -> collections.OrderedDict 

1829 """Creates an OrderedDict containing selected attribute of the METADATA object that will later be included in 

1830 ENVI file headers in the same order. 

1831 

1832 """ 

1833 # FIXME orbit params are missing 

1834 # descr_dic = dict( # FillZeroSaturated von HLP_F ausgeben lassen 

1835 # ALOS_Rad= 

1836 # "(1) GEOCODED Level1B2 Product; '" + self.Dataname + "'\n " 

1837 # "(2) Int16 RadianceData in [W * m-2 * sr-1 * micrometer-1]*10; " 

1838 # "radiance scale factor: 10 (fillPixels: -99, zeroPixels:0, saturatedPixels: 32767 (Max of Int16))'", 

1839 # ALOS_Ref= 

1840 # "(1) Orthorectified JAXA or GFZ; '" + self.Dataname + "'\n " 

1841 # "(2) Int16 TOA_Reflectance in [0-10000]; " 

1842 # "radiance scale factor: 10 (fillPixels: -99, zeroPixels:0, saturatedPixels: 32767 (Max of Int16))'", 

1843 # Terra_Rad= 

1844 # "(1) Orthorectified JAXA or GFZ; '" + self.Dataname + "'\n " 

1845 # "(2) Int16 RadianceData in [W * m-2 * sr-1 * micrometer-1]*10; " 

1846 # "radiance scale factor: 10 (fillPixels: -99, zeroPixels:0, saturatedPixels: 32767 (Max of Int16))'", 

1847 # Terra_Ref= 

1848 # "(1) Orthorectified JAXA or GFZ; '" + self.Dataname + "'\n " 

1849 # "(2) Int16 TOA_Reflectance in [0-10000]; " 

1850 # "radiance scale factor: 10 (fillPixels: -99, zeroPixels:0, saturatedPixels: 32767 (Max of Int16))'", 

1851 # Landsat_Rad= 

1852 # "(1) Landsat DN: "+self.ProcLCode+ " Product; '"+self.Dataname+"'\n " 

1853 # "(2) Int16 RadianceData in [W * m-2 * sr-1 * micrometer-1]*10; " 

1854 # "radiance scale factor: 10 (fillPixels: -99, zeroPixels:0, saturatedPixels: 32767 (Max of Int16))'", 

1855 # Landsat_Ref= 

1856 # "(1) Landsat DN: "+self.ProcLCode + " Product; '" + self.Dataname + "'\n " 

1857 # "(2) Int16 TOA_Reflectance in [0-10000] " 

1858 # "(fillPixels: -99, zeroPixels:0, saturatedPixels: 32767 (Max of Int16))'", 

1859 # RapidEye_Rad= 

1860 # "(1) Ortho Level3A01 Product; '"+self.Dataname+"'\n " 

1861 # "(2) Int16 RadianceData in [W * m-2 * sr-1 * micrometer-1]*10; " 

1862 # "radiance scale factor: 10 (fillPixels: -99, zeroPixels:0, saturatedPixels: 32767 (Max of Int16))'", 

1863 # RapidEye_Ref= 

1864 # "(1) Ortho Level3A01 Product; '" + self.Dataname + "'\n " 

1865 # "(2) Int16 TOA_Reflectance in [0-10000] " 

1866 # "(fillPixels: -99, zeroPixels:0, saturatedPixels: 32767 (Max of Int16))'", 

1867 # SPOT_Rad= 

1868 # "(1) Ortho Level3A01 Product; '" + self.Dataname + "'\n " 

1869 # "(2) Int16 RadianceData in [W * m-2 * sr-1 * micrometer-1]*10; " 

1870 # "radiance scale factor: 10 (fillPixels: -99, zeroPixels:0, saturatedPixels: 32767 (Max of Int16))'", 

1871 # SPOT_Ref= 

1872 # "(1) Ortho Level3A01 Product; '" + self.Dataname + "'\n " 

1873 # "(2) Int16 TOA_Reflectance in [0-10000] " 

1874 # "(fillPixels: -99, zeroPixels:0, saturatedPixels: 32767 (Max of Int16))'") 

1875 

1876 # copy directly compatible keys 

1877 Meta = collections.OrderedDict() 

1878 # Meta['description'] = descr_dic[self.Satellite + '_' + CFG.target_radunit_optical] 

1879 

1880 for odictKey in enviHdr_keyOrder: 

1881 if odictKey in map_odictKeys_objAttrnames: 

1882 attrVal = getattr(self, map_odictKeys_objAttrnames[odictKey]) 

1883 if attrVal and map_odictKeys_objAttrnames[odictKey] in layerdependent_metadata: 

1884 # convert band specific metadata dicts to lists in the order of LayerBandsAssignment 

1885 Meta.update({odictKey: [attrVal[band] for band in self.LayerBandsAssignment if band in attrVal]}) 

1886 else: 

1887 Meta.update({odictKey: attrVal}) 

1888 elif odictKey == 'map info' and self.map_info: 

1889 Meta['map info'] = self.map_info 

1890 elif odictKey == 'coordinate system string' and self.projection: 

1891 Meta['coordinate system string'] = self.projection 

1892 elif odictKey == 'data ignore value': 

1893 Meta['data ignore value'] = self.spec_vals['fill'] if 'fill' in self.spec_vals else None 

1894 elif odictKey == 'corner coordinates lonlat': 

1895 Meta['corner coordinates lonlat'] = str(self.CornerTieP_LonLat).replace('(', '[').replace(')', ']') 

1896 elif odictKey == 'wavelength units': 

1897 Meta['wavelength units'] = "Nanometers" 

1898 elif odictKey == 'band names': 

1899 Meta['band names'] = self.bandnames 

1900 

1901 # add keys that will not be included into ENVI header 

1902 if self.ViewingAngle_arrProv is not None: 

1903 Meta['ViewingAngle_arrProv'] = {k: v.tolist() for k, v in self.ViewingAngle_arrProv.items()} 

1904 Meta['IncidenceAngle'] = self.IncidenceAngle 

1905 if self.IncidenceAngle_arrProv is not None: 

1906 Meta['IncidenceAngle_arrProv'] = {k: v.tolist() for k, v in self.IncidenceAngle_arrProv.items()} 

1907 

1908 return Meta 

1909 

1910 def from_odict(self, odict): 

1911 # type: (collections.OrderedDict) -> METADATA 

1912 

1913 # copy directly compatible keys 

1914 # [setattr(self, attrN, odict[odictKey]) for odictKey, attrN in map_odictKeys_objAttrnames.items()] 

1915 

1916 for odictKey, attrN in map_odictKeys_objAttrnames.items(): 

1917 if attrN not in layerdependent_metadata: 

1918 setattr(self, attrN, odict[odictKey]) 

1919 else: 

1920 # convert band specific metadata to dicts 

1921 setattr(self, attrN, dict(zip(odict['LayerBandsAssignment'], odict[odictKey]))) 

1922 

1923 # set the remaining attributes 

1924 if 'map info' in odict: 

1925 self.map_info = odict['map info'] 

1926 if 'coordinate system string' in odict: 

1927 self.projection = odict['coordinate system string'] 

1928 if 'data ignore value' in odict: 

1929 self.spec_vals['fill'] = odict['data ignore value'] 

1930 if 'ViewingAngle_arrProv' in odict: 

1931 self.ViewingAngle_arrProv = {bN: np.array(odict['ViewingAngle_arrProv'][bN]) 

1932 for bN in self.LayerBandsAssignment if bN in odict['ViewingAngle_arrProv']} 

1933 if 'IncidenceAngle_arrProv' in odict: 

1934 self.IncidenceAngle_arrProv = {bN: np.array(odict['IncidenceAngle_arrProv'][bN]) 

1935 for bN in self.LayerBandsAssignment if bN in odict['IncidenceAngle_arrProv']} 

1936 

1937 return self 

1938 

1939 def calc_solar_irradiance_CWL_FWHM_per_band(self): 

1940 # type: () -> (dict, dict, dict) 

1941 sensorcode = get_GMS_sensorcode(self.GMS_identifier) 

1942 # ms_pan = ('multi' if self.nBands > 1 else 'pan') 

1943 

1944 irr_bands, cwl_bands, fwhm_bands = {}, {}, {} 

1945 

1946 if not sensorcode: 

1947 self.logger.warning('GMS-sensorcode missing. Provider values are used for solar irradiation, CWL and FWHM.') 

1948 else: 

1949 self.logger.info('Calculating solar irradiance, central wavelengths and full width half maxima...') 

1950 

1951 sol_irr = Solar_Irradiance_reader(wvl_min_nm=350, wvl_max_nm=2500) 

1952 srf_dict = SRF_Reader(self.GMS_identifier, no_pan=False, no_thermal=False) 

1953 

1954 for band in srf_dict.keys(): 

1955 if srf_dict[band] is None: 

1956 irr_bands[band] = None 

1957 cwl_bands[band] = None 

1958 fwhm_bands[band] = None 

1959 else: 

1960 WVL_band = (srf_dict[band][:, 0] if 300 < np.max(srf_dict[band][:, 0]) < 15000 else 

1961 srf_dict[band][:, 0] * 1000) # reads wavelengths given in nm and µm 

1962 RSP_band = srf_dict[band][:, 1] 

1963 # sol_irr_at_WVL = \ 

1964 # scipy.interpolate.interp1d(sol_irr[:, 0], sol_irr[:, 1], kind='linear')(WVL_band) # ASTER 

1965 # band 8: ValueError: A value in x_new is above the interpolation range. 

1966 sol_irr_at_WVL = np.interp(WVL_band, sol_irr[:, 0], sol_irr[:, 1], left=0, right=0) 

1967 

1968 irr_bands[band] = round(np.sum(sol_irr_at_WVL * RSP_band) / np.sum(RSP_band), 2) 

1969 cwl_bands[band] = round(np.sum(WVL_band * RSP_band) / np.sum(RSP_band), 2) 

1970 fwhm_bands[band] = round(np.max(WVL_band[RSP_band >= (np.max(RSP_band) / 2.)]) - 

1971 np.min(WVL_band[RSP_band >= (np.max(RSP_band) / 2.)]), 2) 

1972 

1973 return irr_bands, cwl_bands, fwhm_bands 

1974 

1975 def get_EarthSunDistance(self, acqDate): 

1976 """Get earth sun distance (requires file of pre calculated earth sun distance per day) 

1977 

1978 :param acqDate: 

1979 """ 

1980 

1981 if not os.path.exists(CFG.path_earthSunDist): 

1982 self.logger.warning("\n\t WARNING: Earth Sun Distance is assumed to be " 

1983 "1.0 because no database can be found at %s.""" % CFG.path_earthSunDist) 

1984 return 1.0 

1985 if not acqDate: 

1986 self.logger.warning("\n\t WARNING: Earth Sun Distance is assumed to be 1.0 because " 

1987 "acquisition date could not be read from metadata.") 

1988 return 1.0 

1989 

1990 with open(CFG.path_earthSunDist, "r") as EA_dist_f: 

1991 EA_dist_dict = {} 

1992 for line in EA_dist_f: 

1993 date, EA = [item.strip() for item in line.split(",")] 

1994 EA_dist_dict[date] = EA 

1995 

1996 return float(EA_dist_dict[acqDate]) 

1997 

1998 def calc_center_acquisition_time(self, fullSceneCornerLonLat, logger): 

1999 """Calculates a missing center acquistion time using acquisition date, full scene corner coordinates and 

2000 solar azimuth. 

2001 

2002 :param fullSceneCornerLonLat: 

2003 :param logger: 

2004 """ 

2005 assert is_dataset_provided_as_fullScene(self.GMS_identifier) and len(fullSceneCornerLonLat) == 4, \ 

2006 'Center acquisition time can only be computed for datasets provided as full scenes, not for tiles.' 

2007 

2008 ul, lr = fullSceneCornerLonLat[0], fullSceneCornerLonLat[3] 

2009 center_coord = [np.mean([ul[0], lr[0]]), np.mean([ul[1], lr[1]])] 

2010 time0_ord = mdates.date2num( 

2011 datetime.datetime.strptime('%s %s' % (self.AcqDate, '00:00:00'), '%Y-%m-%d %H:%M:%S')) 

2012 time1_ord = mdates.date2num( 

2013 datetime.datetime.strptime('%s %s' % (self.AcqDate, '23:59:59'), '%Y-%m-%d %H:%M:%S')) 

2014 time_stamps_ord = np.linspace(time0_ord, time1_ord, 12 * 60 * 60) 

2015 time_stamps_with_tzinfo = mdates.num2date(time_stamps_ord) 

2016 time_stamps = np.array([time_stamps_with_tzinfo[i].replace(tzinfo=None) 

2017 for i in range(len(time_stamps_with_tzinfo))]) 

2018 sols_az_rad = astronomy.get_alt_az(time_stamps, [center_coord[0]] * time_stamps.size, 

2019 [center_coord[1]] * time_stamps.size)[1] 

2020 sol_azs = 180 * sols_az_rad / math.pi 

2021 diff_az = np.abs(float(self.SunAzimuth) - sol_azs) 

2022 acq_datetime = time_stamps[np.where(diff_az == np.min(diff_az))][0] 

2023 AcqTime = acq_datetime.strftime(format='%H:%M:%S') 

2024 logger.info('Center acquisition time has been calculated: %s' % AcqTime) 

2025 

2026 # update self. 

2027 self.AcqDateTime = datetime.datetime.strptime( 

2028 '%s %s%s' % (self.AcqDate, AcqTime, '.000000+0000'), '%Y-%m-%d %H:%M:%S.%f%z') 

2029 return self.AcqTime 

2030 

2031 def get_overpassDuration_SceneLength(self, fullSceneCornerLonLat, fullSceneCornerPos, shape_fullArr, logger): 

2032 """Calculates duration of image acquisition in seconds. 

2033 

2034 :param fullSceneCornerLonLat: 

2035 :param fullSceneCornerPos: 

2036 :param shape_fullArr: 

2037 :param logger: 

2038 """ 

2039 assert is_dataset_provided_as_fullScene(self.GMS_identifier) and len(fullSceneCornerLonLat) == 4, \ 

2040 'Overpass duration and scene length can only be computed for datasets provided as full scenes, not for ' \ 

2041 'tiles.' 

2042 

2043 # check if current scene is a subset 

2044 assert fullSceneCornerPos != list(([0, 0], [0, shape_fullArr[1] - 1], 

2045 [shape_fullArr[0] - 1, 0], [shape_fullArr[0] - 1, shape_fullArr[1] - 1])),\ 

2046 'Overpass duration and scene length cannot be calculated because the given data represents a subset of ' \ 

2047 'the original scene.' 

2048 

2049 # compute scene length 

2050 orbitAltitudeKm, orbitPeriodMin = self.orbitParams[0], self.orbitParams[2] 

2051 UL, UR, LL, LR = fullSceneCornerLonLat 

2052 geod = pyproj.Geod(ellps='WGS84') 

2053 scene_length = np.mean([geod.inv(UL[0], UL[1], LL[0], LL[1])[2], 

2054 geod.inv(UR[0], UR[1], LR[0], LR[1])[2]]) / 1000 # along-track distance [km] 

2055 logger.info('Calculation of scene length...: %s km' % round(float(scene_length), 2)) 

2056 

2057 # compute overpass duration 

2058 orbitPeriodLength = 2 * math.pi * (6371 + orbitAltitudeKm) 

2059 overpassDurationSec = (scene_length / orbitPeriodLength) * orbitPeriodMin * 60. 

2060 logger.info('Calculation of overpass duration...: %s sec' % round(overpassDurationSec, 2)) 

2061 

2062 return overpassDurationSec, scene_length 

2063 

2064 def filter_layerdependent_metadata(self): 

2065 for attrname in layerdependent_metadata: 

2066 attrVal = getattr(self, attrname) 

2067 

2068 if not attrVal: 

2069 continue 

2070 

2071 if isinstance(attrVal, dict): 

2072 setattr(self, attrname, {bN: attrVal[bN] for bN in self.LayerBandsAssignment}) 

2073 

2074 elif isinstance(attrVal, collections.OrderedDict): 

2075 setattr(self, attrname, collections.OrderedDict((bN, attrVal[bN]) for bN in self.LayerBandsAssignment)) 

2076 

2077 else: 

2078 raise ValueError 

2079 

2080 if attrVal: 

2081 assert len(getattr(self, attrname)) == len(self.LayerBandsAssignment) 

2082 

2083 

2084layerdependent_metadata = ['SolIrradiance', 'CWL', 'FWHM', 'Offsets', 'OffsetsRef', 'Gains', 'GainsRef', 

2085 'ThermalConstK1', 'ThermalConstK2', 'ViewingAngle_arrProv', 'IncidenceAngle_arrProv'] 

2086 

2087 

2088map_odictKeys_objAttrnames = { 

2089 'samples': 'cols', 

2090 'lines': 'rows', 

2091 'bands': 'bands', 

2092 'version_gms_preprocessing': 'version_gms_preprocessing', 

2093 'versionalias_gms_preprocessing': 'versionalias_gms_preprocessing', 

2094 'CS_EPSG': 'CS_EPSG', 

2095 'CS_TYPE': 'CS_TYPE', 

2096 'CS_DATUM': 'CS_DATUM', 

2097 'CS_UTM_ZONE': 'CS_UTM_ZONE', 

2098 'scene length': 'scene_length', 

2099 'wavelength': 'CWL', 

2100 'bandwidths': 'FWHM', 

2101 'LayerBandsAssignment': 'LayerBandsAssignment', 

2102 'data gain values': 'Gains', 

2103 'data offset values': 'Offsets', 

2104 'reflectance gain values': 'GainsRef', 

2105 'reflectance offset values': 'OffsetsRef', 

2106 'Metafile': 'Metafile', 

2107 'Satellite': 'Satellite', 

2108 'Sensor': 'Sensor', 

2109 'Subsystem': 'Subsystem', 

2110 'EntityID': 'EntityID', 

2111 'SceneID': 'SceneID', 

2112 'gResolution': 'gResolution', 

2113 'AcqDate': 'AcqDate', 

2114 'AcqTime': 'AcqTime', 

2115 'overpass duraction sec': 'overpassDurationSec', 

2116 'ProcLCode': 'ProcLCode', 

2117 'SunElevation': 'SunElevation', 

2118 'SunAzimuth': 'SunAzimuth', 

2119 'SolIrradiance': 'SolIrradiance', 

2120 'ThermalConstK1': 'ThermalConstK1', 

2121 'ThermalConstK2': 'ThermalConstK2', 

2122 'EarthSunDist': 'EarthSunDist', 

2123 'ViewingAngle': 'ViewingAngle', 

2124 'IncidenceAngle': 'IncidenceAngle', 

2125 'FieldOfView': 'FOV', 

2126 'PhysUnit': 'PhysUnit', 

2127 'ScaleFactor': 'ScaleFactor', 

2128 'Quality': 'Quality', 

2129 'Additional': 'additional' 

2130} 

2131 

2132 

2133def get_LayerBandsAssignment(GMS_id, nBands=None, sort_by_cwl=None, no_thermal=None, no_pan=None, 

2134 return_fullLBA=False, proc_level=''): 

2135 # type: (GMS_identifier, int, bool, bool, bool, bool, str) -> list 

2136 """Returns LayerBandsAssignment corresponding to given satellite, sensor and subsystem and with respect to 

2137 CFG.sort_bands_by_cwl, CFG.skip_thermal and CFG.skip_pan. 

2138 

2139 :param GMS_id: <dict>, derived from self.get_GMS_identifier() 

2140 NOTE: only if there is an additional key 'proc_level', the processing level will be 

2141 respected. This is needed to get the correct LBA after atm. correction 

2142 :param nBands: should be specified if number of bands differs from standard 

2143 (e.g. SPOT data containing no PAN) 

2144 :param sort_by_cwl: whether to sort the returned bands list by central wavelength position 

2145 (default: CFG.sort_bands_by_cwl) 

2146 :param no_thermal: whether to exclude thermal bands from the returned bands list 

2147 (default: CFG.skip_thermal) 

2148 :param no_pan: whether to exclude panchromatic bands from the returned bands list 

2149 (default: CFG.skip_pan) 

2150 :param return_fullLBA: in case there is a subsystem: 

2151 whether to return LayerBandsAssignment for all bands or for the current subsystem 

2152 :param proc_level: processing level for which the LayerBandsAssignment is returned 

2153 (overrides the proc_level given with GMS_id) 

2154 """ 

2155 # set defaults 

2156 # NOTE: these values cannot be set in function signature because CFG is not present at library import time 

2157 sort_by_cwl = sort_by_cwl if sort_by_cwl is not None else CFG.sort_bands_by_cwl 

2158 no_thermal = no_thermal if no_thermal is not None else CFG.skip_thermal 

2159 no_pan = no_pan if no_pan is not None else CFG.skip_pan 

2160 proc_level = proc_level or GMS_id.proc_level 

2161 

2162 if GMS_id.image_type == 'RSD': 

2163 GMS_sensorcode = get_GMS_sensorcode(GMS_id) 

2164 assert GMS_sensorcode, 'Unable to get Layer Bands Assignment. No valid sensorcode privided (got >None<). ' 

2165 

2166 if return_fullLBA: 

2167 GMS_sensorcode = \ 

2168 'AST_full' if GMS_sensorcode.startswith('AST') else \ 

2169 'S2A_full' if GMS_sensorcode.startswith('S2A') else \ 

2170 'S2B_full' if GMS_sensorcode.startswith('S2B') else GMS_sensorcode 

2171 

2172 dict_LayerBandsAssignment = { 

2173 'AVNIR-2': ['1', '2', '3', '4'], 

2174 'AST_full': ['1', '2', '3N', '3B', '4', '5', '6', '7', '8', '9', '10', '11', '12', '13', '14'], 

2175 'AST_V1': ['1', '2', '3N'], 

2176 'AST_V2': ['3B'], 

2177 'AST_S': ['4', '5', '6', '7', '8', '9'], 

2178 'AST_T': ['10', '11', '12', '13', '14'], 

2179 'TM4': ['1', '2', '3', '4', '5', '6', '7'], 

2180 'TM5': ['1', '2', '3', '4', '5', '6', '7'], 

2181 'TM7': ['1', '2', '3', '4', '5', '6L', '6H', '7', '8'], 

2182 'LDCM': ['1', '2', '3', '4', '5', '6', '7', '8', '9', '10', '11'], 

2183 'SPOT1a': ['1', '2', '3', '4'], 

2184 'SPOT2a': ['1', '2', '3', '4'], 

2185 'SPOT3a': ['1', '2', '3', '4'], 

2186 'SPOT4a': ['1', '2', '3', '4', '5'], 

2187 'SPOT5a': ['1', '2', '3', '4', '5'], 

2188 'SPOT1b': ['1', '2', '3', '4'], 

2189 'SPOT2b': ['1', '2', '3', '4'], 

2190 'SPOT3b': ['1', '2', '3', '4'], 

2191 'SPOT4b': ['1', '2', '3', '4', '5'], 

2192 'SPOT5b': ['1', '2', '3', '4', '5'], 

2193 'RE5': ['1', '2', '3', '4', '5'], 

2194 'S2A_full': ['1', '2', '3', '4', '5', '6', '7', '8', '8A', '9', '10', '11', '12'], 

2195 'S2B_full': ['1', '2', '3', '4', '5', '6', '7', '8', '8A', '9', '10', '11', '12'], 

2196 'S2A10': ['2', '3', '4', '8'], 

2197 'S2A20': ['5', '6', '7', '8A', '11', '12'], 

2198 'S2A60': ['1', '9', '10'], 

2199 'S2B10': ['2', '3', '4', '8'], 

2200 'S2B20': ['5', '6', '7', '8A', '11', '12'], 

2201 'S2B60': ['1', '9', '10'], } 

2202 

2203 dict_cwlSorted_LayerBandsAssignment = { 

2204 'TM4': ['1', '2', '3', '4', '5', '7', '6'], 

2205 'TM5': ['1', '2', '3', '4', '5', '7', '6'], 

2206 'TM7': ['1', '2', '3', '8', '4', '5', '7', '6L', '6H'], 

2207 'LDCM': ['1', '2', '3', '8', '4', '5', '9', '6', '7', '10', '11'], 

2208 'SPOT1a': ['1', '4', '2', '3'], 

2209 'SPOT2a': ['1', '4', '2', '3'], 

2210 'SPOT3a': ['1', '4', '2', '3'], 

2211 'SPOT4a': ['1', '5', '2', '3', '4'], 

2212 'SPOT5a': ['1', '5', '2', '3', '4'], 

2213 'SPOT1b': ['1', '4', '2', '3'], 

2214 'SPOT2b': ['1', '4', '2', '3'], 

2215 'SPOT3b': ['1', '4', '2', '3'], 

2216 'SPOT4b': ['1', '5', '2', '3', '4'], 

2217 'SPOT5b': ['1', '5', '2', '3', '4'], 

2218 } 

2219 

2220 if nBands is None or nBands == len(dict_LayerBandsAssignment[GMS_sensorcode]): 

2221 assert GMS_sensorcode in dict_LayerBandsAssignment, \ 

2222 'Unable to get Layer Bands Assignment. No valid sensorcode privided (got >%s<).' % GMS_sensorcode 

2223 LayerBandsAssignment = dict_LayerBandsAssignment[GMS_sensorcode] 

2224 if sort_by_cwl and GMS_sensorcode in ['TM4', 'TM5', 'TM7', 'LDCM']: 

2225 LayerBandsAssignment = dict_cwlSorted_LayerBandsAssignment[GMS_sensorcode] 

2226 

2227 else: # special case SPOT MSI containing no PAN or SPOT PAN containing only PAN 

2228 assert re.match(r'SPOT', GMS_id.satellite, re.I) and \ 

2229 nBands in [len(dict_LayerBandsAssignment[GMS_sensorcode]) - 1, 1], \ 

2230 "Unable to get Layer Bands Assignment. Provided number of bands doesn´t match known layer band " \ 

2231 "assignments." 

2232 LayerBandsAssignment = [dict_LayerBandsAssignment[GMS_sensorcode][-1]] if nBands == 1 \ 

2233 else dict_LayerBandsAssignment[GMS_sensorcode][:-1] 

2234 

2235 if no_thermal: 

2236 LayerBandsAssignment = [i for i in LayerBandsAssignment if not isTHERMAL(GMS_id, i)] 

2237 if no_pan: 

2238 LayerBandsAssignment = [i for i in LayerBandsAssignment if not isPAN(GMS_id, i)] 

2239 else: 

2240 LayerBandsAssignment = ['1'] 

2241 

2242 # remove those bands that are excluded by atmospheric corrections if proc_level >= L1C 

2243 if proc_level not in [None, 'L1A', 'L1B']: # TODO replace with enum procL 

2244 

2245 if CFG.target_radunit_optical == 'BOA_Ref': 

2246 # return LBA after AC 

2247 try: 

2248 bands_after_ac = get_bands_after_AC(GMS_id) 

2249 LayerBandsAssignment = [i for i in LayerBandsAssignment if i in bands_after_ac] 

2250 except ACNotSupportedError: 

2251 # atmospheric correction is not yet supported -> LBA will be the same after L1C 

2252 pass 

2253 

2254 if proc_level in ['L2B', 'L2C']: 

2255 # handle different number of bands after spectral homogenization to target sensor 

2256 

2257 if GMS_id.dataset_ID == CFG.datasetid_spectral_ref: 

2258 pass # return normal LBA from above 

2259 

2260 elif CFG.datasetid_spectral_ref is not None: 

2261 # the target sensor is NOT a custom sensor but has the spectral characteristics of a known sensor 

2262 # => use the LBA of the target sensor after AC as the new LBA for the requested sensor 

2263 

2264 # find out how the spectral characteristics of this known target sensor look like after AC 

2265 from ..model.gms_object import GMS_identifier # noqa F811 # redefinition of unused 'GMS_identifier' 

2266 tgt_sat, tgt_sen = datasetid_to_sat_sen(CFG.datasetid_spectral_ref) 

2267 tgt_GMSid = GMS_identifier(image_type='RSD', satellite=tgt_sat, sensor=tgt_sen, subsystem='', 

2268 proc_level='L2A', dataset_ID=-9999, logger=None) 

2269 try: 

2270 tgt_sen_LBA = get_bands_after_AC(tgt_GMSid) 

2271 

2272 except ACNotSupportedError: 

2273 # use the target sensor LBA before AC (because target sensor could not be atmospherically corrected) 

2274 tgt_GMSid.proc_level = 'L1B' 

2275 

2276 tgt_sen_LBA = get_LayerBandsAssignment(tgt_GMSid) 

2277 

2278 LayerBandsAssignment = tgt_sen_LBA 

2279 

2280 else: 

2281 # fallback: return a LBA matching the number of bands after spectral homogenization 

2282 LayerBandsAssignment = [str(i + 1) for i in range(len(CFG.target_CWL))] 

2283 

2284 return LayerBandsAssignment 

2285 

2286 

2287def get_bands_after_AC(GMS_id): 

2288 # type: (GMS_identifier) -> List[str] 

2289 """Returns a list of bands that are not removed by atmospheric correction. 

2290 

2291 :param GMS_id: <dict>, derived from self.get_GMS_identifier() 

2292 :return: e.g. ['1', '2', '3', '4', '5', '6', '7', '9'] for Landsat-8 

2293 """ 

2294 path_ac_options = get_path_ac_options(GMS_id) 

2295 

2296 if not path_ac_options or not os.path.exists(path_ac_options): 

2297 raise ACNotSupportedError('Atmospheric correction is not yet supported for %s %s.' 

2298 % (GMS_id.satellite, GMS_id.sensor)) 

2299 

2300 # FIXME this does not work for L7 

2301 # NOTE: don't validate because options contain pathes that do not exist on another server 

2302 ac_bandNs = get_ac_options(path_ac_options, validation=False)['AC']['bands'] 

2303 ac_out_bands = [bN.split('B0')[1] if bN.startswith('B0') else bN.split('B')[1] for bN in ac_bandNs] # sorted 

2304 

2305 return ac_out_bands 

2306 

2307 

2308def get_center_wavelengths_by_LBA(satellite, sensor, LBA, subsystem=None): 

2309 # type: (str, str, list, str) -> List[float] 

2310 """Returns a list of center wavelengths of spectral bands for the given satellite/sensor/LayerBandsAss. combination. 

2311 

2312 :param satellite: target satellite (e.g., 'Sentinel-2A') 

2313 :param sensor: target sensor (e.g., 'MSI') 

2314 :param LBA: LayerBandsAssignment 

2315 :param subsystem: target sensor subsystem (e.g., 'VNIR') 

2316 """ 

2317 srf = RSR(satellite=satellite, sensor=sensor, subsystem=subsystem, LayerBandsAssignment=LBA) 

2318 

2319 return list(srf.wvl) 

2320 

2321 

2322def get_dict_LayerOptTherm(GMS_id, LayerBandsAssignment): 

2323 dict_out = collections.OrderedDict() 

2324 [dict_out.update({lr: 'thermal' if isTHERMAL(GMS_id, lr) else 'optical'}) for lr in LayerBandsAssignment] 

2325 return dict_out 

2326 

2327 

2328def isPAN(GMS_id, LayerNr): 

2329 GMS_sensorcode = get_GMS_sensorcode(GMS_id) 

2330 dict_isPAN = {'TM7': ['8'], 'LDCM': ['8'], 

2331 'SPOT1a': ['4'], 'SPOT2a': ['4'], 'SPOT3a': ['4'], 'SPOT4a': ['5'], 'SPOT5a': ['5'], 

2332 'SPOT1b': ['4'], 'SPOT2b': ['4'], 'SPOT3b': ['4'], 'SPOT4b': ['5'], 'SPOT5b': ['5']} 

2333 return True if GMS_sensorcode in dict_isPAN and LayerNr in dict_isPAN[GMS_sensorcode] else False 

2334 

2335 

2336def isTHERMAL(GMS_id, LayerNr): 

2337 GMS_sensorcode = get_GMS_sensorcode(GMS_id) 

2338 dict_isTHERMAL = {'TM4': ['6'], 'TM5': ['6'], 'TM7': ['6L', '6H'], 'LDCM': ['10', '11'], 

2339 'AST_T': ['10', '11', '12', '13', '14']} 

2340 return True if GMS_sensorcode in dict_isTHERMAL and LayerNr in dict_isTHERMAL[GMS_sensorcode] else False 

2341 

2342 

2343def get_FieldOfView(GMS_id): 

2344 GMS_sensorcode = get_GMS_sensorcode(GMS_id) 

2345 dict_FOV = {'AVNIR-2': 5.79, 

2346 'AST_V1': 6.09, 'AST_V2': 5.19, 'AST_S': 4.9, 'AST_T': 4.9, 

2347 # http://eospso.gsfc.nasa.gov/sites/default/files/mission_handbooks/Terra.pdf 

2348 'TM4': 14.92, 'TM5': 14.92, 'TM7': 14.92, 'LDCM': 14.92, 

2349 'SPOT1a': 4.18, 'SPOT2a': 4.18, 'SPOT3a': 4.18, 'SPOT4a': 4.18, 'SPOT5a': 4.18, 

2350 'SPOT1b': 4.18, 'SPOT2b': 4.18, 'SPOT3b': 4.18, 'SPOT4b': 4.18, 'SPOT5b': 4.18, 

2351 'RE5': 6.99, 

2352 'S2A10': 20.6, 'S2A20': 20.6, 'S2A60': 20.6, 

2353 'S2B10': 20.6, 'S2B20': 20.6, 'S2B60': 20.6} 

2354 return dict_FOV[GMS_sensorcode] 

2355 

2356 

2357def get_orbit_params(GMS_id): 

2358 GMS_sensorcode = get_GMS_sensorcode(GMS_id) 

2359 

2360 # sensor altitude above ground [kilometers] 

2361 dict_altitude = {'AVNIR-2': 691.65, 

2362 'AST_V1': 705, 'AST_V2': 705, 'AST_S': 705, 'AST_T': 705, 

2363 'TM4': 705, 'TM5': 705, 'TM7': 705, 'LDCM': 705, 

2364 'SPOT1a': 822, 'SPOT2a': 822, 'SPOT3a': 822, 'SPOT4a': 822, 'SPOT5a': 822, 

2365 'SPOT1b': 822, 'SPOT2b': 822, 'SPOT3b': 822, 'SPOT4b': 822, 'SPOT5b': 822, 

2366 'RE5': 630, 

2367 'S2A10': 786, 'S2A20': 786, 'S2A60': 786, 

2368 'S2B10': 786, 'S2B20': 786, 'S2B60': 786} 

2369 

2370 # sensor inclination [degrees] 

2371 dict_inclination = {'AVNIR-2': 98.16, 

2372 'AST_V1': 98.3, 'AST_V2': 98.3, 'AST_S': 98.3, 'AST_T': 98.3, 

2373 'TM4': 98.2, 'TM5': 98.2, 'TM7': 98.2, 'LDCM': 98.2, 

2374 'SPOT1a': 98.7, 'SPOT2a': 98.7, 'SPOT3a': 98.7, 'SPOT4a': 98.7, 'SPOT5a': 98.7, 

2375 'SPOT1b': 98.7, 'SPOT2b': 98.7, 'SPOT3b': 98.7, 'SPOT4b': 98.7, 'SPOT5b': 98.7, 

2376 'RE5': 98, 

2377 'S2A10': 98.62, 'S2A20': 98.62, 'S2A60': 98.62, 

2378 'S2B10': 98.62, 'S2B20': 98.62, 'S2B60': 98.62} 

2379 

2380 # time needed for one complete earth revolution [minutes] 

2381 dict_period = {'AVNIR-2': 98.7, 

2382 'AST_V1': 98.88, 'AST_V2': 98.88, 'AST_S': 98.88, 'AST_T': 98.88, 

2383 'TM4': 98.9, 'TM5': 98.9, 'TM7': 98.9, 'LDCM': 98.9, 

2384 'SPOT1a': 101.4, 'SPOT2a': 101.4, 'SPOT3a': 101.4, 'SPOT4a': 101.4, 'SPOT5a': 101.4, 

2385 'SPOT1b': 101.4, 'SPOT2b': 101.4, 'SPOT3b': 101.4, 'SPOT4b': 101.4, 'SPOT5b': 101.4, 

2386 'RE5': 96.7, 

2387 'S2A10': 100.6, 'S2A20': 100.6, 'S2A60': 100.6, 

2388 'S2B10': 100.6, 'S2B20': 100.6, 'S2B60': 100.6} 

2389 

2390 return [dict_altitude[GMS_sensorcode], dict_inclination[GMS_sensorcode], dict_period[GMS_sensorcode]] 

2391 

2392 

2393def get_special_values(GMS_id): 

2394 GMS_sensorcode = get_GMS_sensorcode(GMS_id) # type: str 

2395 dict_fill = {'AVNIR-2': 0, 

2396 'AST_V1': 0, 'AST_V2': 0, 'AST_S': 0, 'AST_T': 0, 

2397 'TM4': 0, 'TM5': 0, 'TM7': 0, 'LDCM': 0, 

2398 'SPOT1a': 0, 'SPOT2a': 0, 'SPOT3a': 0, 'SPOT4a': 0, 'SPOT5a': 0, 

2399 'SPOT1b': 0, 'SPOT2b': 0, 'SPOT3b': 0, 'SPOT4b': 0, 'SPOT5b': 0, 

2400 'RE5': 0, 

2401 'S2A10': 0, 'S2A20': 0, 'S2A60': 0, 

2402 'S2B10': 0, 'S2B20': 0, 'S2B60': 0, 

2403 } 

2404 dict_zero = {'AVNIR-2': None, 

2405 'AST_V1': 1, 'AST_V2': 1, 'AST_S': 1, 'AST_T': 1, 

2406 'TM4': None, 'TM5': None, 'TM7': None, 'LDCM': None, 

2407 'SPOT1a': None, 'SPOT2a': None, 'SPOT3a': None, 'SPOT4a': None, 'SPOT5a': None, 

2408 'SPOT1b': None, 'SPOT2b': None, 'SPOT3b': None, 'SPOT4b': None, 'SPOT5b': None, 

2409 'RE5': None, 

2410 'S2A10': None, 'S2A20': None, 'S2A60': None, 

2411 'S2B10': None, 'S2B20': None, 'S2B60': None, 

2412 } 

2413 dict_saturated = {'AVNIR-2': None, 

2414 'AST_V1': 255, 'AST_V2': 255, 'AST_S': 255, 'AST_T': 65535, 

2415 'TM4': None, 'TM5': None, 'TM7': None, 'LDCM': 65535, 

2416 'SPOT1a': 255, 'SPOT2a': 255, 'SPOT3a': 255, 'SPOT4a': 255, 'SPOT5a': 255, 

2417 'SPOT1b': 255, 'SPOT2b': 255, 'SPOT3b': 255, 'SPOT4b': 255, 'SPOT5b': 255, 

2418 'RE5': None, 

2419 'S2A10': 65535, 'S2A20': 65535, 'S2A60': 65535, 

2420 'S2B10': 65535, 'S2B20': 65535, 'S2B60': 65535 

2421 } 

2422 return {'fill': dict_fill[GMS_sensorcode], 

2423 'zero': dict_zero[GMS_sensorcode], 

2424 'saturated': dict_saturated[GMS_sensorcode]} 

2425 

2426 

2427def metaDict_to_metaODict(metaDict, logger=None): 

2428 """Converts a GMS metadata dictionary to an ordered dictionary according to the sorting given in 

2429 Output_writer.enviHdr_keyOrder. 

2430 

2431 :param metaDict: <dict> GMS metadata dictionary 

2432 :param logger: <logging.logger> if given, warnings will be logged. Otherwise they are raised. 

2433 """ 

2434 

2435 from ..io.output_writer import enviHdr_keyOrder 

2436 expected_keys = [k for k in enviHdr_keyOrder if k in metaDict] 

2437 only_gmsFile_keys = ['ViewingAngle_arrProv', 'IncidenceAngle_arrProv', 'projection'] 

2438 unexpected_keys = [k for k in metaDict.keys() if k not in expected_keys and k not in only_gmsFile_keys] 

2439 

2440 if unexpected_keys: 

2441 msg = 'Got unexpected keys in metadata dictionary: %s. Adding them at the end of output header files.' \ 

2442 % ', '.join(unexpected_keys) 

2443 if logger: 

2444 logger.warning(msg) 

2445 else: 

2446 warnings.warn(msg) 

2447 

2448 meta_vals = [metaDict[k] for k in expected_keys] + [metaDict[k] for k in unexpected_keys] 

2449 return collections.OrderedDict(zip(expected_keys + unexpected_keys, meta_vals)) 

2450 

2451 

2452def LandsatID2dataset(ID_list): 

2453 dataset_list = [] 

2454 for ID in ID_list: 

2455 dataset = dict(image_type='RSD', satellite=None, sensor=None, subsystem=None, acquisition_date=None, 

2456 entity_ID=ID) 

2457 dataset['satellite'] = 'Landsat-5' if ID[:3] == 'LT5' else 'Landsat-7' if ID[:3] == 'LE7' \ 

2458 else 'Landsat-8' if ID[:3] == 'LC8' else dataset['satellite'] 

2459 dataset['sensor'] = 'TM' if ID[:3] == 'LT5' else 'ETM+' if ID[:3] == 'LE7' else \ 

2460 'OLI_TIRS' if ID[:3] == 'LC8' else dataset['satellite'] 

2461 dataset['subsystem'] = None 

2462 dataset['acquisition_date'] = \ 

2463 (datetime.datetime(int(ID[9:13]), 1, 1) + datetime.timedelta(int(ID[13:16]) - 1)).strftime('%Y-%m-%d') 

2464 dataset_list.append(dataset) 

2465 return dataset_list 

2466 

2467 

2468def get_sensormode(dataset): 

2469 if re.search(r'SPOT', dataset['satellite']): 

2470 path_archive = path_generator(dataset).get_local_archive_path_baseN() 

2471 dim_ = open_specific_file_within_archive(path_archive, '*/scene01/metadata.dim')[0] 

2472 SPOT_mode = re.search(r"<SENSOR_CODE>([a-zA-Z0-9]*)</SENSOR_CODE>", dim_, re.I).group(1) 

2473 assert SPOT_mode in ['J', 'X', 'XS', 'A', 'P', 'M'], 'Unknown SPOT sensor mode: %s' % SPOT_mode 

2474 return 'M' if SPOT_mode in ['J', 'X', 'XS'] else 'P' 

2475 else: 

2476 return 'M'