crefl.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

# Copyright (c) 2013, 2014, 2015 Martin Raspaud

# Author(s):

#   Martin Raspaud <martin.raspaud@smhi.se>

# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.

# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.

# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

"""Apply reflectance correction (rayleigh effect removal)

MODIS only at the moment.

Shamelessly translated from the crefl1.7.1.c.

TODO:
 - Interpolation of heights and angle data maybe ?
 - cleanup, optimize
"""

import numpy as np
import logging

logger = logging.getLogger(__name__)

aH2O = np.array(
    [-5.60723, -5.25251, 0, 0, -6.29824, -7.70944, -3.91877, 0, 0, 0, 0, 0, 0, 0, 0, 0])
bH2O = np.array([0.820175, 0.725159, 0, 0, 0.865732, 0.966947,
                 0.745342, 0, 0, 0, 0, 0, 0, 0, 0, 0])
#/*const float aO3[Nbands]={ 0.0711,    0.00313, 0.0104,     0.0930,   0, 0, 0, 0.00244, 0.00383, 0.0225, 0.0663, 0.0836, 0.0485, 0.0395, 0.0119, 0.00263};*/
aO3 = np.array([0.0715289, 0,       0.00743232, 0.089691, 0, 0, 0, 0.001,
                0.00383, 0.0225, 0.0663, 0.0836, 0.0485, 0.0395, 0.0119, 0.00263])
#/*const float taur0[Nbands] = { 0.0507,  0.0164,  0.1915,  0.0948,  0.0036,  0.0012,  0.0004,  0.3109, 0.2375, 0.1596, 0.1131, 0.0994, 0.0446, 0.0416, 0.0286, 0.0155};*/
taur0 = np.array([0.05100, 0.01631, 0.19325, 0.09536, 0.00366, 0.00123, 0.00043,
                  0.3139, 0.2375, 0.1596, 0.1131, 0.0994, 0.0446, 0.0416, 0.0286, 0.0155])

UO3 = 0.319
UH2O = 2.93


def show(data):
    """Show the stetched data.
    """
    import Image as pil
    img = pil.fromarray(np.array((data - data.min()) * 255.0 /
                                 (data.max() - data.min()), np.uint8))
    img.show()


def fintexp1(tau):
    a = [-.57721566, 0.99999193, -0.24991055,
         0.05519968, -0.00976004, 0.00107857]
    a.reverse()
    return np.polyval(a, tau) - np.log(tau)


def fintexp3(tau):
    return (np.exp(-tau) * (1 - tau) + tau * tau * fintexp1(tau)) / 2.0


def csalbr(tau):
    return ((3 * tau - fintexp3(tau) * (4 + 2 * tau) + 2 * np.exp(-tau)) /
            (4 + 3 * tau))


def init_lut():
    vals = np.arange(4000) * 0.0001
    vals[0] = 0.0001
    sphalb = csalbr(vals)
    sphalb[0] = 0
    return sphalb


def chand(phi, muv, mus, taur):
    xfd = 0.958725775
    xbeta2 = 0.5
    musn = mus[:, :, np.newaxis]
    muvn = muv[:, :, np.newaxis]
    phios = phi + np.pi
    xcos1 = 1.0
    xcos2 = np.cos(phios)[:, :, np.newaxis]
    xcos3 = np.cos(2 * phios)[:, :, np.newaxis]
    xph1 = 1 + (3 * musn * musn - 1) * (3 * muvn * muvn - 1) * xfd / 8
    xph2 = -xfd * xbeta2 * 1.5 * musn * muvn * \
        np.sqrt(1 - musn * musn) * np.sqrt(1 - muvn * muvn)
    xph3 = xfd * xbeta2 * 0.375 * (1 - musn * musn) * (1 - muvn * muvn)
    p = np.array([np.ones_like(mus), mus + muv, muv * mus,
                  mus * mus + muv * muv, mus * mus * muv * muv])
    as01 = np.array(
        [0.33243832, 0.16285370, -0.30924818, -0.10324388, 0.11493334])
    as02 = np.array(
        [-6.777104e-02, 1.577425e-03, -1.240906e-02, 3.241678e-02, -3.503695e-02])
    as01 = as01[:, np.newaxis, np.newaxis]
    as02 = as02[:, np.newaxis, np.newaxis]
    fs01 = np.sum(p * as01, axis=0)[:, :, np.newaxis]
    fs02 = np.sum(p * as02, axis=0)[:, :, np.newaxis]
    xlntaur = np.log(taur)
    fs0 = fs01 + fs02 * xlntaur

    as1 = np.array([0.19666292, -5.439061e-02])
    as2 = np.array([0.14545937, -2.910845e-02])

    fs1 = as1[0] + xlntaur * as1[1]
    fs2 = as2[0] + xlntaur * as2[1]

    trdown = np.exp(-taur / musn)
    trup = np.exp(-taur / muvn)
    xitm1 = (1 - trdown * trup) / 4 / (musn + muvn)
    xitm2 = (1 - trdown) * (1 - trup)
    xitot1 = xph1 * (xitm1 + xitm2 * fs0)
    xitot2 = xph2 * (xitm1 + xitm2 * fs1)
    xitot3 = xph3 * (xitm1 + xitm2 * fs2)
    rhoray = xitot1 * xcos1 + xitot2 * xcos2 * 2 + xitot3 * xcos3 * 2
    return rhoray, trdown, trup

