Coverage for enpt/model/images/image_baseclasses.py: 89%

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

2 

3# EnPT, EnMAP Processing Tool - A Python package for pre-processing of EnMAP Level-1B data 

4# 

5# Copyright (C) 2018-2024 Karl Segl (GFZ Potsdam, segl@gfz-potsdam.de), Daniel Scheffler 

6# (GFZ Potsdam, danschef@gfz-potsdam.de), Niklas Bohn (GFZ Potsdam, nbohn@gfz-potsdam.de), 

7# Stéphane Guillaso (GFZ Potsdam, stephane.guillaso@gfz-potsdam.de) 

8# 

9# This software was developed within the context of the EnMAP project supported 

10# by the DLR Space Administration with funds of the German Federal Ministry of 

11# Economic Affairs and Energy (on the basis of a decision by the German Bundestag: 

12# 50 EE 1529) and contributions from DLR, GFZ and OHB System AG. 

13# 

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

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

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

17# version. Please note the following exception: `EnPT` depends on tqdm, which 

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

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

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

21# 

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

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

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

25# details. 

26# 

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

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

29 

30"""EnPT EnMAP object base classes.""" 

31 

32from types import SimpleNamespace 

33from typing import List, Optional # noqa: F401 

34import numpy as np 

35 

36# noinspection PyPackageRequirements 

37from skimage import exposure # contained in package requirements as scikit-image 

38 

39from geoarray import GeoArray 

40 

41from ...model.metadata import EnMAP_Metadata_L2A_MapGeo # noqa: F401 # only used for type hint 

42 

43__author__ = ['Daniel Scheffler', 'Stéphane Guillaso', 'André Hollstein'] 

44 

45 

46class _EnMAP_Image(object): 

47 """EnPT base class for all kinds of EnMAP images. 

48 

49 NOTE: 

50 - Basic functionality that all EnMAP image objects have in common is to be implemented here. 

51 - All EnMAP image classes should (directly or indirectly) inherit from _EnMAP_image. 

52 

53 Attributes: 

54 - to be listed here. Check help(_EnMAP_image) in the meanwhile! 

55 

56 """ 

57 

58 def __init__(self): 

59 """Load and hold data needed for processing EnMAP Level-1B data to Level-2A. 

60 

61 Intended usage: 

62 - only to be used as a base class to be inherited from! 

63 - Example: 

64 class EnMAP_Image_Subclass(_EnMAP_Image): 

65 def __init__(self): 

66 super(EnMAP_Image_Subclass, self).__init__() 

67 

68 """ 

69 # add EnMAP specific attributes 

70 self.paths = SimpleNamespace() 

71 

72 # private attributes 

73 self._data = None 

74 self._mask_landwater = None 

75 self._mask_clouds = None 

76 self._mask_cloudshadow = None 

77 self._mask_haze = None 

78 self._mask_snow = None 

79 self._mask_cirrus = None 

80 self._dem = None 

81 self._deadpixelmap = None 

82 self._polymer_logchl = None 

83 self._polymer_logfb = None 

84 self._polymer_rgli = None 

85 self._polymer_rnir = None 

86 self._polymer_bitmask = None 

87 self._subset = None # FIXME how is _subset to be set? 

88 

89 # defaults 

90 self.entity_ID = '' 

91 self.basename = '' 

92 

93 @property 

94 def data(self) -> GeoArray: 

