diff --git a/ancillary.pyx b/ancillary.pyx
index 818e0a3edc62a0ebca2e7c930ca7bc434f7f537a..2693c0c92c0bd10320dbc96477ffa22ee69414dd 100644
--- a/ancillary.pyx
+++ b/ancillary.pyx
@@ -66,48 +66,28 @@ def py_get_Olson_eco(np.ndarray[float, ndim=1] lat,
def py_get_GEOS(np.ndarray[float, ndim=1] lat, np.ndarray[float, ndim=1] lon, char *startTime,
char *anc_dir, char *geos1, char *geos2, char *geos_lnd, char *geos_ocn, char *geos_cnst,
- tpw, snowfr, icefr, ocnfr, landicefr, sfct):
-
-
-# for v in geos_data:
-# if not geos_data[v].flags['C_CONTIGUOUS']:
-# geos_data[v] = np.ascontiguousarray(geos_data[v])
- if not tpw.flags['C_CONTIGUOUS']:
- tpw = np.ascontiguousarray(tpw)
- if not snowfr.flags['C_CONTIGUOUS']:
- snowfr = np.ascontiguousarray(snowfr)
- if not icefr.flags['C_CONTIGUOUS']:
- icefr = np.ascontiguousarray(icefr)
- if not ocnfr.flags['C_CONTIGUOUS']:
- ocnfr = np.ascontiguousarray(ocnfr)
- if not landicefr.flags['C_CONTIGUOUS']:
- landicefr = np.ascontiguousarray(landicefr)
- if not sfct.flags['C_CONTIGUOUS']:
- sfct = np.ascontiguousarray(sfct)
-
- cdef float[::1] tpw_mv = tpw # geos_data['tpw']
- cdef float[::1] snowfr_mv = snowfr # geos_data['snowfr']
- cdef float[::1] icefr_mv = icefr # geos_data['icefr']
- cdef float[::1] ocnfr_mv = ocnfr # geos_data['geos_ocnfr']
- cdef float[::1] landicefr_mv = landicefr # geos_data['landicefr']
- cdef float[::1] sfct_mv = sfct # geos_data['sfct']
+ geos_data):
+
+ for v in geos_data:
+ if not geos_data[v].flags['C_CONTIGUOUS']:
+ geos_data[v] = np.ascontiguousarray(geos_data[v])
+
+ cdef float[::1] tpw_mv = geos_data['tpw']
+ cdef float[::1] snowfr_mv = geos_data['snowfr']
+ cdef float[::1] icefr_mv = geos_data['icefr']
+ cdef float[::1] ocnfr_mv = geos_data['ocnfr']
+ cdef float[::1] landicefr_mv = geos_data['landicefr']
+ cdef float[::1] sfct_mv = geos_data['sfct']
get_GEOS(&lat[0], &lon[0], startTime, anc_dir, geos1, geos2, geos_lnd, geos_ocn, geos_cnst,
&tpw_mv[0], &snowfr_mv[0], &icefr_mv[0], &ocnfr_mv[0], &landicefr_mv[0], &sfct_mv[0])
-# geos_dict = {'tpw': geos_data['tpw'],
-# 'snowfr': geos_data['snowfr'],
-# 'icefr': geos_data['icefr'],
-# 'geos_ocnfr': geos_data['geos_ocnfr'],
-# 'landicefr': geos_data['landicefr'],
-# 'sfct': geos_data['sfct']
-# }
-# geos_dict = {'tpw': tpw,
-# 'snowfr': snowfr,
-# 'icefr': icefr,
-# 'ocnfr': ocnfr,
-# 'landicefr': landicefr,
-# 'sfct': sfct
-# }
- print(tpw[:10])
- return tpw
+ geos_dict = {'tpw': geos_data['tpw'],
+ 'snowfr': geos_data['snowfr'],
+ 'icefr': geos_data['icefr'],
+ 'ocnfr': geos_data['ocnfr'],
+ 'landicefr': geos_data['landicefr'],
+ 'sfct': geos_data['sfct']
+ }
+
+ return geos_dict
diff --git a/main.py b/main.py
index 00f9bb7321ff62b56189409a874461702942384a..beb40d1e3122d822ac687a5d4b4557945c33dcbe 100644
--- a/main.py
+++ b/main.py
@@ -5,8 +5,6 @@ import numpy as np
import read_data as rd
import tests
-import ancillary_data as anc
-
def main():
datapath = '/ships19/hercules/pveglio/neige_data/snpp_test_input'
@@ -14,61 +12,9 @@ def main():
fname_geo = 'VNP03MOD.A2014213.1548.001.2017301015705.uwssec.nc'
thresh_file = '/home/pveglio/mvcm_leo/thresholds/new_thresholds.mvcm.snpp.v1.0.0.yaml'
- anc_dir = '/ships19/hercules/pveglio/neige_data/snpp_test_input/ancillary/2014_08_01_213'
- sst_file = 'oisst.20140730'
- ndvi_file = 'NDVI.FM.c004.v2.0.WS.00-04.209.hdf'
- geos_cnst = 'GEOS.fp.asm.const_2d_asm_Nx.00000000_0000.V01.nc4'
- geos_lnd = 'GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20140801_1530.V01.nc4'
- geos_ocn = 'GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20140801_1530.V01.nc4'
- geos1 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20140801_1500.V01.nc4'
- geos2 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20140801_1800.V01.nc4'
-
- start_time = '2014-08-01 15:48:00.000'
viirs_data = rd.read_data('viirs', f'{datapath}/{fname_l1b}', f'{datapath}/{fname_geo}')
- sst = np.empty((3232*3200, ), dtype=np.float32)
- ndvi = np.empty((3232*3200, ), dtype=np.float32)
- eco = np.empty((3232*3200, ), dtype=np.ubyte)
-
- geos_data = {'tpw': np.empty((3232*3200, ), dtype=np.float32),
- 'snowfr': np.empty((3232*3200, ), dtype=np.float32),
- 'icefr': np.empty((3232*3200, ), dtype=np.float32),
- 'ocnfr': np.empty((3232*3200, ), dtype=np.float32),
- 'landicefr': np.empty((3232*3200, ), dtype=np.float32),
- 'sfct': np.empty((3232*3200, ), dtype=np.float32),
- }
-
- tpw = np.ones((3232*3200, ), dtype=np.float32)
- snowfr = np.empty((3232*3200, ), dtype=np.float32)
- icefr = np.empty((3232*3200, ), dtype=np.float32)
- ocnfr = np.empty((3232*3200, ), dtype=np.float32)
- landicefr = np.empty((3232*3200, ), dtype=np.float32)
- sfct = np.empty((3232*3200, ), dtype=np.float32)
-
- sst = anc.py_get_Reynolds_SST(viirs_data.latitude.values.reshape((3232*3200, )),
- viirs_data.longitude.values.reshape((3232*3200, )),
- anc_dir, sst_file, sst)
-
- ndvi = anc.py_get_NDVI_background(viirs_data.latitude.values.reshape((3232*3200, )),
- viirs_data.longitude.values.reshape((3232*3200, )),
- anc_dir, ndvi_file, ndvi)
-
- eco = anc.py_get_Olson_eco(viirs_data.latitude.values.reshape((3232*3200, )),
- viirs_data.longitude.values.reshape((3232*3200, )),
- anc_dir, eco)
-
- tpw = anc.py_get_GEOS(viirs_data.latitude.values.reshape((3232*3200, )),
- viirs_data.longitude.values.reshape((3232*3200, )),
- start_time, anc_dir, geos1, geos2, geos_lnd, geos_ocn, geos_cnst,
- tpw, snowfr, icefr, ocnfr, landicefr, sfct)
-
- print(tpw[:10])
- #print(geos_data['tpw'][:10])
- #print(geos_data['snowfr'][:10])
- #print(geos_data['icefr'][:10])
- #print(geos_data['ocnfr'][:10])
- #print(geos_data['landicefr'][:10])
- #print(geos_data['sfct'][:10])
+ viirs_data = rd.read_ancillary_data(viirs_data)
with open(thresh_file) as f:
text = f.read()
diff --git a/read_data.py b/read_data.py
index 06ac0306ddbf96de760d1b4fb4ac1b4ca0069e27..62d1e928a347da839e27bee8bb4dc168cce4e7ab 100644
--- a/read_data.py
+++ b/read_data.py
@@ -3,6 +3,8 @@
import xarray as xr
import numpy as np
+import ancillary_data as anc
+
_DTR = np.pi/180.
_RTD = 180./np.pi
@@ -65,3 +67,58 @@ def scattering_angle(solar_zenith, sensor_zenith, relative_azimuth):
def correct_reflectances():
pass
+
+
+def read_ancillary_data(viirs_data):
+
+ # Ancillary files temporarily defined here. Eventually we will find a better way to pass these
+ start_time = '2014-08-01 15:48:00.000'
+ anc_dir = '/ships19/hercules/pveglio/neige_data/snpp_test_input/ancillary/2014_08_01_213'
+ sst_file = 'oisst.20140730'
+ ndvi_file = 'NDVI.FM.c004.v2.0.WS.00-04.209.hdf'
+ geos_cnst = 'GEOS.fp.asm.const_2d_asm_Nx.00000000_0000.V01.nc4'
+ geos_lnd = 'GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20140801_1530.V01.nc4'
+ geos_ocn = 'GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20140801_1530.V01.nc4'
+ geos1 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20140801_1500.V01.nc4'
+ geos2 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20140801_1800.V01.nc4'
+
+ sst = np.empty((3232*3200, ), dtype=np.float32)
+ ndvi = np.empty((3232*3200, ), dtype=np.float32)
+ eco = np.empty((3232*3200, ), dtype=np.ubyte)
+
+ geos_data = {'tpw': np.empty((3232*3200, ), dtype=np.float32),
+ 'snowfr': np.empty((3232*3200, ), dtype=np.float32),
+ 'icefr': np.empty((3232*3200, ), dtype=np.float32),
+ 'ocnfr': np.empty((3232*3200, ), dtype=np.float32),
+ 'landicefr': np.empty((3232*3200, ), dtype=np.float32),
+ 'sfct': np.empty((3232*3200, ), dtype=np.float32),
+ }
+
+ sst = anc.py_get_Reynolds_SST(viirs_data.latitude.values.reshape((3232*3200, )),
+ viirs_data.longitude.values.reshape((3232*3200, )),
+ anc_dir, sst_file, sst)
+
+ ndvi = anc.py_get_NDVI_background(viirs_data.latitude.values.reshape((3232*3200, )),
+ viirs_data.longitude.values.reshape((3232*3200, )),
+ anc_dir, ndvi_file, ndvi)
+
+ eco = anc.py_get_Olson_eco(viirs_data.latitude.values.reshape((3232*3200, )),
+ viirs_data.longitude.values.reshape((3232*3200, )),
+ anc_dir, eco)
+
+ geos_data = anc.py_get_GEOS(viirs_data.latitude.values.reshape((3232*3200, )),
+ viirs_data.longitude.values.reshape((3232*3200, )),
+ start_time, anc_dir, geos1, geos2, geos_lnd, geos_ocn, geos_cnst, geos_data)
+
+ data = viirs_data
+ data['sst'] = (['number_of_lines', 'number_of_pixels'], sst)
+ data['ndvi'] = (['number_of_lines', 'number_of_pixels'], ndvi)
+ data['eco'] = (['number_of_lines', 'number_of_pixels'], eco)
+ data['geos_tpw'] = (['number_of_lines', 'number_of_pixels'], geos_data['tpw'])
+ data['geos_snowfr'] = (['number_of_lines', 'number_of_pixels'], geos_data['snowfr'])
+ data['geos_icefr'] = (['number_of_lines', 'number_of_pixels'], geos_data['icefr'])
+ data['geos_ocnfr'] = (['number_of_lines', 'number_of_pixels'], geos_data['ocnfr'])
+ data['geos_landicefr'] = (['number_of_lines', 'number_of_pixels'], geos_data['landicefr'])
+ data['geos_sfct'] = (['number_of_lines', 'number_of_pixels'], geos_data['sfct'])
+
+ return data
diff --git a/src/get_GEOS.c b/src/get_GEOS.c
index a379f1a0245184ab8ab57de722f3b8bbf5ee96c0..b8a87cc9d81ca209877ce86bb01d0f98a151335c 100644
--- a/src/get_GEOS.c
+++ b/src/get_GEOS.c
@@ -182,19 +182,21 @@ void get_GEOS(float *lat, float *lon, char *granule_start_time, char *anc_dir, c
/******************************************************************************/
- tpw = grn_anc.tpw;
- snowfr = grn_anc.snowfr;
- icefr = grn_anc.icefr;
- geos_ocnfr = grn_anc.geos_ocnfr;
- landicefr = grn_anc.landicefr;
- sfct = grn_anc.sfct;
-
int i;
-
- for (i=0; i<10; i++) {
- printf("i: %d\ttpw: %f\tsnowfr: %f\ticefr: %f\tgeos_ocnfr: %f\tlandicefr: %f\tsfct: %f\n",
- i, tpw[i], snowfr[i], icefr[i], geos_ocnfr[i], landicefr[i], sfct[i]);
+ for (i=0; i<3232*3200; i++) {
+ tpw[i] = grn_anc.tpw[i];
+ snowfr[i] = grn_anc.snowfr[i];
+ icefr[i] = grn_anc.icefr[i];
+ geos_ocnfr[i] = grn_anc.geos_ocnfr[i];
+ landicefr[i] = grn_anc.landicefr[i];
+ sfct[i] = grn_anc.sfct[i];
}
+ //int i;
+ //
+ //for (i=0; i<10; i++) {
+ // printf("i: %d\ttpw: %f\tsnowfr: %f\ticefr: %f\tgeos_ocnfr: %f\tlandicefr: %f\tsfct: %f\n",
+ // i, tpw[i], snowfr[i], icefr[i], geos_ocnfr[i], landicefr[i], sfct[i]);
+ //}
// return ( return_code);
}
diff --git a/tests.py b/tests.py
index 04667c130e8c62df0401448f5f6f36b08e1e8a8d..044b2066e22c4a07ff19b08fdd5ff060acfdd5a1 100644
--- a/tests.py
+++ b/tests.py
@@ -87,14 +87,16 @@ def vir_refl_test(rad, threshold, viirs_data):
cosvza = np.cos(vza*dtr)
coeffs = utils.get_b1_thresholds()
+ coeffs[:, :3] = coeffs[:, :3] * threshold['b1_bias_adj']
+
# this quantity is the return of get_b1_thresholds() in the C code
- # it's defined here to keep a consistent logic with the original source
- irtn = 1
+ # it's defined here to keep a consistent logic with the original source, for now
+ irtn = 0
if irtn != 0:
coeffs = thr
- coeffs[:3] = coeffs[:3] * 1/np.power(cosvza, vzcpow)
+ coeffs[:, :3] = coeffs[:, :3] * 1/np.power(cosvza, vzcpow)
confidence = conf.conf_test(rad, coeffs)
diff --git a/utils.py b/utils.py
index 21fb39b1551a19dbd2c85d45c35814295cd3e44a..2cf07c3a5c13abf364b2d7e70c46925e81801cca 100644
--- a/utils.py
+++ b/utils.py
@@ -295,3 +295,25 @@ def get_b1_threshold(thresholds, ndvi_thresholds):
thr[:, 3] = 2.0
return thr
+
+
+def find_scene(data):
+ scene = {'ocean_day': 0,
+ 'ocean_night': 0,
+ 'land_day': 0,
+ 'land_night': 0,
+ 'snow_day': 0,
+ 'coast_day': 0,
+ 'desert_day': 0,
+ 'antarctic_day': 0,
+ 'polar_day_snow': 0,
+ 'polar_day_desert': 0,
+ 'polar_day_desert_coast': 0,
+ 'polar_day_coast': 0,
+ 'polar_day_land': 0,
+ 'polar_night_snow': 0,
+ 'polar_night_land': 0,
+ 'polar_ocean_night': 0,
+ }
+
+ return scene