Skip to content
Snippets Groups Projects
NAERSTADMO.py 17.04 KiB
#!/usr/bin/python3

import os, shutil
import subprocess,glob
from dotenv import load_dotenv
from datetime import datetime, timezone, timedelta
from jinja2 import Environment, FileSystemLoader

import logging
import pytz
import configparser
import netCDF4 as nc
import numpy as np 

logging.basicConfig(level=logging.INFO)

load_dotenv()
# Get language stuff
config = configparser.ConfigParser()
config.read("NAERSTADMO.cfg")


infile_path  = os.getenv("WEATHER_DATA_DIR")
outfile_path = os.getenv("DATA_DIR")
outtmp_path = "out/"
tmpfile_path = "tmp/"

bg_filename = (f'{tmpfile_path}background_data.nc')
tmp_filename = (f'{tmpfile_path}background_data_tmp.nc')
prepareWHS = (f'{tmpfile_path}prepare_WHS.nc')

utc_offset = "+02:00"
local_timezone = pytz.timezone("Europe/Oslo")

filename = (f'{tmpfile_path}weather_data.nc')

def create_dataset():

    # Find the latest file from previous run to create a start date
    last_final_date = None
    list_of_files = glob.glob(f"{outtmp_path}final_2[0-9][0-9][0-9]-[01][0-9]-[0123][0-9].nc", recursive=True)
    print(list_of_files)
    if list_of_files:
        latest_file = max(list_of_files,key=os.path.getctime)
        file_name = os.path.basename(latest_file)
        file_date = file_name[file_name.index("final_")+6:file_name.index("final_")+16]
        last_final_date=file_date
        print(last_final_date)

    if last_final_date is None or last_final_date < file_date:
        start_date = datetime.strptime(os.getenv("START_DATE"),'%Y-%m-%d')
    if last_final_date is not None:
        start_date = datetime.strptime(last_final_date,'%Y-%m-%d') - timedelta(days=4)
    print(f"Last date of final calculations is {last_final_date}. Start date = {start_date}")

    # Find the set of data to merge and use as input file based on start date
    list_weatherdata_files = glob.glob(f"{infile_path}/met_1_0km_nordic-2[0-9][0-9][0-9]-[01][0-9]-[0123][0-9].nc")
    
    for file in list_weatherdata_files:
        file_name = os.path.basename(file)
        file_date = file_name[file_name.index("nordic-")+7:file_name.index("nordic-")+17]
        
        end_date = None
        end_date = start_date + timedelta(days=5)
       
        if(file_date >= start_date.strftime('%Y-%m-%d') and file_date <= end_date.strftime('%Y-%m-%d')):
            
            if(os.path.exists(f'{tmpfile_path}weather_data.nc') != True):
                subprocess.run(f'cp {file} {tmpfile_path}weather_data.nc', shell=True)
            else:
                subprocess.run(f'mv {tmpfile_path}weather_data.nc {tmpfile_path}tmpfile.nc',shell=True)
                subprocess.run(f'cdo -f nc2 mergetime {file} {tmpfile_path}tmpfile.nc {tmpfile_path}weather_data.nc', shell=True)
                subprocess.run(f'rm {tmpfile_path}tmpfile.nc',shell=True)

    subprocess.run(f'rm {outtmp_path}final_*',shell=True)

    create_saturation()
    create_pressure()
    create_wvd()
    create_BT()
    prepare_WHS()
    create_matrix()
    create_warning_status(start_date)