95 """Return the actual EnMAP image data. 

96 

97 Bundled with all the corresponding metadata. 

98 

99 Attributes and functions (most important; for a full list check help(self.data)!): 

100 

101 - ALL attributes of numpy.ndarray! 

102 - is_inmem(bool): 

103 True if the image data are completely loaded into memory; False if GeoArray only holds a link to a file 

104 on disk. 

105 - arr: np.ndarray holding the pixel values (if is_mem is True) 

106 - rows(int) 

107 - cols(int) 

108 - bands(int) 

109 - shape(tuple) 

110 - gt(list): GDAL geotransform: contains the geocoding 

111 - prj(str): WKT projection string 

112 - show(): plot the image 

113 - show_map(): plot a map of the image (based on cartopy library) 

114 - reproject_to_new_grid() 

115 

116 Usage (there will soon be detailed instructions on usage at https://git.gfz-potsdam.de/danschef/geoarray): 

117 

118 - Use self.data like a normal numpy.ndarray! 

119 - NOTE: Operators like *, /, + , - will soon be implemented. In the meanwhile use: 

120 result = self.data[:] *10 

121 

122 - How to set self.data? 

123 - Link an image file to self.data -> all raster data is read into memory ON DEMAND: 

124 self.data = '/path/to/image.tif' # sets self.data to GeoArray('/path/to/image.tif') 

125 

126 - Link a numpy.ndarray instance with self.data (remaining attributes like geocoding, projection, etc. 

127 are copied from the previous self.data attribute). 

128 self.data = numpy.array([[1,2,3],[4,5,6]]) 

129 

130 - Set self.data to an existing instance of GeoArray 

131 (allows to set specific attributes of GeoArray by yourself) 

132 self.data = GeoArray('/path/to/image.tif', geotransform=[...], projection='WKTString') 

133 

134 - Delete self.data: 

135 del self.data or 'self.data = None' 

136 

137 :return: instance of geoarray.GeoArray 

138 """ 

139 if self._subset is None: 

140 return self._data 

141 

142 return GeoArray(self._data[self._subset[0]:self._subset[1], self._subset[2]:self._subset[3], :], 

143 geotransform=self._data.gt, projection=self._data.prj) 

144 

145 @data.setter 

146 def data(self, *geoArr_initArgs): 

147 if geoArr_initArgs[0] is not None: 

148 # TODO this must be able to handle subset inputs in tiled processing 

149 if isinstance(geoArr_initArgs[0], np.ndarray): 

150 self._data = GeoArray(geoArr_initArgs[0], geotransform=self._data.gt, projection=self._data.prj) 

151 else: 

152 self._data = GeoArray(*geoArr_initArgs) 

153 else: 

154 del self.data 

155 

156 @data.deleter 

157 def data(self): 

158 self._data = None 

159 

160 @property 

161 def mask_landwater(self) -> GeoArray: 

162 """Return the land/water mask. 

163 

164 pixel values: 

165 - 0: background within scene dimensions, e.g. due to missing values/errors 

166 - 1: no water 

167 - 2: water 

168 - 3: background outside the scene dimensions (artifact from resampling between map and sensor geometry) 

169 

170 :return: geoarray.GeoArray 

171 """ 

172 return self._mask_landwater 

173 

174 @mask_landwater.setter 

175 def mask_landwater(self, *geoArr_initArgs): 

176 self._mask_landwater = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'mask_landwater', nodataVal=3) 

177 

178 @mask_landwater.deleter 

179 def mask_landwater(self): 

180 self._mask_landwater = None 

181 

182 @property 

183 def mask_clouds(self) -> GeoArray: 

184 """Return the cloud mask. 

185 

186 :return: geoarray.GeoArray 

187 """ 

188 return self._mask_clouds 

189 

190 @mask_clouds.setter 

191 def mask_clouds(self, *geoArr_initArgs): 

192 self._mask_clouds = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'mask_clouds', nodataVal=0) 

193 

194 @mask_clouds.deleter 

195 def mask_clouds(self): 

196 self._mask_clouds = None 

197 

198 @property 

199 def mask_cloudshadow(self) -> GeoArray: 

200 """Return the cloud shadow mask (0=no cloud shadow, 1=cloud shadow). 

201 

202 :return: geoarray.GeoArray 

203 """ 

204 return self._mask_cloudshadow 

205 

206 @mask_cloudshadow.setter 

