...
 
Commits (13)
......@@ -2,19 +2,28 @@
History
=======
0.12.? (2020-04-??)
-------------------
coming soon
-----------
* Added config parameters to run EnPT in 3 AC modes: 'land', 'water', 'combined'.
* Added some boilerplate code in atmospheric_correction.py which is to be replaced by separate AC calls for water and
land surfaces later.
0.12.x (2020-05-??)
0.13.0 (2020-05-18)
-------------------
* Renamed DEM for Arcachon test dataset.
* Fixed typo.
* Added config parameter 'output_format'.
* Implementated ENVI output format.
* Fixed log message.
* The output interleave is now also configurable via the parameter 'output_interleave'.
* Implemented 3 new config parameters: 'target_projection_type', 'target_epsg' and 'target_coord_grid'. This allows
the user to choose between UTM and geographic L2A projection, to specify a custom L2A projection by providing an
EPSG code or to specify a certain L2A coordinate grid.
* Added 'grid_res' and 'epsg' attributes to EnMAP_Metadata_L2A_MapGeo object.
* The L2A projection metadata is now correctly written to the XML file.
0.12.8 (2020-05-13)
......
......@@ -80,6 +80,10 @@ def get_enpt_argparser():
'ignored if DEM is given')
add('-od', '--output_dir', type=str, default=None,
help='output directory where processed data and log files are saved')
add('-of', '--output_format', type=str, default='GTiff',
help="file format of all raster output files ('GTiff': GeoTIFF, 'ENVI': ENVI BSQ; default: 'ENVI')")
add('-ointlv', '--output_interleave', type=str, default='pixel',
help="raster data interleaving type ('band', 'line', 'pixel'; default: 'pixel')")
add('-wd', '--working_dir', type=str, default=None,
help='directory to be used for temporary files')
add('-nla', '--n_lines_to_append', type=int, default=None,
......@@ -129,6 +133,12 @@ def get_enpt_argparser():
add('--vswir_overlap_algorithm', type=str, default='swir_only',
help="Algorithm specifying how to deal with the spectral bands in the VNIR/SWIR spectral overlap region "
"('order_by_wvl', 'average', 'vnir_only', 'swir_only')")
add('-tgtprj', '--target_projection_type', type=str, default='UTM',
help="Projection type of the raster output files ('UTM', 'Geographic') (default: 'UTM')")
add('-tgtepsg', '--target_epsg', type=int, default=None,
help="Custom EPSG code of the target projection (overrides target_projection_type)")
add('-tgtgrid', '--target_coord_grid', nargs=4, type=float, default=None,
help="Custom target coordinate grid where is output is resampled to ([x0, x1, y0, y1], e.g., [0, 30, 0, 30])")
# link parser to run function
parser.set_defaults(func=run_job)
......
......@@ -239,7 +239,16 @@ class EnMAPL2Product_MapGeo(_EnMAP_Image):
makedirs(product_dir, exist_ok=True)
# save raster data
kwargs_save = dict(fmt='GTiff', creationOptions=["COMPRESS=LZW"])
kwargs_save = \
dict(fmt='GTiff',
creationOptions=["COMPRESS=LZW",
"NUM_THREADS=%d" % self.cfg.CPUs,
"INTERLEAVE=%s" % ('BAND' if self.cfg.output_interleave == 'band' else 'PIXEL')]
) if self.cfg.output_format == 'GTiff' else \
dict(fmt='ENVI',
creationOptions=["INTERLEAVE=%s" % ("BSQ" if self.cfg.output_interleave == 'band' else
"BIL" if self.cfg.output_interleave == 'line' else
"BIP")])
outpaths = dict(metaxml=path.join(product_dir, self.meta.filename_metaxml))
for attrName in ['data', 'mask_landwater', 'mask_clouds', 'mask_cloudshadow', 'mask_haze', 'mask_snow',
......
......@@ -120,7 +120,7 @@ class EnMAP_Detector_SensorGeo(_EnMAP_Image):
method_spectral, method_spatial = self.cfg.deadpix_P_interp_spectral, self.cfg.deadpix_P_interp_spatial
self.logger.info("Correcting dead pixels of %s detector...\n"
"Used algorithm: %s interpolation in the %s domain"
% (self.detector_name, algo, method_spectral if algo == 'spectral' else method_spatial))
% (self.detector_name, method_spectral, algo if algo == 'spectral' else method_spatial))
self.data = \
Dead_Pixel_Corrector(algorithm=algo,
......
......@@ -41,7 +41,7 @@ from py_tools_ds.geo.vector.topology import Polygon, get_footprint_polygon # no
from geoarray import GeoArray
from .metadata_sensorgeo import EnMAP_Metadata_L1B_SensorGeo
from ...options.config import EnPTConfig, enmap_xres
from ...options.config import EnPTConfig
from ..srf import SRF
from ...version import __version__
......@@ -54,6 +54,7 @@ class EnMAP_Metadata_L2A_MapGeo(object):
meta_l1b: EnMAP_Metadata_L1B_SensorGeo,
wvls_l2a: Union[List, np.ndarray],
dims_mapgeo: Tuple[int, int, int],
grid_res_l2a: Tuple[float, float],
logger=None):
"""EnMAP Metadata class for the metadata of the complete EnMAP L2A product in map geometry incl. VNIR and SWIR.
......@@ -61,10 +62,12 @@ class EnMAP_Metadata_L2A_MapGeo(object):
:param meta_l1b: metadata object of the L1B dataset in sensor geometry
:param wvls_l2a: list of center wavelengths included in the L2A product
:param dims_mapgeo: dimensions of the EnMAP raster data in map geometry, e.g., (1024, 1000, 218)
:param grid_res_l2a Coordinate grid resolution of the L2A product (x, y)
:param logger: instance of logging.logger or subclassed
"""
self.cfg = config
self._meta_l1b = meta_l1b
self.grid_res = grid_res_l2a
self.logger = logger or logging.getLogger()
# defaults
......@@ -88,22 +91,32 @@ class EnMAP_Metadata_L2A_MapGeo(object):
self.earthSunDist: float = meta_l1b.earthSunDist # earth-sun distance
# generate file names for L2A output
file_ext_l1b = os.path.splitext(meta_l1b.vnir.filename_data)[1]
file_ext_l2a = \
'.TIF' if self.cfg.output_format == 'GTiff' else \
'.bsq' if self.cfg.output_interleave == 'band' else \
'.bil' if self.cfg.output_interleave == 'line' else \
'.bip'
def convert_fn(fn):
return fn.replace('L1B-', 'L2A-').replace(file_ext_l1b, file_ext_l2a)
if not self.cfg.is_dummy_dataformat:
self.scene_basename = meta_l1b.vnir.filename_data.split('-SPECTRAL_IMAGE')[0].replace('L1B-', 'L2A-')
self.scene_basename = convert_fn(meta_l1b.vnir.filename_data.split('-SPECTRAL_IMAGE')[0])
else:
self.scene_basename = os.path.splitext(meta_l1b.vnir.filename_data)[0]
self.filename_data = meta_l1b.vnir.filename_data.replace('L1B-', 'L2A-').replace('_VNIR', '')
self.filename_deadpixelmap_vnir = meta_l1b.vnir.filename_deadpixelmap.replace('L1B-', 'L2A-')
self.filename_deadpixelmap_swir = meta_l1b.swir.filename_deadpixelmap.replace('L1B-', 'L2A-')
self.filename_quicklook_vnir = meta_l1b.vnir.filename_quicklook.replace('L1B-', 'L2A-')
self.filename_quicklook_swir = meta_l1b.swir.filename_quicklook.replace('L1B-', 'L2A-')
self.filename_mask_landwater = meta_l1b.vnir.filename_mask_landwater.replace('L1B-', 'L2A-')
self.filename_mask_clouds = meta_l1b.vnir.filename_mask_clouds.replace('L1B-', 'L2A-')
self.filename_mask_cloudshadow = meta_l1b.vnir.filename_mask_cloudshadow.replace('L1B-', 'L2A-')
self.filename_mask_haze = meta_l1b.vnir.filename_mask_haze.replace('L1B-', 'L2A-')
self.filename_mask_snow = meta_l1b.vnir.filename_mask_snow.replace('L1B-', 'L2A-')
self.filename_mask_cirrus = meta_l1b.vnir.filename_mask_cirrus.replace('L1B-', 'L2A-')
self.filename_metaxml = meta_l1b.filename_metaxml.replace('L1B-', 'L2A-')
self.filename_data = convert_fn(meta_l1b.vnir.filename_data).replace('_VNIR', '')
self.filename_deadpixelmap_vnir = convert_fn(meta_l1b.vnir.filename_deadpixelmap)
self.filename_deadpixelmap_swir = convert_fn(meta_l1b.swir.filename_deadpixelmap)
self.filename_quicklook_vnir = convert_fn(meta_l1b.vnir.filename_quicklook)
self.filename_quicklook_swir = convert_fn(meta_l1b.swir.filename_quicklook)
self.filename_mask_landwater = convert_fn(meta_l1b.vnir.filename_mask_landwater)
self.filename_mask_clouds = convert_fn(meta_l1b.vnir.filename_mask_clouds)
self.filename_mask_cloudshadow = convert_fn(meta_l1b.vnir.filename_mask_cloudshadow)
self.filename_mask_haze = convert_fn(meta_l1b.vnir.filename_mask_haze)
self.filename_mask_snow = convert_fn(meta_l1b.vnir.filename_mask_snow)
self.filename_mask_cirrus = convert_fn(meta_l1b.vnir.filename_mask_cirrus)
self.filename_metaxml = convert_fn(meta_l1b.filename_metaxml)
# fuse band-wise metadata (sort all band-wise metadata by wavelengths but band number keeps as it is)
# get band index order
......@@ -144,6 +157,7 @@ class EnMAP_Metadata_L2A_MapGeo(object):
self.lat_UL_UR_LL_LR = [lat for lon, lat in common_UL_UR_LL_LR]
self.ll_mapPoly = get_footprint_polygon(tuple(zip(self.lon_UL_UR_LL_LR,
self.lat_UL_UR_LL_LR)), fix_invalid=True)
self.epsg = self._meta_l1b.vnir.epsg_ortho
if meta_l1b.vnir.unit != meta_l1b.swir.unit or meta_l1b.vnir.unitcode != meta_l1b.swir.unitcode:
raise RuntimeError('L2A data should have the same radiometric unit for VNIR and SWIR. '
......@@ -198,9 +212,9 @@ class EnMAP_Metadata_L2A_MapGeo(object):
'.TIF',
'.TIFF',
'.GTIFF',
'.BSQ',
'.BIL',
'.BIP',
'.BSQ', '.bsq',
'.BIL', '.bil',
'.BIP', '.bip',
'.JPEG2000'] \
else 'xml' if ext == '.XML' \
else 'NA'
......@@ -407,9 +421,20 @@ class EnMAP_Metadata_L2A_MapGeo(object):
# ortho #
#########
# TODO update product/ortho/projection
# {UTM_ZoneX_North, UTM_ZoneX_South (where X in {1..60}), Geographic, LAEA-ETRS89, NA}
uk('product/ortho/resolution', enmap_xres)
if self.epsg == 4326:
proj_str = 'Geographic'
elif self.epsg == 3035:
proj_str = 'LAEA-ETRS89'
elif len(str(self.epsg)) == 5 and str(self.epsg)[:3] in ['326', '327']:
zone = int(str(self.epsg)[-2:])
if not 0 <= zone <= 60:
raise RuntimeError('Invalid L2A UTM zone: %d.' % zone)
proj_str = 'UTM_Zone%d_North' % zone if str(self.epsg).startswith('326') else 'UTM_Zone%d_South' % zone
else:
proj_str = 'NA'
uk('product/ortho/projection', proj_str) # {UTM_ZoneX_North, UTM_ZoneX_South (where X in {1..60}), Geographic, LAEA-ETRS89, NA} # noqa
uk('product/ortho/resolution', self.grid_res[0])
uk('product/ortho/resampling', self.cfg.ortho_resampAlg)
# band statistics
......
......@@ -271,7 +271,10 @@ class EnMAP_Metadata_L1B_Detector_SensorGeo(object):
self.ll_mapPoly = get_footprint_polygon(tuple(zip(self.lon_UL_UR_LL_LR,
self.lat_UL_UR_LL_LR)), fix_invalid=True)
from ...processors.spatial_transform import get_UTMEPSG_from_LonLat_cornersXY
self.epsg_ortho = get_UTMEPSG_from_LonLat_cornersXY(lons=self.lon_UL_UR_LL_LR, lats=self.lat_UL_UR_LL_LR)
# NOTE: self.cfg.target_epsg is set if user-provided or in case of Lon/Lat coordinates
self.epsg_ortho = \
self.cfg.target_epsg or \
get_UTMEPSG_from_LonLat_cornersXY(lons=self.lon_UL_UR_LL_LR, lats=self.lat_UL_UR_LL_LR)
def calc_smile(self):
"""Compute smile for each EnMAP column.
......
......@@ -125,6 +125,11 @@ config_for_testing_dlr = dict(
n_lines_to_append=50,
disable_progress_bars=False,
is_dummy_dataformat=False,
# output_format='ENVI',
# output_interleave='band',
# target_projection_type='Geographic',
# target_epsg=32632,
# target_coord_grid=[-1.37950, -1.37923, 44.60710, 44.60737],
enable_ac=True,
mode_ac='land',
enable_ice_retrieval=False,
......@@ -134,10 +139,9 @@ config_for_testing_dlr = dict(
)
enmap_coordinate_grid = dict(x=np.array([0, 30]),
y=np.array([0, 30]))
enmap_xres, enmap_yres = np.ptp(enmap_coordinate_grid['x']), np.ptp(enmap_coordinate_grid['y'])
assert enmap_xres == enmap_yres, 'Unequal X/Y resolution of the output grid!'
enmap_coordinate_grid_utm = dict(x=np.array([0, 30]),
y=np.array([0, 30]))
enmap_xres, enmap_yres = np.ptp(enmap_coordinate_grid_utm['x']), np.ptp(enmap_coordinate_grid_utm['y'])
class EnPTConfig(object):
......@@ -168,6 +172,15 @@ class EnPTConfig(object):
:key output_dir:
output directory where processed data and log files are saved
:key output_format:
file format of all raster output files ('GTiff': GeoTIFF, 'ENVI': ENVI BSQ; default: 'ENVI')
:key output_interleave:
raster data interleaving type (default: 'pixel')
- 'band': band-sequential (BSQ),
- 'line': data interleaved-by-line (BIL; only usable for ENVI output format),
- 'pixel' data interleaved-by-pixel (BIP)
:key working_dir:
directory to be used for temporary files
......@@ -236,8 +249,18 @@ class EnPTConfig(object):
:key deadpix_P_interp_spatial:
Spatial interpolation algorithm to be used during dead pixel correction
('linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic')
:key ortho_resampAlg:
Ortho-rectification resampling algorithm ('nearest', 'bilinear', 'gauss')
:key target_projection_type:
Projection type of the raster output files ('UTM', 'Geographic') (default: 'UTM')
:key target_epsg:
Custom EPSG code of the target projection (overrides target_projection_type)
:key target_coord_grid:
Custom target coordinate grid where the output is resampled to ([x0, x1, y0, y1], e.g., [0, 30, 0, 30])
"""
# fixed attributes
......@@ -263,7 +286,7 @@ class EnPTConfig(object):
self.is_dummy_dataformat = gp('is_dummy_dataformat')
if 'is_dlr_dataformat' in user_opts:
warnings.warn("The 'is_dlr_dataformat' flag is deprectated and will not exist in future. "
warnings.warn("The 'is_dlr_dataformat' flag is deprecated and will not exist in future. "
"Please set 'is_dummy_dataformat' to False instead.", DeprecationWarning)
self.is_dummy_dataformat = user_opts['is_dlr_dataformat'] is False
......@@ -284,6 +307,8 @@ class EnPTConfig(object):
##################
self.output_dir = self.absPath(gp('output_dir', fallback=os.path.abspath(os.path.curdir)))
self.output_format = gp('output_format')
self.output_interleave = gp('output_interleave')
###########################
# processor configuration #
......@@ -322,6 +347,10 @@ class EnPTConfig(object):
# orthorectification / VSWIR fusion
self.ortho_resampAlg = gp('ortho_resampAlg')
self.vswir_overlap_algorithm = gp('vswir_overlap_algorithm')
self.target_projection_type = gp('target_projection_type')
self.target_epsg = gp('target_epsg')
grid = gp('target_coord_grid')
self.target_coord_grid = dict(x=np.array(grid[:2]), y=np.array(grid[2:])) if grid else None
#########################
# validate final config #
......@@ -329,6 +358,32 @@ class EnPTConfig(object):
EnPTValidator(allow_unknown=True, schema=enpt_schema_config_output).validate(self.to_dict())
# check invalid interleave
if self.output_interleave == 'line' and self.output_format == 'GTiff':
warnings.warn("The interleaving type 'line' is not supported by the GTiff output format. Using 'pixel'.",
UserWarning)
self.output_interleave = 'pixel'
# override target_projection_type if target_epsg is given
if self.target_epsg:
self.target_projection_type = \
'Geographic' if self.target_epsg == 4326 else \
'UTM' if len(str(self.target_epsg)) == 5 and str(self.target_epsg)[:3] in ['326', '327'] else \
'NA'
if self.target_projection_type == 'Geographic':
self.target_epsg = 4326
# set target coordinate grid to the UTM EnMAP grid if no other grid is provided and target projection is UTM
self.target_coord_grid = \
self.target_coord_grid if self.target_coord_grid else \
enmap_coordinate_grid_utm if self.target_projection_type == 'UTM' else None
# bug warning regarding holes in bilinear resampling output
if self.target_projection_type == 'Geographic' and self.ortho_resampAlg == 'bilinear':
warnings.warn("There is currently a bug that causes holes in the bilinear resampling results if the "
"target projection is 'Geographic'. It is recommended to use 'nearest' or 'gauss' instead.",
UserWarning)
@staticmethod
def absPath(path):
return path if not path or os.path.isabs(path) else os.path.abspath(os.path.join(path_enptlib, path))
......
......@@ -21,7 +21,14 @@
},
"output": {
"output_dir": "" /*output directory where processed data and log files are saved*/
"output_dir": "", /*output directory where processed data and log files are saved*/
"output_format": "GTiff", /*file format of all raster output files
'GTiff': GeoTIFF
'ENVI': ENVI BSQ*/
"output_interleave": "pixel" /*raster data interleaving type (default: 'pixel')
- 'band': band-sequential (BSQ),
- 'line': data interleaved-by-line (BIL; only usable for ENVI output format),
- 'pixel' data interleaved-by-pixel (BIP)*/
},
"processors": {
......@@ -69,12 +76,17 @@
"orthorectification": {
"resamp_alg": "bilinear", /*Ortho-rectification resampling algorithm ('nearest', 'bilinear', 'gauss')*/
"vswir_overlap_algorithm": "swir_only" /*Algorithm how to output the spectral bands in the VNIR/SWIR spectral overlap region
'order_by_wvl': keep spectral bands unchanged, order bands by wavelength
'average': average the spectral information within the overlap
'vnir_only': only use the VNIR bands (cut overlapping SWIR bands)
'swir_only': only use the SWIR bands (cut overlapping VNIR bands)
*/
"vswir_overlap_algorithm": "swir_only", /*Algorithm how to output the spectral bands in the VNIR/SWIR spectral overlap region
'order_by_wvl': keep spectral bands unchanged, order bands by wavelength
'average': average the spectral information within the overlap
'vnir_only': only use the VNIR bands (cut overlapping SWIR bands)
'swir_only': only use the SWIR bands (cut overlapping VNIR bands)
*/
"target_projection_type": "UTM", /*projection type of the raster output files (default: 'UTM')
('UTM': Universal Transverse Mercator, 'Geographic': Longitude/Latitude)*/
"target_epsg": "None", /*custom EPSG code of the target projection (overrides target_projection_type)*/
"target_coord_grid": "None" /*custom target coordinate grid where is output is resampled to
([x0, x1, y0, y1], e.g., [0, 30, 0, 30])*/
}
}
}
\ No newline at end of file
......@@ -45,7 +45,7 @@ enpt_schema_input = dict(
average_elevation=dict(type='integer', required=False),
path_l1b_snr_model=dict(type='string', required=False),
working_dir=dict(type='string', required=False, nullable=True),
n_lines_to_append=dict(type='integer', required=False, nullable=True),
n_lines_to_append=dict(type='integer', required=False, nullable=True, min=0),
disable_progress_bars=dict(type='boolean', required=False, nullable=True),
)),
......@@ -53,6 +53,8 @@ enpt_schema_input = dict(
type='dict', required=False,
schema=dict(
output_dir=dict(type='string', required=False),
output_format=dict(type='string', required=False, allowed=['GTiff', 'ENVI']),
output_interleave=dict(type='string', required=False, allowed=['band', 'line', 'pixel'])
)),
processors=dict(
......@@ -64,7 +66,7 @@ enpt_schema_input = dict(
schema=dict(
path_earthSunDist=dict(type='string', required=False),
path_solar_irr=dict(type='string', required=False),
scale_factor_toa_ref=dict(type='integer', required=False),
scale_factor_toa_ref=dict(type='integer', required=False, min=1),
)),
geometry=dict(
......@@ -83,7 +85,7 @@ enpt_schema_input = dict(
auto_download_ecmwf=dict(type='boolean', required=False),
enable_ice_retrieval=dict(type='boolean', required=False),
enable_cloud_screening=dict(type='boolean', required=False),
scale_factor_boa_ref=dict(type='integer', required=False),
scale_factor_boa_ref=dict(type='integer', required=False, min=1),
)),
smile=dict(
......@@ -109,6 +111,9 @@ enpt_schema_input = dict(
resamp_alg=dict(type='string', required=False, allowed=['nearest', 'bilinear', 'gauss']),
vswir_overlap_algorithm=dict(type='string', required=False,
allowed=['order_by_wvl', 'average', 'vnir_only', 'swir_only']),
target_projection_type=dict(type='string', required=False, allowed=['UTM', 'Geographic']),
target_epsg=dict(type='integer', required=False, nullable=True, min=0, forbidden=[0]),
target_coord_grid=dict(type='list', required=False, nullable=True, minlength=4, maxlength=4)
))
))
)
......@@ -131,6 +136,8 @@ parameter_mapping = dict(
# output
output_dir=('output', 'output_dir'),
output_format=('output', 'output_format'),
output_interleave=('output', 'output_interleave'),
# processors > toa_ref
path_earthSunDist=('processors', 'toa_ref', 'path_earthSunDist'),
......@@ -162,6 +169,9 @@ parameter_mapping = dict(
# processors > orthorectification
ortho_resampAlg=('processors', 'orthorectification', 'resamp_alg'),
vswir_overlap_algorithm=('processors', 'orthorectification', 'vswir_overlap_algorithm'),
target_projection_type=('processors', 'orthorectification', 'target_projection_type'),
target_epsg=('processors', 'orthorectification', 'target_epsg'),
target_coord_grid=('processors', 'orthorectification', 'target_coord_grid'),
)
......@@ -179,6 +189,10 @@ def get_updated_schema(source_schema, key2update, new_value):
from copy import deepcopy
tgt_schema = deepcopy(source_schema)
tgt_schema['processors']['schema']['orthorectification']['schema']['target_coord_grid'] = \
dict(type='dict', required=False, nullable=True)
return deep_update(tgt_schema, key2update, new_value)
......
......@@ -47,8 +47,7 @@ from ...model.metadata import EnMAP_Metadata_L2A_MapGeo
from ..spatial_transform import \
Geometry_Transformer, \
Geometry_Transformer_3D, \
move_extent_to_EnMAP_grid, \
get_UTMEPSG_from_LonLat_cornersXY
move_extent_to_coord_grid
__author__ = 'Daniel Scheffler'
......@@ -87,8 +86,8 @@ class Orthorectifier(object):
lons_vnir, lats_vnir = enmap_ImageL1.vnir.detector_meta.lons, enmap_ImageL1.vnir.detector_meta.lats
lons_swir, lats_swir = enmap_ImageL1.swir.detector_meta.lons, enmap_ImageL1.swir.detector_meta.lats
# get target UTM zone and common extent # TODO add optionally user defined UTM zone?
tgt_epsg = self._get_tgt_UTMepsg(enmap_ImageL1)
# get target EPSG code and common extent
tgt_epsg = enmap_ImageL1.meta.vnir.epsg_ortho
tgt_extent = self._get_common_extent(enmap_ImageL1, tgt_epsg, enmap_grid=True)
kw_init = dict(resamp_alg=self.cfg.ortho_resampAlg,
nprocs=self.cfg.CPUs,
......@@ -98,7 +97,10 @@ class Orthorectifier(object):
# increase that if the resampling result contains gaps (default is 32 but this is quite slow)
kw_init['neighbours'] = 8
kw_trafo = dict(tgt_prj=tgt_epsg, tgt_extent=tgt_extent)
kw_trafo = dict(tgt_prj=tgt_epsg, tgt_extent=tgt_extent,
tgt_coordgrid=((self.cfg.target_coord_grid['x'], self.cfg.target_coord_grid['y'])
if self.cfg.target_coord_grid else None)
)
# transform VNIR and SWIR to map geometry
GeoTransformer = Geometry_Transformer if lons_vnir.ndim == 2 else Geometry_Transformer_3D
......@@ -163,6 +165,7 @@ class Orthorectifier(object):
meta_l1b=enmap_ImageL1.meta,
wvls_l2a=L2_obj.data.meta.band_meta['wavelength'],
dims_mapgeo=L2_obj.data.shape,
grid_res_l2a=(L2_obj.data.gt[1], abs(L2_obj.data.gt[5])),
logger=L2_obj.logger)
L2_obj.meta.add_band_statistics(L2_obj.data)
......@@ -174,13 +177,8 @@ class Orthorectifier(object):
return L2_obj
@staticmethod
def _get_tgt_UTMepsg(enmap_ImageL1: EnMAPL1Product_SensorGeo) -> int:
return get_UTMEPSG_from_LonLat_cornersXY(lons=enmap_ImageL1.vnir.detector_meta.lon_UL_UR_LL_LR,
lats=enmap_ImageL1.vnir.detector_meta.lat_UL_UR_LL_LR)
@staticmethod
def _get_common_extent(enmap_ImageL1: EnMAPL1Product_SensorGeo,
def _get_common_extent(self,
enmap_ImageL1: EnMAPL1Product_SensorGeo,
tgt_epsg: int,
enmap_grid: bool = True) -> Tuple[float, float, float, float]:
"""Get common target extent for VNIR and SWIR.
......@@ -199,37 +197,39 @@ class Orthorectifier(object):
S_UL_UR_LL_LR_ll = [(S_lons[y, x], S_lats[y, x]) for y, x in [(0, 0), (0, -1), (-1, 0), (-1, -1)]]
# transform them to UTM
V_UL_UR_LL_LR_utm = [transform_any_prj(EPSG2WKT(4326), EPSG2WKT(tgt_epsg), x, y) for x, y in V_UL_UR_LL_LR_ll]
S_UL_UR_LL_LR_utm = [transform_any_prj(EPSG2WKT(4326), EPSG2WKT(tgt_epsg), x, y) for x, y in S_UL_UR_LL_LR_ll]
V_UL_UR_LL_LR_prj = [transform_any_prj(EPSG2WKT(4326), EPSG2WKT(tgt_epsg), x, y) for x, y in V_UL_UR_LL_LR_ll]
S_UL_UR_LL_LR_prj = [transform_any_prj(EPSG2WKT(4326), EPSG2WKT(tgt_epsg), x, y) for x, y in S_UL_UR_LL_LR_ll]
# separate X and Y
V_X_utm, V_Y_utm = zip(*V_UL_UR_LL_LR_utm)
S_X_utm, S_Y_utm = zip(*S_UL_UR_LL_LR_utm)
V_X_prj, V_Y_prj = zip(*V_UL_UR_LL_LR_prj)
S_X_prj, S_Y_prj = zip(*S_UL_UR_LL_LR_prj)
# in case of 3D geolayers, the corner coordinates have multiple values for multiple bands
# -> use the innermost coordinates to avoid pixels with VNIR-only/SWIR-only values due to keystone
# (these pixels would be set to nodata later anyways, so we don need to increase the extent for them)
# (these pixels would be set to nodata later anyways, so we don't need to increase the extent for them)
if V_lons.ndim == 3:
V_X_utm = (V_X_utm[0].max(), V_X_utm[1].min(), V_X_utm[2].max(), V_X_utm[3].min())
V_Y_utm = (V_Y_utm[0].min(), V_Y_utm[1].min(), V_Y_utm[2].max(), V_Y_utm[3].max())
S_X_utm = (S_X_utm[0].max(), S_X_utm[1].min(), S_X_utm[2].max(), S_X_utm[3].min())
S_Y_utm = (S_Y_utm[0].min(), S_Y_utm[1].min(), S_Y_utm[2].max(), S_Y_utm[3].max())
V_X_prj = (V_X_prj[0].max(), V_X_prj[1].min(), V_X_prj[2].max(), V_X_prj[3].min())
V_Y_prj = (V_Y_prj[0].min(), V_Y_prj[1].min(), V_Y_prj[2].max(), V_Y_prj[3].max())
S_X_prj = (S_X_prj[0].max(), S_X_prj[1].min(), S_X_prj[2].max(), S_X_prj[3].min())
S_Y_prj = (S_Y_prj[0].min(), S_Y_prj[1].min(), S_Y_prj[2].max(), S_Y_prj[3].max())
# use inner coordinates of VNIR and SWIR as common extent
xmin_utm = max([min(V_X_utm), min(S_X_utm)])
ymin_utm = max([min(V_Y_utm), min(S_Y_utm)])
xmax_utm = min([max(V_X_utm), max(S_X_utm)])
ymax_utm = min([max(V_Y_utm), max(S_Y_utm)])
common_extent_utm = (xmin_utm, ymin_utm, xmax_utm, ymax_utm)
xmin_prj = max([min(V_X_prj), min(S_X_prj)])
ymin_prj = max([min(V_Y_prj), min(S_Y_prj)])
xmax_prj = min([max(V_X_prj), max(S_X_prj)])
ymax_prj = min([max(V_Y_prj), max(S_Y_prj)])
common_extent_prj = (xmin_prj, ymin_prj, xmax_prj, ymax_prj)
# move the extent to the EnMAP coordinate grid
if enmap_grid:
common_extent_utm = move_extent_to_EnMAP_grid(common_extent_utm)
if enmap_grid and self.cfg.target_coord_grid:
common_extent_prj = move_extent_to_coord_grid(common_extent_prj,
self.cfg.target_coord_grid['x'],
self.cfg.target_coord_grid['y'],)
enmap_ImageL1.logger.info('Computed common target extent of orthorectified image (xmin, ymin, xmax, ymax in '
'EPSG %s): %s' % (tgt_epsg, str(common_extent_utm)))
'EPSG %s): %s' % (tgt_epsg, str(common_extent_prj)))
return common_extent_utm
return common_extent_prj
class VNIR_SWIR_Stacker(object):
......
......@@ -29,7 +29,7 @@
"""EnPT module 'spatial transform', containing everything related to spatial transformations."""
from typing import Union, Tuple, List, Optional # noqa: F401
from typing import Union, Tuple, List, Optional, Sequence # noqa: F401
from multiprocessing import Pool, cpu_count
from collections import OrderedDict
import numpy as np
......@@ -41,11 +41,11 @@ import numpy_indexed as npi
from sensormapgeo import SensorMapGeometryTransformer, SensorMapGeometryTransformer3D
from sensormapgeo.transformer_2d import AreaDefinition
from py_tools_ds.geo.projection import get_proj4info, prj_equal, EPSG2WKT, WKT2EPSG, proj4_to_WKT
from py_tools_ds.geo.projection import prj_equal, EPSG2WKT
from py_tools_ds.geo.coord_grid import find_nearest
from py_tools_ds.geo.coord_trafo import transform_any_prj, transform_coordArray
from ...options.config import enmap_coordinate_grid
from ...options.config import enmap_coordinate_grid_utm
__author__ = 'Daniel Scheffler'
......@@ -80,10 +80,6 @@ class Geometry_Transformer(SensorMapGeometryTransformer):
if data_sensorgeo.is_map_geo:
raise RuntimeError('The dataset to be transformed into map geometry already represents map geometry.')
# get EnMAP grid
tgt_epsg = WKT2EPSG(proj4_to_WKT(get_proj4info(proj=tgt_prj)))
tgt_coordgrid = (enmap_coordinate_grid['x'], enmap_coordinate_grid['y']) if tgt_epsg != 4326 else None
# run transformation (output extent/area definition etc. is internally computed from the geolayers if not given)
out_data, out_gt, out_prj = \
super(Geometry_Transformer, self).to_map_geometry(data_sensorgeo[:],
......@@ -126,10 +122,6 @@ class Geometry_Transformer_3D(SensorMapGeometryTransformer3D):
if data_sensorgeo.is_map_geo:
raise RuntimeError('The dataset to be transformed into map geometry already represents map geometry.')
# get EnMAP grid
tgt_epsg = WKT2EPSG(proj4_to_WKT(get_proj4info(proj=tgt_prj)))
tgt_coordgrid = (enmap_coordinate_grid['x'], enmap_coordinate_grid['y']) if tgt_epsg != 4326 else None
# run transformation (output extent/area definition etc. is internally computed from the geolayers if not given)
out_data, out_gt, out_prj = \
super(Geometry_Transformer_3D, self).to_map_geometry(data_sensorgeo[:],
......@@ -234,14 +226,17 @@ class VNIR_SWIR_SensorGeometryTransformer(object):
return tgt_data_sensorgeo
def move_extent_to_EnMAP_grid(extent_utm: Tuple[float, float, float, float]) -> Tuple[float, float, float, float]:
"""Move a UTM coordinate extent to the EnMAP coordinate grid (30m x 30m, origin at 0/0).
def move_extent_to_coord_grid(extent_utm: Tuple[float, float, float, float],
tgt_xgrid: Sequence[float],
tgt_ygrid: Sequence[float],
) -> Tuple[float, float, float, float]:
"""Move the given coordinate extent to a coordinate grid.
:param extent_utm: xmin, ymin, xmax, ymax coordinates
:param tgt_xgrid: target X coordinate grid, e.g. [0, 30]
:param tgt_ygrid: target Y coordinate grid, e.g. [0, 30]
"""
xmin, ymin, xmax, ymax = extent_utm
tgt_xgrid = enmap_coordinate_grid['x']
tgt_ygrid = enmap_coordinate_grid['y']
tgt_xmin = find_nearest(tgt_xgrid, xmin, roundAlg='off', extrapolate=True)
tgt_xmax = find_nearest(tgt_xgrid, xmax, roundAlg='on', extrapolate=True)
tgt_ymin = find_nearest(tgt_ygrid, ymin, roundAlg='off', extrapolate=True)
......@@ -443,7 +438,9 @@ class RPC_Geolayer_Generator(object):
ymin, ymax = cornerCoordsUTM[:, 1].min(), cornerCoordsUTM[:, 1].max()
# get UTM bounding box and move it to the EnMAP grid
xmin, ymin, xmax, ymax = move_extent_to_EnMAP_grid((xmin, ymin, xmax, ymax))
xmin, ymin, xmax, ymax = \
move_extent_to_coord_grid((xmin, ymin, xmax, ymax),
enmap_coordinate_grid_utm['x'], enmap_coordinate_grid_utm['y'])
# set up a regular grid of UTM points with a specific mesh width
meshwidth = 300 # 10 EnMAP pixels
......
......@@ -27,6 +27,6 @@
# You should have received a copy of the GNU Lesser General Public License along
# with this program. If not, see <http://www.gnu.org/licenses/>.
__version__ = '0.12.8'
__versionalias__ = '20200513.01'
__version__ = '0.13.0'
__versionalias__ = '20200518.01'
__author__ = 'Daniel Scheffler'
......@@ -69,6 +69,18 @@ class Test_CLIParser(TestCase):
self.assertNotIsInstance(config.CPUs, str)
self.assertEqual(config.CPUs, cpu_count())
def test_param_list(self):
parsed_args = self.parser_run.parse_args(self.baseargs +
['--target_coord_grid', '0', '30', '0', '30'])
self.get_config(parsed_args) # we don't check the result here as EnPT_Config generates a dict from it
try:
self.parser_run.parse_args(self.baseargs + ['--target_coord_grid', '0', '30'])
except SystemExit as e:
assert isinstance(e.__context__, ArgumentError)
else:
raise ValueError("Exception not raised")
def test_param_boolean(self):
parsed_args = self.parser_run.parse_args(self.baseargs +
['--enable_ac', 'True'])
......
......@@ -49,7 +49,7 @@ from enpt.processors.spatial_transform import \
Geometry_Transformer, RPC_Geolayer_Generator, RPC_3D_Geolayer_Generator, VNIR_SWIR_SensorGeometryTransformer
from enpt.options.config import config_for_testing, EnPTConfig, config_for_testing_dlr
from enpt.io.reader import L1B_Reader
from enpt.options.config import enmap_coordinate_grid
from enpt.options.config import enmap_coordinate_grid_utm
__author__ = 'Daniel Scheffler'
......@@ -83,11 +83,14 @@ class Test_Geometry_Transformer(TestCase):
GT.to_map_geometry(self.gA2transform_mapgeo, tgt_prj=32632)
# test transformation to UTM zone 32
data_mapgeo, gt, prj = GT.to_map_geometry(self.gA2transform_sensorgeo, tgt_prj=32632)
self.assertEqual((gt[1], -gt[5]), (np.ptp(enmap_coordinate_grid['x']),
np.ptp(enmap_coordinate_grid['x']))) # 30m output
data_mapgeo, gt, prj = GT.to_map_geometry(self.gA2transform_sensorgeo, tgt_prj=32632,
tgt_coordgrid=(enmap_coordinate_grid_utm['x'],
enmap_coordinate_grid_utm['y']))
self.assertEqual((gt[1], -gt[5]), (np.ptp(enmap_coordinate_grid_utm['x']),
np.ptp(enmap_coordinate_grid_utm['y']))) # 30m output
self.assertTrue(is_point_on_grid((gt[0], gt[3]),
xgrid=enmap_coordinate_grid['x'], ygrid=enmap_coordinate_grid['y']))
xgrid=enmap_coordinate_grid_utm['x'],
ygrid=enmap_coordinate_grid_utm['y']))
# test transformation to LonLat
GT.to_map_geometry(self.gA2transform_sensorgeo, tgt_prj=4326)
......