ssr_algore.py

#!/usr/bin/env python
############################################################################
#
# MODULE:       ssr_algore.py
# AUTHOR:       Collin Bode, UC Berkeley
#
# PURPOSE:
#       Al Gore Rhythm combines solar radiation model with Light Penetration
#       Index (LPI). Merges all solar radiation runs into a single estimate
#       of Total Solar Radiation in watt-hours per meter squared per day.
#
# Modified:     Collin Bode, October, 2012
#               Migrated to unified parameter set.
#
# COPYRIGHT:    (c) 2011 Collin Bode
#               This program is free software under the GNU General Public
#               License (>=v2). Read the file COPYING that comes with GRASS
#               for details.
#
#############################################################################

import os
import sys
import traceback
import datetime as dt
import numpy as np

from ssr_params import *
from ssr_utilities import (printout, ensure_dir, raster_path, raster_exists,
                            read_raster, save_raster, setup_workspace)


def main():
    setup_workspace()

    tlog = dt.datetime.strftime(dt.datetime.now(), "%Y-%m-%d_h%Hm%M")
    log_path = os.path.join(workspace, dir_logs, 'ssr_' + tlog + '_algore.log')
    lf = open(log_path, 'a')

    ow = int(algore_run - 1)  # 0 = no overwrite, 1 = overwrite

    printout('---------------------------------------', lf)
    printout('-- ALGORITHM FOR CLEAR SKY RADIATION --', lf)
    printout('          LPI year: ' + year, lf)
    printout('          LPI pref: ' + lpipref, lf)
    printout('            region: ' + bregion, lf)
    printout('        sun subdir: ' + dir_sun, lf)
    printout(' SSR output subdir: ' + dir_ssr, lf)
    printout('    max veg height: ' + maxheight, lf)
    printout('    Algorithm code: ' + algore, lf)
    printout('keep intermediates: ' + str(keeptemp), lf)
    printout('   overwrite files: ' + str(ow), lf)
    printout('---------------------------------------', lf)

    ssr_out_dir = os.path.join(workspace, dir_ssr)
    ensure_dir(ssr_out_dir)

    r1start = dt.datetime.now()
    printout("Starting Al Gore Rhythm at " + str(r1start), lf)

    # Vegetation height raster (used for open-canopy masking)
    veg_path = raster_path(vegheight, dir_sun)
    veg_data, veg_profile = read_raster(veg_path)
    maxheight_f = float(maxheight)

    for doyn in range(5, 366, 7):
        doy = str(doyn).zfill(3)
        month = dt.datetime.strftime(
            dt.datetime(2011, 1, 1) + dt.timedelta(doyn - 1), "%m")
        printout("Processing DOY " + doy + " (month " + month + ")", lf)

        # Input solar radiation rasters
        sundem = bregion + C + 'mdem'
        suncan = bregion + C + 'mcan'
        dembeam_path = raster_path(sundem + doy + 'beam', dir_sun)
        demdiff_path = raster_path(sundem + doy + 'diff', dir_sun)
        canbeam_path = raster_path(suncan + doy + 'beam', dir_sun)
        candiff_path = raster_path(suncan + doy + 'diff', dir_sun)

        # LPI raster for this month
        lpi_month = lpipref + 'm' + month if not lpivsjune else lpipref + 'm06'
        lpi_path = raster_path(lpi_month, dir_lpi)

        # Output raster names
        lpipart = C + 'm' + year + 's' + boxsize + 'm' + algore
        if lpivsjune:
            lpipart = C + 'm' + year + 's' + boxsize + 'mjune' + algore
        ssr_name = 'ssr_' + lpipart + doy
        ssr_path = raster_path(ssr_name, dir_ssr)

        if not ow and os.path.exists(ssr_path):
            printout("Skipping (exists): " + ssr_name, lf)
            continue

        # Check required inputs exist
        missing = [p for p in [dembeam_path, demdiff_path, canbeam_path,
                                candiff_path, lpi_path] if not os.path.exists(p)]
        if missing:
            printout("SKIPPING DOY " + doy + " - missing inputs: " +
                     str([os.path.basename(p) for p in missing]), lf)
            continue

        dembeam, profile = read_raster(dembeam_path)
        demdiff, _ = read_raster(demdiff_path)
        canbeam, _ = read_raster(canbeam_path)
        candiff, _ = read_raster(candiff_path)
        lpi, _ = read_raster(lpi_path)

        ###################################################################
        # 1. SUBCANOPY: merge LPI with bare-earth solar radiation
        printout("DOY " + doy + " 1. merging lpi+dem using: " + algore, lf)
        if algore == 'cl':
            lpibeam = dembeam * lpi
            lpidiff = 0.94 * lpi * demdiff
        elif algore == 'cn':
            lpibeam = 1.428 * dembeam * lpi
            lpidiff = 0.94 * (1.428 * lpi) * demdiff
        elif algore == 'gn':
            lpibeam = (1.428 * lpi) * dembeam
            lpidiff = 0.01719 + 1.024 * (1.428 * lpi) * demdiff
        elif algore == 'gl':
            lpibeam = lpi * dembeam
            lpidiff = 0.01719 + 1.024 * lpi * demdiff
        else:  # 'pl' power law
            lpibeam = 1.2567 * dembeam * lpi
            lpidiff = 1.1224 * (lpi ** 0.3157) * demdiff

        subcanopy = lpibeam + lpidiff

        ###################################################################
        # 2. OPEN CANOPY: canopy solar radiation, mask tall vegetation
        printout("DOY " + doy + " 2. masking areas under tall trees", lf)
        canglob = canbeam + candiff
        opencanopy = np.where(veg_data < maxheight_f, canglob, -88.0)

        ###################################################################
        # 3. Merge: keep whichever is higher (open or subcanopy)
        printout("DOY " + doy + " 3. merging -> " + ssr_name, lf)
        ssr_data = np.where(opencanopy > subcanopy, opencanopy, subcanopy)
        save_raster(ssr_data.astype(np.float32), profile, ssr_path)

        # 4. Save intermediate rasters if keeptemp
        if keeptemp:
            lpipart_name = 'subcanopy_' + lpipart + doy
            save_raster(subcanopy.astype(np.float32), profile,
                        raster_path(lpipart_name, dir_ssr))
            opencanopy_name = 'opencanopy_' + lpipart + doy
            save_raster(opencanopy.astype(np.float32), profile,
                        raster_path(opencanopy_name, dir_ssr))

    r1end = dt.datetime.now()
    printout("DONE! Al Gore Rhythm finished at " + str(r1end), lf)
    printout("Total time: " + str(r1end - r1start), lf)
    printout("--------------------------------------", lf)
    lf.close()
    sys.exit("FINISHED.")


if __name__ == "__main__":
    try:
        main()
    except Exception:
        traceback.print_exc()
        sys.exit(1)