def create_warning_status(start_date):
    
    risk_files = glob.glob(f"{outtmp_path}/final_*.nc")
    
    # Classifying warning status for the WHS model
    # 0  == RISK            --> Grey
    # 0  <  RISK <  1       --> Yellow
    # 1 <= RISK < 2.5       --> Orange
    # 2.5 <= RISK           --> Red

    timestep_dates = [] # Used in the mapfile template

    for file in risk_files:
      
        subprocess.run(f'cdo -aexpr,"WARNING_STATUS = RISK < 0 ? 0 : -1; WARNING_STATUS = RISK <= 1 && WARNING_STATUS == -1 ? 2 : WARNING_STATUS; WARNING_STATUS = RISK <= 2.5 && WARNING_STATUS == -1 ? 3 : WARNING_STATUS;WARNING_STATUS = RISK >= 2.5 && WARNING_STATUS == -1 ? 4 : WARNING_STATUS" -timmax {file} {tmpfile_path}result_unmasked.nc', shell=True)
        file_name = os.path.basename(file)
        file_date = file_name[file_name.index("final_")+6:file_name.index("final_")+16]
    
        timestep_dates.append(file_date)
        # 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} {tmpfile_path}result_unmasked.nc {tmpfile_path}result_{file_date}.nc', shell=True)
        else:
            os.rename(f"{tmpfile_path}result_unmasked.nc", f"{tmpfile_path}result{file_date}.nc")

        if(file_date > start_date.strftime('%Y-%m-%d')):
            
            with open("/dev/null", "w") as devnull:
                subprocess.run(f'gdal_translate -ot Int16 -of GTiff  NETCDF:"{tmpfile_path}result_{file_date}.nc":WARNING_STATUS {tmpfile_path}result_WARNING_STATUS_{file_date}.tif', shell=True)
                subprocess.run(f'gdal_translate -ot Float32 -of GTiff  NETCDF:"{tmpfile_path}result_{file_date}.nc":RISK {tmpfile_path}RISK_result_{file_date}.tif', shell=True)
                subprocess.run(f'gdal_translate -ot Float32 -of GTiff  NETCDF:"{tmpfile_path}result_{file_date}.nc":IR {tmpfile_path}IR_result_{file_date}.tif', shell=True)
        
                # 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 {tmpfile_path}result_WARNING_STATUS_{file_date}.tif {outfile_path}result_WARNING_STATUS_{file_date}.tif', shell=True, stdout=devnull)
                subprocess.run(f'gdalwarp -t_srs EPSG:4326 {tmpfile_path}RISK_result_{file_date}.tif {outfile_path}RISK_result_{file_date}.tif', shell=True, stdout=devnull)
                subprocess.run(f'gdalwarp -t_srs EPSG:4326 {tmpfile_path}IR_result_{file_date}.tif {outfile_path}IR_result_{file_date}.tif', shell=True, stdout=devnull)


    # 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}/NAERSTADMO.map", 'w') as f:
        f.write(output)

    subprocess.run(f"rm {tmpfile_path}*", shell=True)



