Revision eaf6c10b
Added by Jim Regetz about 12 years ago
- ID eaf6c10b5bd677d33c8d02e35c6b128e39cb9aa0
| climate/extra/aggregate-daily-lst.py | ||
|---|---|---|
| 7 | 7 |
# |
| 8 | 8 |
# TODO: |
| 9 | 9 |
# - functionalize to encapsulate high level procedural steps |
| 10 |
# - as specific example of above, write a function that: |
|
| 11 |
# - takes hdf path |
|
| 12 |
# - reads in both lst and qc |
|
| 13 |
# - nulls any lst values with bad qc scores |
|
| 14 |
# - removes qc |
|
| 15 |
# - returns name of stored qc-adjusted lst raster |
|
| 16 | 10 |
# - create a 'main' function to allow script to be run as a command |
| 17 | 11 |
# that can accept arguments? |
| 18 | 12 |
# - put downloaded data somewhere other than working directory? |
| ... | ... | |
| 74 | 68 |
gs.mapcalc('mean_%s = %s/nobs_%s' % (name, numerator, name),
|
| 75 | 69 |
overwrite=overwrite) |
| 76 | 70 |
|
| 71 |
def load_qc_adjusted_lst(hdf): |
|
| 72 |
"""Load LST_Day_1km into GRASS from the specified hdf file, and |
|
| 73 |
nullify any cells for which QA flags indicate anything other than |
|
| 74 |
high quality. |
|
| 75 |
|
|
| 76 |
Argument: |
|
| 77 |
hdf -- local path to the 11A1 HDF file |
|
| 78 |
|
|
| 79 |
Returns the name of the QC-adjusted LST map created in GRASS. |
|
| 80 |
""" |
|
| 81 |
lstname = 'LST_' + '_'.join(os.path.basename(hdf).split('.')[1:3])
|
|
| 82 |
# read in lst grid |
|
| 83 |
gs.run_command('r.in.gdal',
|
|
| 84 |
input='HDF4_EOS:EOS_GRID:%s:MODIS_Grid_Daily_1km_LST:LST_Day_1km' % hdf, |
|
| 85 |
output=lstname) |
|
| 86 |
# read in qc grid |
|
| 87 |
gs.run_command('r.in.gdal',
|
|
| 88 |
input = 'HDF4_EOS:EOS_GRID:%s:MODIS_Grid_Daily_1km_LST:QC_Day' % hdf, |
|
| 89 |
output = 'qc_this_day') |
|
| 90 |
gs.run_command('g.region', rast=lstname)
|
|
| 91 |
# null out any LST cells for which QC is not high [00] |
|
| 92 |
gs.mapcalc('${lst} = if((qc_this_day & 0x03)==0, ${lst}, null())',
|
|
| 93 |
lst=lstname, overwrite=True) |
|
| 94 |
# clean up |
|
| 95 |
gs.run_command('g.remove', rast='qc_this_day')
|
|
| 96 |
# return name of qc-adjusted LST raster in GRASS |
|
| 97 |
return lstname |
|
| 98 |
|
|
| 77 | 99 |
#---------------------------------------- |
| 78 | 100 |
# set up a (temporary?) GRASS data store |
| 79 | 101 |
#---------------------------------------- |
| ... | ... | |
| 130 | 152 |
|
| 131 | 153 |
#----------------------------------------------------------------------- |
| 132 | 154 |
gs.os.environ['GRASS_OVERWRITE'] = '1' |
| 133 |
|
|
| 134 |
LST = list() |
|
| 135 |
for hdf in hdfs: |
|
| 136 |
hdflabel = '_'.join(hdf.split('.')[1:3])
|
|
| 137 |
LST.append('LST_' + hdflabel)
|
|
| 138 |
# NOTE: this gives 'datum not recognized' warning, ok to ignore? |
|
| 139 |
# read in lst grid |
|
| 140 |
gs.run_command('r.in.gdal',
|
|
| 141 |
input='HDF4_EOS:EOS_GRID:%s:MODIS_Grid_Daily_1km_LST:LST_Day_1km' % hdf, |
|
| 142 |
output=LST[-1]) |
|
| 143 |
# read in qc grid |
|
| 144 |
gs.run_command('r.in.gdal',
|
|
| 145 |
input = 'HDF4_EOS:EOS_GRID:%s:MODIS_Grid_Daily_1km_LST:QC_Day' % hdf, |
|
| 146 |
output = 'qc_this_day') |
|
| 147 |
gs.run_command('g.region', rast=LST[-1])
|
|
| 148 |
# null out any LST cells for which QC is not high [00] |
|
| 149 |
gs.mapcalc('${lst} = if((qc_this_day & 0x03)==0, ${lst}, null())',
|
|
| 150 |
lst=LST[-1]) |
|
| 155 |
# load and qc-adjust all daily LSTs, return list of the GRASS map names |
|
| 156 |
LST = [load_qc_adjusted_lst(hdf) for hdf in hdfs] |
|
| 151 | 157 |
# calculate mean and nobs rasters |
| 152 | 158 |
calc_clim(LST, 'LST_%s_%d_%03d_%03d' % (tile, year, start_doy, end_doy)) |
| 153 | 159 |
# clean up |
| 154 | 160 |
gs.run_command('g.remove', rast=','.join(LST))
|
| 155 |
gs.run_command('g.remove', rast='qc_this_day')
|
|
| 156 | 161 |
gs.os.environ['GRASS_OVERWRITE'] = '0' |
| 157 | 162 |
#----------------------------------------------------------------------- |
| 158 | 163 |
|
Also available in: Unified diff
wrote function to encapsulate LST loading and QC-adjustment