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

revised calculation of acquisition geometry arrays for Landsat

algorithms.GEOPROCESSING:
- calc_VAA_using_fullSceneCornerLonLat():
    - implemented keyword 'orbit_params' as fallback; added docstring
    -added 90 degrees to VAA
- calc_VZA_array(): added some notes
- calc_RAA_array(): added docstring

algorithms.L1A_P.L1A_object():
- calc_mean_VAA(): updated calc_VAA_using_fullSceneCornerLonLat call

algorithms.L1C_P.L1C_object():
- VZA_arr.getter: added parameter 'meshwidth'
- revised VAA_arr.getter

algorithms.METADATA:
- get_orbit_params(): added some notes

- updated __version__


Former-commit-id: b8a055bf
parent 93bc4d81
......@@ -15,7 +15,7 @@ from . import config
from .processing.process_controller import process_controller
__version__ = '20170126.01'
__version__ = '20170126.02'
__author__ = 'Daniel Scheffler'
__all__ = ['algorithms',
'io',
......
......@@ -23,6 +23,7 @@ import re
import subprocess
import sys
import time
import warnings
import numpy as np
import scipy
......@@ -3055,12 +3056,28 @@ def get_lonlat_coord_array(shape_fullArr, arr_pos, geotransform, projection, mes
return xcoord_ycoord_arr.astype(np.float32), [UL_lonlat, UR_lonlat, LL_lonlat, LR_lonlat]
def calc_VAA_using_fullSceneCornerLonLat(fullSceneCornerLonLat):
# type: (list) -> float
def calc_VAA_using_fullSceneCornerLonLat(fullSceneCornerLonLat, orbit_params):
# type: (list, list) -> float
"""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 orbit_params: list of [altitude, inclination, period] => inclination is used as fallback
"""
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]
VAA_mean = float(np.mean([forward_az_left, forward_az_right]))
VAA_mean = float(np.mean([forward_az_left, forward_az_right])) + 90
# validation:
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).'
'Using sensor inclination (%s degrees) as fallback.' %(VAA_mean, orbit_params[1]))
VAA_mean = float(orbit_params[1]) # inclination # FIXME is this correct?
return VAA_mean
......@@ -3080,6 +3097,11 @@ def calc_VZA_array(shape_fullArr, arr_pos, fullSceneCornerPos, viewing_angle, FO
:param nodata_mask: <numpy array>, used for declaring nodata values in the output VZA array
:param outFill: the value that is assigned to nodata area in the output VZA array"""
# FIXME in case of Sentinel-2 the viewing_angle corresponds to the center point of the image footprint (trueDataCornerPos)
# FIXME => the algorithm must use the center viewing angle + orbit inclination and must calculate the FOV to be used
# FIXME via the widths of the footprint at the center position
# FIXME since S2 brings its own VZA array this is only relevant other scenes provided as granules (e.g. RapidEye)
if nodata_mask is not None: assert isinstance(nodata_mask, (GeoArray,np.ndarray)), \
"'nodata_mask' must be a numpy array or an instance of GeoArray. Got %s" % type(nodata_mask)
colsFullArr, rowsFullArr, bandsFullArr = shape_fullArr
......@@ -3100,10 +3122,10 @@ def calc_VZA_array(shape_fullArr, arr_pos, fullSceneCornerPos, viewing_angle, FO
[1, ur[1], ur[0], ur[1] * ur[0]],
[1, ll[1], ll[0], ll[1] * ll[0]],
[1, lr[1], lr[0], lr[1] * lr[0]]])
const_matrix = np.array([viewing_angle - FOV / 2., # VZA@UL
viewing_angle + FOV / 2., # VZA@UR
viewing_angle - FOV / 2., # VZA@LL
viewing_angle + FOV / 2.]) # VZA@LR
const_matrix = np.array([viewing_angle + FOV / 2., # VZA@UL # +- aligning seems to match with Sentinel-2 # TODO correct?
viewing_angle - FOV / 2., # VZA@UR
viewing_angle + FOV / 2., # VZA@LL
viewing_angle - FOV / 2.]) # VZA@LR
factors = np.linalg.solve(coeff_matrix, const_matrix)
VZA_array = (factors[0] + factors[1] * cols_arr + factors[2] * rows_arr + factors[3] *
......@@ -3255,6 +3277,14 @@ def calc_SZA_SAA_array(shape_fullArr, arr_pos, AcqDate, CenterAcqTime, fullScene
def calc_RAA_array(SAA_array, VAA_array, nodata_mask=None, outFill=None):
"""Calculate relative azimuth angle between solar azimuth and viewing azimuth in degrees.
:param SAA_array:
:param VAA_array:
:param nodata_mask:
:param outFill: the value to be used to fill areas outside the actual image bounds
:return:
"""
if nodata_mask is not None: assert isinstance(nodata_mask, (GeoArray, np.ndarray)), \
"'nodata_mask' must be a numpy array or an instance of GeoArray. Got %s" % type(nodata_mask)
......
......@@ -927,7 +927,7 @@ class L1A_object(GMS_object):
def calc_mean_VAA(self):
"""Calculates mean viewing azimuth angle using sensor flight line derived from corner coordinates."""
"""Calculates mean viewing azimuth angle using sensor flight line derived from full scene corner coordinates."""
if re.search('Sentinel-2',self.satellite,re.I):
self.logger.warning('No precise calculation of mean viewing azimuth angle possible because orbit track '
......@@ -936,6 +936,7 @@ class L1A_object(GMS_object):
'in metadata.')
self.VAA_mean = self.MetaObj.IncidenceAngle
else:
self.VAA_mean = GEOP.calc_VAA_using_fullSceneCornerLonLat(self.fullSceneCornerLonLat)
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))
......@@ -77,7 +77,7 @@ class L1C_object(L1B_object):
GEOP.get_lonlat_coord_array(self.shape_fullArr, self.arr_pos,
mapinfo2geotransform(self.meta_odict['map info']),
self.meta_odict['coordinate system string'],
meshwidth=10,
meshwidth=10, # for faster processing
nodata_mask=None, # dont overwrite areas outside the image with nodata
outFill=get_outFillZeroSaturated(np.float32)[0])[0]
return self._lonlat_arr
......@@ -105,11 +105,12 @@ class L1C_object(L1B_object):
bandwise=0)
else:
self._VZA_arr = GEOP.calc_VZA_array(self.shape_fullArr, self.arr_pos, self.fullSceneCornerPos,
float(self.meta_odict['ViewingAngle']),
float(self.meta_odict['FieldOfView']),
self.logger,
nodata_mask=None, # dont overwrite areas outside the image with nodata
outFill=get_outFillZeroSaturated(np.float32)[0])
float(self.meta_odict['ViewingAngle']),
float(self.meta_odict['FieldOfView']),
self.logger,
nodata_mask=None, # dont overwrite areas outside the image with nodata
outFill=get_outFillZeroSaturated(np.float32)[0],
meshwidth=10) # for faster processing
return self._VZA_arr
......@@ -136,10 +137,11 @@ class L1C_object(L1B_object):
else:
# only a mean VAA is available
if self.VAA_mean is None:
self.VAA_mean = GEOP.calc_VAA_using_fullSceneCornerLonLat(self.fullSceneCornerLonLat)
self.VAA_mean = \
GEOP.calc_VAA_using_fullSceneCornerLonLat(self.fullSceneCornerLonLat, self.MetaObj.orbitParams)
assert isinstance(self.VAA_mean, float)
self._VAA_arr = np.full(self.shape_fullArr, self.VAA_mean, np.float32)
self._VAA_arr = np.full(self.VZA_arr.shape, self.VAA_mean, np.float32)
return self._VAA_arr
......
......@@ -1431,6 +1431,7 @@ class METADATA(object):
ENVI file headers in the same order.
"""
# FIXME orbit params are missing
# descr_dic = { ### FillZeroSaturated von HLP_F ausgeben lassen
# 'ALOS_Rad' :"(1) GEOCODED Level1B2 Product; '"+self.Dataname+"'\n (2) Int16 RadianceData in [W * m-2 * sr-1 * micrometer-1]*10; radiance scale factor: 10 (fillPixels: -99, zeroPixels:0, saturatedPixels: 32767 (Max of Int16))'",
# 'ALOS_Ref' :"(1) Orthorectified JAXA or GFZ; '"+self.Dataname+"'\n (2) Int16 TOA_Reflectance in [0-10000]; radiance scale factor: 10 (fillPixels: -99, zeroPixels:0, saturatedPixels: 32767 (Max of Int16))'",
......@@ -1814,6 +1815,8 @@ def get_FieldOfView(GMS_identifier):
def get_orbit_params(GMS_identifier):
GMS_sensorcode = get_GMS_sensorcode(GMS_identifier)
# sensor altitude above ground [kilometers]
dict_altitude = {'AVNIR-2':691.65,
'AST_V1':705, 'AST_V2':705, 'AST_S':705, 'AST_T':705,
'TM4':705,'TM5':705,'TM7':705,'LDCM':705,
......@@ -1822,6 +1825,8 @@ def get_orbit_params(GMS_identifier):
'RE5':630,
'S2A10':786,'S2A20':786,'S2A60':786,
'S2B10':786,'S2B20':786,'S2B60':786 }
# sensor inclination [degrees]
dict_inclination = {'AVNIR-2':98.16,
'AST_V1':98.3, 'AST_V2':98.3, 'AST_S':98.3, 'AST_T':98.3,
'TM4':98.2,'TM5':98.2,'TM7':98.2,'LDCM':98.2,
......@@ -1830,6 +1835,8 @@ def get_orbit_params(GMS_identifier):
'RE5':98,
'S2A10':98.62,'S2A20':98.62,'S2A60':98.62,
'S2B10':98.62,'S2B20':98.62,'S2B60':98.62 }
# time needed for one complete earth revolution [minutes]
dict_period = {'AVNIR-2':98.7,
'AST_V1':98.88, 'AST_V2':98.88, 'AST_S':98.88, 'AST_T':98.88,
'TM4':98.9,'TM5':98.9,'TM7':98.9,'LDCM':98.9,
......@@ -1838,6 +1845,7 @@ def get_orbit_params(GMS_identifier):
'RE5':96.7,
'S2A10':100.6,'S2A20':100.6,'S2A60':100.6,
'S2B10':100.6,'S2B20':100.6,'S2B60':100.6 }
return [dict_altitude[GMS_sensorcode],dict_inclination[GMS_sensorcode],dict_period[GMS_sensorcode]]
def get_special_values(GMS_identifier):
......
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