diff --git a/read_data.py b/read_data.py index 9f27a6d12b25b49fb5bf63529d44383194946184..f6ccb8f816b7d348b56801cf09a13037297fd085 100644 --- a/read_data.py +++ b/read_data.py @@ -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) + + + + +