diff --git a/COMMAND.back b/COMMAND.back new file mode 100644 index 0000000..e69de29 diff --git a/RELEASES.back b/RELEASES.back new file mode 100644 index 0000000..34492a4 --- /dev/null +++ b/RELEASES.back @@ -0,0 +1,30 @@ +*************************************************************************************************************** +* * +* * +* * +* Input file for the Lagrangian particle dispersion model FLEXPART * +* Please select your options * +* * +* * +* * +*************************************************************************************************************** +&RELEASES_CTRL + NSPEC = 1, ! Total number of species + SPECNUM_REL= 16, ! Species numbers in directory SPECIES + / +&RELEASE +IDATE1 = start_date, +ITIME1 = start_time, +IDATE2 = end_date, +ITIME2 = end_time, +LON1 = lon_1, +LON2 = 30.10, +LAT1 = lat_1, +LAT2 = 50.90, +Z1 = 0.00, +Z2 = 200.00, +ZKIND = 1, +MASS = 1.00e-04, +PARTS = 10000, +COMMENT = release_comment, +/ diff --git a/csv_parser.py b/csv_parser.py new file mode 100644 index 0000000..3a83ab6 --- /dev/null +++ b/csv_parser.py @@ -0,0 +1,46 @@ +import logging,os,csv + +logging.basicConfig(filename="main.log", level=logging.INFO, + format="%(asctime)s %(message)s") + +path_to_data_files = '/data/calculations/' +basename = os.path.basename(os.getcwd()) + +with open('measurem.csv', newline='') as csvfile: + csv_reader = csv.reader(csvfile, delimiter=';') + csv_header = next(csv_reader) + for row in csv_reader: + calc_folder_id = basename + '_' + row[1] + path_to_calc_dir = path_to_data_files + calc_folder_id + lat = row[5] + lon = row[6] + start_date = row[7].replace('.', '') + start_time = row[8].replace(':', '') + end_date = row[9].replace('.', '') + end_time = row[10].replace(':', '') + comment = "RELEASE " + row[1] + print('comment ', comment) + print('comment ', calc_folder_id) + print('path_to_data_files', path_to_calc_dir) + + if not os.path.exists(path_to_calc_dir): + os.makedirs(path_to_calc_dir) + logging.info('Folder ' + calc_folder_id + ' created') + logging.info('Parsing ' + comment + ' file') + + # Read in the default RELEASES file + with open('/data/opt/RELEASES.back', 'r') as file: + filedata = file.read() + + # Replace the target string + filedata = filedata.replace('start_date', start_date) + filedata = filedata.replace('start_time', start_time) + filedata = filedata.replace('end_date', end_date) + filedata = filedata.replace('end_time', end_time) + filedata = filedata.replace('lat_1', lat) + filedata = filedata.replace('lon_1', lon) + filedata = filedata.replace('release_comment', comment) + + # Write the file out to the calc dir + with open(path_to_calc_dir+'/'+'RELEASES', 'w') as file: + file.write(filedata) diff --git a/parse_csv_files.py b/parse_csv_files.py new file mode 100644 index 0000000..b1a371e --- /dev/null +++ b/parse_csv_files.py @@ -0,0 +1,81 @@ +''' +description: Download weather forecasts from GFS for flexpart +license: APACHE 2.0 +author: Synkevych Roman, (synkevych.roman@gmail.com) +https://www.ncei.noaa.gov/data/global-forecast-system/access/ +grid-003-1.0-degree/forecast +/202205 +/20220511 +/gfs_3_20220511_0000_000.grb2 + +//www.ncei.noaa.gov/data/global-forecast-system/access/ +grid-003-1.0-degree/analysis/202004/20200418/gfs_3_20200418_0000_000.grb2 +''' + +from pygfs.utils import * +import shutil +import os +import requests +path_to_csv_file='' +path_to_default_files='' + +number_of_calculations=0 +calc_folder_name='' + +# Create simlink for simflex and flexpart + +# provide start_data and end_date for download forecast +# and create AVAILABLE +# this script should test if file available before downloading + +# create a series_id.log that saves info about each calculations + + +class gfs: + site_url = "https://www.ncei.noaa.gov/data/global-forecast-system/access/" + + def download(self, date, fcsthours, resolution): + self.date = date + self.fcsthours = fcsthours + self.resolution = resolution + self.check_fcsthours() + self.build_filelist() + self.download_files() + + def build_filelist(self): + self.filelist = [self.get_filename(fcsthour) + for fcsthour in self.fcsthours] + + def get_filename(self, fcsthour): + int_, dec = (str(float(self.resolution))).split('.') + yr = str(self.date.year).zfill(4) + mnth = str(self.date.month).zfill(2) + day = str(self.date.day).zfill(2) + baseurl = 'https://rda.ucar.edu/data/ds084.1/' + fpath = os.path.join(yr, yr + mnth + day, 'gfs.' + + int_ + 'p' + dec.zfill(2) + '.' + + yr + mnth + day + '00.f' + + str(fcsthour).zfill(3) + '.grib2') + return (os.path.join(baseurl, fpath)) + + def download_files(self): + ''' + Download the actual files + ''' + [self.download_file(fl) for fl in self.filelist] + + def download_file(self, url): + ''' + Stream download url to localfile + ''' + local_filename = url.split('/')[-1] + with self.session.get(url, stream=True) as r: + with open(local_filename, 'wb') as f: + shutil.copyfileobj(r.raw, f) + #check_gribfile(local_filename) + return local_filename + +def parse_releases(): + # csv file name + # row line + # releases template file name and location diff --git a/template/1.txt b/template/1.txt new file mode 100644 index 0000000..6c93ff9 --- /dev/null +++ b/template/1.txt @@ -0,0 +1,31 @@ +id_calc use_no id_measu id_station station country s_lat s_lng id_nuclide name_nuclide date_start time_start date_end time_end val sigma_or_ldl backgr +1 1 1 1 node_291_137 Ukraine 50.85 30.05 3 Ru-106 2020-04-18 05:00:00 2020-04-18 05:00:10 0.00524513 8.3922E-4 6.0E-6 +1 1 3 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-18 07:00:00 2020-04-18 07:00:10 0.01818158 0.00290905 6.0E-6 +1 1 4 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-18 08:00:00 2020-04-18 08:00:10 0.03814321 0.00610291 6.0E-6 +1 1 5 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-18 09:00:00 2020-04-18 09:00:10 0.04771746 0.00763479 6.0E-6 +1 1 6 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-18 10:00:00 2020-04-18 10:00:10 0.03077517 0.00492403 6.0E-6 +1 1 7 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-18 11:00:00 2020-04-18 11:00:10 0.00891589 0.00142654 6.0E-6 +1 1 8 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-18 12:00:00 2020-04-18 12:00:10 0.00124652 1.9944E-4 6.0E-6 +1 1 9 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-18 13:00:00 2020-04-18 13:00:10 1.3097E-4 2.1E-5 6.0E-6 +1 1 10 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-18 14:00:00 2020-04-18 14:00:10 9.6E-7 1.5E-7 6.0E-6 +1 1 11 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-18 15:00:00 2020-04-18 15:00:10 0.0 0.0 6.0E-6 +1 1 12 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-18 16:00:00 2020-04-18 16:00:10 1.1426E-4 1.83E-5 6.0E-6 +1 1 13 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-18 17:00:00 2020-04-18 17:00:10 0.00414009 6.6241E-4 6.0E-6 +1 1 14 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-18 18:00:00 2020-04-18 18:00:10 0.02838107 0.00454097 6.0E-6 +1 1 15 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-18 19:00:00 2020-04-18 19:00:10 0.0516992 0.00827187 6.0E-6 +1 1 16 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-18 20:00:00 2020-04-18 20:00:10 0.03917277 0.00626764 6.0E-6 +1 1 17 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-18 21:00:00 2020-04-18 21:00:10 0.01503785 0.00240606 6.0E-6 +1 1 18 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-18 22:00:00 2020-04-18 22:00:10 0.00435961 6.9754E-4 6.0E-6 +1 1 19 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-18 23:00:00 2020-04-18 23:00:10 0.0026768 4.2829E-4 6.0E-6 +1 1 20 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-19 00:00:00 2020-04-19 00:00:10 0.01201856 0.00192297 6.0E-6 +1 1 21 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-19 01:00:00 2020-04-19 01:00:10 0.0566713 0.00906741 6.0E-6 +1 1 22 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-19 02:00:00 2020-04-19 02:00:10 0.07849871 0.01255979 6.0E-6 +1 1 23 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-19 03:00:00 2020-04-19 03:00:10 0.06059366 0.00969499 6.0E-6 +1 1 24 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-19 04:00:00 2020-04-19 04:00:10 0.05526348 0.00884216 6.0E-6 +1 1 25 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-19 05:00:00 2020-04-19 05:00:10 0.02873416 0.00459747 6.0E-6 +1 1 26 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-19 06:00:00 2020-04-19 06:00:10 0.00174476 2.7916E-4 6.0E-6 +1 1 27 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-19 07:00:00 2020-04-19 07:00:10 0.0 1.0E-6 6.0E-6 +1 1 28 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-19 08:00:00 2020-04-19 08:00:10 0.0 1.0E-6 6.0E-6 +1 1 29 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-19 09:00:00 2020-04-19 09:00:10 0.0 1.0E-6 6.0E-6 +1 1 30 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-19 11:00:00 2020-04-19 11:00:10 0.0 1.0E-6 6.0E-6 +1 1 31 1 node_291_138 Ukraine 50.85 30.05 3 Ru-106 2020-04-19 12:00:00 2020-04-19 12:00:10 0.0 1.0E-6 6.0E-6 diff --git a/template/1.xml b/template/1.xml new file mode 100644 index 0000000..5f1beed --- /dev/null +++ b/template/1.xml @@ -0,0 +1,17 @@ + + + 1 + 1 + 0 + 37.85 + 0.05 + 400.0 + 322.0 + 0.1 + 0.1 + 2020-04-18 00:00:00 + 2020-04-21 00:00:00 + 0 + 34 + 1 + diff --git a/template/calc_releases.sh b/template/calc_releases.sh new file mode 100644 index 0000000..dfd004b --- /dev/null +++ b/template/calc_releases.sh @@ -0,0 +1,25 @@ +#!/bin/bash + +sed_length=15 +startline=15 +path=$(pwd) +for i in `seq 1 31`; do + # replace RELEASES file with begining of the file + cp options/RELEASES.tmp options/RELEASES + start=$(( startline * i )) + if [ "$i" -eq 2 ] ; then + start=$((start+1)) + elif (( $i > 2 )) ; then + start=$((start+i-1)) + fi + end=$(( start + sed_length )) + # add to RELEASES one RELEASE from + sed -n "${start},${end}p" options/RELEASES.back >> options/RELEASES + ./FLEXPART>>allout.txt 2>&1 + mv output/grid_time_20200421000000.nc output/grid_time_20200421000000_$i.nc + if [ "$i" -eq 1 ] ; then + echo "#id_obs; path_to_file; pointspec_ind;" >> table_srs_paths.txt + else + echo "$i; $(pwd)/output/grid_time_20200421000000_$i.nc; 1">>table_srs_paths.txt + fi +done diff --git a/template/download_grib.py b/template/download_grib.py new file mode 100644 index 0000000..dcdfd25 --- /dev/null +++ b/template/download_grib.py @@ -0,0 +1,88 @@ +from datetime import datetime +from datetime import timedelta, date +import os, urllib.request + +# 2 month equal to 44 Gb / calc speed and time needed to downloads this data + +HHMMSS = ['030000', '060000', '090000', '120000', + '150000', '180000', '210000', '000000'] +FILE_HOURS = ['0000', '0000', '0600', '0600', '1200', '1200', '1800', '1800'] +FILE_SUFFIX = ['003', '006'] +NCEI_URL = "https://www.ncei.noaa.gov/data/global-forecast-system/access/historical/" +FILE_PREFIX = "gfs_4_" +DOMAIN = NCEI_URL + "forecast/grid-004-0.5-degree/" +basename = os.getcwd() +DATA_FOLDER = '/data/grib_data/' +available_template_header = """XXXXXX EMPTY LINES XXXXXXXXX +XXXXXX EMPTY LINES XXXXXXXX +YYYYMMDD HHMMSS name of the file(up to 80 characters) +""" + +def write_to_file(file_name, contents, mode='w'): + file = open(basename + '/' + file_name, mode) + file.write(contents) + file.close() + +def create_folder(directory=None): + if not os.path.exists(directory): + os.makedirs(directory) + +def parse_available_file(date=None, file_name=None): + if date is not None and file_name is not None: + available_template_body = """{yyyymmdd} {hhmmss} {file_name} ON DISC +""".format(yyyymmdd=date.strftime('%Y%m%d'), + hhmmss=date.strftime('%H%M%S'), + file_name=file_name) + write_to_file('AVAILABLE', available_template_body, 'a') + +def download_grib(date_start=None, date_end=None): + if date_start and date_end is not None: + # dates should ends with hours that divides to 3 + start_date = date_start - timedelta(hours = date_start.hour % 3) + end_date = date_end + + # download last dataset using end date + if end_date.hour % 3 == 1: + end_date = date_end + timedelta(hours=2) + elif end_date.hour% 3 == 2: + end_date = date_end + timedelta(hours=1) + else: + end_date = date_end + timedelta(hours=3) + + days, seconds = ( + end_date - start_date).days, (end_date - start_date).seconds + hours = (days * 24 + seconds / 3600) // 8 + print('amount of datasets: ', hours) + if hours <= 0: + print('Error, invalid START or END date') + return + + create_folder(DATA_FOLDER) + write_to_file('AVAILABLE', available_template_header) + + end_forecast_date = start_forecast_date = start_date + while(end_forecast_date < end_date): + forecast_suffix = '' + if end_forecast_date.hour % 6 == 0: + # start_forecast_date - in Available + start_forecast_date = end_forecast_date - timedelta(hours = 6) + forecast_suffix = FILE_SUFFIX[1] + elif end_forecast_date.hour % 6 == 3: + start_forecast_date = end_forecast_date - timedelta(hours = 3) + forecast_suffix = FILE_SUFFIX[0] + file_name = FILE_PREFIX + \ + start_forecast_date.strftime('%Y%m%d_%H%M_') + forecast_suffix + ".grb2" + path_to_file = os.path.join(DATA_FOLDER + '/', file_name) + # test if file exist + print('file_name', file_name) + if os.path.isfile(path_to_file): + print("File", file_name, " exist.") + parse_available_file(end_forecast_date, file_name) + end_forecast_date = start_forecast_date = end_forecast_date + timedelta(hours=3) + continue + else: + URL = DOMAIN + start_forecast_date.strftime('%Y%m/%Y%m%d/') + file_name + urllib.request.urlretrieve(URL, path_to_file) + parse_available_file(end_forecast_date, file_name) + end_forecast_date = start_forecast_date = end_forecast_date + \ + timedelta(hours=3) diff --git a/template/parse_input.py b/template/parse_input.py new file mode 100644 index 0000000..ec56c61 --- /dev/null +++ b/template/parse_input.py @@ -0,0 +1,222 @@ +import os +import logging +import csv +import xml.etree.ElementTree as ET +from subprocess import call + +logging.basicConfig(filename="parsing.log", level=logging.INFO, + format="%(asctime)s %(message)s") + +basename = os.getcwd() +template_header = """*************************************************************************************************************** +* * +* Input file for the Lagrangian particle dispersion model FLEXPART * +* Please select your options * +* * +*************************************************************************************************************** +""" + + +def write_to_file(file_name, contents, mode='w'): + file = open(basename + file_name, mode) + file.write(contents) + file.close() + + +def parse_command_file(start_date, start_time, end_date, end_time): + command_body = """&COMMAND + LDIRECT= -1, ! Simulation direction in time ; 1 (forward) or -1 (backward) + IBDATE= {date_1}, ! Start date of the simulation ; YYYYMMDD: YYYY=year, MM=month, DD=day + IBTIME= {time_1}, ! Start time of the simulation ; HHMISS: HH=hours, MI=min, SS=sec; UTC + IEDATE= {date_2}, ! End date of the simulation ; same format as IBDATE + IETIME= {time_2}, ! End time of the simulation ; same format as IBTIME + LOUTSTEP= 3600, ! Interval of model output; average concentrations calculated every LOUTSTEP (s) + LOUTAVER= 3600, ! Interval of output averaging (s) + LOUTSAMPLE= 3600, ! Interval of output sampling (s), higher stat. accuracy with shorter intervals + ITSPLIT= 99999999, ! Interval of particle splitting (s) + LSYNCTIME= 60, ! All processes are synchronized to this time interval (s) + CTL= -5.0000000, ! CTL>1, ABL time step = (Lagrangian timescale (TL))/CTL, uses LSYNCTIME if CTL<0 + IFINE= 4, ! Reduction for time step in vertical transport, used only if CTL>1 + IOUT= 9, ! Output type: [1]mass 2]pptv 3]1&2 4]plume 5]1&4, +8 for NetCDF output + IPOUT= 0, ! Particle position output: 0]no 1]every output 2]only at end 3]time averaged + LSUBGRID= 0, ! Increase of ABL heights due to sub-grid scale orographic variations;[0]off 1]on + LCONVECTION= 1, ! Switch for convection parameterization;0]off [1]on + LAGESPECTRA= 0, ! Switch for calculation of age spectra (needs AGECLASSES);[0]off 1]on + IPIN= 0, ! Warm start from particle dump (needs previous partposit_end file); [0]no 1]yes + IOUTPUTFOREACHRELEASE= 1, ! Separate output fields for each location in the RELEASE file; [0]no 1]yes + IFLUX= 0, ! Output of mass fluxes through output grid box boundaries + MDOMAINFILL= 0, ! Switch for domain-filling, if limited-area particles generated at boundary + IND_SOURCE= 1, ! Unit to be used at the source ; [1]mass 2]mass mixing ratio + IND_RECEPTOR= 1, ! Unit to be used at the receptor; [1]mass 2]mass mixing ratio 3]wet depo. 4]dry depo. + MQUASILAG= 0, ! Quasi-Lagrangian mode to track individual numbered particles + NESTED_OUTPUT= 0, ! Output also for a nested domain + LINIT_COND= 0, ! Output sensitivity to initial conditions (bkw mode only) [0]off 1]conc 2]mmr + SURF_ONLY= 0, ! Output only for the lowest model layer, used w/ LINIT_COND=1 or 2 + CBLFLAG= 0, ! Skewed, not Gaussian turbulence in the convective ABL, need large CTL and IFINE + OHFIELDS_PATH= "../../flexin/", ! Default path for OH file + / +""".format(date_1=start_date, time_1=start_time, date_2=end_date, time_2=end_time) + write_to_file('/options/COMMAND', template_header + command_body) + + +def parse_outgrid_file(args): + outgrid_template = """!******************************************************************************* +! * +! Input file for the Lagrangian particle dispersion model FLEXPART * +! Please specify your output grid * +! * +! OUTLON0 = GEOGRAPHYICAL LONGITUDE OF LOWER LEFT CORNER OF OUTPUT GRID * +! OUTLAT0 = GEOGRAPHYICAL LATITUDE OF LOWER LEFT CORNER OF OUTPUT GRID * +! NUMXGRID = NUMBER OF GRID POINTS IN X DIRECTION (= No. of cells + 1) * +! NUMYGRID = NUMBER OF GRID POINTS IN Y DIRECTION (= No. of cells + 1) * +! DXOUT = GRID DISTANCE IN X DIRECTION * +! DYOUN = GRID DISTANCE IN Y DIRECTION * +! OUTHEIGHTS = HEIGHT OF LEVELS (UPPER BOUNDARY) * +!******************************************************************************* +&OUTGRID + OUTLON0= {out_lon}, + OUTLAT0= {out_lat}, + NUMXGRID= {x_grid}, + NUMYGRID= {y_grid}, + DXOUT= {dx_out}, + DYOUT= {dy_out}, + OUTHEIGHTS= 200.0, + / +""".format(out_lon=args['out_longitude'], + out_lat=args['out_latitude'], + x_grid=args['num_x_grid'].split('.')[0], + y_grid=args['num_y_grid'].split('.')[0], + dx_out=args['dx_out'], + dy_out=args['dy_out']) + write_to_file('/options/OUTGRID', outgrid_template) + + +def parse_simflex_input(end_date, end_time): + dir_name = 'simflex' + if not os.path.exists(dir_name): + os.makedirs(dir_name) + simflex_paths = """#id_obs; path_to_file; srs_ind; +{obs_id}; {path_to_file}; {srs_id} +""".format(obs_id=1, path_to_file=basename+'/grid_time_' + end_date + end_time + '.nc', srs_id=1) + write_to_file('/' + dir_name + '/table_srs_paths.txt', simflex_paths) + + +def get_xml_data(): + def remove_hyphens_and_colon(str): + # example 2020-04-18 00:00:00 + return str.replace('-', '').replace(':', '') + + xml_tree = ET.parse(os.path.basename(basename) + '.xml') + xml_root = xml_tree.getroot() + + start_date_time_str = xml_root.find('imin').text + end_date_time_str = xml_root.find('imax').text + start_date, start_time = remove_hyphens_and_colon( + start_date_time_str).split() + end_date, end_time = remove_hyphens_and_colon(end_date_time_str).split() + args = { + 'out_longitude': xml_root.find('se_lon').text, + 'out_latitude': xml_root.find('se_lat').text, + 'num_x_grid': xml_root.find('nx').text, + 'num_y_grid': xml_root.find('ny').text, + 'dx_out': xml_root.find('dlat').text, + 'dy_out': xml_root.find('dlon').text + } + + parse_command_file(start_date, start_time, end_date, end_time) + logging.info('Parsing COMMAND file compleated.') + parse_outgrid_file(args) + logging.info('Parsing OUTGRID file compleated.') + parse_simflex_input(end_date, end_time) + logging.info('Parsing table_srs_paths.txt file compleated.') + + # Download grib_data using received dates + rc = call("./download_grib.sh" + + ' ' + start_date + + ' ' + end_date, shell=True) + + +def parse_release_file(args): + SPECIES_BY_ID = {"O3": '002', "NO": '003', "NO2": '004', + "HNO3": '005', "HNO2": '006', "H2O2": '007', + "NO2": '008', "HCHO": '009', "PAN": '010', + "NH3": '011', "SO4-aero": '012', "NO3-aero": '013', + "I2-131": '014', "I-131": '015', "Cs-137": '016', + "Y-91": '017', "Ru-106": '018', "Kr-85": '019', + "Sr-90": '020', "Xe-133": '021', "CO": '022', + "SO2": '023', "AIRTRACER": '024', "AERO-TRACE": '025', + "CH4": '026', "C2H6": '027', "C3H8": '028', + "PCB28": '031', "G-HCH": '034', "BC": '040'} + species_id = int(SPECIES_BY_ID.get(args['species_name'])) + release_header = """&RELEASES_CTRL + NSPEC = 1, ! Total number of species + SPECNUM_REL= {sn}, ! Species numbers in directory SPECIES + / +""".format(sn=species_id) + + release_body = """&RELEASE + IDATE1 = {date_1}, + ITIME1 = {time_1}, + IDATE2 = {date_2}, + ITIME2 = {time_2}, + LON1 = {lon_1}, + LON2 = {lon_2}, + LAT1 = {lat_1}, + LAT2 = {lat_2}, + Z1 = 0.00, + Z2 = 200.00, + ZKIND = 1, + MASS = {mass}, + PARTS = 10000, ! >= 10000, <=600000, + COMMENT = "{comment}", + / +""".format(date_1=args['start_date'], + time_1=args['start_time'], + date_2=args['end_date'], + time_2=args['end_time'], + lon_1=args['longitude_1'], + lon_2=args['longitude_2'], + lat_1=args['latitude_1'], + lat_2=args['latitude_2'], + mass=args['mass'], + comment=args['comment']) + + if args['id'] == '1': + write_to_file('/options/RELEASES', template_header + + release_header + release_body) + else: + write_to_file('/options/RELEASES', release_body, 'a') + + +def get_data_from_txt_file(): + with open(os.path.basename(basename) + '.txt', newline='') as csvfile: + csv_reader = csv.reader(csvfile, delimiter='\t') + csv_header = next(csv_reader) + for row in csv_reader: + # calc_id(0), use(1), m_id(2), s_id(3), station(4), country(5), s_lat(6),s_lng(7), + # id_nuclide(8), name_nuclide(9), date_start(10), time_start(11), date_end(12), time_end(13), val(14), sigma + measurement_id = row[2] + latitude_1 = float(row[6]) + longitude_1 = float(row[7]) + species_mass = "{:e}".format(int(row[14])) + args = { + 'id': measurement_id, + 'latitude_1': latitude_1, + 'longitude_1': longitude_1, + 'latitude_2': round(latitude_1+0.10, 2), + 'longitude_2': round(longitude_1+0.10, 2), + 'species_name': row[9], + 'start_date': row[10].replace('-', ''), + 'start_time': row[11].replace(':', ''), + 'end_date': row[12].replace('-', ''), + 'end_time': row[13].replace(':', ''), + 'mass': species_mass, + 'comment': "RELEASE " + measurement_id + } + + parse_release_file(args) + logging.info('Parsing RELEASE file compleated.') + + +get_xml_data() +get_data_from_txt_file() diff --git a/template/parser.py b/template/parser.py new file mode 100644 index 0000000..6bf8b43 --- /dev/null +++ b/template/parser.py @@ -0,0 +1,302 @@ +#!/usr/bin/env python3 + +import os +import logging +import csv +import xml.etree.ElementTree as ET +from subprocess import run +from download_grib import * +from datetime import datetime, timedelta + +logging.basicConfig(filename="parsing.log", level=logging.INFO, + format="%(asctime)s %(message)s") + +basename = os.getcwd() +simflex_dir_name = basename + '/simflex' +template_header = """*************************************************************************************************************** +* * +* Input file for the Lagrangian particle dispersion model FLEXPART * +* Please select your options * +* * +*************************************************************************************************************** +""" +user_params = {} +releases_params = [] + +if not os.path.exists('simflex'): + os.makedirs('simflex') + +if not os.path.exists('output'): + os.makedirs('output') + +def write_to_file(file_name, contents, mode='w'): + file = open(basename + file_name, mode) + file.write(contents) + file.close() + + +def get_xml_params(): + xml_tree = ET.parse(os.path.basename(basename) + '.xml') + xml_root = xml_tree.getroot() + start_date_time_str = xml_root.find('imin').text + end_date_time_str = xml_root.find('imax').text + start_date_time = datetime.strptime(start_date_time_str, '%Y-%m-%d %H:%M:%S') + end_date_time = datetime.strptime(end_date_time_str, '%Y-%m-%d %H:%M:%S') + nx = (xml_root.find('nx').text).split('.')[0] + ny = (xml_root.find('ny').text).split('.')[0] + loutstep = xml_root.find('loutstep').text if xml_root.find( + 'loutstep') is not None else '3600' + return { + 'start_date_time': start_date_time, + 'end_date_time': end_date_time, + 'out_longitude': xml_root.find('se_lon').text, + 'out_latitude': xml_root.find('se_lat').text, + 'num_x_grid': nx, + 'num_y_grid': ny, + 'dx_out': xml_root.find('dlat').text, + 'dy_out': xml_root.find('dlon').text, + 'minheight': xml_root.find('minheight').text, + 'maxheight': xml_root.find('maxheight').text, + 'loutstep': loutstep + } + +with open(os.path.basename(basename) + '.txt', newline='') as csvfile: + csv_reader = csv.reader(csvfile, delimiter='\t') + csv_header = next(csv_reader) + for row in csv_reader: + # name of params in row and it's order: + # calc_id(0), use(1), m_id(2), s_id(3), station(4), country(5), s_lat(6),s_lng(7), + # id_nuclide(8), name_nuclide(9), date_start(10), time_start(11), date_end(12), time_end(13), val(14), sigma + measurement_id = row[2] + latitude_1 = float(row[6]) - 0.001 + latitude_2 = float(row[6]) + 0.001 + longitude_1 = float(row[7]) - 0.001 + longitude_2 = float(row[7]) + 0.001 + start_date_time = datetime.strptime( + row[10] + row[11], '%Y-%m-%d%H:%M:%S') + end_date_time = datetime.strptime( + row[12] + row[13], '%Y-%m-%d%H:%M:%S') + species_mass = "{:e}".format(float(row[14])) + releases_params.append({ + 'id': measurement_id, + 'latitude_1': round(latitude_1,3), + 'longitude_1': round(longitude_1,3), + 'latitude_2': round(latitude_2, 3), + 'longitude_2': round(longitude_2, 3), + 'species_name': row[9], + 'start_date_time': start_date_time, + 'end_date_time': end_date_time, + 'mass': species_mass, + 'comment': "RELEASE " + measurement_id + }) + +def parse_command_file(): + start_date_time = user_params['start_date_time'] + end_date_time = releases_params[-1]['end_date_time'] + timedelta(hours=1) + command_body = """&COMMAND + LDIRECT= -1, ! Simulation direction in time ; 1 (forward) or -1 (backward) + IBDATE= {date_1}, ! Start date of the simulation ; YYYYMMDD: YYYY=year, MM=month, DD=day + IBTIME= {time_1}, ! Start time of the simulation ; HHMISS: HH=hours, MI=min, SS=sec; UTC + IEDATE= {date_2}, ! End date of the simulation ; same format as IBDATE + IETIME= {time_2}, ! End time of the simulation ; same format as IBTIME + LOUTSTEP= {loutstep}, ! Interval of model output; average concentrations calculated every LOUTSTEP (s) + LOUTAVER= {loutstep}, ! Interval of output averaging (s) + LOUTSAMPLE= {loutstep}, ! Interval of output sampling (s), higher stat. accuracy with shorter intervals + ITSPLIT= 99999999, ! Interval of particle splitting (s) + LSYNCTIME= 60, ! All processes are synchronized to this time interval (s) + CTL= -5.0000000, ! CTL>1, ABL time step = (Lagrangian timescale (TL))/CTL, uses LSYNCTIME if CTL<0 + IFINE= 4, ! Reduction for time step in vertical transport, used only if CTL>1 + IOUT= 9, ! Output type: [1]mass 2]pptv 3]1&2 4]plume 5]1&4, +8 for NetCDF output + IPOUT= 0, ! Particle position output: 0]no 1]every output 2]only at end 3]time averaged + LSUBGRID= 1, ! Increase of ABL heights due to sub-grid scale orographic variations;[0]off 1]on + LCONVECTION= 1, ! Switch for convection parameterization;0]off [1]on + LAGESPECTRA= 0, ! Switch for calculation of age spectra (needs AGECLASSES);[0]off 1]on + IPIN= 0, ! Warm start from particle dump (needs previous partposit_end file); [0]no 1]yes + IOUTPUTFOREACHRELEASE= 1, ! Separate output fields for each location in the RELEASE file; [0]no 1]yes + IFLUX= 0, ! Output of mass fluxes through output grid box boundaries + MDOMAINFILL= 0, ! Switch for domain-filling, if limited-area particles generated at boundary + IND_SOURCE= 1, ! Unit to be used at the source ; [1]mass 2]mass mixing ratio + IND_RECEPTOR= 1, ! Unit to be used at the receptor; [1]mass 2]mass mixing ratio 3]wet depo. 4]dry depo. + MQUASILAG= 0, ! Quasi-Lagrangian mode to track individual numbered particles + NESTED_OUTPUT= 0, ! Output also for a nested domain + LINIT_COND= 0, ! Output sensitivity to initial conditions (bkw mode only) [0]off 1]conc 2]mmr + SURF_ONLY= 0, ! Output only for the lowest model layer, used w/ LINIT_COND=1 or 2 + CBLFLAG= 0, ! Skewed, not Gaussian turbulence in the convective ABL, need large CTL and IFINE + OHFIELDS_PATH= "../../flexin/", ! Default path for OH file + / +""".format(date_1=start_date_time.strftime('%Y%m%d'), + time_1=start_date_time.strftime('%H%M%S'), + date_2=end_date_time.strftime('%Y%m%d'), + time_2=end_date_time.strftime('%H%M%S'), + loutstep=user_params['loutstep']) + write_to_file('/options/COMMAND', template_header + command_body) + logging.info('Parsing COMMAND file compleated.') + +def parse_outgrid_file(): + outgrid_template = """!******************************************************************************* +! * +! Input file for the Lagrangian particle dispersion model FLEXPART * +! Please specify your output grid * +! * +! OUTLON0 = GEOGRAPHYICAL LONGITUDE OF LOWER LEFT CORNER OF OUTPUT GRID * +! OUTLAT0 = GEOGRAPHYICAL LATITUDE OF LOWER LEFT CORNER OF OUTPUT GRID * +! NUMXGRID = NUMBER OF GRID POINTS IN X DIRECTION (= No. of cells + 1) * +! NUMYGRID = NUMBER OF GRID POINTS IN Y DIRECTION (= No. of cells + 1) * +! DXOUT = GRID DISTANCE IN X DIRECTION * +! DYOUN = GRID DISTANCE IN Y DIRECTION * +! OUTHEIGHTS = HEIGHT OF LEVELS (UPPER BOUNDARY) * +!******************************************************************************* +&OUTGRID + OUTLON0= {out_lon}, + OUTLAT0= {out_lat}, + NUMXGRID= {x_grid}, + NUMYGRID= {y_grid}, + DXOUT= {dx_out}, + DYOUT= {dy_out}, + OUTHEIGHTS= {maxheight}, + / +""".format(out_lon=user_params['out_longitude'], + out_lat=user_params['out_latitude'], + x_grid=user_params['num_x_grid'], + y_grid=user_params['num_y_grid'], + dx_out=user_params['dx_out'], + dy_out=user_params['dy_out'], + maxheight=user_params['maxheight']) + write_to_file('/options/OUTGRID', outgrid_template) + logging.info('Parsing OUTGRID file compleated.') + + +def parse_simflex_input(id, new_output_file_path): + flexpart_output_header = """#id_obs; path_to_file; srs_ind; +""" + flexpart_output_template = """{obs_id}; {path_to_file}; {srs_id} +""".format(obs_id=id, path_to_file=new_output_file_path, srs_id=1) + table_srs_file_name = '/simflex' + '/table_srs_paths.txt' + if id > 1: + write_to_file(table_srs_file_name, flexpart_output_template, 'a') + else: + write_to_file(table_srs_file_name, flexpart_output_header + flexpart_output_template) + +def parse_simflexinp(): + start_date_time = user_params['start_date_time'] + simflexinp_template = """$simflexinp +redirect_console=.false., +Niso_=11, +Isolines_(1:11) = 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.95 0.99, +Threshprob_=0.9, +syear_={start_year}, +smon_={start_month}, +sday_={start_day}, +shr_={start_hour}, +sminut_={start_min}, +loutstep_=3600, ! in flexpart COMMAND +tstart_max_=-9999.0, +thresh_start_=1.0, +min_duration_=3600, +dlon_={dlon}, +dlat_={dlat}, +outlon_={se_lon}, ! 0. default +outlat_={se_lat}, +nlon_={nx}, +nlat_={ny}, +nhgt_={minheight}, +DHgt_={maxheight} +$end +""".format(start_year=start_date_time.strftime('%Y'), + start_month=start_date_time.strftime('%m'), + start_day=start_date_time.strftime('%d'), + start_hour=start_date_time.strftime('%H'), + start_min=start_date_time.strftime('%M'), + dlon=user_params['dy_out'], + dlat=user_params['dx_out'], + se_lon=user_params['out_longitude'], + se_lat=user_params['out_latitude'], + nx=user_params['num_x_grid'], + ny=user_params['num_y_grid'], + minheight=user_params['minheight'], + maxheight=user_params['maxheight']) + write_to_file('/simflex/simflexinp.nml', simflexinp_template) + logging.info('Parsing simflexinp file compleated.') + +def parse_releases_file(releases_params): + SPECIES_BY_ID = {"O3": '002', "NO": '003', "NO2": '004', + "HNO3": '005', "HNO2": '006', "H2O2": '007', + "NO2": '008', "HCHO": '009', "PAN": '010', + "NH3": '011', "SO4-aero": '012', "NO3-aero": '013', + "I2-131": '014', "I-131": '015', "Cs-137": '016', + "Y-91": '017', "Ru-106": '018', "Kr-85": '019', + "Sr-90": '020', "Xe-133": '021', "CO": '022', + "SO2": '023', "AIRTRACER": '024', "AERO-TRACE": '025', + "CH4": '026', "C2H6": '027', "C3H8": '028', + "PCB28": '031', "G-HCH": '034', "BC": '040'} + species_id = int(SPECIES_BY_ID.get(releases_params['species_name'])) + release_header = """&RELEASES_CTRL + NSPEC = 1, ! Total number of species + SPECNUM_REL= {species_number}, ! Species numbers in directory SPECIES + / +""".format(species_number=species_id) + + release_body = """&RELEASE + IDATE1 = {date_1}, + ITIME1 = {time_1}, + IDATE2 = {date_2}, + ITIME2 = {time_2}, + LON1 = {lon_1}, + LON2 = {lon_2}, + LAT1 = {lat_1}, + LAT2 = {lat_2}, + Z1 = {minheight}, + Z2 = {maxheight}, + ZKIND = 1, + MASS = {mass}, + PARTS = 10000, + COMMENT = "{comment}", + / +""".format(date_1=releases_params['start_date_time'].strftime('%Y%m%d'), + time_1=releases_params['start_date_time'].strftime('%H%M%S'), + date_2=releases_params['end_date_time'].strftime('%Y%m%d'), + time_2=releases_params['end_date_time'].strftime('%H%M%S'), + lon_1=releases_params['longitude_1'], + lon_2=releases_params['longitude_2'], + lat_1=releases_params['latitude_1'], + lat_2=releases_params['latitude_2'], + mass=releases_params['mass'], + comment=releases_params['comment'], + minheight=user_params['minheight'], + maxheight=user_params['maxheight']) + + write_to_file('/options/RELEASES', template_header + + release_header + release_body) + logging.info('Parsing RELEASE file compleated.') + +user_params = get_xml_params() + +# First date from user last is the last release date + 1 hour +# It also creates AVAILABLE file and fill it +logging.info('Started loading grib data.') +download_grib(user_params['start_date_time'], releases_params[-1]['end_date_time']) +logging.info('Finished loading grib data and filling AVAILABLE file.\n') + +parse_command_file() +parse_outgrid_file() +parse_simflexinp() + +output_file_id = 1 +end_date_time_str = ( + releases_params[-1]['end_date_time'] + timedelta(hours=1)).strftime('%Y%m%d%H%M%S') +for param in releases_params: + logging.info('Running FLEXPART ' + str(output_file_id) + + ' iteration in ' + str(len(releases_params)) + '.') + parse_releases_file(param) + rc = run("time ./FLEXPART", shell=True) + # move output prognose to simflex folder and rename it according to the release comments + output_file_name = '/grid_time_' + end_date_time_str + old_output_file_path = basename + '/output' + output_file_name + '.nc' + new_output_file_path = simflex_dir_name + output_file_name + '_' + str(output_file_id) + '.nc' + os.rename(old_output_file_path, new_output_file_path) + + parse_simflex_input(output_file_id, new_output_file_path) + output_file_id = output_file_id + 1 + logging.info('FLEXPART completed calculation.\n') + +# nohup /data/calculations/33_14/parse_input.py &