1 |
7526fb1c
|
Jim Regetz
|
#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 |
|
|
|