MAXAIRMASS = 18
SCALEHEIGHT = 8000


def compute_atm_variables(mus, muv, phi, height):
    sphalb0 = init_lut()
    print sphalb0

    musn = mus[:, :, np.newaxis]
    muvn = muv[:, :, np.newaxis]
    m = np.ma.masked_greater(1 / musn + 1 / muvn, MAXAIRMASS)
    psurfratio = np.exp(-height / SCALEHEIGHT)
    taur = taur0[np.newaxis, np.newaxis, :] * psurfratio[:, :, np.newaxis]
    rhoray, trdown, trup = chand(phi, muv, mus, taur)
    logger.debug("rhoray: " + str((rhoray.max(), rhoray.min())))

    sphalb = sphalb0[(taur / 0.0001 + 0.5).astype(int)]
    sphalb = np.ma.masked_greater_equal(taur, 4000 * 0.0001)
    Ttotrayu = ((2 / 3.0 + muvn) + (2 / 3.0 - muvn) * trup) / (4 / 3.0 + taur)
    Ttotrayd = ((2 / 3.0 + musn) + (2 / 3.0 - musn) * trdown) / \
        (4 / 3.0 + taur)

    tO3 = np.exp(-m * UO3 * aO3)
    tO3[:, :, aO3 == 0] = 1
    tH2O = np.exp(-np.exp(aH2O + bH2O * np.log(m * UH2O)))
    tH2O[:, :, bH2O == 0] = 1

    TtotraytH2O = np.ma.masked_invalid(Ttotrayu * Ttotrayd * tH2O)

    # tO2 = np.exp(-m *aO2)
    tO2 = 1

    tOG = np.ma.masked_invalid(tO3 * tO2)
    return sphalb, rhoray, TtotraytH2O, tOG


def correct_refl(refl, sphalb, rhoray, TtotraytH2O, tOG):
    corr_refl = (refl / tOG - rhoray) / TtotraytH2O
    corr_refl /= (1 + corr_refl * sphalb)
    return corr_refl

if __name__ == '__main__':

    from mpop.satellites import PolarFactory
    from datetime import datetime
    from mpop.utils import debug_on

    debug_on()

    t = datetime(2012, 12, 10, 10, 29, 35)

    g = PolarFactory.create_scene("terra", "", "modis", t)

    g.load(["1", "2", "3", "4"], resolution=1000)

    l = g.project("spain")

    img1 = l.image.truecolor()
    img1.save("before.png")

    height = g.hdfeos_l1b_reader.get_height()[:]
    sunz = (g.hdfeos_l1b_reader.get_sunz()[:] *
            g.hdfeos_l1b_reader.get_sunz().attributes()["scale_factor"])
    satz = g.hdfeos_l1b_reader.get_satz(
    )[:] * g.hdfeos_l1b_reader.get_satz().attributes()["scale_factor"]
    suna = g.hdfeos_l1b_reader.get_suna(
    )[:] * g.hdfeos_l1b_reader.get_suna().attributes()["scale_factor"]
    sata = g.hdfeos_l1b_reader.get_sata(
    )[:] * g.hdfeos_l1b_reader.get_sata().attributes()["scale_factor"]

    sunz = np.ma.masked_outside(sunz, -89, 89)

    phi = np.deg2rad(suna - sata)
    sphalb, rhoray, TtotraytH2O, tOG = compute_atm_variables(np.cos(np.deg2rad(sunz)),
                                                             np.cos(
                                                                 np.deg2rad(satz)),
                                                             phi,
                                                             height)

    #refl = np.dstack([g["1"].data, g["2"].data, g["3"].data, g["4"].data])[2::5, 2::5, :] / 400
    #corr_refl = (refl / tOG[:, :, :4] - rhoray[:, :, :4]) / TtotraytH2O[:, :, :4]
    #corr_refl /= (1 + corr_refl * sphalb[:, :, :4])

    from geotiepoints.interpolator import generic_modis5kmto1km
    chans = ["1", "2", "3", "4"]
    refl = np.dstack(
        [g["1"].data, g["2"].data, g["3"].data, g["4"].data]) / 400
    mus = generic_modis5kmto1km(np.cos(np.deg2rad(sunz)))[0][:, :, np.newaxis]
    refl /= mus
    print "interpolation"
    for num, chan in enumerate(chans):
        print chan
        sph, tog, rho, th2o = generic_modis5kmto1km(sphalb[:, :, num], rhoray[:, :, num],
                                                    tOG[:, :, num], TtotraytH2O[:, :, num])
        g[chan].data = correct_refl(refl[:, :, num], sph, tog, rho, th2o)

    l = g.project("spain")

    img = l.image.truecolor()

    img.save("after.png")
    img.show()