Computing Weather Data with Python: NetCDF4 and Xarray

When I’m building a system that’s able to compute the daily Canadian Forest Fire Weather Indices, I began with the Google Earth Engine platform. It is capable of doing raster computation at a scale by harnessing the capabilities of the Google Cloud computing platform. Global datasets are updated daily on the Google Earth Engine platform.

In many cases, you can also download the files directly from the data source providers, like ECMWF with their ERA5 product or the GloH2O with their MSWX and MSWEP products. MSWX and MSWEP near real-time data is provided for free (time delay < 1 day) and ERA5T products can be accessible with a lag of 5 days. MSWX and MSWEP can be accessed through Google Drive link, while ERA5T data can be accessed through the Climate Data Store (CDS) API.

To prepare a development environment, I am using Python 3.9 that is provided by an Anaconda environment. I am using Miniconda as a Package Manager / Environment Manager.

$ conda create -n fwi-xarray python=3.9
$ conda activate fwi-xarray
$ conda install -c conda-forge jupyterlab numpy xarray[complete] netcdf4

Once I have my environment, I began writing code for Fine Fuel Moisture Code. I created a FineFuelMoistureCode class to help calculate the FFMC.

import xarray as xr
import numpy as np

class FineFuelMoistureCode:
    """
    Fine Fuel Moisture Code Calculation based on the
    Canadian Forest Fire Weather Index System using xarray
    """
    def __init__(self, temp, rhum, wind, rain, ffmc_prev):
        self.temp = temp
        self.rhum = rhum
        self.wind = wind
        self.rain = rain

        if type(ffmc_prev) == float or type(ffmc_prev) == int:
            self.ffmc_prev = rain.where(rain == ffmc_prev, \
                other = ffmc_prev)
        else:
            self.ffmc_prev = ffmc_prev

    def compute(self):
        """
        Computes the FFMC through Raining Phase 
        and Drying Phase
        
        Returns
        -------
        ffmc : xarray.DataArray
            today's FFMC
        """
        self.raining_phase()
        self.drying_phase()
        return self.ffmc
    
    def __moisture_content(self, ffmc):
        return 147.2 * (101.0 - ffmc) / (59.5 + ffmc)

    def __moisture_code(self, moisture):
        return 59.5 * (250.0 - moisture) / (147.2 + moisture)

    def __rain_normal(self, r_f, m_o):
        delta_m = (42.5 * np.exp(-100.0 / (251.0 - m_o)) * \
            (1 - np.exp(-6.93 / r_f))) * r_f
        return m_o + delta_m

    def __rain_compensation(self, r_f, m_o):    
        mr = self.__rain_normal(r_f, m_o)
        corrective = 0.0015 * (m_o - 150.0) ** 2 * r_f **0.5
        return mr + corrective

    def __no_rain(self, rain, m_o):
        return m_o + 0.0 * rain

    def raining_phase(self):
        """
        Moisture change due to rain from past 24 hours
        """
        # Moisture content before rain
        m_o = xr.apply_ufunc(self.__moisture_content, self.ffmc_prev)

        # Calculate effective rain due to canopy
        no_rain = self.rain.where(self.rain <= 0.5)
        effective_rain = xr.apply_ufunc(lambda x: x - 0.5, \
            self.rain.where(self.rain > 0.5))

        # Use corrective equation for high moisture content
        compensation = m_o.where(m_o > 150.0)
        normal = m_o.where(m_o <= 150.0)

        mo_rc = xr.apply_ufunc(self.__rain_compensation, \
            effective_rain, compensation)
        mo_r = xr.apply_ufunc(self.__rain_normal, \
            effective_rain, normal)
        mo = xr.apply_ufunc(self.__no_rain, no_rain, m_o)

        self.mr = mo_r.fillna(0) + mo_rc.fillna(0) + mo.fillna(0)
        self.mr.rename('moisture_after_rain')

    def __drying(self, mr, E_d, k_d):
        return E_d + (mr - E_d) / 10 ** k_d

    def __wetting(self, mr, E_w, k_w):
        return E_w - (E_w - mr) / 10 ** k_w

    def drying_phase(self):
        """
        Moisture change due to drying phase from noon to afternoon
        """        
        # Equilibrium moisture content for drying and wetting phase
        E_d = 0.942 * self.rhum ** 0.679 + \
            11.0 * np.exp((self.rhum - 100.0) / 10) + \
            0.18 * (21.1 - self.temp) * (1 - np.exp(-0.115 * self.rhum))
        E_w = 0.618 * self.rhum ** 0.753 + \
            10.0 * np.exp((self.rhum - 100.0) / 10) + \
            0.18 * (21.1 - self.temp) * (1 - np.exp(-0.115 * self.rhum))
        
        # Calculate the log drying/wetting rate
        k_1 = 0.424 * (1 - ((100.0 - self.rhum) / 100.0) ** 1.7) + \
            0.0694 * self.wind ** 0.5 * \
            (1 - ((100.0 - self.rhum) / 100.0) ** 8)
        k_0 = 0.424 * (1 - ((100.0 - self.rhum) / 100.0) ** 1.7) + \
            0.0694 * self.wind ** 0.5 * \
            (1 - (self.rhum/100) ** 8)
        k_d = k_0 * 0.581 * np.exp(0.0365 * self.temp)
        k_w = k_1 * 0.581 * np.exp(0.0365 * self.temp)

        # Wetting and drying conditions
        drying = self.mr.where(self.mr > E_d)
        wetting = self.mr.where(self.mr < E_w)
        no_change = self.mr.where((self.mr >= E_w) & (self.mr <= E_d))

        # Moisture content after drying
        m_d = xr.apply_ufunc(self.__drying, drying, E_d, k_d)
        m_w = xr.apply_ufunc(self.__wetting, wetting, E_w, k_w)
        m = m_d.fillna(0) + m_w.fillna(0) + no_change.fillna(0)
                                        
        # Calculate Fine Fuel Moisture Code
        self.ffmc = xr.apply_ufunc(self.__moisture_code, m)
        self.ffmc = self.ffmc.rename('fine_fuel_moisture_code')

By separating the raining phase and the drying phase, I was able to create a more readable code, at least for me. I’ll share another snippet from my final code once I have it all cleaned up and fixed. But, when I had finished coding, I found out that there’s already a submodule in xclim to calculate Fire Weather Indices!

# temp, rhum, wind, rain is a DataArray
# for daily weather measurements

ffmc_prev = 85.0
ffmc_calc = FineFuelMoistureCode(temp, rhum, wind, rain, ffmc_prev)
ffmc = ffmc_calc.compute()

However, I think that this code is much simpler for equatorial area and may be easier to implement for daily FWI calculation.

Diterbitkan oleh josefmtd

Electronics Engineer

Tinggalkan Balasan

Isikan data di bawah atau klik salah satu ikon untuk log in:

Logo WordPress.com

You are commenting using your WordPress.com account. Logout /  Ubah )

Foto Google

You are commenting using your Google account. Logout /  Ubah )

Gambar Twitter

You are commenting using your Twitter account. Logout /  Ubah )

Foto Facebook

You are commenting using your Facebook account. Logout /  Ubah )

Connecting to %s

%d blogger menyukai ini: