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', 'snow', 'ice', 'new_zealand']
# temp value, need to verify what the actual bad_data value is in the C code
_bad_data = -999.0
# 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
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
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
# Force consistency between lsf and ecosystem type for water
idx = np.nonzero((lsf == 0) | (lsf >= 5) & (lsf < 7))
eco[idx] = 14
# start by defining anythings as land
scene_flag['land'] = 1
scene_flag['water'] = 0
# 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) & (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))
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))
scene_flag['coast'][idx] = 1
idx = np.nonzero(lsf == 2)
scene_flag['coast'][idx] = 1
scene_flag['land'][idx] = 1
idx = np.nonzero(lsf == 3)
scene_flag['land'][idx] = 1
scene_flag['sh_lake'][idx] = 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'] = 1
idx = np.nonzero((lsf == 255) & (b065 != _bad_data) & (b086 != _bad_data) & (b086/b065 <= 0.9))
# Check surface elevation
# First, define "Greenland".
idx = np.nonzero((scene_flag['land'] == 1) &
(lat >= 60.0) & (lat < 67.0) & (lon >= -60.0) & (lon < -30.0))
scene_flag['greenland'][idx] = 1
idx = np.nonzero((scene_flag['land'] == 1) &
(lat >= 67.0) & (lat < 75.0) & (lon >= -60.0) & (lon < -10.0))
scene_flag['greenland'][idx] = 1
idx = np.nonzero((scene_flag['land'] == 1) &
(lat >= 75.0) & (lon >= -70.0) & (lon < -10.0))
scene_flag['greenland'][idx] = 1
scene_flag['high_elevation'][elev > 2000] = 1
idx = np.nonzero((elev > 200) & (scene_flag['greenland'] == 1) & (scene_flag['land'] == 1))
scene_flag['high_elevation'][idx] = 1
idx = np.nonzero((lat >= 75.0) & (lat <= 79.0) & (lon >= -73.0) & (lon <= -50.0) &
(scene_flag['land'] == 1))
scene_flag['high_elevation'][idx] = 1
scene_flag['antarctica'][lat < -60.0] = 1
##################################
# somewhere here I need to add #
# the 11um elevation correction #
##################################
# this is a temporary variable for the 11um elevation correction
elev_correction = elev/1000.0 * 5.0
## Get surface temperature from NWP and SST fields
## if it's land get it from GDAS/GEOS5
#sfctmp[scene_flag['land'] == 1] = sfct
## otherwise use the ReynoldsSST
#sfctmp[scene_flag['land'] == 0] = reynSST
# Use background NDVI to define "desert"
idx = np.nonzero((scene_flag['land'] == 1) & (ndvibk < 0.3))
scene_flag['desert'][idx] = 1
idx = np.nonzero((scene_flag['land'] == 1) & (lat < -69.0))
scene_flag['desert'][idx] = 1
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'][idx] = 0
snow_fraction = geos_data['snowfr']
perm_ice_fraction = geos_data['landicefr']
ice_fraction = geos_data['icefr']
idx = np.nonzero((snow_fraction > 0.10) & (snow_fraction <= 1.0))
scene_flag['map_snow'] = 1
idx = np.nonzero((perm_ice_fraction > 0.10) & (perm_ice_fraction <= 1.0))
scene_flag['map_snow'][idx] = 1
idx = np.nonzero((ice_fraction > 0.10) & (ice_fraction <= 1.0))
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()
idx = np.nonzero((day == 1) & (water == 1) & (lat >= -60.0) & (lat <= 25.0) &
(scene_flag['map_snow'] == 1) & (scene_flag['ndsi_snow'] == 1))
scene_flag['ice'][idx] = 1
idx = np.nonzero((day == 1) & (water == 1) & (lat < -60.0) &
(scene_flag['ndsi_snow'] == 1))
scene_flag['ice'][idx] = 1
idx = np.nonzero((day == 1) & (water == 1) & (lsf == 3) | (lsf == 5) &
(scene_flag['ndsi_snow'] == 1))
scene_flag['ice'][idx] = 1
idx = np.nonzero((day == 1) & (water == 1) &
(scene_flag['map_ice'] == 1) & (scene_flag['ndsi_snow'] == 1))
scene_flag['ice'][idx] = 1
# 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
idx = np.nonzero((day == 1) & (land == 1) & (lat >= -60.0) & (lat <= 25.0) &
(scene_flag['map_snow'] == 1) & (scene_flag['ndsi_snow'] == 1) |
(scene_flag['new_zealand'] == 1))
scene_flag['snow'][idx] = 1
idx = np.nonzero((day == 1) & (land == 1) & (lat < -60.0))
scne_flag['snow'][idx] = 1
idx = np.nonzero((day == 1) & (land == 1) & (scene_flag['ndsi_snow'] == 1))
scne_flag['snow'][idx] = 1
idx = np.nonzero((day == 0) & (scene_flag['map_snow'] == 1) &
(sfctmp > 280.0) & (elev < 500.0))
scene_flag['snow'][idx] = 0
idx = np.nonzero((day == 0) & (scene_flag['map_snow'] == 1) &
(sfctmp > 280.0) & (elev < 500.0))
scene_flag['ice'][idx] = 0
idx = np.nonzero((day == 0) & (lat > 86.0))
scene_flag['ice'] = 1
##################################
# CONTINUE FROM HERE #
##################################
# Check regional uniformity
# Check for data border pixels
# NEED TO UNDERSTAND WHAT THIS PART DOES
return scene