-
Brita Linnestad authoredBrita Linnestad authored
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()