diff --git a/ancillary.pyx b/ancillary.pyx
index 2693c0c92c0bd10320dbc96477ffa22ee69414dd..cde00422b7c3d449d79e366734f88fc1b1992164 100644
--- a/ancillary.pyx
+++ b/ancillary.pyx
@@ -6,6 +6,7 @@ cdef extern void get_NDVI_background(float *, float *, char *, char *, float *)
cdef extern void get_Olson_eco(float *, float *, char *, unsigned char *)
cdef extern void get_GEOS(float *, float *, char *, char *, char *, char *, char *, char *, char *,
float *, float *, float *, float *, float *, float *)
+cdef extern void snow_mask(char *, unsigned char)
import cython
from cython.view cimport array as cvarray
@@ -91,3 +92,9 @@ def py_get_GEOS(np.ndarray[float, ndim=1] lat, np.ndarray[float, ndim=1] lon, ch
}
return geos_dict
+
+
+def py_snow_mask(char *satname, unsigned char lsf):
+ # need to have for loop here to compute all the pixels since the function, as with everything else,
+ # is run per pixel.
+ pass
diff --git a/main.py b/main.py
index beb40d1e3122d822ac687a5d4b4557945c33dcbe..6c7c0ed7fd983859a53b08332ff833df3c1e2a08 100644
--- a/main.py
+++ b/main.py
@@ -5,18 +5,55 @@ import numpy as np
import read_data as rd
import tests
-
-def main():
- datapath = '/ships19/hercules/pveglio/neige_data/snpp_test_input'
- fname_l1b = 'VNP02MOD.A2014213.1548.001.2017301015346.uwssec.bowtie_restored_scaled.nc'
- fname_geo = 'VNP03MOD.A2014213.1548.001.2017301015705.uwssec.nc'
- thresh_file = '/home/pveglio/mvcm_leo/thresholds/new_thresholds.mvcm.snpp.v1.0.0.yaml'
-
- viirs_data = rd.read_data('viirs', f'{datapath}/{fname_l1b}', f'{datapath}/{fname_geo}')
-
- viirs_data = rd.read_ancillary_data(viirs_data)
-
- with open(thresh_file) as f:
+# #################################################################### #
+# TEST CASE
+# data:
+_datapath = '/ships19/hercules/pveglio/mvcm_viirs_hires'
+_fname_mod02 = 'VNP02MOD.A2022173.0336.001.2022173154147.uwssec.nc'
+_fname_mod03 = 'VNP03MOD.A2022173.0336.001.2022173155356.uwssec.nc'
+_fname_img02 = 'VNP02IMG.A2022173.0336.001.2022173154147.uwssec.nc'
+_fname_img03 = 'VNP03IMG.A2022173.0336.001.2022173155356.uwssec.nc'
+
+# thresholds:
+_threshold_file = '/home/pveglio/mvcm_leo/thresholds/new_thresholds.mvcm.snpp.v1.0.0.yaml'
+
+# ancillary files:
+_geos_atm_1 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_0300.V01.nc4'
+_geos_atm_2 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_0600.V01.nc4'
+_geos_land = 'GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20220622_0330.V01.nc4'
+_geos_ocean = 'GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20220622_0330.V01.nc4'
+_geos_constants = 'GEOS.fp.asm.const_2d_asm_Nx.00000000_0000.V01.nc4'
+_ndvi_file = 'NDVI.FM.c004.v2.0.WS.00-04.177.hdf'
+_sst_file = 'oisst.20220622'
+
+# #################################################################### #
+
+
+def main(*, data_path=_datapath, mod02=_fname_mod02, mod03=_fname_mod03,
+ img02=_fname_img02, img03=_fname_img03, threshold_file=_threshold_file,
+ geos_atm_1=_geos_atm_1, geos_atm_2=_geos_atm_2, geos_land=_geos_land,
+ geos_ocean=_geos_ocean, geos_constants=_geos_constants, ndvi_file=_ndvi_file, sst_file=_sst_file):
+
+ # datapath = '/ships19/hercules/pveglio/neige_data/snpp_test_input'
+ # fname_l1b = 'VNP02MOD.A2014213.1548.001.2017301015346.uwssec.bowtie_restored_scaled.nc'
+ # fname_geo = 'VNP03MOD.A2014213.1548.001.2017301015705.uwssec.nc'
+ # thresh_file = '/home/pveglio/mvcm_leo/thresholds/new_thresholds.mvcm.snpp.v1.0.0.yaml'
+
+ ancillary_file_names = {'GEOS_atm_1': f'{data_path}/ancillary/{geos_atm_1}',
+ 'GEOS_atm_2': f'{data_path}/ancillary/{geos_atm_2}',
+ 'GEOS_land': f'{data_path}/ancillary/{geos_land}',
+ 'GEOS_ocean': f'{data_path}/ancillary/{geos_ocean}',
+ 'GEOS_constants': f'{data_path}/ancillary/{geos_constants}',
+ 'NDVI': f'{data_path}/ancillary/{ndvi_file}',
+ 'SST': f'{data_path}/ancillary/{sst_file}',
+ 'ANC_DIR': f'{data_path}/ancillary'
+ }
+
+ viirs_data = rd.read_data('viirs', f'{data_path}/{mod02}', f'{data_path}/{mod03}')
+
+ viirs_data = rd.read_ancillary_data(ancillary_file_names, viirs_data)
+
+ with open(threshold_file) as f:
text = f.read()
thresholds = yml.safe_load(text)
@@ -30,7 +67,7 @@ def main():
confidence[2, :, :] = tests.nir_refl_test(viirs_data.M07.values, thresholds['Daytime_Ocean'],
thresholds['Sun_Glint'], viirs_data)
- # Note that here I'm using M05/M07 but the corresponding hi-res channelsa re I1/I2
+ # Note that here I'm using M05/M07 but the corresponding hi-res channels are I1/I2
# IMPORTANT: conf_test_dble() needs to be verified. I don't think it's working as intended at the moment
confidence[3, :, :] = tests.vis_nir_ratio_test(viirs_data.M05.values, viirs_data.M07.values,
thresholds['Daytime_Ocean'], thresholds['Sun_Glint'])
diff --git a/read_data.py b/read_data.py
index 62d1e928a347da839e27bee8bb4dc168cce4e7ab..e0b70886f313fc30f57182b726ccbe4e9d3d530a 100644
--- a/read_data.py
+++ b/read_data.py
@@ -69,56 +69,58 @@ def correct_reflectances():
pass
-def read_ancillary_data(viirs_data):
+def read_ancillary_data(file_names, 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),
+ start_time = '2022-06-22 03:36:00.000'
+ anc_dir = file_names['ANC_DIR']
+ sst_file = file_names['SST']
+ ndvi_file = file_names['NDVI']
+ geos_cnst = file_names['GEOS_constants']
+ geos_lnd = file_names['GEOS_land']
+ geos_ocn = file_names['GEOS_ocean']
+ geos1 = file_names['GEOS_atm_1']
+ geos2 = file_names['GEOS_atm_2']
+
+ dim1 = viirs_data.latitude.shape[0]
+ dim2 = viirs_data.latitude.shape[1]
+ sst = np.empty((dim1*dim2, ), dtype=np.float32)
+ ndvi = np.empty((dim1*dim2, ), dtype=np.float32)
+ eco = np.empty((dim1*dim2, ), 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),
}
- sst = anc.py_get_Reynolds_SST(viirs_data.latitude.values.reshape((3232*3200, )),
- viirs_data.longitude.values.reshape((3232*3200, )),
+ sst = anc.py_get_Reynolds_SST(viirs_data.latitude.values.reshape((dim1*dim2, )),
+ viirs_data.longitude.values.reshape((dim1*dim2, )),
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, )),
+ ndvi = anc.py_get_NDVI_background(viirs_data.latitude.values.reshape((dim1*dim2, )),
+ viirs_data.longitude.values.reshape((dim1*dim2, )),
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, )),
+ eco = anc.py_get_Olson_eco(viirs_data.latitude.values.reshape((dim1*dim2, )),
+ viirs_data.longitude.values.reshape((dim1*dim2, )),
anc_dir, eco)
- geos_data = anc.py_get_GEOS(viirs_data.latitude.values.reshape((3232*3200, )),
- viirs_data.longitude.values.reshape((3232*3200, )),
+ geos_data = anc.py_get_GEOS(viirs_data.latitude.values.reshape((dim1*dim2, )),
+ viirs_data.longitude.values.reshape((dim1*dim2, )),
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'])
+ 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)))
return data
diff --git a/scene.py b/scene.py
index 2413767020df0c44faf53549a3062ef4e41e44d7..7e9af40ac75436be840727672e3e0cd0e13b2a03 100644
--- a/scene.py
+++ b/scene.py
@@ -1,4 +1,9 @@
import numpy as np
+import ruamel_yaml as yml
+
+from glob import glob
+
+import read_data as rd
# lsf: land sea flag
_scene_list = ['ocean_day', 'ocean_night', 'land_day', 'land_night', 'snow_day', 'coast_day',
@@ -6,15 +11,18 @@ _scene_list = ['ocean_day', 'ocean_night', 'land_day', 'land_night', 'snow_day',
'polar_day_desert_coast', 'polar_day_coast', 'polar_day_land', 'polar_night_snow',
'polar_night_land', 'polar_ocean_night']
_flags = ['day', 'night', 'land', 'coast', 'sh_lake', 'sh_ocean', 'water', 'polar', 'sunglint',
- 'greenland', 'high_elevation', 'antarctica', 'desert', 'vrused', 'map_snow', 'map_ice',
+ 'greenland', 'high_elevation', 'antarctica', 'desert', 'visusd', 'vrused', 'map_snow', 'map_ice',
'ndsi_snow', 'snow', 'ice', 'new_zealand']
# temp value, need to verify what the actual bad_data value is in the C code
_bad_data = -999.0
+_dtr = np.pi/180.
+_rtd = 180./np.pi
+
# I'm defining here the flags for difference scenes. Eventually I want to find a better way of doing this
land = 1
-coast = .2
+#coast = .2
sh_lake = .3
sh_ocean = .4
water = 5
@@ -23,32 +31,92 @@ sunglint = 70
day = 100
night = 200
+# #################################################################### #
+# TEST CASE
+# data:
+datapath = '/ships19/hercules/pveglio/mvcm_viirs_hires'
+fname_mod02 = glob(f'{datapath}/VNP02MOD.A2022173.1448.001.*.uwssec.nc')[0]
+fname_mod03 = glob(f'{datapath}/VNP03MOD.A2022173.1448.001.*.uwssec.nc')[0]
+fname_img02 = glob(f'{datapath}/VNP02IMG.A2022173.1448.001.*.uwssec.nc')[0]
+fname_img03 = glob(f'{datapath}/VNP03IMG.A2022173.1448.001.*.uwssec.nc')[0]
+
+# thresholds:
+threshold_file = '/home/pveglio/mvcm_leo/thresholds/new_thresholds.mvcm.snpp.v1.0.0.yaml'
+
+# ancillary files:
+geos_atm_1 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1500.V01.nc4'
+geos_atm_2 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1800.V01.nc4'
+geos_land = 'GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20220622_1630.V01.nc4'
+geos_ocean = 'GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20220622_1630.V01.nc4'
+geos_constants = 'GEOS.fp.asm.const_2d_asm_Nx.00000000_0000.V01.nc4'
+ndvi_file = 'NDVI.FM.c004.v2.0.WS.00-04.177.hdf'
+sst_file = 'oisst.20220622'
+eco_file = 'goge1_2_img.v1'
+
+# #################################################################### #
+
+
+def test_scene():
+
+ ancillary_file_names = {'GEOS_atm_1': geos_atm_1,
+ 'GEOS_atm_2': geos_atm_2,
+ 'GEOS_land': geos_land,
+ 'GEOS_ocean': geos_ocean,
+ 'GEOS_constants': geos_constants,
+ 'NDVI': ndvi_file,
+ 'SST': sst_file,
+ 'ANC_DIR': f'{datapath}/ancillary'
+ }
+
+ viirs_data = rd.read_data('viirs', f'{fname_mod02}', f'{fname_mod03}')
+
+ viirs_data = rd.read_ancillary_data(ancillary_file_names, viirs_data)
+
+ with open(threshold_file) as f:
+ text = f.read()
+ thresholds = yml.safe_load(text)
+
+ sunglint_angle = thresholds['Sun_Glint']['bounds'][3]
+
+ scene = find_scene(viirs_data, sunglint_angle)
+
+ return scene
+
-def find_scene(data):
+def find_scene(data, sunglint_angle):
eco = data['eco'].values
lsf = data['land_water_mask'].values
lat = data['latitude'].values
lon = data['longitude'].values
sza = data['solar_zenith'].values
+ vza = data['sensor_zenith'].values
+ raz = data['relative_azimuth'].values
b065 = data['M05'].values
b086 = data['M07'].values
- elev = data[].values # !!!!!!!!!!! THIS NEEDS TO BE DEFINED IN read_data()
- day = np.zeros((3232, 3200))
- day[sza <= 85] = 1
+ elev = data['height'].values
+ ndvibk = data['ndvi'].values
+ sfctmp = data['geos_sfct'].values
+
+ dim1 = data.latitude.shape[0]
+ dim2 = data.latitude.shape[1]
- tmp = np.ones((3232, 3200))
+ day = np.zeros((dim1, dim2))
+ day[sza <= 85] = 1
+ cos_refang = np.sin(vza*_dtr) * np.sin(sza*_dtr) * np.cos(raz*_dtr) + np.cos(vza*_dtr) * np.cos(sza*_dtr)
+ refang = np.arccos(cos_refang) * _rtd
- tmp[day == 1] = day
- tmp[day == 0] = night
+ # tmp = np.ones((dim1, dim2))
+ # tmp[day == 1] = day
+ # tmp[day == 0] = night
- scene = {scn: np.zeros((3232, 3200)) for scn in _scene_list}
- scene_flag = {flg: np.zeros((3232, 3200)) for flg in _flags}
+ scene = {scn: np.zeros((dim1, dim2)) for scn in _scene_list}
+ scene_flag = {flg: np.zeros((dim1, dim2)) for flg in _flags}
scene_flag['day'][sza <= 85] = 1
scene_flag['visusd'][sza <= 85] = 1
scene_flag['night'][sza > 85] = 1
- scene_flag[np.abs(lat) > 60]['polar'] = 1
+ scene_flag['polar'][np.abs(lat) > 60] = 1
# ################# need to pass refang (once I figure out what it is) and sunglint_angle. The latter
# comes from the thresholds file. In the C code is snglnt_bounds[3]
idx = np.nonzero((scene_flag['day'] == 1) & (refang <= sunglint_angle))
@@ -59,25 +127,31 @@ def find_scene(data):
eco[idx] = 14
# start by defining anythings as land
- scene_flag['land'] = 1
- scene_flag['water'] = 0
+ scene_flag['land'] = np.ones((dim1, dim2))
+ scene_flag['water'] = np.zeros((dim1, dim2))
# Fix-up for missing ecosystem data in eastern Greenland and north-eastern Siberia.
# Without this, these regions become completely "coast".
idx = np.nonzero((lsf != 255) & (lsf == 1) | (lsf == 4))
scene_flag['land'][idx] = 1
- idx = np.nonzero((lsf != 255) & (eco == 14))
+ # idx = np.nonzero((lsf != 255) & (eco == 14))
- idx = np.nonzero((lsf != 255) & (eco == 14) & (lat < 64.0))
+ idx = np.nonzero((lsf != 255) & (eco == 14) & (lsf == 1) | (lsf == 4) & (lat < 64.0))
scene_flag['coast'][idx] = 1
- idx = np.nonzero((lsf != 255) & (eco == 14) &
- (lat >= 67.5) & (lon < -40.0) & (lon > -168.6) | (lon > -12.5))
+ idx = np.nonzero((lsf != 255) & (eco == 14) & (lsf == 1) | (lsf == 4) &
+ (lat >= 67.5) & (lon < -40.0) & (lon > -168.6))
+ scene_flag['coast'][idx] = 1
+ idx = np.nonzero((lsf != 255) & (eco == 14) & (lsf == 1) | (lsf == 4) &
+ (lat >= 67.5) & (lon > -12.5))
scene_flag['coast'][idx] = 1
- idx = np.nonzero((lsf != 255) & (eco == 14) &
- (lat >= 64.0) & (lat < 67.5) & (lon < -40.0) & (lon > -168.5) | (lon > -30.0))
+ idx = np.nonzero((lsf != 255) & (eco == 14) & (lsf == 1) | (lsf == 4) &
+ (lat >= 64.0) & (lat < 67.5) & (lon < -40.0) & (lon > -168.5))
+ scene_flag['coast'][idx] = 1
+ idx = np.nonzero((lsf != 255) & (eco == 14) & (lsf == 1) | (lsf == 4) &
+ (lat >= 64.0) & (lat < 67.5) & (lon > -30.0))
scene_flag['coast'][idx] = 1
idx = np.nonzero(lsf == 2)
@@ -88,18 +162,23 @@ def find_scene(data):
scene_flag['land'][idx] = 1
scene_flag['sh_lake'][idx] = 1
+ idx = np.nonzero((lsf == 0) | (lsf >= 5) & (lsf <= 7))
+ scene_flag['water'][idx] = 1
+ scene_flag['land'][idx] = 0
+
+ scene_flag['sh_ocean'][lsf == 0] = 1
+
# Need shallow lakes to be processed as "coast" for day, but not night
idx = np.nonzero((lsf == 3) & (day == 1))
scene_flag['coast'][idx] = 1
- idx = np.nonzero((lsf == 3) & (day == 0))
- scene_flag['sh_ocean'][idx] = 1
-
# if land/sea flag is missing, then calculate visible ratio to determine if land or water.
idx = np.nonzero((lsf == 255) & (b065 != _bad_data) & (b086 != _bad_data) & (b086/b065 > 0.9))
+ scene_flag['land'][idx] = 1
- scene_flag['land'] = 1
idx = np.nonzero((lsf == 255) & (b065 != _bad_data) & (b086 != _bad_data) & (b086/b065 <= 0.9))
+ scene_flag['land'][idx] = 0
+ scene_flag['water'][idx] = 1
# Check surface elevation
# First, define "Greenland".
@@ -147,15 +226,15 @@ def find_scene(data):
idx = np.nonzero((eco == 2) | (eco == 8) | (eco == 11) | (eco == 40) | (eco == 41) | (eco == 46) |
(eco == 51) | (eco == 52) | (eco == 59) | (eco == 71) | (eco == 50))
- scene_flag['vrused'] = 1
+ scene_flag['vrused'] = np.ones((dim1, dim2))
scene_flag['vrused'][idx] = 0
- snow_fraction = geos_data['snowfr']
- perm_ice_fraction = geos_data['landicefr']
- ice_fraction = geos_data['icefr']
+ snow_fraction = data['geos_snowfr']
+ perm_ice_fraction = data['geos_landicefr']
+ ice_fraction = data['geos_icefr']
idx = np.nonzero((snow_fraction > 0.10) & (snow_fraction <= 1.0))
- scene_flag['map_snow'] = 1
+ scene_flag['map_snow'][idx] = 1
idx = np.nonzero((perm_ice_fraction > 0.10) & (perm_ice_fraction <= 1.0))
scene_flag['map_snow'][idx] = 1
@@ -164,9 +243,11 @@ def find_scene(data):
scene_flag['map_ice'][idx] = 1
# need to define this function and write this block better
- if day == 1:
- # Run quick version of D. Hall's snow detection algorithm
- scene_flag['ndsi_snow'] = run_snow_mask()
+ # if day == 1:
+ # # Run quick version of D. Hall's snow detection algorithm
+ # # Temporarily disabled because this function does not exist and I need to decide how to implement it
+ # # scene_flag['ndsi_snow'] = run_snow_mask()
+ # scene_flag['ndsi_snow'] = np.zeros((dim1, dim2))
idx = np.nonzero((day == 1) & (water == 1) & (lat >= -60.0) & (lat <= 25.0) &
(scene_flag['map_snow'] == 1) & (scene_flag['ndsi_snow'] == 1))
@@ -187,7 +268,7 @@ def find_scene(data):
# Define New Zealand region which receives snow but snow map does not show it.
idx = np.nonzero((day == 1) & (land == 1) &
(lat >= 48.0) & (lat <= -34.0) & (lon >= 165.0) & (lon <= 180.0))
- scene_flag['new_zealand'] = 1
+ scene_flag['new_zealand'][idx] = 1
idx = np.nonzero((day == 1) & (land == 1) & (lat >= -60.0) & (lat <= 25.0) &
(scene_flag['map_snow'] == 1) & (scene_flag['ndsi_snow'] == 1) |
@@ -195,10 +276,10 @@ def find_scene(data):
scene_flag['snow'][idx] = 1
idx = np.nonzero((day == 1) & (land == 1) & (lat < -60.0))
- scne_flag['snow'][idx] = 1
+ scene_flag['snow'][idx] = 1
idx = np.nonzero((day == 1) & (land == 1) & (scene_flag['ndsi_snow'] == 1))
- scne_flag['snow'][idx] = 1
+ scene_flag['snow'][idx] = 1
idx = np.nonzero((day == 0) & (scene_flag['map_snow'] == 1) &
(sfctmp > 280.0) & (elev < 500.0))
@@ -209,7 +290,7 @@ def find_scene(data):
scene_flag['ice'][idx] = 0
idx = np.nonzero((day == 0) & (lat > 86.0))
- scene_flag['ice'] = 1
+ scene_flag['ice'][idx] = 1
##################################
# CONTINUE FROM HERE #
@@ -218,7 +299,12 @@ def find_scene(data):
# Check regional uniformity
# Check for data border pixels
# NEED TO UNDERSTAND WHAT THIS PART DOES
- return scene
+
+ # TEMP VALUES FOR DEBUGGING
+ scene_flag['lat'] = lat
+ scene_flag['lon'] = lon
+
+ return scene_flag
diff --git a/utils.py b/utils.py
index 2cf07c3a5c13abf364b2d7e70c46925e81801cca..7f90e49341d66e61b1270ad13d2afbb22bee02b4 100644
--- a/utils.py
+++ b/utils.py
@@ -317,3 +317,13 @@ def find_scene(data):
}
return scene
+
+
+def snow_mask(viirs_data):
+ # (I1 - I3) / (I1 + I3)
+ # ndsi = (b04 - b10) / (b04 + b10)
+ # ndvi = (b07 - b05) / (b05 + b07)
+ # pred_ndvi = prd_ndvi_const[0] + 3.0*ndsi
+
+ pass
+