def create_matrix():
    
    dataset=nc.Dataset(f'{tmpfile_path}background_data.nc', 'r')
    data = dataset.variables['time'][:]
    ntimesteps = np.size(data)

    for time in range(ntimesteps-5): # Need to cut the last day before the end du to WH calculations

        # For the first hour, all LastHour values are set to 0
        if time==0:
            # Some values requires a 'LastHourValue'. For the first hour of the dataset, these are set to 0 
            subprocess.run(f'cdo -seltimestep,1 -expr,"WVDLastHour=((TM >= 0)?0:0);Q0LastHour=((TM >= 0)?0:0);TSSHLastHour=((TM >= 0)?0:0);VASLastHour=((TM >= 0)?0:0);VRSLastHour=((TM >= 0)?0:0);WHSLastHour=((TM >= 0)?0:0);WDLastHour=((TM >= 0)?0:0);" -shifttime,1hour {tmpfile_path}background_data.nc {tmpfile_path}LastHour.nc', shell=True)

        
        if time>=1:
            
            # Create file with only one hour, add last hour values 
            create_HH1_HH2(time)
            
            #subprocess.run(f'cdo -seltimestep,{str(time)} background_data.nc bg_tmp.nc', shell=True) #denne virker med skifter WD;WH med 1 
            subprocess.run(f'cdo -seltimestep,{str(time+1)} {tmpfile_path}background_data.nc {tmpfile_path}bg_tmp.nc', shell=True) 
            subprocess.run(f'cdo -O -merge {tmpfile_path}bg_tmp.nc {tmpfile_path}LastHour.nc {tmpfile_path}hh1_hh2.nc {tmpfile_path}this_hour.nc ', shell=True)
            subprocess.run(f'rm {tmpfile_path}bg_tmp.nc', shell=True)


            # Create TSSF, VAS, VRS
            create_TSSH_VAS(time)
            create_VRS(time)
            subprocess.run(f'cdo -selname,TM,UM,Q0,RR,FF,BT,TSSH,VAS,SPH,WVD {tmpfile_path}this_hour.nc {tmpfile_path}background_tmp.nc ',shell=True)
            subprocess.run(f'cdo -merge {tmpfile_path}background_tmp.nc {tmpfile_path}VRS.nc {tmpfile_path}background_data_{str(time)}.nc',shell=True)
           
            # Create WHS            
            create_WHS_WH(time)
            file=nc.Dataset(f'{tmpfile_path}bg_final_{str(time)}.nc', 'r')
            timestep = file.variables["time"][:]
            date = datetime.fromtimestamp(timestep[0]).astimezone(local_timezone).strftime('%Y-%m-%d')
            hour = datetime.fromtimestamp(timestep[0]).astimezone(local_timezone).strftime('%H')
            
            subprocess.run(f'cdo -expr,"TSWH=(WHS==0)?0:sumWHTM;IR=(TSWH>=40?1:0);RISK=(IR==1?(TSWH/40)*VRS:0);VRS=(VRS==0?0:VRS);" {tmpfile_path}bg_final_{str(time)}.nc {tmpfile_path}final_{str(date)}_{str(hour)}.nc' , shell=True)
            
            if (str(hour) == '00' or time == 1):
                subprocess.run(f'cp {tmpfile_path}final_{str(date)}_{str(hour)}.nc {outtmp_path}final_{str(date)}.nc', shell=True)
            else:
                subprocess.run(f'cp {outtmp_path}final_{str(date)}.nc {tmpfile_path}tmp_{str(date)}.nc', shell=True)
                subprocess.run(f'cdo -O -mergetime {tmpfile_path}tmp_{str(date)}.nc {tmpfile_path}final_{str(date)}_{str(hour)}.nc {outtmp_path}final_{str(date)}.nc', shell=True)
                subprocess.run(f'rm {tmpfile_path}tmp_{str(date)}.nc {tmpfile_path}final_{str(date)}_{str(hour)}.nc', shell=True)
               
            # Create last hour values
            subprocess.run(f'rm {tmpfile_path}LastHour.nc', shell=True)
            subprocess.run(f'cdo -selname,WHS,TSSH,VAS,VRS,WD,Q0,WVD {tmpfile_path}bg_final_{str(time)}.nc {tmpfile_path}LastHour_tmp.nc', shell=True)
            subprocess.run(f'cdo -chname,WD,WDLastHour -chname,WHS,WHSLastHour -chname,TSSH,TSSHLastHour -chname,VAS,VASLastHour -chname,VRS,VRSLastHour -chname,WVD,WVDLastHour -chname,Q0,Q0LastHour -shifttime,1hour {tmpfile_path}LastHour_tmp.nc {tmpfile_path}LastHour.nc', shell=True)
            subprocess.run(f'rm {tmpfile_path}bg_final_{str(time)}.nc', shell=True)
            

