Newer
Older
Paolo Veglio
committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
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