diff --git a/main.py b/main.py index db6d6e21dfcd39cfdcd50132ac327952aadb7f07..08ca5a4a02046427c275b8e6c7bb9f866084e269 100644 --- a/main.py +++ b/main.py @@ -38,12 +38,13 @@ _eco_file = 'goge1_2_img.v1' # #################################################################### # -def main(*, +def main(satellite: str = 'snpp', + sensor: str = 'viirs', data_path: str = _datapath, mod02: str = _fname_mod02, mod03: str = _fname_mod03, - img02: str = _fname_img02, - img03: str = _fname_img03, + img02: str = None, + img03: str = None, threshold_file: str = _threshold_file, geos_atm_1: str = _geos_atm_1, geos_atm_2: str = _geos_atm_2, @@ -56,15 +57,19 @@ def main(*, Parameters ---------- + satellite: str + satellite name, not case-sensitive. Available options: [snpp, ] + sensor: str + sensor name, not case-sensitive. Available options: [viirs, ] data_path: str path where the data is stored mod02: str VIIRS MOD02 file name mod03: str VIIRS MOD03 file name - img02: str + img02: [optional] str VIIRS IMG02 file name - img03: str + img03: [optional] str VIIRS IMG03 file name threshold_file: str thresholds file name @@ -108,7 +113,7 @@ def main(*, sunglint_angle = thresholds['Sun_Glint']['bounds'][3] - viirs_data = rd.get_data(file_names, sunglint_angle) + viirs_data = rd.get_data(satellite, sensor, file_names, sunglint_angle) # Initialize the confidence arrays for the various test groups cmin_G1 = np.ones(viirs_data.M01.shape) diff --git a/read_data.py b/read_data.py index 2a12b609fc58f479b05823c926f9a3d62f97283e..9f27a6d12b25b49fb5bf63529d44383194946184 100644 --- a/read_data.py +++ b/read_data.py @@ -5,6 +5,7 @@ import ancillary_data as anc import scene as scn import os +import logging from datetime import datetime as dt from typing import Dict @@ -12,30 +13,47 @@ _DTR = np.pi/180. _RTD = 180./np.pi _bad_data = -999.0 +logging.basicConfig(level=logging.INFO, format='%(name)s - %(levelname)s - %(message)s') +# logging.basicConfig(level=logging.INFO, filename='logfile.log', 'filemode='w', +# format='%(name)s %(levelname)s %(message)s') -def read_data(sensor: str, - l1b_filename: str, - geo_filename: str) -> xr.Dataset: + +def read_viirs_data(l1b_filename: str, + geo_filename: str) -> xr.Dataset: + """Read VIIRS MOD or IMG data + + 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) + """ in_data = xr.open_dataset(geo_filename, group='geolocation_data') data = xr.open_dataset(l1b_filename, group='observation_data', decode_cf=False) - if sensor.lower() == 'viirs': - 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 - 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)) - - 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]) + 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 else: - pass + 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)) + + 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 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) @@ -50,6 +68,19 @@ def read_data(sensor: str, def relative_azimuth_angle(sensor_azimuth: np.ndarray, solar_azimuth: np.ndarray) -> np.ndarray: + """Computation of the relative azimuth angle + + 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)) return rel_azimuth @@ -57,6 +88,21 @@ def relative_azimuth_angle(sensor_azimuth: np.ndarray, def sun_glint_angle(sensor_zenith: np.ndarray, solar_zenith: np.ndarray, rel_azimuth: np.ndarray) -> np.ndarray: + """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() + + 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 @@ -67,6 +113,21 @@ def sun_glint_angle(sensor_zenith: np.ndarray, def scattering_angle(solar_zenith: np.ndarray, sensor_zenith: np.ndarray, relative_azimuth: np.ndarray) -> np.ndarray: + """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() + + 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)) @@ -79,9 +140,22 @@ def correct_reflectances(): def read_ancillary_data(file_names: str, - viirs_data: xr.Dataset, + latitude: np.ndarray, + longitude: np.ndarray, resolution: int = 1) -> xr.Dataset: - + """Read ancillary data using C functions from original MVCM + + Parameters + ---------- + file_names: str + latitude: np.ndarray + longitude: np.ndarray + resolution: int + + Returns + ------- + data: xarray.Dataset + """ # Ancillary files temporarily defined here. Eventually we will find a better way to pass these anc_dir = file_names['ANC_DIR'] sst_file = file_names['SST'] @@ -93,58 +167,51 @@ def read_ancillary_data(file_names: str, geos2 = file_names['GEOS_atm_2'] vnptime = '.'.join(os.path.basename(file_names['MOD02']).split('.')[1:3]) start_time = dt.strftime(dt.strptime(vnptime, 'A%Y%j.%H%M'), '%Y-%m-%d %H:%M:%S.000') - # start_time = '2022-06-22 14:54:00.000' - dim1 = viirs_data.latitude.shape[0] - dim2 = viirs_data.latitude.shape[1] + out_shape = latitude.shape - sst = np.empty((dim1*dim2, ), dtype=np.float32) - ndvi = np.empty((dim1*dim2, ), dtype=np.float32) - eco = np.empty((dim1*dim2, ), dtype=np.ubyte) + sst = np.empty((np.prod(out_shape), ), dtype=np.float32) + ndvi = np.empty((np.prod(out_shape), ), dtype=np.float32) + eco = np.empty((np.prod(out_shape), ), dtype=np.ubyte) - geos_data = {'tpw': np.empty((dim1*dim2, ), dtype=np.float32), - 'snowfr': np.empty((dim1*dim2, ), dtype=np.float32), - 'icefr': np.empty((dim1*dim2, ), dtype=np.float32), - 'ocnfr': np.empty((dim1*dim2, ), dtype=np.float32), - 'landicefr': np.empty((dim1*dim2, ), dtype=np.float32), - 'sfct': np.empty((dim1*dim2, ), dtype=np.float32), + geos_data = {'tpw': np.empty((np.prod(out_shape), ), dtype=np.float32), + 'snowfr': np.empty((np.prod(out_shape), ), dtype=np.float32), + 'icefr': np.empty((np.prod(out_shape), ), dtype=np.float32), + 'ocnfr': np.empty((np.prod(out_shape), ), dtype=np.float32), + 'landicefr': np.empty((np.prod(out_shape), ), dtype=np.float32), + 'sfct': np.empty((np.prod(out_shape), ), dtype=np.float32), } - sst = anc.py_get_Reynolds_SST(viirs_data.latitude.values.reshape((dim1*dim2, )), - viirs_data.longitude.values.reshape((dim1*dim2, )), - resolution, anc_dir, sst_file, sst) - print('SST read successfully') + sst = anc.py_get_Reynolds_SST(latitude.ravel(), longitude.ravel(), resolution, anc_dir, sst_file, sst) + logging.info('SST read successfully') - ndvi = anc.py_get_NDVI_background(viirs_data.latitude.values.reshape((dim1*dim2, )), - viirs_data.longitude.values.reshape((dim1*dim2, )), + ndvi = anc.py_get_NDVI_background(latitude.ravel(), longitude.ravel(), resolution, anc_dir, ndvi_file, ndvi) - print('NDVI read successfully') - - eco = anc.py_get_Olson_eco(viirs_data.latitude.values.reshape((dim1*dim2, )), - viirs_data.longitude.values.reshape((dim1*dim2, )), - resolution, anc_dir, eco) - print('Olson eco read successfully') - - geos_data = anc.py_get_GEOS(viirs_data.latitude.values.reshape((dim1*dim2, )), - viirs_data.longitude.values.reshape((dim1*dim2, )), - resolution, start_time, anc_dir, geos1, geos2, geos_lnd, geos_ocn, geos_cnst, geos_data) - print('GEOS5 data read successfully') - - data = viirs_data - data['sst'] = (('number_of_lines', 'number_of_pixels'), sst.reshape((dim1, dim2))) - data['ndvi'] = (('number_of_lines', 'number_of_pixels'), ndvi.reshape((dim1, dim2))) - data['eco'] = (('number_of_lines', 'number_of_pixels'), eco.reshape((dim1, dim2))) - data['geos_tpw'] = (('number_of_lines', 'number_of_pixels'), geos_data['tpw'].reshape((dim1, dim2))) - data['geos_snowfr'] = (('number_of_lines', 'number_of_pixels'), geos_data['snowfr'].reshape((dim1, dim2))) - data['geos_icefr'] = (('number_of_lines', 'number_of_pixels'), geos_data['icefr'].reshape((dim1, dim2))) - data['geos_ocnfr'] = (('number_of_lines', 'number_of_pixels'), geos_data['ocnfr'].reshape((dim1, dim2))) - data['geos_landicefr'] = (('number_of_lines', 'number_of_pixels'), geos_data['landicefr'].reshape((dim1, dim2))) - data['geos_sfct'] = (('number_of_lines', 'number_of_pixels'), geos_data['sfct'].reshape((dim1, dim2))) + logging.info('NDVI read successfully') + + eco = anc.py_get_Olson_eco(latitude.ravel(), longitude.ravel(), resolution, anc_dir, eco) + logging.info('Olson eco read successfully') + + geos = anc.py_get_GEOS(latitude.ravel(), longitude.ravel(), resolution, start_time, anc_dir, + geos1, geos2, geos_lnd, geos_ocn, geos_cnst, geos_data) + logging.info('GEOS5 data read successfully') + + ancillary = {'sst': sst.reshape(out_shape), + 'ndvi': ndvi.reshape(out_shape), + 'eco': eco.reshape(out_shape) + } + for var in list(geos): + ancillary[f'geos_{var}'] = geos[var].reshape(out_shape) + + dims = ('number_of_lines', 'number_of_pixels') + data = xr.Dataset.from_dict({var: {'dims': dims, 'data': ancillary[var]} for var in list(ancillary)}) return data -def get_data(file_names: Dict[str, str], +def get_data(satellite: str, + sensor: str, + file_names: Dict[str, str], sunglint_angle: float, hires: bool = False) -> xr.Dataset: @@ -156,7 +223,7 @@ def get_data(file_names: Dict[str, str], img03 = file_names['IMG03'] if hires is False: - viirs_data = read_data('viirs', f'{mod02}', f'{mod03}') + viirs_data = read_viirs_data(sensor, f'{mod02}', f'{mod03}') viirs_data = read_ancillary_data(file_names, viirs_data) if (('M05' in viirs_data) and ('M07' in viirs_data)): diff --git a/setup.py b/setup.py index 192f920fcd412fe3cfd636130a63a61790d64452..18f61ecb8153d4e8c74898ebc43eaafcf3426252 100644 --- a/setup.py +++ b/setup.py @@ -26,14 +26,14 @@ sourcefiles = ['src/get_Reynolds_SST.c', # SST ] include_dirs = ['include', - '/opt/hdf4/4.2.9-gcc-4.9.2/include', - '/opt/hdfeos2/2.19-gcc-4.9.2/include', - '/opt/netcdf4/4.3.3-gcc-4.9.2/include', + '/opt/hdf4/4.2.14-gcc-8.3/include', + '/opt/hdfeos2/2.20-gcc-8.3/include', + '/opt/netcdf4/4.7.0-gcc-8.3/include', numpy.get_include(), ] -library_dirs = ['/opt/hdf4/4.2.9-gcc-4.9.2/lib', - '/opt/hdfeos2/2.19-gcc-4.9.2/lib', - '/opt/netcdf4/4.3.3-gcc-4.9.2/lib', +library_dirs = ['/opt/hdf4/4.2.14-gcc-8.3/lib', + '/opt/hdfeos2/2.20-gcc-8.3/lib', + '/opt/netcdf4/4.7.0-gcc-8.3/lib', ] extensions = [Extension('ancillary_data', sourcefiles,