207 def mask_cloudshadow(self, *geoArr_initArgs): 

208 self._mask_cloudshadow = \ 

209 self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'mask_cloudshadow', nodataVal=0) 

210 

211 @mask_cloudshadow.deleter 

212 def mask_cloudshadow(self): 

213 self._mask_cloudshadow = None 

214 

215 @property 

216 def mask_haze(self) -> GeoArray: 

217 """Return the haze mask (0=no haze, 1=haze).. 

218 

219 :return: geoarray.GeoArray 

220 """ 

221 return self._mask_haze 

222 

223 @mask_haze.setter 

224 def mask_haze(self, *geoArr_initArgs): 

225 self._mask_haze = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'mask_haze', nodataVal=0) 

226 

227 @mask_haze.deleter 

228 def mask_haze(self): 

229 self._mask_haze = None 

230 

231 @property 

232 def mask_snow(self) -> GeoArray: 

233 """Return the snow mask (0=no snow, 1=snow). 

234 

235 :return: geoarray.GeoArray 

236 """ 

237 return self._mask_snow 

238 

239 @mask_snow.setter 

240 def mask_snow(self, *geoArr_initArgs): 

241 self._mask_snow = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'mask_snow', nodataVal=0) 

242 

243 @mask_snow.deleter 

244 def mask_snow(self): 

245 self._mask_snow = None 

246 

247 @property 

248 def mask_cirrus(self) -> GeoArray: 

249 """Return the cirrus mask (0=none, 1=thin, 2=medium, 3=thick). 

250 

251 :return: geoarray.GeoArray 

252 """ 

253 return self._mask_cirrus 

254 

255 @mask_cirrus.setter 

256 def mask_cirrus(self, *geoArr_initArgs): 

257 self._mask_cirrus = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'mask_cirrus', nodataVal=0) 

258 

259 @mask_cirrus.deleter 

260 def mask_cirrus(self): 

261 self._mask_cirrus = None 

262 

263 @property 

264 def dem(self) -> GeoArray: 

265 """Return a DEM in the exact dimension and pixel grid of self.data. 

266 

267 :return: geoarray.GeoArray 

268 """ 

269 if self._dem is None: 

270 raise NotImplementedError('DEM is missing. An automatic DEM getter is currently not implemented.') 

271 return self._dem 

272 

273 @dem.setter 

274 def dem(self, *geoArr_initArgs): 

275 self._dem = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'dem', nodataVal=0) # FIXME 0? 

276 

277 @dem.deleter 

278 def dem(self): 

279 self._dem = None 

280 

281 @property 

282 def deadpixelmap(self) -> GeoArray: 

283 """Return the dead pixel map. 

284 

285 :return: geoarray.GeoArray 

286 """ 

287 if self._deadpixelmap is not None: 

288 self._deadpixelmap.arr = self._deadpixelmap[:].astype(bool) # ensure boolean map 

289 

290 return self._deadpixelmap 

291 

292 @deadpixelmap.setter 

293 def deadpixelmap(self, *geoArr_initArgs): 

294 if geoArr_initArgs[0] is not None: 

295 dpm = GeoArray(*geoArr_initArgs) 

296 

297 if dpm.ndim == 3 and dpm.shape != self.data.shape: 

298 raise ValueError("The 'deadpixelmap' GeoArray can only be instanced with a 3D array with the same size " 

299 "like %s.data, i.e.: %s Received %s." 

300 % (self.__class__.__name__, str(self.data.shape), str(dpm.shape))) 

301 elif dpm.ndim == 2 and dpm.shape != (self.data.bands, self.data.cols): 

302 raise ValueError("The 'deadpixelmap' GeoArray can only be instanced with an array with the size " 

303 "'bands x columns' of the GeoArray %s.data. Received %s. Expected %s" 

304 % (self.__class__.__name__, str(dpm.shape), str((self.data.bands, self.data.cols)))) 

