Project

General

Profile

« Previous | Next » 

Revision 934d6ab7

Added by Jim Regetz over 12 years ago

  • ID 934d6ab78e67fe65dd10eaa0acdf44b9cdb06618

added script to reproduce calculation of 8-day LST from daily LST

View differences:

climate/extra/test-lst-8day-avg.py
1
#!/usr/bin/python
2
#
3
# Python test script using GRASS to load 8 daily LST grids and calculate
4
# pixel-wise averages excluding missing values, for comparison against
5
# the corresponding 8-day LST product obtained from LP DAAC. Script is
6
# hard-coded to use the 8-day LST grid for tile h09v04 on 2000.03.05,
7
# for no particular reason.
8
#
9
# 8-day target:
10
#  MOD11A2.A2000065.h09v04.005.2007176162644.hdf
11
#
12
# Note that it appears Version 5 8-day means were rounded to the nearest
13
# integer, as implemented below. However, reproducing the Version 4
14
# product instead appears to require truncating the mean. This
15
# difference can be expressed as follows in terms of GRASS mapcalc,
16
# where 'sum' and 'n' are integer rasters:
17
#  * v004: sum/n                <-- integer division
18
#  * v005: round(float(sum)/n)  <-- floating point division, then round
19
#
20
# Script can be run from within an existing GRASS 7 session; all work is
21
# done within a temporary LOCATION that gets removed at the end.
22
# Untested in earlier versions of GRASS.
23
#
24
# Jim Regetz
25
# NCEAS
26
# 16-May-2012
27

  
28
import os, glob, shutil
29
import numpy as np
30
import grass.script as gs
31
import grass.script.array as garray
32

  
33
tmp_location = 'tmp' + str(os.getpid())
34
orig_location = gs.gisenv()['LOCATION_NAME']
35
orig_mapset = gs.gisenv()['MAPSET']
36

  
37
gs.os.environ['GRASS_OVERWRITE'] = '1'
38

  
39
# load 11A2 (8-day) tile of interest
40
name_of_8day_hdf = 'MOD11A2.A2000065.h09v04.005.2007176162644.hdf'
41
lst8day = 'HDF4_EOS:EOS_GRID:%s/%s:%s' % (
42
     '/home/layers/data/climate/MOD11A2.005-8day-1km-L3-LST',
43
     name_of_8day_hdf,
44
     'MODIS_Grid_8Day_1km_LST:LST_Day_1km')
45
gs.run_command('r.in.gdal', input=lst8day, output='LST_8day_h09v04',
46
    location=tmp_location)
47

  
48
# switch to new location
49
gs.run_command('g.gisenv', set='LOCATION_NAME=%s' % tmp_location)
50
gs.run_command('g.gisenv', set='MAPSET=PERMANENT')
51

  
52
# set GRASS region to match
53
gs.run_command('g.proj', flags='c',
54
    proj4='+proj=sinu +a=6371007.181 +b=6371007.181 +ellps=sphere')
55
gs.run_command('g.region', rast='LST_8day_h09v04')
56

  
57
# load the corresponding dailies
58
ordir = '/home/layers/data/climate/MOD11A1.004-OR-orig-2000'
59
LST = list()
60
for doy in range(65, 72+1):
61
    pathglob = os.path.join(ordir, 'MOD11A1.A2000%03d.h09v04*hdf' % doy)
62
    try:
63
        hdfpath = glob.glob(pathglob)[0]
64
    except:
65
        print '*** problem finding', pathglob, '***'
66
        continue
67
    hdfname = os.path.basename(hdfpath)
68
    LST.append('LST_' + '_'.join(hdfname.split('.')[1:3]))
69
    gs.run_command('r.in.gdal',
70
        input='HDF4_EOS:EOS_GRID:%s:MODIS_Grid_Daily_1km_LST:LST_Day_1km' % hdfpath,
71
        output=LST[-1])
72

  
73
# generate the 8-day average
74
numerator = '(%s)' % '+'.join(['if(isnull(%s), 0, %s)' % (m, m)
75
    for m in LST])
76
denominator = '(%s)' % '+'.join(['if(!isnull(%s))' % m
77
    for m in LST])
78
gs.mapcalc('LST_mean_h09v04 = round(float(%s)/%s)' % (numerator, denominator))
79

  
80
gs.os.environ['GRASS_OVERWRITE'] = '0'
81

  
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()
86

  
87
ours = garray.array()
88
ours.read('LST_mean_h09v04')
89
print 'Average of our 8-day product:', ours.mean()
90

  
91
if (ours==theirs).all():
92
    print 'Successfully reproduced', name_of_8day_hdf
93
else:
94
    print 'Failed to reproduce', name_of_8day_hdf
95

  
96
# clean up
97
gs.run_command('g.gisenv', set='LOCATION_NAME=%s' % orig_location)
98
gs.run_command('g.gisenv', set='MAPSET=%s' % orig_mapset)
99
shutil.rmtree(os.path.join(gs.gisenv()['GISDBASE'], tmp_location))

Also available in: Unified diff