PSILARTEMP.py

#!/usr/bin/python3
"""
    Copyright (C) 2023 NIBIO <https://www.nibio.no/>. 

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU Affero 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 Affero General Public License for more details.

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

import subprocess
import os
from dotenv import load_dotenv
from datetime import datetime
import pytz
import netCDF4 as nc
import configparser
from jinja2 import Environment, FileSystemLoader

# Load deployment specific settings
load_dotenv()
# Get language stuff
config = configparser.ConfigParser()
config.read("PSILARTEMP.cfg")

year = os.getenv("YEAR")
recurring_start_date = os.getenv("RECURRING_START_DATE")
weatherdata_path = os.getenv("WEATHER_DATA_DIR")
tmp_path = "tmp/"
out_path = os.getenv("DATA_DIR")

# TODO: Put timezone out in .env
local_timezone = pytz.timezone("Europe/Oslo")

# If we are before the recurring start date, exit nicely
start_time = datetime.strptime(f"{year}-{recurring_start_date}","%Y-%m-%d")
if datetime.now() <= start_time:
    print(f"Today is before the configured start time of {start_time}. Exiting.")
    exit(0)
    
# """
# Calculate cumulated degree days above 5 degrees after 1st of April
# Remove all values before April 1st
subprocess.run(
    f"cdo -selname,TM -seldate,{year}-{recurring_start_date}T00:00:00,{year}-12-31T00:00:00 {weatherdata_path}met_1_0km_nordic-{year}.nc {tmp_path}TM_from_april.nc",
    shell=True,
)
# Subtracting 5 deg C from all cells
subprocess.run(
    f"cdo -subc,5 {tmp_path}TM_from_april.nc {tmp_path}TM_minus_5.nc", shell=True
)
# Create an .nc file with all cells containing 1 if cell in tgminus5.nc is >=0, 0 otherwise
subprocess.run(
    f"cdo -gtc,0 {tmp_path}TM_minus_5.nc {tmp_path}TM_minus_5_gtc.nc", shell=True
)

# Multiplying tgminus5.nc with tgminus5gtc.nc, so that all sub zero values are set to 0
subprocess.run(
    f"cdo -mul {tmp_path}TM_minus_5.nc {tmp_path}TM_minus_5_gtc.nc {tmp_path}TM_minus_5_nozero.nc",
    shell=True,
)
# Accumulate day degrees with the corrected values
subprocess.run(
    f"cdo -timcumsum {tmp_path}TM_minus_5_nozero.nc {tmp_path}DD_unmasked.nc",
    shell=True,
)


# Mask results using a CSV file with polygons
# Env variable MASK_FILE must be set
if os.getenv("MASK_FILE") is not None:
    mask_file = os.getenv("MASK_FILE")
    print(f"Applying mask file {mask_file} to result.nc")
    subprocess.run(
        f"cdo -maskregion,{mask_file} {tmp_path}DD_unmasked.nc {tmp_path}DD.nc",
        shell=True,
    )
else:
    os.rename(f"{tmp_path}DD_unmasked.nc", f"{tmp_path}DD.nc")

# Add the DD threshold classification => warning status
# 0: DD == 0            # Forecast not started
# 2: DD <= 260          # Flight period not started
# 3: 260 < DD <= 360    # Flight period starting
# 4: 360 < DD <= 560    # Peak flight period
# 5: 560 < DD           # 1st generation flight period ended
# (A>0)*(A<=260)*2 + (A>260)*(A<=360)*3 + (A>360)*(A<=560)*4 + (A>560)*5
subprocess.run(
    f'cdo -aexpr,"WARNING_STATUS=(TM>0)*(TM<=260)*2 + (TM>260)*(TM<=360)*3 + (TM>360)*(TM<=560)*4 + (TM>560)*5" {tmp_path}DD.nc {tmp_path}result.nc',
    shell=True,
)

# """

# Split the combined file into daily .nc files again, with YYYY-MM-DD in the filename. Convert to corresponding GeoTIFF files
# Variables that needs discrete classification, must be integers in order for mapserver to work properly (Don't ask!)
# Since we need WARNING_STATUS to be discretely classified, we need to create a separate GeoTIFF file for it
result = nc.Dataset(f"{tmp_path}result.nc", "r")
timesteps = result.variables["time"][:]
# print(timesteps)
timestep_index = 1
timestep_dates = []
for timestep in timesteps:
    timestep_date = datetime.fromtimestamp(timestep).astimezone(local_timezone)
    file_date = timestep_date.astimezone(local_timezone).strftime("%Y-%m-%d")
    # print(f"{timestep_index}: {file_date}")
    timestep_dates.append(file_date)
    # We must delete the GeoTIFF file before merging
    # """
    subprocess.run(f"rm {tmp_path}result_*{file_date}_lcc.tif", shell=True)
    subprocess.run(f"rm {out_path}result_*{file_date}.tif", shell=True)
    with open("/dev/null", "w") as devnull:
        subprocess.run(
            f'gdal_translate -ot Int16 -of GTiff  -b {timestep_index} NETCDF:"{tmp_path}result.nc":WARNING_STATUS {tmp_path}result_WARNING_STATUS_{file_date}_lcc.tif',
            shell=True,
            stdout=devnull,
        )
        subprocess.run(
            f'gdal_translate -ot Float32 -of GTiff -b {timestep_index} NETCDF:"{tmp_path}result.nc":TM {tmp_path}result_{file_date}_lcc.tif',
            shell=True,
            stdout=devnull,
        )
        # Need to reproject the files, to ensure we have the projection given in the generted mapfile. We always use EPSG:4326 for this
        subprocess.run(
            f"gdalwarp -t_srs EPSG:4326 {tmp_path}result_WARNING_STATUS_{file_date}_lcc.tif {out_path}result_WARNING_STATUS_{file_date}.tif",
            shell=True,
            stdout=devnull,
        )
        subprocess.run(
            f"gdalwarp -t_srs EPSG:4326 {tmp_path}result_{file_date}_lcc.tif {out_path}result_{file_date}.tif",
            shell=True,
            stdout=devnull,
        )
    # """
    timestep_index = timestep_index + 1

if len(timestep_dates) != len(set(timestep_dates)):
    print("ERROR: Something is wrong with your timesteps")
    # print(timestep_dates)

# print(len(timestep_dates))

# Generate mapfile
# Building data sets for language specific legends
languages = []
language_codes = config["i18n"]["languages"].split(",")
for language_code in language_codes:
    language = {"language_code": language_code}
    if ("i18n.%s" % language_code) in config:
        for keyword in config["i18n.%s" % language_code]:
            language[keyword] = config["i18n.%s" % language_code][keyword]
        languages.append(language)

# The paths should be set in a .env file
env = Environment(loader=FileSystemLoader("."))
template = env.get_template("mapfile/template.j2")
output = template.render(
    {
        "timestep_dates": timestep_dates,
        "mapserver_data_dir": os.getenv("MAPSERVER_DATA_DIR"),
        "mapserver_mapfile_dir": os.getenv("MAPSERVER_MAPFILE_DIR"),
        "mapserver_log_file": os.getenv("MAPSERVER_LOG_FILE"),
        "mapserver_image_path": os.getenv("MAPSERVER_IMAGE_PATH"),
        "mapserver_extent": os.getenv("MAPSERVER_EXTENT"),
        "languages": languages,
        "language_codes": language_codes,
    }
)
mapfile_outdir = os.getenv("MAPFILE_DIR")
with open(f"{mapfile_outdir}/PSILARTEMP.map", "w") as f:
    f.write(output)