def create_WHS_WH(time):
    
    for j in range(5):
        
        subprocess.run(f'cdo -O -chname,WVD,WVDLastHour -selname,WVD -seltimestep,{str(time+j+1)} {tmpfile_path}prepare_WHS.nc {tmpfile_path}WVD_LastHourtmp.nc', shell=True)
        subprocess.run(f'cdo -O -shifttime,1hour {tmpfile_path}WVD_LastHourtmp.nc {tmpfile_path}WVD_LastHour.nc', shell=True)

        subprocess.run(f'cdo -O -merge -seltimestep,{str(time+j+1)} {tmpfile_path}prepare_WHS.nc {tmpfile_path}WVD_LastHour.nc {tmpfile_path}WHStmp.nc', shell=True)
        os.remove(f'{tmpfile_path}WVD_LastHour.nc')

        subprocess.run(f'cdo -O -chname,BT,BT1 -selname,BT -seltimestep,{str(time+j+2)} {tmpfile_path}prepare_WHS.nc {tmpfile_path}BT1.nc', shell=True) 
        subprocess.run(f'cdo -O -chname,BT,BT2 -selname,BT -seltimestep,{str(time+j+3)} {tmpfile_path}prepare_WHS.nc {tmpfile_path}BT2.nc', shell=True)
        subprocess.run(f'cdo -O -merge {tmpfile_path}WHStmp.nc {tmpfile_path}BT1.nc {tmpfile_path}BT2.nc {tmpfile_path}WHS_background.nc', shell=True)
        os.remove(f'{tmpfile_path}WHStmp.nc')
       
        # TODO WHS is not equal to 0 r 1 for all cases
        subprocess.run(f'cdo -aexpr,"WHC=((WVD < 340) || (RR > 0.1)) ? 1 : 0; WHS=(((WVDLastHour + WVD) < 180) || (RR > 0.1) || (BT > 42) || ((BT+BT1+BT2)>120)) ? 1 : 0;" {tmpfile_path}WHS_background.nc {tmpfile_path}WHS_{str(time+j)}.nc', shell=True)
        os.remove(f'{tmpfile_path}WHS_background.nc')

        if(j==0):
            subprocess.run(f'cdo -O -selname,WDLastHour {tmpfile_path}LastHour.nc {tmpfile_path}LastHourWD.nc',shell=True)
          
        subprocess.run(f'cdo -O -merge {tmpfile_path}LastHourWD.nc {tmpfile_path}WHS_{str(time+j)}.nc {tmpfile_path}WHS_final.nc', shell=True)
        
        subprocess.run(f'cdo -aexpr,"WD=(WHC*(WDLastHour+WHS))" {tmpfile_path}WHS_final.nc {tmpfile_path}WD_{str(time+j)}.nc', shell=True)
        subprocess.run(f'cdo -aexpr,"WH=(WD > 0 ? 1 : 0)" {tmpfile_path}WD_{str(time+j)}.nc {tmpfile_path}WH_{str(time+j)}.nc', shell=True)

        if(j==0):
            subprocess.run(f'cdo -selname,WHS,WD,WH {tmpfile_path}WH_{str(time+j)}.nc {tmpfile_path}WHS.nc', shell=True)

        subprocess.run(f'cdo -O chname,WD,WDLastHour -selname,WD -shifttime,1hour {tmpfile_path}WD_{str(time+j)}.nc {tmpfile_path}LastHourWD.nc', shell=True)
        
        if(j==4):
            first_file=nc.Dataset(f'{tmpfile_path}background_data_{str(time)}.nc', 'r')
            time_var = first_file.variables["time"]
            dtime = nc.num2date(time_var[:],time_var.units)
            date = dtime[0].isoformat().split('T')

            subprocess.run(f'cdo -O mergetime {tmpfile_path}WH_*.nc {tmpfile_path}WH{str(time)}.nc', shell=True)
            
            # Må summere til 1 verdi
            subprocess.run(f'cdo -setdate,{date[0]} -settime,{date[1]} -timsum -expr,"sumWHTM = (WH == 1) ? (TM) : 0" {tmpfile_path}WH{str(time)}.nc {tmpfile_path}WHTM.nc', shell=True)
            subprocess.run(f'cdo -merge {tmpfile_path}background_data_{str(time)}.nc {tmpfile_path}WHTM.nc {tmpfile_path}WHS.nc {tmpfile_path}bg_final_{str(time)}.nc', shell=True )
            subprocess.run(f'rm {tmpfile_path}background_data_{str(time)}.nc', shell=True)
            subprocess.run(f'rm {tmpfile_path}W*.nc', shell=True)
        
      


