Revision 90913004
Added by Jim Regetz over 12 years ago
- ID 90913004954434385f3614b97ee21b8a855f249e
- Parent e84c3d48
| climate/extra/test-lst-8day-avg.py | ||
|---|---|---|
| 4 | 4 |
# pixel-wise averages excluding missing values, for comparison against |
| 5 | 5 |
# the corresponding 8-day LST product obtained from LP DAAC. Script is |
| 6 | 6 |
# hard-coded to use the 8-day LST grid for tile h09v04 on 2000.03.05, |
| 7 |
# for no particular reason. |
|
| 7 |
# for no particular reason. Also compares the observed count of non-null |
|
| 8 |
# daily LST values (0-8) to the number of clear sky days as indicated in |
|
| 9 |
# the 8-day Clear_sky_days layer. |
|
| 8 | 10 |
# |
| 9 | 11 |
# 8-day target: |
| 10 | 12 |
# MOD11A2.A2000065.h09v04.005.2007176162644.hdf |
| ... | ... | |
| 38 | 40 |
|
| 39 | 41 |
# load 11A2 (8-day) tile of interest |
| 40 | 42 |
name_of_8day_hdf = 'MOD11A2.A2000065.h09v04.005.2007176162644.hdf' |
| 43 |
# first load daytime LST |
|
| 41 | 44 |
lst8day = 'HDF4_EOS:EOS_GRID:%s/%s:%s' % ( |
| 42 | 45 |
'/home/layers/data/climate/MOD11A2.005-8day-1km-L3-LST', |
| 43 | 46 |
name_of_8day_hdf, |
| ... | ... | |
| 54 | 57 |
proj4='+proj=sinu +a=6371007.181 +b=6371007.181 +ellps=sphere') |
| 55 | 58 |
gs.run_command('g.region', rast='LST_8day_h09v04')
|
| 56 | 59 |
|
| 60 |
# now load clear sky days and generate count raster |
|
| 61 |
csd8day = 'HDF4_EOS:EOS_GRID:%s/%s:%s' % ( |
|
| 62 |
'/home/layers/data/climate/MOD11A2.005-8day-1km-L3-LST', |
|
| 63 |
name_of_8day_hdf, |
|
| 64 |
'MODIS_Grid_8Day_1km_LST:Clear_sky_days') |
|
| 65 |
gs.run_command('r.in.gdal', input=csd8day, output='CSD_8day_h09v04')
|
|
| 66 |
gs.mapcalc('nCSD_8day_h09v04 = %s' %
|
|
| 67 |
'+'.join(['((CSD_8day_h09v04 >> %d) & 1)' % pos for pos in range(8)])) |
|
| 68 |
|
|
| 57 | 69 |
# load the corresponding dailies |
| 58 | 70 |
ordir = '/home/layers/data/climate/MOD11A1.004-OR-orig-2000' |
| 59 | 71 |
LST = list() |
| 72 |
CDC = list() |
|
| 60 | 73 |
for doy in range(65, 72+1): |
| 61 | 74 |
pathglob = os.path.join(ordir, 'MOD11A1.A2000%03d.h09v04*hdf' % doy) |
| 62 | 75 |
try: |
| ... | ... | |
| 69 | 82 |
gs.run_command('r.in.gdal',
|
| 70 | 83 |
input='HDF4_EOS:EOS_GRID:%s:MODIS_Grid_Daily_1km_LST:LST_Day_1km' % hdfpath, |
| 71 | 84 |
output=LST[-1]) |
| 85 |
CDC.append('CDC_' + '_'.join(hdfname.split('.')[1:3]))
|
|
| 86 |
gs.run_command('r.in.gdal',
|
|
| 87 |
input='HDF4_EOS:EOS_GRID:%s:MODIS_Grid_Daily_1km_LST:Clear_day_cov' % hdfpath, |
|
| 88 |
output=CDC[-1]) |
|
| 72 | 89 |
|
| 73 |
# generate the 8-day average
|
|
| 90 |
# generate the 8-day mean and number of observations
|
|
| 74 | 91 |
numerator = '(%s)' % '+'.join(['if(isnull(%s), 0, %s)' % (m, m) |
| 75 | 92 |
for m in LST]) |
| 76 | 93 |
denominator = '(%s)' % '+'.join(['if(!isnull(%s))' % m |
| 77 | 94 |
for m in LST]) |
| 78 |
gs.mapcalc('LST_mean_h09v04 = round(float(%s)/%s)' % (numerator, denominator))
|
|
| 95 |
gs.mapcalc('nobs_LST_h09v04 = %s' % denominator)
|
|
| 96 |
gs.mapcalc('mean_LST_h09v04 = round(float(%s)/nobs_LST_h09v04)' % numerator)
|
|
| 79 | 97 |
|
| 80 |
gs.os.environ['GRASS_OVERWRITE'] = '0' |
|
| 98 |
# generate the 8-day sum of Clear_day_cov |
|
| 99 |
gs.mapcalc('sum_CDC_h09v04 = %s' % '+'.join(['if(isnull(%s), 0, %s)' % (m, m)
|
|
| 100 |
for m in CDC])) |
|
| 81 | 101 |
|
| 82 |
# compare grid values and report out... |
|
| 83 |
theirs = garray.array() |
|
| 84 |
theirs.read('LST_8day_h09v04')
|
|
| 85 |
print 'Average of their 8-day product:', theirs.mean() |
|
| 102 |
gs.os.environ['GRASS_OVERWRITE'] = '0' |
|
| 86 | 103 |
|
| 87 |
ours = garray.array() |
|
| 88 |
ours.read('LST_mean_h09v04')
|
|
| 89 |
print 'Average of our 8-day product:', ours.mean() |
|
| 104 |
# compare LST grid values and report out... |
|
| 105 |
mean_theirs = garray.array() |
|
| 106 |
mean_theirs.read('LST_8day_h09v04')
|
|
| 107 |
print 'Average of their 8-day product:', mean_theirs.mean() |
|
| 108 |
mean_ours = garray.array() |
|
| 109 |
mean_ours.read('mean_LST_h09v04')
|
|
| 110 |
print 'Average of our 8-day product:', mean_ours.mean() |
|
| 111 |
if (mean_ours==mean_theirs).all(): |
|
| 112 |
print 'Successfully reproduced LST_Day_1km from', name_of_8day_hdf |
|
| 113 |
else: |
|
| 114 |
print 'Failed to reproduce LST_Day_1km from', name_of_8day_hdf |
|
| 90 | 115 |
|
| 91 |
if (ours==theirs).all(): |
|
| 92 |
print 'Successfully reproduced', name_of_8day_hdf |
|
| 116 |
# compare nObs grid values and report out... |
|
| 117 |
nobs_theirs = garray.array() |
|
| 118 |
nobs_theirs.read('nCSD_8day_h09v04')
|
|
| 119 |
print 'Average of counts based on Clear_sky_days:', nobs_theirs.mean() |
|
| 120 |
nobs_ours = garray.array() |
|
| 121 |
nobs_ours.read('nobs_LST_h09v04')
|
|
| 122 |
print 'Average of our N_obs:', nobs_ours.mean() |
|
| 123 |
if (nobs_ours==nobs_theirs).all(): |
|
| 124 |
print 'Successfully reproduced Clear_sky_days counts from', name_of_8day_hdf |
|
| 93 | 125 |
else: |
| 94 |
print 'Failed to reproduce', name_of_8day_hdf |
|
| 126 |
print 'Failed to reproduce Clear_sky_days counts from', name_of_8day_hdf
|
|
| 95 | 127 |
|
| 96 | 128 |
# clean up |
| 97 | 129 |
gs.run_command('g.gisenv', set='LOCATION_NAME=%s' % orig_location)
|
Also available in: Unified diff
augmented test script to reproduce 8-day LST counts from daily LST