1
|
#This is the python script for checking the quality control and SDS fields of MODIS_Oregon hdf4 tiles.
|
2
|
#File contains scripting for all checks and annotation re: what I was looking for, what I found, etc.
|
3
|
#Created by Natalie Robinson- Oct. 2011
|
4
|
|
5
|
#-------------------------------------------------------------------------------------
|
6
|
#Need to run through original .hdf files.
|
7
|
#Looking for string match: 'AUTOMATICQUALITYFLAG': 'PASSED'
|
8
|
# 'CLOUDCONTAMINATED_LST_SCREENED': 'YES'
|
9
|
|
10
|
import gdal
|
11
|
import gdalnumeric
|
12
|
import sys
|
13
|
import string
|
14
|
import os
|
15
|
import re
|
16
|
import numpy
|
17
|
from pyhdf.SD import *
|
18
|
from numpy import *
|
19
|
import warnings
|
20
|
warnings.simplefilter('ignore', DeprecationWarning)
|
21
|
import glob #probably don't need
|
22
|
import optparse #probably don't need
|
23
|
import fnmatch #probably don't need
|
24
|
from pyhdf import SD
|
25
|
import scipy
|
26
|
from scipy import *
|
27
|
|
28
|
|
29
|
sourceDir = "/data/project/organisms/MODIS_LST_Oregon"
|
30
|
|
31
|
#Find out how many files in directory for indexing
|
32
|
len= len(os.listdir(sourceDir)) #15274
|
33
|
|
34
|
#Find other QC tags to check- test run on single file
|
35
|
|
36
|
testset= gdal.Open ('/data/project/organisms/MODIS_LST_Oregon/MOD11A1.A2000065.h08v04.005.2007176143143.hdf')
|
37
|
allMeta= testset.GetMetadata ()
|
38
|
dict.keys (allMeta)
|
39
|
|
40
|
#Keys of potential interest- what are OK values?
|
41
|
|
42
|
allMeta['QAPERCENTNOTPRODUCEDCLOUD'] #Here= 9
|
43
|
allMeta['QAPERCENTNOTPRODUCEDOTHER'] #Here= 90
|
44
|
allMeta['QAPERCENTMISSINGDATA'] #Here= 0
|
45
|
allMeta['QAPERCENTINTERPOLATEDDATA'] #Here= 0
|
46
|
allMeta['QAPERCENTGOODQUALITY'] #Here= 0
|
47
|
allMeta['QAPERCENTCLOUDCOVER'] #Here= 9
|
48
|
allMeta['DAYNIGHTFLAG'] #Here= Both
|
49
|
allMeta['SCIENCEQUALITYFLAG'] #Here= Not Investigated
|
50
|
allMeta['QAPERCENTOTHERQUALITY'] #Here= 0
|
51
|
allMeta['AUTOMATICQUALITYFLAG'] #Here= Passed
|
52
|
allMeta['CLOUD_CONTAMINATED_LST_SCREENED'] #Here= YES
|
53
|
#-----------------------------------------------------------------------------------------
|
54
|
#Loop through folder and grab values for each quality control tag- evaluate values
|
55
|
|
56
|
Cloud_contam=[]
|
57
|
Day_Night=[]
|
58
|
Sci_Qual=[]
|
59
|
Auto_flag=[]
|
60
|
for dirname, dirnames, filenames in os.walk(sourceDir):
|
61
|
for filename in filenames:
|
62
|
ext = os.path.splitext(filename)[1]
|
63
|
if ext == ".hdf":
|
64
|
infile= sourceDir + "/" + filename
|
65
|
ds= gdal.Open(infile)
|
66
|
hdf_md = ds.GetMetadata()
|
67
|
for k,v in hdf_md.iteritems():
|
68
|
if k=='CLOUD_CONTAMINATED_LST_SCREENED':
|
69
|
if v== 'YES':
|
70
|
Cloud_contam.append ("ok")
|
71
|
else:
|
72
|
Cloud_contam.append ("Bad")
|
73
|
if k=='DAYNIGHTFLAG':
|
74
|
if v== 'Both':
|
75
|
Day_Night.append ("ok")
|
76
|
else:
|
77
|
Day_Night.append ("Check")
|
78
|
if k=='SCIENCEQUALITYFLAG':
|
79
|
if v== 'Not Investigated':
|
80
|
Sci_Qual.append ("?")
|
81
|
else:
|
82
|
Sci_Qual.append ("Check")
|
83
|
if k== 'AUTOMATICQUALITYFLAG':
|
84
|
if v== 'Passed':
|
85
|
Auto_flag.append("ok")
|
86
|
else:
|
87
|
Auto_flag.append("Bad")
|
88
|
|
89
|
|
90
|
if "Bad" in Cloud_contam: print ("There's a bad Cloud Contaminated LST Screen")
|
91
|
else: print ("All good on the cloud contaminated LST screen front")
|
92
|
|
93
|
|
94
|
if "Bad" in Auto_flag: print ("There's a bad auto flag")
|
95
|
else: print ("All good on the auto flag front")
|
96
|
|
97
|
|
98
|
if "Check" in Day_Night: print ("Check the day/night flags")
|
99
|
else: print ("All good on the day/night flag front")
|
100
|
|
101
|
|
102
|
if "Check" in Sci_Qual: print ("Check the science quality flags")
|
103
|
else: print ("All good on the science quality flag front?")
|
104
|
|
105
|
#------------------------------------------------------------------------------------------
|
106
|
#Loop through the folder and find ranges for some of the metadata tags
|
107
|
|
108
|
#Percent not produced other
|
109
|
Pct_Other= [] #Create empty list to hold values of specific key
|
110
|
Pct_Cloud=[]
|
111
|
Pct_Missing= [] #Create empty list to hold values of specific key
|
112
|
Pct_Interp=[]
|
113
|
Pct_Good=[]
|
114
|
Pct_CC=[]
|
115
|
Pct_OtherQual=[]
|
116
|
|
117
|
for dirname, dirnames, filenames in os.walk(sourceDir):
|
118
|
for filename in filenames:
|
119
|
ext = os.path.splitext(filename)[1]
|
120
|
if ext == ".hdf":
|
121
|
infile= sourceDir + "/" + filename
|
122
|
ds= gdal.Open(infile)
|
123
|
hdf_md = ds.GetMetadata()
|
124
|
for k,v in hdf_md.iteritems():
|
125
|
if k== 'QAPERCENTNOTPRODUCEDOTHER': Pct_Other.append(v)
|
126
|
if k== 'QAPERCENTNOTPRODUCEDCLOUD': Pct_Cloud.append(v)
|
127
|
if k== 'QAPERCENTMISSINGDATA': Pct_Missing.append(v)
|
128
|
if k== 'QAPERCENTINTERPOLATEDDATA': Pct_Interp.append(v)
|
129
|
if k== 'QAPERCENTGOODQUALITY': Pct_Good.append(v)
|
130
|
if k== 'QAPERCENTCLOUDCOVER': Pct_CC.append(v)
|
131
|
if k== 'QAPERCENTOTHERQUALITY': Pct_OtherQual.append(v)
|
132
|
|
133
|
|
134
|
print ("range Pct_Other =" + min (Pct_Other) + " to " + max(Pct_Other))
|
135
|
print ("range Pct_Cloud =" + min (Pct_Cloud) + " to " + max(Pct_Cloud))
|
136
|
print ("range Pct_Missing =" + min (Pct_Missing) + " to " + max(Pct_Missing))
|
137
|
print ("range Pct_Interp =" + min (Pct_Interp) + " to " + max(Pct_Interp))
|
138
|
print ("range Pct_Good =" + min (Pct_Good) + " to " + max(Pct_Good))
|
139
|
print ("range Pct_CC =" + min (Pct_CC) + " to " + max(Pct_CC))
|
140
|
print ("range Pct_OtherQual =" + min (Pct_OtherQual) + " to " + max(Pct_OtherQual))
|
141
|
#-----------------------------------------------------------------------------------------------
|
142
|
#Potential problem file re: Science quality flag explanation (http://landweb.nascom.nasa.gov/cgi-bin/QA_WWW/detailInfo.cgi?prod_id=MOD11A1&ver=C4)
|
143
|
|
144
|
ohno= gdal.Open('/data/project/organisms/MODIS_LST_Oregon/MOD11A1.A2003359.h09v04.005.2008163212412.hdf')
|
145
|
ohnoMeta= ohno.GetMetadata()
|
146
|
ohnoMeta['SCIENCEQUALITYFLAG'] #Comes back "Not Investigated" but this tile SHOULD BE "Suspect"
|
147
|
#--------------------------------------------------------------------------------------------
|
148
|
#After checking through known issues for MOD11A1 data (from http://landweb.nascom.nasa.gov/cgi-bin/QA_WWW/getSummary.cgi?esdt=MOD11), there are some versions that are affected and others that aren't. Check PGEversion of tiles to see if they are affected:
|
149
|
|
150
|
PGE=[]
|
151
|
for dirname, dirnames, filenames in os.walk(sourceDir):
|
152
|
for filename in filenames:
|
153
|
ext = os.path.splitext(filename)[1]
|
154
|
if ext == ".hdf":
|
155
|
infile= sourceDir + "/" + filename
|
156
|
ds= gdal.Open(infile)
|
157
|
hdf_md = ds.GetMetadata()
|
158
|
for k,v in hdf_md.iteritems():
|
159
|
if k== 'PGEVERSION': PGE.append(v)
|
160
|
|
161
|
|
162
|
print ("range PGE =" + min (PGE) + " to " + max(PGE)) #Range= 5.3.6 to 5.5.6
|
163
|
|
164
|
|
165
|
#THIS IS IMPORTANT:
|
166
|
# MOre infromation on QA/QC tags and how to check them available at: https://lpdaac.usgs.gov/products/modis_products_table/land_surface_temperature_emissivity/daily_l3_global_1km/mod11a1
|
167
|
#---------------------------------------------------------------------------------
|
168
|
#Check out SDS's
|
169
|
|
170
|
SDS's: Desc. Acceptable range: Fill Value:
|
171
|
# LST_Day_1km Daily day 1km grid LST 7500-65535 0
|
172
|
# QC_Day QC for day LST emiss. 0-255 0 (overlap)
|
173
|
# Day_view_time Time: day temp obs. 0-240 255
|
174
|
# Day_view_angle Zenith ang. for AM LST 0-130 255
|
175
|
# LST_Night_1km Daily PM 1km grid LST 7500-65535 0
|
176
|
# QC_Night QC for night LST emiss. 0-255 0 (overlap)
|
177
|
# Night_view_time Time: night temp obs. 0-240 255
|
178
|
# Night_view_angle Zenith ang. for PM LST 0-130 255
|
179
|
# Emis_31 Band 31 emiss. 1-255 0
|
180
|
# Emis_32 Band 32 emiss. 1-255 0
|
181
|
# Clear_day_cov AM clear-sky cov. 0-65535 0 (overlap)
|
182
|
# Clear_night_cov PM clear-sky cov. 0-65535 0 (overlap)
|
183
|
|
184
|
|
185
|
#For one file- tester-----------------------------------------------------------------------
|
186
|
file='/data/project/organisms/MODIS_LST_Oregon/MOD11A1.A2000065.h08v04.005.2007176143143.hdf' # The hdf file to read
|
187
|
|
188
|
#Open hdf and read SDS data
|
189
|
hdf=SD.SD(file)
|
190
|
qcday=hdf.select('QC_Day')
|
191
|
QC_Day=qcday.get()
|
192
|
qcnight=hdf.select('QC_Night')
|
193
|
QC_Night=qcnight.get()
|
194
|
lstday1km=hdf.select('LST_Day_1km')
|
195
|
LSTDay1km=lstday1km.get()
|
196
|
dvtime=hdf.select('Day_view_time')
|
197
|
DayViewTime=dvtime.get()
|
198
|
nvtime=hdf.select('Night_view_time')
|
199
|
NightViewTime=nvtime.get()
|
200
|
dvang=hdf.select('Day_view_angl')
|
201
|
DayViewAngle=dvang.get()
|
202
|
nvang=hdf.select('Night_view_angl')
|
203
|
NightViewAngle=nvang.get()
|
204
|
lstnight1km=hdf.select('LST_Night_1km')
|
205
|
LSTNight1km=lstnight1km.get()
|
206
|
clearday=hdf.select('Clear_day_cov')
|
207
|
ClearDayCov=clearday.get()
|
208
|
clearnight=hdf.select('Clear_night_cov')
|
209
|
ClearNightCov=clearnight.get()
|
210
|
emis31=hdf.select('Emis_31')
|
211
|
Emiss_Band31=emis31.get()
|
212
|
emis32=hdf.select('Emis_31')
|
213
|
Emiss_Band32=emis32.get()
|
214
|
|
215
|
print "QC_Day range:", QC_Day.min(), " to ",QC_Day.max()
|
216
|
print "QC_Night range:", QC_Night.min(), " to ",QC_Night.max()
|
217
|
print "LST Day 1km range:",LSTDay1km.min(), " to ",LSTDay1km.max() #Need to remove 0 results to see minimum real value (0=fill, don't know if values <7500 exist b/c masked by 0
|
218
|
print "LST Night 1km range:",LSTNight1km.min(), " to ",LSTNight1km.max()
|
219
|
print "Day View Time range:",DayViewTime.min(), " to ",DayViewTime.max()
|
220
|
print "Night View Time range:",NightViewTime.min(), " to ",NightViewTime.max()
|
221
|
print "Day View Angle range:",DayViewAngle.min(), " to ",DayViewAngle.max()
|
222
|
print "Night View Angle range:",NightViewAngle.min(), " to ",NightViewAngle.max()
|
223
|
print "Clear Day Coverage range:",ClearDayCov.min(), " to ",ClearDayCov.max()
|
224
|
print "Clear Night Coverage range:",ClearNightCov.min(), " to ",ClearNightCov.max()
|
225
|
print "Emissivity Band 31 range:",Emiss_Band31.min(), " to ",Emiss_Band31.max()
|
226
|
print "Emissivity Band 32 range:",Emiss_Band32.min(), " to ",Emiss_Band32.max()
|
227
|
|
228
|
#For some of these the fill value is the highest or lowest value of the array, so pixels with non-fill values need to still be checked to make sure they are within the valid range
|
229
|
|
230
|
scipy.unique(LSTDay1km) #Need to see 0 (fill), and all others >=7500
|
231
|
scipy.unique(LSTNight1km) #Need to see 0 (fill), and all others >=7500
|
232
|
scipy.unique(DayViewTime) #Need to see 255 (fill), and all others <=240
|
233
|
scipy.unique(NightViewTime) #Need to see 255 (fill), and all others <=240
|
234
|
scipy.unique(DayViewAngle) #Need to see 255 (fill), and all others <=130
|
235
|
scipy.unique(NightViewAngle) #Need to see 255 (fill), and all others <=130
|
236
|
scipy.unique(Emiss_Band31)
|
237
|
scipy.unique(Emiss_Band32)
|
238
|
|
239
|
#Passes check: LSTDay1km, LSTNight1km, DayViewTime, NightViewTime, DayViewAngle, NightViewAngle (barely, highest non-fill value=130), both Emiss Bands.
|
240
|
|
241
|
#Loop through files----------------------------------------------------------------------
|
242
|
#Bad= ANY - value
|
243
|
#Bad QC_Day and QC_Night= max > 255
|
244
|
#Bad ClearDay/ClearNight= max > 65535
|
245
|
#Bad LSTDay1km/LSTNight1km= max > 65535 or scipy.unique(LSTXXX1km)[1] < 7500
|
246
|
#Bad DayViewTime/NightViewTime= max > 255 or scipy.unique(XXXViewTime)[-2] > 240
|
247
|
#Bad DayViewAngle/NightViewAngle= max > 255 or scipy.unique(XXXViewAngle)[-2] > 130
|
248
|
#Bad Emiss_Band31/Emiss_Band32= max > 255
|
249
|
|
250
|
QCDayBad=[]
|
251
|
QCNightBad=[]
|
252
|
ClearDayBad=[]
|
253
|
ClearNightBad=[]
|
254
|
Emiss31Bad=[]
|
255
|
Emiss32Bad=[]
|
256
|
LSTDayBad=[]
|
257
|
LSTNightBad=[]
|
258
|
DayViewTimeBad=[]
|
259
|
NightViewTimeBad=[]
|
260
|
DayViewAngleBad=[]
|
261
|
NightViewAngleBad=[]
|
262
|
|
263
|
for dirname, dirnames, filenames in os.walk(sourceDir):
|
264
|
for filename in filenames:
|
265
|
ext = os.path.splitext(filename)[1]
|
266
|
if ext == ".hdf":
|
267
|
infile= sourceDir + "/" + filename
|
268
|
hdf4SDS= SD.SD(infile)
|
269
|
qcday=hdf4SDS.select('QC_Day')
|
270
|
QC_Day=qcday.get()
|
271
|
qcnight=hdf4SDS.select('QC_Night')
|
272
|
QC_Night=qcnight.get()
|
273
|
lstday1km=hdf4SDS.select('LST_Day_1km')
|
274
|
LSTDay1km=lstday1km.get()
|
275
|
dvtime=hdf4SDS.select('Day_view_time')
|
276
|
DayViewTime=dvtime.get()
|
277
|
nvtime=hdf4SDS.select('Night_view_time')
|
278
|
NightViewTime=nvtime.get()
|
279
|
dvang=hdf4SDS.select('Day_view_angl')
|
280
|
DayViewAngle=dvang.get()
|
281
|
nvang=hdf4SDS.select('Night_view_angl')
|
282
|
NightViewAngle=nvang.get()
|
283
|
lstnight1km=hdf4SDS.select('LST_Night_1km')
|
284
|
LSTNight1km=lstnight1km.get()
|
285
|
clearday=hdf4SDS.select('Clear_day_cov')
|
286
|
ClearDayCov=clearday.get()
|
287
|
clearnight=hdf4SDS.select('Clear_night_cov')
|
288
|
ClearNightCov=clearnight.get()
|
289
|
emis31=hdf4SDS.select('Emis_31')
|
290
|
Emiss_Band31=emis31.get()
|
291
|
emis32=hdf4SDS.select('Emis_31')
|
292
|
Emiss_Band32=emis32.get()
|
293
|
if QC_Day.min()<0 or QC_Day.max()>255:
|
294
|
QCDayBad.append (filename)
|
295
|
if QC_Night.min()<0 or QC_Night.max()>255:
|
296
|
QCNightBad.append(filename)
|
297
|
if ClearDayCov.min()<0 or ClearDayCov.max()>65535:
|
298
|
ClearDayBad.append (filename)
|
299
|
if ClearNightCov.min()<0 or ClearNightCov.max()>65535:
|
300
|
ClearNightBad.append (filename)
|
301
|
if Emiss_Band31.min()<0 or Emiss_Band31.max()>255:
|
302
|
Emiss31Bad.append (filename)
|
303
|
if Emiss_Band32.min()<0 or Emiss_Band32.max()>255:
|
304
|
Emiss32Bad.append (filename)
|
305
|
if LSTDay1km.min()<0 or LSTDay1km.max()>65535:
|
306
|
LSTDayBad.append (filename)
|
307
|
if LSTNight1km.min()<0 or LSTNight1km.max()>65535:
|
308
|
LSTNightBad.append (filename)
|
309
|
if DayViewTime.min()<0 or DayViewTime.max()>255:
|
310
|
DayViewTimeBad.append (filename)
|
311
|
if NightViewTime.min()<0 or NightViewTime.max()>255:
|
312
|
NightViewTimeBad.append (filename)
|
313
|
if DayViewAngle.min()<0 or DayViewAngle.max()>255:
|
314
|
DayViewAngleBad.append (filename)
|
315
|
if NightViewAngle.min()<0 or NightViewAngle.max()>255:
|
316
|
NightViewAngleBad.append (filename)
|
317
|
|
318
|
|
319
|
print QCDayBad #None found
|
320
|
print QCNightBad #None found
|
321
|
print ClearDayBad #None found
|
322
|
print ClearNightBad #None found
|
323
|
print Emiss31Bad #None found
|
324
|
print Emiss32Bad #None found
|
325
|
print LSTDayBad #None found
|
326
|
print LSTNightBad #None found
|
327
|
print DayViewTimeBad #None found
|
328
|
print NightViewTimeBad #None found
|
329
|
print DayViewAngleBad #None found
|
330
|
print NightViewAngleBad #None found
|
331
|
|
332
|
#Do searches for masked (by fill value) min or max values separately------------------------
|
333
|
LSTDayRealMin=[]
|
334
|
LSTNightRealMin=[]
|
335
|
DayViewTimeRealMax=[]
|
336
|
NightViewTimeRealMax=[]
|
337
|
DayViewAngleRealMax=[]
|
338
|
NightViewAngleRealMax=[]
|
339
|
|
340
|
for dirname, dirnames, filenames in os.walk(sourceDir):
|
341
|
for filename in filenames:
|
342
|
ext = os.path.splitext(filename)[1]
|
343
|
if ext == ".hdf":
|
344
|
infile= sourceDir + "/" + filename
|
345
|
hdf4SDS= SD.SD(infile)
|
346
|
lstday1km=hdf4SDS.select('LST_Day_1km')
|
347
|
LSTDay1km=lstday1km.get()
|
348
|
dvtime=hdf4SDS.select('Day_view_time')
|
349
|
DayViewTime=dvtime.get()
|
350
|
nvtime=hdf4SDS.select('Night_view_time')
|
351
|
NightViewTime=nvtime.get()
|
352
|
dvang=hdf4SDS.select('Day_view_angl')
|
353
|
DayViewAngle=dvang.get()
|
354
|
nvang=hdf4SDS.select('Night_view_angl')
|
355
|
NightViewAngle=nvang.get()
|
356
|
lstnight1km=hdf4SDS.select('LST_Night_1km')
|
357
|
LSTNight1km=lstnight1km.get()
|
358
|
if scipy.unique(LSTDay1km)[1]<100:
|
359
|
LSTDayRealMin.append (filename)
|
360
|
#if scipy.unique(LSTNight1km)[1]<7500: #Only value for pixels=0
|
361
|
# LSTNightRealMin.append (filename)
|
362
|
if scipy.unique(DayViewTime)[-2]>240:
|
363
|
DayViewTimeRealMax.append (filename)
|
364
|
#if scipy.unique(NightViewTime)[-2]>240: #Only value for pixels=255
|
365
|
# NightViewTimeRealMax.append (filename)
|
366
|
if scipy.unique(DayViewAngle)[-2]>130:
|
367
|
DayViewAngleRealMax.append (filename)
|
368
|
#if scipy.unique(NightViewAngle)[-2]>130: #Only value for pixels=255
|
369
|
# NightViewAngleRealMax.append (filename)
|
370
|
|
371
|
print LSTDayRealMin #Found no bad tiles
|
372
|
print DayViewTimeRealMax #Found no bad tiles
|
373
|
print DayViewAngleRealMax #Found no bad tiles
|
374
|
#---------------------------------------------------------------------------------------------
|
375
|
#Find % fill values in each file (this only makes sense for SDS's where fill value
|
376
|
# is NOT within acceptabe range.
|
377
|
|
378
|
sourceDir = "/data/project/organisms/MODIS_LST_Oregon"
|
379
|
|
380
|
#Create empty lists to hold data
|
381
|
Emiss_Band31PCTFill=[]
|
382
|
DayViewTimePCTFill=[]
|
383
|
NightViewTimePCTFill=[]
|
384
|
DayViewAnglePCTFill=[]
|
385
|
NightViewAnglePCTFill=[]
|
386
|
LSTDayPCTFill=[]
|
387
|
LSTNightPCTFill=[]
|
388
|
Emiss_Band32PCTFill=[]
|
389
|
|
390
|
for dirname, dirnames, filenames in os.walk(sourceDir):
|
391
|
for filename in filenames:
|
392
|
ext = os.path.splitext(filename)[1]
|
393
|
if ext == ".hdf":
|
394
|
infile= sourceDir + "/" + filename
|
395
|
hdf4SDS= SD.SD(infile) #Get SDS's
|
396
|
lstday1km=hdf4SDS.select('LST_Day_1km') #Select specific SDS
|
397
|
LSTDay1km=lstday1km.get() #Get data (this is an array)
|
398
|
LSTDayFlat= LSTDay1km.flatten ([LSTDay1km]) #Flatten array for searchability
|
399
|
LSTDaySort=sort(LSTDayFlat)
|
400
|
LSTDayFill= LSTDaySort.tolist() #Coerce flattened array into list
|
401
|
LSTDayPCTFill.append("%.5f" % (LSTDayFill.count (0)/1440000.0)) #Count #fills/tot pixels
|
402
|
lstnight1km=hdf4SDS.select('LST_Night_1km')
|
403
|
LSTNight1km=lstnight1km.get()
|
404
|
LSTNightFlat= LSTNight1km.flatten ([LSTNight1km])
|
405
|
LSTNightSort=sort(LSTNightFlat)
|
406
|
LSTNightFill= LSTNightSort.tolist()
|
407
|
LSTNightPCTFill.append("%.5f" % (LSTNightFill.count (0)/1440000.0))
|
408
|
dvtime=hdf4SDS.select('Day_view_time')
|
409
|
DayViewTime=dvtime.get()
|
410
|
DayViewTimeFlat= DayViewTime.flatten ([DayViewTime])
|
411
|
DayViewTimeSort=sort(DayViewTimeFlat)
|
412
|
DayViewTimeFill= DayViewTimeSort.tolist()
|
413
|
DayViewTimePCTFill.append("%.5f" % (DayViewTimeFill.count (255)/1440000.0))
|
414
|
nvtime=hdf4SDS.select('Night_view_time')
|
415
|
NightViewTime=nvtime.get()
|
416
|
NightViewTimeFlat= NightViewTime.flatten ([NightViewTime])
|
417
|
NightViewTimeSort=sort(NightViewTimeFlat)
|
418
|
NightViewTimeFill= NightViewTimeSort.tolist()
|
419
|
NightViewTimePCTFill.append("%.5f" % (NightViewTimeFill.count (255)/1440000.0))
|
420
|
dvang=hdf4SDS.select('Day_view_angl')
|
421
|
DayViewAngle=dvang.get()
|
422
|
DayViewAngleFlat= DayViewAngle.flatten ([DayViewAngle])
|
423
|
DayViewAngleSort=sort(DayViewAngleFlat)
|
424
|
DayViewAngleFill= DayViewAngleSort.tolist()
|
425
|
DayViewAnglePCTFill.append("%.5f" % (DayViewAngleFill.count (255)/1440000.0))
|
426
|
nvang=hdf4SDS.select('Night_view_angl')
|
427
|
NightViewAngle=nvang.get()
|
428
|
NightViewAngleFlat= NightViewAngle.flatten ([NightViewAngle])
|
429
|
NightViewAngleSort=sort(NightViewAngleFlat)
|
430
|
NightViewAngleFill= NightViewAngleSort.tolist()
|
431
|
NightViewAnglePCTFill.append("%.5f" % (NightViewAngleFill.count (255)/1440000.0))
|
432
|
emis31=hdf4SDS.select('Emis_31')
|
433
|
Emiss_Band31=emis31.get()
|
434
|
Emiss_Band31Flat= Emiss_Band31.flatten ([Emiss_Band31])
|
435
|
Emiss_Band31Sort=sort(Emiss_Band31Flat)
|
436
|
Emiss_Band31Fill= Emiss_Band31Sort.tolist()
|
437
|
Emiss_Band31PCTFill.append("%.5f" % (Emiss_Band31Fill.count (0)/1440000.0))
|
438
|
emis32=hdf4SDS.select('Emis_31')
|
439
|
Emiss_Band32=emis32.get()
|
440
|
Emiss_Band32Flat= Emiss_Band32.flatten ([Emiss_Band32])
|
441
|
Emiss_Band32Sort=sort(Emiss_Band32Flat)
|
442
|
Emiss_Band32Fill= Emiss_Band32Sort.tolist()
|
443
|
Emiss_Band32PCTFill.append("%.5f" % (Emiss_Band32Fill.count (0)/1440000.0))
|
444
|
|
445
|
|
446
|
|
447
|
SDSFile= open('/data/project/organisms/MODIS_LST_Oregon/Test/SDS_PctFills.txt','w')
|
448
|
#SDSFile.write(" ".join(["%s" % el for el in filenames])) #Doesn't work- filenames is a list
|
449
|
SDSFile.write("%s\n" % LSTDayPCTFill)
|
450
|
SDSFile.write("%s\n" % LSTNightPCTFill)
|
451
|
SDSFile.write("%s\n" % DayViewTimePCTFill)
|
452
|
SDSFile.write("%s\n" % NightViewTimePCTFill)
|
453
|
SDSFile.write("%s\n" % DayViewAnglePCTFill)
|
454
|
SDSFile.write("%s\n" % NightViewAnglePCTFill)
|
455
|
SDSFile.write("%s\n" % Emiss_Band31PCTFill)
|
456
|
SDSFile.write("%s\n" % Emiss_Band32PCTFill)
|
457
|
SDSFile.close()
|
458
|
|
459
|
|