Skip to content
Snippets Groups Projects
Commit eae30170 authored by Paolo Veglio's avatar Paolo Veglio
Browse files

reworked read_data with annotations and attrs

parent e699059e
No related branches found
No related tags found
No related merge requests found
......@@ -8,6 +8,7 @@ import os
import logging
from datetime import datetime as dt
from typing import Dict
from attrs import define, field, validators
_DTR = np.pi/180.
_RTD = 180./np.pi
......@@ -18,121 +19,264 @@ logging.basicConfig(level=logging.INFO, format='%(name)s - %(levelname)s - %(mes
# format='%(name)s %(levelname)s %(message)s')
def read_viirs_data(l1b_filename: str,
geo_filename: str) -> xr.Dataset:
"""Read VIIRS MOD or IMG data
@define(kw_only=True, slots=True)
class CollectInputs(object):
"""Class that collects all input files
Parameters
----------
l1b_filename: str
L1b VIIRS filename
geo_filename: str
geoocation VIIRS filename
Returns
-------
in_data: xarray.Dataset
dataset containing VIIRS bands, geolocation data and additional info (i.e. relative azimuth,
sunglint angle and scattering angle)
file_name_geo: str
name of VIIRS geolocation file
file_name_l1b: str
name of VIIRS L1b file
ancillary_dir: str
path for the ancillary files. (this is necessary because the C functions
still require path and file name separately)
sst_file: str
file name for the Reynolds SST
ndvi_file: str
file name for the NDVI hdf4 file
geos_file_before: str
file name of the GEOS-5 containing atmospheric profiles before the VIIRS timestamp
geos_file_after: str
file name of the GEOS-5 containing atmospheric profiles after the VIIRS timestamp
geos_land: str
file name of the GEOS-5 file that has the land ice/snow information
geos_ocean: str
file name of the GEOS-5 file that has the ocean information
geos_constants: str
file name for the GEOS-5 constans
dims [optional]: str
name of the dimensions for the arrays in the xarray.Dataset output
"""
file_name_geo: str = field([validators.of_instance(str), ])
file_name_l1b: str = field([validators.of_instance(str), ])
ancillary_dir: str = field([validators.of_instance(str), ])
sst_file: str = field([validators.of_instance(str), ])
ndvi_file: str = field([validators.of_instance(str), ])
geos_file_1: str = field([validators.of_instance(str), ])
geos_file_2: str = field([validators.of_instance(str), ])
geos_land: str = field([validators.of_instance(str), ])
geos_ocean: str = field([validators.of_instance(str), ])
geos_constants: str = field([validators.of_instance(str), ])
dims: tuple = field(default=('number_of_lines', 'number_of_pixels'),
validator=[validators.instance_of(tuple), ])
in_data = xr.open_dataset(geo_filename, group='geolocation_data')
data = xr.open_dataset(l1b_filename, group='observation_data', decode_cf=False)
@define(slots=True, kw_only=True)
class ReadData(CollectInputs):
"""Class that reads the L1b/geolocation data from VIIRS. Inherits file names from CollectInputs
for band in list(data.variables):
if 'reflectance' in data[band].long_name:
if hasattr(data[band], 'VCST_scale_factor'):
scale_factor = data[band].VCST_scale_factor * data[band].bias_correction
Parameters
----------
satellite: str
satellite name.
sensor: str
sensor name
"""
satellite: str = field(validator=[validators.instance_of(str),
validators.in_(['snpp'])])
sensor: str = field(validator=[validators.instance_of(str),
validators.in_(['viirs'])])
logging.debug('Instance of ReadData created')
def read_viirs_geo(self) -> xr.Dataset:
"""Read VIIRS geolocation data and generate additional angles
Parameters
----------
self.file_name_geo: str
file name of the geolocation file
Returns
-------
geo_data xarray.Dataset
dataset containing all geolocation data
"""
logging.debug(f'Reading {self.file_name_geo}')
geo_data = xr.open_dataset(self.file_name_geo, group='geolocation_data')
relazi = self.relative_azimuth_angle(geo_data.sensor_azimuth.values, geo_data.solar_azimuth.values)
sunglint = self.sun_glint_angle(geo_data.sensor_zenith.values, geo_data.solar_zenith.values, relazi)
scatt_angle = self.scattering_angle(geo_data.solar_zenith.values, geo_data.sensor_zenith.values,
relazi)
geo_data['relative_azimuth'] = (self.dims, relazi)
geo_data['sunglint_angle'] = (self.dims, sunglint)
geo_data['scattering_angle'] = (self.dims, scatt_angle)
logging.debug('Geolocation file read correctly')
return geo_data
def read_viirs_l1b(self,
solar_zenith: np.ndarray(float)) -> xr.Dataset:
"""Read VIIRS L1b data
Parameters
----------
self.file_name_l1b: str
file name of the L1b file
solar_zenith: np.ndarray
solar zenith angle derived from the geolocation file
"""
logging.debug(f'Reading {self.file_name_l1b}')
l1b_data = xr.open_dataset(self.file_name_l1b, group='observation_data', decode_cf=False)
rad_data = xr.Dataset()
for band in list(l1b_data.variables):
if 'reflectance' in l1b_data[band].long_name:
if hasattr(l1b_data[band], 'VCST_scale_factor'):
scale_factor = l1b_data[band].VCST_scale_factor * l1b_data[band].bias_correction
else:
scale_factor = l1b_data[band].scale_factor
rad_data[band] = (self.dims, l1b_data[band].values * scale_factor / np.cos(solar_zenith*_DTR))
elif 'radiance' in l1b_data[band].long_name:
bt_lut = f'{band}_brightness_temperature_lut'
rad_data[band] = (self.dims,
l1b_data[bt_lut].values[l1b_data[band].values])
else:
scale_factor = data[band].scale_factor
in_data[band] = (('number_of_lines', 'number_of_pixels'),
data[band].values * scale_factor / np.cos(in_data.solar_zenith.values*_DTR))
pass
elif 'radiance' in data[band].long_name:
in_data[band] = (('number_of_lines', 'number_of_pixels'),
data[f'{band}_brightness_temperature_lut'].values[data[band].values])
else:
pass
logging.debug('L1b file read correctly')
relazi = relative_azimuth_angle(in_data.sensor_azimuth.values, in_data.solar_azimuth.values)
sunglint = sun_glint_angle(in_data.sensor_zenith.values, in_data.solar_zenith.values, relazi)
scatt_angle = scattering_angle(in_data.solar_zenith.values, in_data.sensor_zenith.values, relazi)
return rad_data
in_data['relative_azimuth'] = (('number_of_lines', 'number_of_pixels'), relazi)
in_data['sunglint_angle'] = (('number_of_lines', 'number_of_pixels'), sunglint)
in_data['scattering_angle'] = (('number_of_lines', 'number_of_pixels'), scatt_angle)
def relative_azimuth_angle(sensor_azimuth: np.ndarray(float),
solar_azimuth: np.ndarray(float)) -> np.ndarray(float):
"""Computation of the relative azimuth angle
return in_data
Parameters
----------
sensor_azimuth: np.ndarray
sensor azimuth angle from the geolocation file
solar_azimuth: np.ndarray
solar azimuth angle from the geolocation file
Returns
-------
relative_azimuth: np.ndarray
"""
rel_azimuth = np.abs(180. - np.abs(sensor_azimuth - solar_azimuth))
def relative_azimuth_angle(sensor_azimuth: np.ndarray,
solar_azimuth: np.ndarray) -> np.ndarray:
"""Computation of the relative azimuth angle
logging.debug('Relative azimuth calculated successfully.')
Parameters
----------
sensor_azimuth: np.ndarray
sensor azimuth angle from the geolocation file
solar_azimuth: np.ndarray
solar azimuth angle from the geolocation file
return rel_azimuth
Returns
-------
relative_azimuth: np.ndarray
"""
rel_azimuth = np.abs(180. - np.abs(sensor_azimuth - solar_azimuth))
return rel_azimuth
def sun_glint_angle(sensor_zenith: np.ndarray(float),
solar_zenith: np.ndarray(float),
rel_azimuth: np.ndarray(float)) -> np.ndarray(float):
"""Computation of the sun glint angle
Parameters
----------
sensor_zenith: np.ndarray
sensor zenith angle from the geolocation file
solar_zenith: np.ndarray
solar zenith angle from the geolocation file
relative_azimuth: np.ndarray
relative azimuth computed from function relative_azimuth_angle()
def sun_glint_angle(sensor_zenith: np.ndarray,
solar_zenith: np.ndarray,
rel_azimuth: np.ndarray) -> np.ndarray:
"""Computation of the sun glint angle
Returns
-------
sunglint_angle: np.ndarray
"""
cossna = (np.sin(sensor_zenith*_DTR) * np.sin(solar_zenith*_DTR) * np.cos(rel_azimuth*_DTR) +
np.cos(sensor_zenith*_DTR) * np.cos(solar_zenith*_DTR))
cossna[cossna > 1] = 1
sunglint_angle = np.arccos(cossna) * _RTD
Parameters
----------
sensor_zenith: np.ndarray
sensor zenith angle from the geolocation file
solar_zenith: np.ndarray
solar zenith angle from the geolocation file
relative_azimuth: np.ndarray
relative azimuth computed from function relative_azimuth_angle()
logging.debug('Sunglint generated')
Returns
-------
sunglint_angle: np.ndarray
"""
cossna = (np.sin(sensor_zenith*_DTR) * np.sin(solar_zenith*_DTR) * np.cos(rel_azimuth*_DTR) +
np.cos(sensor_zenith*_DTR) * np.cos(solar_zenith*_DTR))
cossna[cossna > 1] = 1
sunglint_angle = np.arccos(cossna) * _RTD
return sunglint_angle
return sunglint_angle
def scattering_angle(solar_zenith: np.ndarray(float),
sensor_zenith: np.ndarray(float),
relative_azimuth: np.ndarray(float)) -> np.ndarray(float):
"""Computation of the scattering angle
Parameters
----------
solar_zenith: np.ndarray
solar zenith angle from the geolocation file
sensor_zenith: np.ndarray
sensor zenith angle angle from the geolocation file
relative_azimuth: np.ndarray
relative azimuth computed from function relative_azimuth_angle()
def scattering_angle(solar_zenith: np.ndarray,
sensor_zenith: np.ndarray,
relative_azimuth: np.ndarray) -> np.ndarray:
"""Computation of the scattering angle
Returns
-------
scattering_angle: np.ndarray
"""
cos_scatt_angle = -1. * (np.cos(solar_zenith*_DTR) * np.cos(sensor_zenith*_DTR) -
np.sin(solar_zenith*_DTR) * np.sin(sensor_zenith*_DTR) *
np.cos(relative_azimuth*_DTR))
scatt_angle = np.arccos(cos_scatt_angle) * _RTD
logging.debug('Scattering angle calculated successfully')
return scatt_angle
@define(kw_only=True, slots=True)
class ReadAncillary(CollectInputs):
"""Class tha processes all the ancillary files and makes them available for the MVCM"
Parameters
----------
solar_zenith: np.ndarray
solar zenith angle from the geolocation file
sensor_zenith: np.ndarray
sensor zenith angle angle from the geolocation file
relative_azimuth: np.ndarray
relative azimuth computed from function relative_azimuth_angle()
Returns
-------
scattering_angle: np.ndarray
latitude: np.ndarray(float)
array of latitudes for the granule that is being processed
longitude: np.ndarray(float)
array of longitudes for the granule that is being processed
resolution: int
flag to switch between MOD (1) and IMG (2) resolution
"""
cos_scatt_angle = -1. * (np.cos(solar_zenith*_DTR) * np.cos(sensor_zenith*_DTR) -
np.sin(solar_zenith*_DTR) * np.sin(sensor_zenith*_DTR) *
np.cos(relative_azimuth*_DTR))
scatt_angle = np.arccos(cos_scatt_angle) * _RTD
return scatt_angle
latitude: np.ndarray(float) = field([validators.of_instance(np.ndarray), ])
longitude: np.ndarray(float) = field([validators.of_instance(np.ndarray), ])
resolution: int = field([validators.of_instance(int),
validators.in_([1, 2]), ])
vnp_time = '.'.join(os.path.basename(CollectInputs.file_name_l1b).split('.')[1:3])
start_time = dt.strftime(dt.strptime(vnp_time, 'A%Y%j.%H%M'), '%Y-%m-%d %H:%M:%S.000')
out_shape = latitude.shape
def get_sst(self):
sst = np.empty(self.out_shape, dtype=np.float32)
sst = anc.py.get_Reynolds_SST(self.latitude.ravel(), self.longitude.ravel(), self.resolution,
self.ancillary_dir, self.sst_file, sst)
logging.debug('SST read successfully')
return sst.reshape(self.out_shape)
def get_ndvi(self):
ndvi = np.empty(self.out_shape, dtype=np.float32)
ndvi = anc.py_get_NDVI_background(self.latitude.ravel(), self.longitude.ravel(), self.resolution,
self.ancillary_dir, self.ndvi_file, ndvi)
logging.debug('NDVI read successfully')
return ndvi.reshape(self.out_shape)
def get_eco(self):
pass
def get_geos(self):
geos_variables = ['tpw', 'snow_fraction', 'ice_fraction', 'ocean_fraction',
'land_ice_fraction', 'surface_temperature']
geos_data = {var: np.empty(self.out_shape, dtype=np.float32) for var in geos_variables}
geos_data = anc.py_get_GEOS(self.latitude.ravel(), self.longitude.ravel(), self.resolution,
self.start_time, self.ancillary_dir, self.geos_file_before,
self.geos_file_after, self.geos_land, self.geos_ocean,
self.geos_constants, geos_data)
for var in list(geos_variables):
geos_data[var] = geos_data[var].reshape(self.out_shape)
logging.debug('GEOS data read successfully')
return geos_data
def correct_reflectances():
......@@ -279,7 +423,7 @@ def get_data(satellite: str,
viirs_data['M128'] = (('number_of_lines', 'number_of_pixels'), np.zeros(viirs_data.M15.shape))
else:
viirs_data = read_data('viirs', f'{img02}', f'{img03}')
viirs_data = read_viirs_data('viirs', f'{img02}', f'{img03}')
viirs_data['M05'] = viirs_data.I01
viirs_data['M07'] = viirs_data.I02
......@@ -320,3 +464,22 @@ def get_data(satellite: str,
data.drop_vars(['latitude', 'longitude'])
return data
def get_data_new(satellite: str,
sensor: str,
file_names: Dict[str, str],
sunglint_angle: float,
hires: bool = False) -> xr.Dataset:
"""Reads satellite and ancillary data and packs everything in an xarray dataset for the main function
"""
mod02 = file_names['MOD02']
mod03 = file_names['MOD03']
viirs_data = read_viirs_data(sensor, f'{mod02}', f'{mod03}')
ancillary_data = read_ancillary_data(file_names, viirs_data.latitude.values, viirs_data.longitude.values)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment