Commit ff2f340f authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Fixed AssertionError 'exactly 4 image corners must be present within the...

Fixed AssertionError 'exactly 4 image corners must be present within the dataset'. Updated version info.
Former-commit-id: e75b3615
Former-commit-id: e4bf232e
parent 0d15094c
......@@ -13,8 +13,8 @@ from .processing.process_controller import process_controller # noqa: E402
__author__ = """Daniel Scheffler"""
__email__ = 'daniel.scheffler@gfz-potsdam.de'
__version__ = '0.8.9'
__versionalias__ = '20171018.04'
__version__ = '0.8.10'
__versionalias__ = '20171020.01'
__all__ = ['algorithms',
'io',
'misc',
......
......@@ -203,6 +203,7 @@ class L1A_object(GMS_object):
i, list_matching_dsIdx = 0, []
while True: # Get dataset indices within HDF file
# noinspection PyBroadException
try:
ds = hdfFile.select(i)
if subsystem_identifier in str(ds.dimensions()) and 'ImagePixel' in str(ds.dimensions()):
......@@ -496,10 +497,11 @@ class L1A_object(GMS_object):
'mask_nodata') and self.mask_nodata is not None, "The L1A object needs to have a nodata mask."
self.logger.info('Calculating true data corner positions (image and world coordinates)...')
if re.search('ETM+', self.sensor) and self.acq_datetime > datetime.datetime(year=2003, month=5, day=31,
tzinfo=datetime.timezone.utc):
# if re.search('ETM+', self.sensor) and self.acq_datetime > datetime.datetime(year=2003, month=5, day=31,
# tzinfo=datetime.timezone.utc):
if is_dataset_provided_as_fullScene(self.GMS_identifier):
self.trueDataCornerPos = calc_FullDataset_corner_positions(self.mask_nodata, algorithm='numpy',
assert_four_corners=False)
assert_four_corners=True)
else:
self.trueDataCornerPos = calc_FullDataset_corner_positions(self.mask_nodata, assert_four_corners=False)
......@@ -546,9 +548,10 @@ class L1A_object(GMS_object):
def calc_orbit_overpassParams(self):
"""Calculate orbit parameters."""
self.MetaObj.overpassDurationSec, self.MetaObj.scene_length = \
self.MetaObj.get_overpassDuration_SceneLength(
self.fullSceneCornerLonLat, self.fullSceneCornerPos, self.shape_fullArr, self.logger)
if is_dataset_provided_as_fullScene(self.GMS_identifier):
self.MetaObj.overpassDurationSec, self.MetaObj.scene_length = \
self.MetaObj.get_overpassDuration_SceneLength(
self.fullSceneCornerLonLat, self.fullSceneCornerPos, self.shape_fullArr, self.logger)
def add_rasterInfo_to_MetaObj(self, custom_rasObj=None):
"""
......@@ -602,13 +605,15 @@ class L1A_object(GMS_object):
def calc_mean_VAA(self):
"""Calculates mean viewing azimuth angle using sensor flight line derived from full scene corner coordinates."""
if re.search('Sentinel-2', self.satellite, re.I):
if is_dataset_provided_as_fullScene(self.GMS_identifier):
self.VAA_mean = \
GEOP.calc_VAA_using_fullSceneCornerLonLat(self.fullSceneCornerLonLat, self.MetaObj.orbitParams)
else:
# e.g. Sentinel-2 / RapidEye
self.logger.warning('No precise calculation of mean viewing azimuth angle possible because orbit track '
'cannot be reconstructed from dataset since full scene corner positions are unknown. '
'Mean VAA angle is filled with the mean value of the viewing azimuth array provided '
'in metadata.')
self.VAA_mean = self.MetaObj.IncidenceAngle
else:
self.VAA_mean = \
GEOP.calc_VAA_using_fullSceneCornerLonLat(self.fullSceneCornerLonLat, self.MetaObj.orbitParams)
self.logger.info('Calculation of mean VAA...: %s' % round(self.VAA_mean, 2))
......@@ -1163,9 +1163,12 @@ def calc_VAA_using_fullSceneCornerLonLat(fullSceneCornerLonLat, orbit_params):
"""Calculates the Viewing azimuth angle (defined as 90 degrees from the flight line),
e.g. if flight line is 8 degrees from North -> VAA will be 98 degrees.
:param fullSceneCornerLonLat:
:param fullSceneCornerLonLat: UL, UR, LL, LR
:param orbit_params: list of [altitude, inclination, period] => inclination is used as fallback
"""
assert len(fullSceneCornerLonLat) == 4, \
'VAA can only be calculated with fullSceneCornerLonLat representing 4 coordinates (UL, UR, LL, LR).'
UL_LonLat, UR_LonLat, LL_LonLat, LR_LonLat = fullSceneCornerLonLat
forward_az_left = pyproj.Geod(ellps='WGS84').inv(*LL_LonLat, *UL_LonLat)[0]
forward_az_right = pyproj.Geod(ellps='WGS84').inv(*LR_LonLat, *UR_LonLat)[0]
......@@ -1175,7 +1178,7 @@ def calc_VAA_using_fullSceneCornerLonLat(fullSceneCornerLonLat, orbit_params):
if abs(VAA_mean - 90) < 1:
# fullSceneCornerLonLat obviously don't belong to a full scene but a granule
assert orbit_params
warnings.warn('Derivation of mean VAA angle from flight line delivered a non reasonable value (%s degress).'
warnings.warn('Derivation of mean VAA angle from flight line delivered a non reasonable value (%s degrees).'
'Using sensor inclination (%s degrees) as fallback.' % (VAA_mean, orbit_params[1]))
VAA_mean = float(orbit_params[1]) # inclination # FIXME is this correct?
......@@ -1304,7 +1307,7 @@ def calc_SZA_SAA_array(shape_fullArr, arr_pos, AcqDate, CenterAcqTime, fullScene
:param AcqDate:
:param CenterAcqTime:
:param fullSceneCornerPos:
:param fullSceneCornerLonLat:
:param fullSceneCornerLonLat: UL, UR, LL, LR
:param overpassDurationSec:
:param logger:
:param meshwidth: <int> defines the density of the mesh used for generating the output
......
......@@ -153,7 +153,7 @@ def is_dataset_provided_as_fullScene(GMS_identifier):
sensorcode = get_GMS_sensorcode(GMS_identifier)
dict_fullScene_or_tiles = {
'AVNIR-2': True,
'AST_full': False,
'AST_full': True,
'AST_V1': True,
'AST_V2': True,
'AST_S': True,
......@@ -173,7 +173,7 @@ def is_dataset_provided_as_fullScene(GMS_identifier):
'SPOT4b': True,
'SPOT5b': True,
'RE5': False,
'S2A_full': False, # FIXME this changed for S2 in 08/2016
'S2A_full': False,
'S2A10': False,
'S2A20': False,
'S2A60': False,
......
......@@ -187,6 +187,7 @@ class Dataset(object):
@property
def arr(self):
# type: () -> GeoArray
# TODO this must return a subset if self.subset is not None
return self._arr
......
......@@ -29,7 +29,7 @@ from gms_preprocessing.algorithms import geoprocessing as GEOP
from gms_preprocessing.misc import helper_functions as HLP_F
from gms_preprocessing.misc import database_tools as DB_T
from gms_preprocessing.misc.path_generator import path_generator, get_path_ac_options
from gms_preprocessing.misc.definition_dicts import get_GMS_sensorcode
from gms_preprocessing.misc.definition_dicts import get_GMS_sensorcode, is_dataset_provided_as_fullScene
from sicor.options import get_options as get_ac_options
......@@ -254,6 +254,7 @@ class METADATA(object):
# Solar irradiance, central wavelengths , full width half maximum per band
self.wvlUnit = 'Nanometers'
# derive number of bands from number of given gains/offsets in header file or from raster data
# noinspection PyBroadException
try:
self.nBands = (np.mean([len(self.Gains), len(self.Offsets)])
if np.std([len(self.Gains), len(self.Offsets)]) == 0 and len(self.Gains) != 0 else
......@@ -743,6 +744,7 @@ class METADATA(object):
# Solar irradiance, central wavelengths , full width half maximum per band
self.wvlUnit = 'Nanometers'
# derive number of bands from number of given gains/offsets in header file or from raster data
# noinspection PyBroadException
try:
self.nBands = (np.mean([len(self.Gains), len(self.Offsets)])
if np.std([len(self.Gains), len(self.Offsets)]) == 0 and len(self.Gains) != 0 else
......@@ -1109,6 +1111,7 @@ class METADATA(object):
# Solar irradiance, central wavelengths, full width half maximum per band
self.wvlUnit = 'Nanometers'
# derive number of bands from number of given gains/offsets in header file or from raster data
# noinspection PyBroadException
try:
self.nBands = (np.mean([len(self.Gains), len(self.Offsets)])
if np.std([len(self.Gains), len(self.Offsets)]) == 0 and len(self.Gains) != 0 else
......@@ -1597,6 +1600,8 @@ class METADATA(object):
:param fullSceneCornerLonLat:
:param logger:
"""
assert is_dataset_provided_as_fullScene(self.GMS_identifier) and len(fullSceneCornerLonLat) == 4, \
'Center acquisition time can only be computed for datasets provided as full scenes, not for tiles.'
ul, lr = fullSceneCornerLonLat[0], fullSceneCornerLonLat[3]
center_coord = [np.mean([ul[0], lr[0]]), np.mean([ul[1], lr[1]])]
......@@ -1629,22 +1634,30 @@ class METADATA(object):
:param shape_fullArr:
:param logger:
"""
if fullSceneCornerPos != list(([0, 0], [0, shape_fullArr[1] - 1],
[shape_fullArr[0] - 1, 0], [shape_fullArr[0] - 1, shape_fullArr[1] - 1])):
orbitAltitudeKm, orbitPeriodMin = self.orbitParams[0], self.orbitParams[2]
UL, UR, LL, LR = fullSceneCornerLonLat
geod = pyproj.Geod(ellps='WGS84')
scene_length = np.mean([geod.inv(UL[0], UL[1], LL[0], LL[1])[2],
geod.inv(UR[0], UR[1], LR[0], LR[1])[2]]) / 1000 # along-track distance [km]
logger.info('Calculation of scene length...: %s km' % round(scene_length, 2))
orbitPeriodLength = 2 * math.pi * (6371 + orbitAltitudeKm)
overpassDurationSec = (scene_length / orbitPeriodLength) * orbitPeriodMin * 60.
logger.info('Calculation of overpass duration...: %s sec' % round(overpassDurationSec, 2))
return overpassDurationSec, scene_length
else:
logger.warning('Overpass duration and scene length cannot be calculated because the given data represents'
'a subset of the original scene.')
assert is_dataset_provided_as_fullScene(self.GMS_identifier) and len(fullSceneCornerLonLat) == 4, \
'Overpass duration and scene length can only be computed for datasets provided as full scenes, not for ' \
'tiles.'
# check if current scene is a subset
assert fullSceneCornerPos != list(([0, 0], [0, shape_fullArr[1] - 1],
[shape_fullArr[0] - 1, 0], [shape_fullArr[0] - 1, shape_fullArr[1] - 1])),\
'Overpass duration and scene length cannot be calculated because the given data represents a subset of ' \
'the original scene.'
# compute scene length
orbitAltitudeKm, orbitPeriodMin = self.orbitParams[0], self.orbitParams[2]
UL, UR, LL, LR = fullSceneCornerLonLat
geod = pyproj.Geod(ellps='WGS84')
scene_length = np.mean([geod.inv(UL[0], UL[1], LL[0], LL[1])[2],
geod.inv(UR[0], UR[1], LR[0], LR[1])[2]]) / 1000 # along-track distance [km]
logger.info('Calculation of scene length...: %s km' % round(float(scene_length), 2))
# compute overpass duration
orbitPeriodLength = 2 * math.pi * (6371 + orbitAltitudeKm)
overpassDurationSec = (scene_length / orbitPeriodLength) * orbitPeriodMin * 60.
logger.info('Calculation of overpass duration...: %s sec' % round(overpassDurationSec, 2))
return overpassDurationSec, scene_length
def filter_layerdependent_metadata(self):
FULL_LayerBandsAssignment = \
......
......@@ -25,7 +25,7 @@ test_requirements = requirements + ['coverage', 'nose', 'nose-htmloutput', 'redn
setup(
name='gms_preprocessing',
version='0.8.9',
version='0.8.10',
description="GeoMultiSens - Scalable Multi-Sensor Analysis of Remote Sensing Data",
long_description=readme + '\n\n' + history,
author="Daniel Scheffler",
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment