#!/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. Also compares the observed count of non-null

# daily LST values (0-8) to the number of clear sky days as indicated in

# the 8-day Clear_sky_days layer.

#

# 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'

# first load daytime LST

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')

# now load clear sky days and generate count raster

csd8day = '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:Clear_sky_days')

gs.run_command('r.in.gdal', input=csd8day, output='CSD_8day_h09v04')

gs.mapcalc('nCSD_8day_h09v04 = %s' %

    '+'.join(['((CSD_8day_h09v04 >> %d) & 1)' % pos for pos in range(8)]))

# load the corresponding dailies

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

LST = list()

CDC = 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])

    CDC.append('CDC_' + '_'.join(hdfname.split('.')[1:3]))

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

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

        output=CDC[-1])

# generate the 8-day mean and number of observations

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('nobs_LST_h09v04 = %s' % denominator)

gs.mapcalc('mean_LST_h09v04 = round(float(%s)/nobs_LST_h09v04)' % numerator)

# generate the 8-day sum of Clear_day_cov

gs.mapcalc('sum_CDC_h09v04 = %s' % '+'.join(['if(isnull(%s), 0, %s)' % (m, m)

    for m in CDC]))

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

# compare LST grid values and report out...

mean_theirs = garray.array()

mean_theirs.read('LST_8day_h09v04')

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

mean_ours = garray.array()

mean_ours.read('mean_LST_h09v04')

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

if (mean_ours==mean_theirs).all():

    print 'Successfully reproduced LST_Day_1km from', name_of_8day_hdf

else:

    print 'Failed to reproduce LST_Day_1km from', name_of_8day_hdf

# compare nObs grid values and report out...

nobs_theirs = garray.array()

nobs_theirs.read('nCSD_8day_h09v04')

print 'Average of counts based on Clear_sky_days:', nobs_theirs.mean()

nobs_ours = garray.array()

nobs_ours.read('nobs_LST_h09v04')

print 'Average of our N_obs:', nobs_ours.mean()

if (nobs_ours==nobs_theirs).all():

    print 'Successfully reproduced Clear_sky_days counts from', name_of_8day_hdf

else:

    print 'Failed to reproduce Clear_sky_days counts from', 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))