scene seems to be indentified correctly. code works (i.e. tested) up to scene identification

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

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'])

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

scene.py

 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
 
 
 

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
+