305 

306 self._deadpixelmap = dpm 

307 else: 

308 del self.deadpixelmap 

309 

310 @deadpixelmap.deleter 

311 def deadpixelmap(self): 

312 self._deadpixelmap = None 

313 

314 @property 

315 def polymer_logchl(self) -> GeoArray: 

316 """Chlorophyll concentration, in mg/m3, in 10-based logarithm (generated by Polymer AC for water areas). 

317 

318 Data range: -3 to 2.3 

319 

320 :return: geoarray.GeoArray 

321 """ 

322 return self._polymer_logchl 

323 

324 @polymer_logchl.setter 

325 def polymer_logchl(self, *geoArr_initArgs): 

326 self._polymer_logchl = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'polymer_logchl', 

327 nodataVal=-9999) 

328 

329 @polymer_logchl.deleter 

330 def polymer_logchl(self): 

331 self._polymer_logchl = None 

332 

333 @property 

334 def polymer_logfb(self) -> GeoArray: 

335 """Particle scattering factor fb in Park & Ruddick (2005) (in 10-based logarithm; generated by Polymer AC). 

336 

337 Data range: positive and negative values, maximum range not defined 

338 

339 :return: geoarray.GeoArray 

340 """ 

341 return self._polymer_logfb 

342 

343 @polymer_logfb.setter 

344 def polymer_logfb(self, *geoArr_initArgs): 

345 self._polymer_logfb = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'polymer_logfb', 

346 nodataVal=-9999) 

347 

348 @polymer_logfb.deleter 

349 def polymer_logfb(self): 

350 self._polymer_logfb = None 

351 

352 @property 

353 def polymer_rgli(self) -> GeoArray: 

354 """Reflectance of the sun glint predicted from ECMWF wind speed (generated by Polymer AC for water areas). 

355 

356 Typical data range: 0-0.2 

357 

358 :return: geoarray.GeoArray 

359 """ 

360 return self._polymer_rgli 

361 

362 @polymer_rgli.setter 

363 def polymer_rgli(self, *geoArr_initArgs): 

364 self._polymer_rgli = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'polymer_rgli', 

365 nodataVal=-9999) 

366 

367 @polymer_rgli.deleter 

368 def polymer_rgli(self): 

369 self._polymer_rgli = None 

370 

371 @property 

372 def polymer_rnir(self) -> GeoArray: 

373 """TOA reflectance at 865 nm corrected for Rayleigh scattering (generated by Polymer AC for water areas). 

374 

375 Typical data range: 0-0.2 

376 

377 :return: geoarray.GeoArray 

378 """ 

379 return self._polymer_rnir 

380 

381 @polymer_rnir.setter 

382 def polymer_rnir(self, *geoArr_initArgs): 

383 self._polymer_rnir = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'polymer_rnir', 

384 nodataVal=-9999) 

385 

386 @polymer_rnir.deleter 

387 def polymer_rnir(self): 

388 self._polymer_rnir = None 

389 

390 @property 

391 def polymer_bitmask(self) -> GeoArray: 

