Newer
Older
Paolo Veglio
committed
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
Paolo Veglio
committed
# 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()
Paolo Veglio
committed
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
Paolo Veglio
committed
# 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
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
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
# 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
Paolo Veglio
committed
##################################
# somewhere here I need to add #
# the 11um elevation correction #
Paolo Veglio
committed
##################################
# this is a temporary variable for the 11um elevation correction
elev_correction = elev/1000.0 * 5.0
Paolo Veglio
committed
## 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
Paolo Veglio
committed
# 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
Paolo Veglio
committed
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
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 #
##################################
Paolo Veglio
committed
# Check regional uniformity
# Check for data border pixels
# NEED TO UNDERSTAND WHAT THIS PART DOES