def create_VRS(time):
    subprocess.run(f'cdo -aexpr,"RTA=(((Q0-Q0LastHour)>7)?1:0)+((WVD-WVDLastHour)>=15?1:0);IRTA=(1-(BT/80));SFRS=((1-(Q0-270)/540)/1.5);" {tmpfile_path}this_hour.nc {tmpfile_path}this_hr.nc', shell=True)
    subprocess.run(f'cdo -expr,"VRS=(((VAS*RTA*IRTA*SFRS + (0.95*VRSLastHour*SFRS))/(1+(0.1*RR)))/(1+(0.1*WHSLastHour*VRSLastHour)))" {tmpfile_path}this_hr.nc {tmpfile_path}VRS.nc', shell=True)
    subprocess.run(f'rm {tmpfile_path}this_hr.nc', shell=True)

def create_TSSH_VAS(time):
    subprocess.run(f'cdo -O -aexpr,"TSSH=((HH1==1)?HH2*(TSSHLastHour+TM):HH2*0.75*TSSHLastHour);SPH=((TSSH>87)?1:0);VAS=((VASLastHour*0.95*(1-((WVD-220)/6000))+SPH)/(1+0.3*RR))" {tmpfile_path}this_hour.nc {tmpfile_path}this_hour_tmp.nc', shell=True)
    os.rename(src='tmp/this_hour_tmp.nc', dst='tmp/this_hour.nc')

def create_HH1_HH2(time):
    # on hour basis
    subprocess.run(f'cdo -O -expr,"HH1=(WVD <= 220) ? 1 : 0; HH2=(WVD <= 520) ? 1 : 0" -seltimestep,{str(time+1)} {tmpfile_path}background_data.nc {tmpfile_path}hh1_hh2.nc', shell=True)


def create_saturation():
    
    #This is fixed for all hours and should be available in the background data
    expr="aexpr,SP=(0.61078*exp(17.269*TM/(TM+237.3)))"
    cdo_command = [
            "cdo",
            expr,
            filename,
            bg_filename,
            ]
    run_cdo(cdo_command=cdo_command)


def create_pressure():
    
    #This is fixed for all hours and should be available in the background data
    expr="aexpr,PP=UM*SP/100"
    cdo_command = [
            "cdo",
            expr,
            bg_filename,
            tmp_filename,
            ]
    run_cdo(cdo_command=cdo_command)
    os.rename(src=tmp_filename, dst=bg_filename)


def create_wvd():
    
    #This is fixed for all hours and should be available in the background data
    expr="aexpr,WVD=(SP-PP)*1000"
    cdo_command = [
            "cdo",
            expr,
            bg_filename,
            tmp_filename,
            ]
    run_cdo(cdo_command=cdo_command)
    os.rename(src=tmp_filename, dst=bg_filename)


def create_BT():
    
    # BT is not available in the dataset and need to be calculted and added to background_data.nc

    expr="aexpr,BT=((RR > 0 || (((100-UM)/100)*6.112*exp(17.67*TM/(TM+243.5))) < 2)) ? 60 : 0"
    cdo_command = [
            "cdo",
            expr,
            bg_filename,
            tmp_filename,
            ]
    run_cdo(cdo_command=cdo_command)
    os.rename(src=tmp_filename, dst=bg_filename)


def prepare_WHS():
    #system("cdo selname,RR,BT,WVD background_data.nc prepareWHS.nc")
    my_variables = ['TM', 'RR', 'BT', 'WVD']
    variable_list=','.join(my_variables)
    cdo_command = [
            "cdo",
            "selname,"+variable_list,
            bg_filename,
            prepareWHS,
            ]
    run_cdo(cdo_command=cdo_command)

    
    
def run_cdo(cdo_command):
    try:
        print(cdo_command)
        subprocess.run(cdo_command, check=True)
        logging.info(f"CDO command {cdo_command[1]} executed successfully.")
    except subprocess.CalledProcessError as e:
        logging.error(f"CDO command {cdo_command[1]} failed:", e)
        quit()

create_dataset()