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

first draft of testing functions

parent 3e6611a0
No related branches found
No related tags found
No related merge requests found
......@@ -3,7 +3,6 @@ import numpy as np
import numpy.typing as npt
import ancillary_data as anc
import scene as scn
import os
import logging
......@@ -14,6 +13,7 @@ from attrs import define, field, validators, Factory
_DTR = np.pi/180.
_RTD = 180./np.pi
_bad_data = -999.0
_datapath = '/ships19/hercules/pveglio/mvcm_viirs_hires'
logging.basicConfig(level=logging.INFO, format='%(name)s - %(levelname)s - %(message)s')
# logging.basicConfig(level=logging.INFO, filename='logfile.log', 'filemode='w',
......@@ -50,16 +50,28 @@ class CollectInputs(object):
dims [optional]: str
name of the dimensions for the arrays in the xarray.Dataset output
"""
file_name_geo: str = field(default='', validator=[validators.instance_of(str), ])
file_name_l1b: str = field(default='', validator=[validators.instance_of(str), ])
ancillary_dir: str = field(default='', validator=[validators.instance_of(str), ])
sst_file: str = field(default='', validator=[validators.instance_of(str), ])
ndvi_file: str = field(default='', validator=[validators.instance_of(str), ])
geos_file_before: str = field(default='', validator=[validators.instance_of(str), ])
geos_file_after: str = field(default='', validator=[validators.instance_of(str), ])
geos_land: str = field(default='', validator=[validators.instance_of(str), ])
geos_ocean: str = field(default='', validator=[validators.instance_of(str), ])
geos_constants: str = field(default='', validator=[validators.instance_of(str), ])
file_name_geo: str = field(default=f'{_datapath}/VNP03MOD.A2022173.1312.001.2022174012746.uwssec.nc',
validator=[validators.instance_of(str), ])
file_name_l1b: str = field(default=f'{_datapath}/VNP02MOD.A2022173.1312.001.2022174011547.uwssec.nc',
validator=[validators.instance_of(str), ])
ancillary_dir: str = field(default=f'{_datapath}/ancillary',
validator=[validators.instance_of(str), ])
sst_file: str = field(default='oisst.20220622',
validator=[validators.instance_of(str), ])
eco_file: str = field(default='goge1_2_img.v1',
validator=[validators.instance_of(str), ])
ndvi_file: str = field(default='NDVI.FM.c004.v2.0.WS.00-04-177.hdf',
validator=[validators.instance_of(str), ])
geos_file_1: str = field(default='GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1200.V01.nc4',
validator=[validators.instance_of(str), ])
geos_file_2: str = field(default='GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1500.V01.nc4',
validator=[validators.instance_of(str), ])
geos_land: str = field(default='GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20220622_1330.V01.nc4',
validator=[validators.instance_of(str), ])
geos_ocean: str = field(default='GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20220622_1330.V01.nc4',
validator=[validators.instance_of(str), ])
geos_constants: str = field(default='GEOS.fp.asm.const_2d_asm_Nx.00000000_0000.V01.nc4',
validator=[validators.instance_of(str), ])
dims: tuple = field(default=('number_of_lines', 'number_of_pixels'),
validator=[validators.instance_of(tuple), ])
......@@ -335,10 +347,12 @@ class ReadAncillary(CollectInputs):
sst: npt.NDArray[float]
array containing the Reynolds SST interpolated at the sensor's resolution
"""
if not os.path.isfile(os.path.join(self.ancillary_dir, self.sst_file)):
logging.error('SST file not found')
sst = np.empty(self.out_shape, dtype=np.float32).ravel()
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')
logging.debug('SST file read successfully')
return sst.reshape(self.out_shape)
def get_ndvi(self) -> npt.NDArray[float]:
......@@ -352,14 +366,30 @@ class ReadAncillary(CollectInputs):
ndvi: npt.NDArray[float]
NDVI interpolated at the sensor's resolution
"""
if not os.path.isfile(os.path.join(self.ancillary_dir, self.ndvi_file)):
logging.error('NDVI file not found')
ndvi = np.empty(self.out_shape, dtype=np.float32).ravel()
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')
logging.debug('NDVI file read successfully')
return ndvi.reshape(self.out_shape)
def get_eco(self):
pass
def get_eco(self) -> npt.NDArray[float]:
"""Read ECO file
Parameters
----------
Returns
-------
eco: npt.NDArray[float]
Olson ecosystem type interpolated at the sensor's resolution
"""
eco = np.empty(self.out_shape, dtype=np.ubyte).ravel()
eco = anc.py_get_Olson_eco(self.latitude.ravel(), self.longitude.ravel(), self.resolution,
self.ancillary_dir, eco)
logging.debug('Olson ecosystem file read successfully')
return eco.reshape(self.out_shape)
def get_geos(self) -> Dict:
"""Read GEOS-5 data and interpolate the fields to the sensor resolution
......@@ -372,13 +402,24 @@ class ReadAncillary(CollectInputs):
geos_data: Dict
dictionary containing all quantities required by MVCM (see geos_variables here below)
"""
if not os.path.isfile(os.path.join(self.ancillary, self.geos_file_1)):
logging.error('GEOS-5 file 1 not found')
if not os.path.isfile(os.path.join(self.ancillary, self.geos_file_2)):
logging.error('GEOS-5 file 2 not found')
if not os.path.isfile(os.path.join(self.ancillary, self.geos_land_file)):
logging.error('GEOS-5 land file not found')
if not os.path.isfile(os.path.join(self.ancillary, self.geos_ocean_file)):
logging.error('GEOS-5 ocean file not found')
if not os.path.isfile(os.path.join(self.ancillary, self.geos_constants_file)):
logging.error('GEOS-5 constants file not found')
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).ravel() for var in geos_variables}
geos_data = anc.py_get_GEOS(self.latitude.ravel(), self.longitude.ravel(), self.resolution,
self.get_granule_time(), self.ancillary_dir, self.geos_file_before,
self.geos_file_after, self.geos_land, self.geos_ocean,
self.get_granule_time(), self.ancillary_dir, self.geos_file_1,
self.geos_file_2, self.geos_land, self.geos_ocean,
self.geos_constants, geos_data)
for var in list(geos_variables):
......@@ -400,6 +441,7 @@ class ReadAncillary(CollectInputs):
"""
ancillary_data = xr.Dataset()
ancillary_data['sst'] = (self.dims, self.get_sst())
ancillary_data['ecosystem'] = (self.dims, self.get_eco())
ancillary_data['ndvi'] = (self.dims, self.get_ndvi())
geos_tmp = self.get_geos()
......@@ -408,6 +450,8 @@ class ReadAncillary(CollectInputs):
return ancillary_data
"""
scene_flags = scn.find_scene(viirs, sunglint_angle)
scene = scn.scene_id(scene_flags)
......@@ -534,3 +578,4 @@ def get_data(satellite: str,
data.drop_vars(['latitude', 'longitude'])
return data
"""
import numpy as np
import xarray as xr
import os
import pytest
import read_data as rd
@pytest.fixture
def fixturepath():
return os.path.join(os.path.dirname(__file__), 'fixtures')
# @pytest.fixture()
# def l1b_file(fixturepath):
# return os.path.join(fixturepath, '')
# @pytest.fixture()
# def geo_file(fixturepath):
# return os.path.join(fixturepath, '')
# @pytest.fixture
# def ancillary_dir(fixturepath):
# return os.path.join(fixturepath, '')
@pytest.fixture
def sst_file():
return 'oisst.20220622'
@pytest.fixture
def ndvi_file():
return 'NDVI.FM.c004.v2.0.WS.00-04.177.hdf'
@pytest.fixture
def geos_file_1():
return 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1200.V01.nc4'
@pytest.fixture
def geos_file_2():
return 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1500.V01.nc4'
@pytest.fixture
def geos_land():
return 'GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20220622_1330.V01.nc4'
@pytest.fixture
def geos_ocean():
return 'GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20220622_1330.V01.nc4'
@pytest.fixture
def geos_constants():
return 'GEOS.fp.asm.const_2d_asm_Nx.00000000_0000.V01.nc4'
@pytest.fixture
def ref_file(fixturepath):
return os.path.join(fixturepath, 'ref_ancillary.nc')
def test_sst(fixturepath, sst_file, ref_file):
viirs = rd.ReadData(satellite='snpp',
sensor='viirs')
geo = viirs.read_viirs_geo()
ancillary = rd.ReadAncillary(latitude=geo.latitude.values,
longitude=geo.longitude.values,
resolution=1,
ancillary_dir=fixturepath,
sst_file=sst_file)
sst = ancillary.get_sst()
check_differences(ref_file, sst, 'sst')
def test_ndvi(fixturepath, ndvi_file, ref_file):
viirs = rd.ReadData(satellite='snpp',
sensor='viirs')
geo = viirs.read_viirs_geo()
ancillary = rd.ReadAncillary(latitude=geo.latitude.values,
longitude=geo.longitude.values,
resolution=1,
ancillary_dir=fixturepath,
ndvi_file=ndvi_file)
ndvi = ancillary.get_ndvi()
check_differences(ref_file, ndvi, 'ndvi')
def test_geos(fixturepath, geos_file_1, geos_file_2, geos_land, geos_ocean,
geos_constants, ref_file):
viirs = rd.ReadData(satellite='snpp',
sensor='viirs')
geo = viirs.read_viirs_geo()
ancillary = rd.ReadAncillary(latitude=geo.latitude.values,
longitude=geo.longitude.values,
resolution=1,
ancillary_dir=fixturepath,
geos_file_1=geos_file_1,
geos_file_2=geos_file_2,
geos_land=geos_land,
geos_ocean=geos_ocean,
geos_constants=geos_constants)
geos_data = ancillary.get_geos()
for var in list(geos_data.keys()):
check_differences(ref_file, geos_data[var], f'geos_{var}')
def check_differences(ref_file, test_data, var_name):
ref_data = xr.open_dataset(ref_file)
check = np.allclose(ref_data[var_name].values, test_data, equal_nan=True)
assert check
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