From d29519e77327d40432b857f75842543467ae5f0c Mon Sep 17 00:00:00 2001
From: Paolo Veglio <paolo.veglio@ssec.wisc.edu>
Date: Tue, 31 Jan 2023 17:22:14 +0000
Subject: [PATCH] first draft of testing functions
---
read_data.py | 79 ++++++++++++++++++++------
tests/test_read_data.py | 120 ++++++++++++++++++++++++++++++++++++++++
2 files changed, 182 insertions(+), 17 deletions(-)
create mode 100644 tests/test_read_data.py
diff --git a/read_data.py b/read_data.py
index a2e64c5..415fc69 100644
--- a/read_data.py
+++ b/read_data.py
@@ -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
+"""
diff --git a/tests/test_read_data.py b/tests/test_read_data.py
new file mode 100644
index 0000000..e070170
--- /dev/null
+++ b/tests/test_read_data.py
@@ -0,0 +1,120 @@
+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
--
GitLab