#!/usr/bin/python

#

# Python test script using GRASS to load 8 daily LST grids and calculate

# pixel-wise averages excluding missing values, for comparison against

# the corresponding 8-day LST product obtained from LP DAAC. Script is

# hard-coded to use the 8-day LST grid for tile h09v04 on 2000.03.05,

# for no particular reason.

#

# 8-day target:

#  MOD11A2.A2000065.h09v04.005.2007176162644.hdf

#

# Note that it appears Version 5 8-day means were rounded to the nearest

# integer, as implemented below. However, reproducing the Version 4

# product instead appears to require truncating the mean. This

# difference can be expressed as follows in terms of GRASS mapcalc,

# where 'sum' and 'n' are integer rasters:

#  * v004: sum/n                <-- integer division

#  * v005: round(float(sum)/n)  <-- floating point division, then round

#

# Script can be run from within an existing GRASS 7 session; all work is

# done within a temporary LOCATION that gets removed at the end.

# Untested in earlier versions of GRASS.

#

# Jim Regetz

# NCEAS

# 16-May-2012

import os, glob, shutil

import numpy as np

import grass.script as gs

import grass.script.array as garray

tmp_location = 'tmp' + str(os.getpid())

orig_location = gs.gisenv()['LOCATION_NAME']

orig_mapset = gs.gisenv()['MAPSET']

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

# load 11A2 (8-day) tile of interest

name_of_8day_hdf = 'MOD11A2.A2000065.h09v04.005.2007176162644.hdf'

lst8day = 'HDF4_EOS:EOS_GRID:%s/%s:%s' % (

     '/home/layers/data/climate/MOD11A2.005-8day-1km-L3-LST',

     name_of_8day_hdf,

     'MODIS_Grid_8Day_1km_LST:LST_Day_1km')

gs.run_command('r.in.gdal', input=lst8day, output='LST_8day_h09v04',

    location=tmp_location)

# switch to new location

gs.run_command('g.gisenv', set='LOCATION_NAME=%s' % tmp_location)

gs.run_command('g.gisenv', set='MAPSET=PERMANENT')

# set GRASS region to match

gs.run_command('g.proj', flags='c',

    proj4='+proj=sinu +a=6371007.181 +b=6371007.181 +ellps=sphere')

gs.run_command('g.region', rast='LST_8day_h09v04')

# load the corresponding dailies

ordir = '/home/layers/data/climate/MOD11A1.004-OR-orig-2000'

LST = list()

for doy in range(65, 72+1):

    pathglob = os.path.join(ordir, 'MOD11A1.A2000%03d.h09v04*hdf' % doy)

    try:

        hdfpath = glob.glob(pathglob)[0]

    except:

        print '*** problem finding', pathglob, '***'

        continue

    hdfname = os.path.basename(hdfpath)

    LST.append('LST_' + '_'.join(hdfname.split('.')[1:3]))

    gs.run_command('r.in.gdal',

        input='HDF4_EOS:EOS_GRID:%s:MODIS_Grid_Daily_1km_LST:LST_Day_1km' % hdfpath,

        output=LST[-1])

# generate the 8-day average

numerator = '(%s)' % '+'.join(['if(isnull(%s), 0, %s)' % (m, m)

    for m in LST])

denominator = '(%s)' % '+'.join(['if(!isnull(%s))' % m

    for m in LST])

gs.mapcalc('LST_mean_h09v04 = round(float(%s)/%s)' % (numerator, denominator))

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

# compare grid values and report out...

theirs = garray.array()

theirs.read('LST_8day_h09v04')

print 'Average of their 8-day product:', theirs.mean()

ours = garray.array()

ours.read('LST_mean_h09v04')

print 'Average of our 8-day product:', ours.mean()

if (ours==theirs).all():

    print 'Successfully reproduced', name_of_8day_hdf

else:

    print 'Failed to reproduce', name_of_8day_hdf

# clean up

gs.run_command('g.gisenv', set='LOCATION_NAME=%s' % orig_location)

gs.run_command('g.gisenv', set='MAPSET=%s' % orig_mapset)

shutil.rmtree(os.path.join(gs.gisenv()['GISDBASE'], tmp_location))