392 """Quality flags as generated by Polymer AC (bit-encoded mask). 

393 

394 Value 0 represents water, all fine, no flags. Other pixel values are explained below: 

395 

396 +--------------------+-------------+--------------------------------------------+ 

397 | Flag name | Flag value | Description | 

398 +====================+=============+============================================+ 

399 | LAND | 1 | Land mask | 

400 +--------------------+-------------+--------------------------------------------+ 

401 | CLOUD_BASE | 2 | Polymer's basic cloud mask | 

402 +--------------------+-------------+--------------------------------------------+ 

403 | L1_INVALID | 4 | Invalid level1 pixel | 

404 +--------------------+-------------+--------------------------------------------+ 

405 | NEGATIVE_BB | 8 | (deprecated flag) | 

406 +--------------------+-------------+--------------------------------------------+ 

407 | OUT_OF_BOUNDS | 16 | Retrieved marine parameters are outside | 

408 | | | valid bounds | 

409 +--------------------+-------------+--------------------------------------------+ 

410 | EXCEPTION | 32 | A processing error was encountered | 

411 +--------------------+-------------+--------------------------------------------+ 

412 | THICK_AEROSOL | 64 | Thick aerosol flag | 

413 +--------------------+-------------+--------------------------------------------+ 

414 | HIGH_AIR_MASS | 128 | Air mass exceeds 5 | 

415 +--------------------+-------------+--------------------------------------------+ 

416 | EXTERNAL_MASK | 512 | Pixel was masked using external mask | 

417 +--------------------+-------------+--------------------------------------------+ 

418 | CASE2 | 1024 | Pixel was processed in "case2" mode | 

419 +--------------------+-------------+--------------------------------------------+ 

420 | INCONSISTENCY | 2048 | Inconsistent result was detected | 

421 | | | (atmospheric reflectance out of bounds) | 

422 +--------------------+-------------+--------------------------------------------+ 

423 | ANOMALY_RWMOD_BLUE | 4096 | Excessive difference was found at 412nm | 

424 | | | between Rw and Rwmod | 

425 +--------------------+-------------+--------------------------------------------+ 

426 

427 :return: geoarray.GeoArray 

428 """ 

429 return self._polymer_bitmask 

430 

431 @polymer_bitmask.setter 

432 def polymer_bitmask(self, *geoArr_initArgs): 

433 self._polymer_bitmask = self._get_geoarray_with_datalike_geometry(geoArr_initArgs, 'polymer_bitmask', 

434 nodataVal=-9999) 

435 

436 @polymer_bitmask.deleter 

437 def polymer_bitmask(self): 

438 self._polymer_bitmask = None 

439 

440 def _get_geoarray_with_datalike_geometry(self, 

441 geoArr_initArgs: tuple, 

442 attrName: str, 

443 nodataVal: int = None, 

444 specialclass=None) -> Optional[GeoArray]: 

445 if geoArr_initArgs[0] is None: 

446 return None 

447 else: 

448 GeoArrayOrSubclass = GeoArray if not specialclass else specialclass 

449 gA = GeoArrayOrSubclass(*geoArr_initArgs) 

450 

451 if gA.shape[:2] != self.data.shape[:2]: 

452 raise ValueError("The '%s' GeoArray can only be instanced with an array with " 

453 "the same X/Y dimensions like %s.data %s. Got %s." % 

454 (attrName, self.__class__.__name__, str(self.data.shape[:2]), str(gA.shape[:2]))) 

455 

456 # noinspection PyProtectedMember 

457 if gA._nodata is None and nodataVal is not None: 

458 gA.nodata = nodataVal 

459 gA.gt = self.data.gt 

460 gA.prj = self.data.prj 

461 

462 return gA 

463 

464 def generate_quicklook(self, bands2use: List[int]) -> GeoArray: 

465 """ 

466 Generate image quicklook and save into a file. 

467 

468 :param bands2use: band indices of self.data to be used as (red, green, blue) for quicklook image, 

469 e.g., [3, 2, 1] 

470 :return: GeoArray 

471 """ 

472 red, green, blue = [self.data[:, :, bandidx] for bandidx in bands2use] 

473 

474 def rescale(bandarray): 

475 pixvals = np.ma.masked_equal(bandarray, self.data.nodata).compressed() \ 

476 if self.data.nodata is not None else bandarray 

477 p2, p98 = np.nanpercentile(pixvals, 2), np.nanpercentile(pixvals, 98) 

478 

479 return exposure.rescale_intensity(bandarray, in_range=(p2, p98), out_range=(0, 255)) 

480 

481 pix = np.dstack((rescale(red), rescale(green), rescale(blue))).astype(np.uint8) 

482 

483 return GeoArray(pix)