initial draft of the function that defines the scene based on ancillary inputs and geolocation

scene.py

+import numpy as np
+
+# lsf: land sea flag
+_scene_list = ['ocean_day', 'ocean_night', 'land_day', 'land_night', 'snow_day', 'coast_day',
+               'desert_day', 'antarctic_day', 'polar_day_snow', 'polar_day_desert',
+               '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',
+          'ndsi_snow', 'ice', 'new_zealand']
+
+# I'm defining here the flags for difference scenes. Eventually I want to find a better way of doing this
+land = 1
+coast = .2
+sh_lake = .3
+sh_ocean = .4
+water = 5
+polar = 60
+sunglint = 70
+day = 100
+night = 200
+
+
+def find_scene(data):
+    eco = data['eco'].values
+    lsf = data['land_water_mask'].values
+    lat = data['latitude'].values
+    lon = data['longitude'].values
+    sza = data['solar_zenith'].values
+    day = np.zeros((3232, 3200))
+    day[sza <= 85] = 1
+
+    tmp = np.ones((3232, 3200))
+
+    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_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
+    ################## 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))
+    scene_flag['sunglint'][idx] = 1
+
+    ##################################
+    #      CONTINUE FROM HERE        #
+    ##################################
+
+    # Force consistency between lsf and ecosystem type for water
+    idx = np.nonzero((lsf == 0) | ((lsf >= 5) & (lsf <= 7)))
+    eco[idx] = 14
+
+    idx = np.nonzero((lsf == 1) | (lsf == 4))
+    tmp[idx] = tmp[idx] + land
+
+    idx = np.nonzero(eco == 14)
+
+    # fix-up for missing ecosystem data in eastern Greenland and north-eastern Siberia.
+    # Without this, these regions become completely "coast".
+    idx1 = np.nonzero((eco == 14) & (lat < 64))
+    idx2 = np.nonzero((eco == 14) & (lat >= 67.5) & (lon > -168.5) & (lon < -40))
+    idx3 = np.nonzero((eco == 14) & (lat >= 67.5) & (lon > -30))
+    idx4 = np.nonzero((eco == 14) & (lat >= 64) & (lat < 67.5) & (lon > -168.5) & (lon < -40))
+    idx5 = np.nonzero((eco == 14) & (lat >= 64) & (lat < 67.5) & (lon > -30))
+    tmp[idx1] = tmp[idx1] + coast
+    tmp[idx2] = tmp[idx2] + coast
+    tmp[idx3] = tmp[idx3] + coast
+    tmp[idx4] = tmp[idx4] + coast
+    tmp[idx5] = tmp[idx5] + coast
+
+    tmp[lsf == 2] = tmp[lsf == 2] + land + coast
+    # need shallow lakes to be processed as "coast" for day, but not night
+    tmp[(lsf == 3) & (day == 1)] = tmp[(lsf == 3) & (day == 1)] + land + sh_lake + coast
+    tmp[(lsf > 4) & (lsf != 255)] = tmp[(lsf > 4) & (lsf != 255)] + water + sh_ocean
+
+    # need to cover the case where eco is 255 and I use te vis ratio between gbnd065 and gbnd086
+    # to determine if land or water
+
+    greenland_flag = 0
+    hi_elevation_flag = 0
+    antarctic_flag = 0
+    desert_flag = 0
+    vrused = 0 # visible ratio test used based on eco_type. Check C code
+
+    map_snow_flag = 0
+    map_ice_flag = 0
+    ndsi_snow_flag = 0
+    ice_flag = 0
+    new_zealand_flag = 0
+
+    return scene
+
+
+