Commit d1919418 authored by Christoph Heim's avatar Christoph Heim
Browse files

Working on ARPEGE still.

parent 0348d228
......@@ -4,7 +4,7 @@
description: Extract lat-lon box of data from model ICON.
author: Christoph Heim
date created: 27.06.2019
date changed: 15.07.2019
date changed: 17.07.2019
usage: arguments:
1st: n jobs for multiprocessing pool
python: 3.5.2
......@@ -26,8 +26,8 @@ def comp_weights_file(target_grid, weights_file, inp_file, grid_def_file,
res, box, options):
"""
"""
if (not os.path.exists(target_grid)) or (options['recompute']):
write_grid_file(box, target_grid, res)
#if (not os.path.exists(target_grid)) or (options['recompute']):
# write_grid_file(box, target_grid, res)
print('Compute weights file')
ofile = cdo.gennn(target_grid, input=(' -setgrid,'+grid_def_file+
......@@ -194,6 +194,7 @@ if __name__ == '__main__':
# target grid on which to interpolate the model output
target_grid = os.path.join('grids','latlon_{}km_dom_{}'.format(
res, domain['code']))
write_grid_file(box, target_grid, res)
# find times and files that should be extracted
# and prepare arguments for function
......
......@@ -4,7 +4,7 @@
description: Extract lat-lon box of data from model MPAS.
author: Christoph Heim
date created: 05.07.2019
date changed: 16.07.2019
date changed: 17.07.2019
usage: arguments:
1st: n jobs for multiprocessing pool
for 2D field: 8 jobs possible
......@@ -27,8 +27,8 @@ def comp_weights_file(target_grid, weights_file, inp_file, grid_def_file,
res, box, options):
"""
"""
if (not os.path.exists(target_grid)) or (options['recompute']):
write_grid_file(box, target_grid, res)
#if (not os.path.exists(target_grid)) or (options['recompute']):
# write_grid_file(box, target_grid, res)
print('Compute weights file')
......@@ -37,7 +37,8 @@ def comp_weights_file(target_grid, weights_file, inp_file, grid_def_file,
# ' '+inp_file),
# output=weights_file)
input = "-sellonlatbox,{},{},{},{} -setgrid,mpas:{} {}".format(box['lon'].start,box['lon'].stop,
input = "-sellonlatbox,{},{},{},{} -setgrid,mpas:{} {}".format(
box['lon'].start,box['lon'].stop,
box['lat'].start,box['lat'].stop,
grid_def_file,
inp_file)
......@@ -72,13 +73,14 @@ def sellatlon_MPAS(inp_file, out_file, dt, box, options, var_dict,
'-selname,'+var_dict['file']+' '+inp_file
),
output=tmp_file, options='-f nc4')
TM.stop('prepr')
TM.start('grid')
if not os.path.exists(weights_file):
comp_weights_file(target_grid, weights_file,
tmp_file, grid_def_file,
res, box, options)
TM.stop('prepr')
TM.stop('grid')
# cdo
TM.start('cdo')
......@@ -131,7 +133,7 @@ if __name__ == '__main__':
# variables to extract
var_names = ['QC', 'T',
'TQC',]
var_names = ['LWUTOA']
var_names = ['QC']
# model resolutions [km] of simulations
ress = [7.5, 3.75]
......@@ -140,7 +142,7 @@ if __name__ == '__main__':
# date range
first_date = datetime(2016,8,10)
last_date = datetime(2016,8,19)
last_date = datetime(2016,8,10)
# options for computation
options = {}
......@@ -218,6 +220,7 @@ if __name__ == '__main__':
# target grid on which to interpolate the model output
target_grid = os.path.join('grids','latlon_{}km_dom_{}'.format(
res, domain['code']))
write_grid_file(box, target_grid, res)
# find times and files that should be extracted
# and prepare arguments for function
......@@ -234,7 +237,6 @@ if __name__ == '__main__':
# compute weights file if necessary
if not os.path.exists(weights_file):
TM.start('prepr')
dt = dt_range[0]
inp_file = glob.glob(os.path.join(
inp_dir,var_dict[var_name]['type']+
......@@ -244,7 +246,6 @@ if __name__ == '__main__':
sellatlon_MPAS(inp_file, out_file, dt, box, options,
var_dict[var_name], weights_file,
target_grid, var_name, res)
TM.stop('prepr')
# run function serial or parallel
if n_tasks > 1:
......
......@@ -4,7 +4,7 @@
description: Extract lat-lon box of data from model ARPEGE-NH.
author: Christoph Heim
date created: 09.07.2019
date changed: 16.07.2019
date changed: 17.07.2019
usage: arguments:
1st: n jobs for multiprocessing pool
python: 3.5.2
......@@ -16,13 +16,44 @@ from datetime import datetime, timedelta
from multiprocessing import Pool
from pathlib import Path
from cdo import Cdo
from package.utilities import Timer, cdo_mergetime
from package.utilities import Timer, cdo_mergetime, write_grid_file
from namelist import domain
from functions import paste_dir_names
###############################################################################
def sellatlon_ARPEGE(inp_file, out_file, dt, box, options, var_name, var_dict,
res):
def get_splf(split_files, var_dict):
"""
Get glob search pattern for grib var split file.
"""
split_file_descr = '{}.{}*{}.gp'
splf = split_file_descr.format(split_files, var_dict['grb_srf'],
var_dict['file'])
return(splf)
def comp_weights_file(target_grid, weights_file, inp_file, grid_def_file,
res, box, options):
"""
"""
print('Compute weights file')
#ofile = cdo.gennn(target_grid,
# input=(' -setgrid,mpas:'+grid_def_file+
# ' '+inp_file),
# output=weights_file)
frmt_string = ("-sellonlatbox,{},{},{},{} -setgrid,{}"+
" -setgridtype,regular {}")
input = frmt_string.format(
box['lon'].start,box['lon'].stop,
box['lat'].start,box['lat'].stop,
grid_def_file,
inp_file)
ofile = cdo.gennn(target_grid,
input=input, output=weights_file)
def sellatlon_ARPEGE(inp_file, out_file, dt, box, options, var_name, var_dicts,
res, weights_file, target_grid):
TM = Timer()
file_code = '{}km_{}_{:%Y%m%d%H%M}'.format(res, var_name, dt)
......@@ -35,27 +66,96 @@ def sellatlon_ARPEGE(inp_file, out_file, dt, box, options, var_name, var_dict,
print('\t\t'+file_code)
TM.start('prep')
split = os.path.split(out_file)
tmp_file = os.path.join(split[0],'tmp_'+split[1][:-3])
tmp_dir = os.path.join(split[0],'dirtmp_{:%Y%m%d%H%M}'.format(dt))
Path(tmp_dir).mkdir(parents=True, exist_ok=True)
split_files = os.path.join(tmp_dir,'split')
# run gribsplit if not already done
search = glob.glob(get_splf(split_files, var_dicts[var_name]))
if len(search) == 0:
# Split original grib files
command = './gribsplit'
subprocess.call([command, inp_file, split_files],
stdout=subprocess.DEVNULL)
# remove all split variables that are irrelevant
keep_files = []
for key,var_dict in var_dicts.items():
search = glob.glob(get_splf(split_files, var_dict))
if var_dict['vdim'] == '3D':
# Remove unnecessary levels
levels = ['l{}00'.format(lev) for lev in range(box['vert0'],box['vert1']+1)]
remove = []
for file in search:
match = False
for lev in levels:
if lev in file:
match = True
if not match:
remove.append(file)
for file in remove:
search.remove(file)
keep_files.extend(search)
search = glob.glob('{}/*'.format(tmp_dir))
for file in search:
if file not in keep_files:
#print('remove ',file)
os.remove(file)
#command = './gribsplit {} {}'.format(inp_file, tmp_file)
command = './gribsplit'
subprocess.call([command, inp_file, tmp_file],
stdout=subprocess.DEVNULL)
tmp_files = glob.glob(get_splf(split_files, var_dicts[var_name]))
TM.stop('prep')
quit()
TM.start('grid')
if not os.path.exists(weights_file):
comp_weights_file(target_grid, weights_file,
tmp_files[0], grid_def_file,
res, box, options)
TM.stop('grid')
# cdo
TM.start('cdo')
#ofile = cdo.sellonlatbox(
# box['lon0'],box['lon1'],
# box['lat0'],box['lat1'],
# input=(' -sellevidx,'+str(box['vert0'])+'/'+
# str(box['vert1'])+
# ' -selname,'+var_dict[var_name]['key']+
# ' '+inp_file),
# output=out_file)
merge_files = []
for tmp_file in tmp_files:
print(tmp_file)
input = ("-sellonlatbox,{},{},{},{} -setgrid,{}"+
" -setgridtype,regular {}").format(
box['lon'].start, box['lon'].stop,
box['lat'].start, box['lat'].stop,
grid_def_file,
tmp_file)
#input = ("-remap,{},{} -setgrid,{}"+
# " -setgridtype,regular {}").format(
# target_grid, weights_file,
# grid_def_file, tmp_file)
if var_dicts[var_name]['vdim'] == '3D':
out_file_use = tmp_file + '.nc'
merge_files.append(out_file_use)
else:
out_file_use = out_file
ofile = cdo.remap(target_grid, weights_file,
input=input, output=out_file_use, options='-f nc4')
#ofile = cdo.sellonlatbox(
# box['lon'].start, box['lon'].stop,
# box['lat'].start, box['lat'].stop,
# input=input, output=out_file_use, options='-f nc4')
# merge vertical levels
if var_dicts[var_name]['vdim'] == '3D':
merge_files.sort()
cdo.merge(input=merge_files, output=out_file)
if options['rm_tmp_files'] and (
var_dicts[var_name]['vdim'] == '3D'):
for file in merge_fies:
os.remove(file)
TM.stop('cdo')
#quit()
return(TM)
......@@ -77,9 +177,13 @@ if __name__ == '__main__':
model_name = 'ARPEGE-NH'
# variables to extract
var_names = ['QC', 'T', 'H']
var_names = ['H']
var_names = ['LWUTOA']
var_names = ['QV', 'QC', 'T', 'W', 'H',
'U10M', 'V10M', 'T2M', 'LWUTOA', 'SWDSFC',
'SLHFLX', 'SSHFLX', 'TQC']
#var_names = ['QC']
#var_names = ['W']
var_names = ['U10M', 'V10M', 'T2M', 'LWUTOA', 'SWDSFC',
'SLHFLX', 'SSHFLX', 'TQC']
# model resolutions [km] of simulations
ress = [2.5]
......@@ -91,6 +195,7 @@ if __name__ == '__main__':
# options for computation
options = {}
options['recompute'] = 0
options['rm_tmp_files'] = 0
options['rm_tmp_folder'] = 0
###########################################################################
......@@ -98,12 +203,23 @@ if __name__ == '__main__':
# ARPEGE SPECIFIC SETTINGS
###########################################################################
var_dict = {
'QC' :{'file':'0.1.83','vdim':'3D',},
'T' :{'file':'0.0.0' ,'vdim':'3D',},
'H' :{'file':'0.3.4' ,'vdim':'3D',},
'LWUTOA':{'file':'0.5.5' ,'vdim':'2D',},
'QV' :{'file':'0.1.0', 'grb_srf':'t119', 'vdim':'3D',},
'QC' :{'file':'0.1.83', 'grb_srf':'t119', 'vdim':'3D',},
'T' :{'file':'0.0.0', 'grb_srf':'t119', 'vdim':'3D',},
'W' :{'file':'0.2.9', 'grb_srf':'t119', 'vdim':'3D',},
'H' :{'file':'0.3.4' , 'grb_srf':'t119', 'vdim':'3D',},
'U10M' :{'file':'0.2.2', 'grb_srf':'t103', 'vdim':'2D',},
'V10M' :{'file':'0.2.3', 'grb_srf':'t103', 'vdim':'2D',},
'T2M' :{'file':'0.0.0', 'grb_srf':'t103', 'vdim':'2D',},
'LWUTOA':{'file':'0.5.5', 'grb_srf':'t8', 'vdim':'2D',},
'SWDSFC':{'file':'0.4.7', 'grb_srf':'t1', 'vdim':'2D',},
'SLHFLX':{'file':'0.0.10', 'grb_srf':'t1', 'vdim':'2D',},
'SSHFLX':{'file':'0.0.11', 'grb_srf':'t1', 'vdim':'2D',},
'TQC' :{'file':'192.128.78', 'grb_srf':'t1', 'vdim':'2D',},
}
inc_min = {'3D':180, '2D':15}
remap_res = 2.5
###########################################################################
## PREPARING STEPS
......@@ -141,23 +257,46 @@ if __name__ == '__main__':
shutil.rmtree(out_tmp_dir)
Path(out_tmp_dir).mkdir(parents=True, exist_ok=True)
# prepare grid interpolation
# ARPGEGE grid definition file
grid_def_file = 'griddes.arpege1'
# weights file to compute once for a grid
weights_file = os.path.join(out_tmp_dir,
'weights_{}km_dom_{}'.format(remap_res,
domain['code']))
# target grid on which to interpolate the model output
target_grid = os.path.join('grids','latlon_{}km_dom_{}'.format(
remap_res, domain['code']))
write_grid_file(box, target_grid, remap_res)
# find times and files that should be extracted
# and prepare arguments for function
# (treat ingenious idea of using 2400 for midnight instead of 0000)
args = []
for dt in dt_range:
for ti,dt in enumerate(dt_range):
if (dt.hour == 0) and (dt.minute == 0):
dt = dt - timedelta(days=1)
dt_tmp = dt - timedelta(days=1)
file_name = 'ARPNH{}{:%Y%m%d}2400'.format(
var_dict[var_name]['vdim'], dt)
var_dict[var_name]['vdim'], dt_tmp)
inp_file = glob.glob(os.path.join(inp_dir,
'{:%Y%m%d}'.format(dt_tmp), file_name))[0]
else:
file_name = 'ARPNH{}{:%Y%m%d%H%M}'.format(
var_dict[var_name]['vdim'], dt)
inp_file = glob.glob(os.path.join(inp_dir,
inp_file = glob.glob(os.path.join(inp_dir,
'{:%Y%m%d}'.format(dt), file_name))[0]
out_file = os.path.join(out_tmp_dir,
var_name+'_{:%Y%m%d%H%M}'.format(dt)+'.nc')
args.append( (inp_file, out_file, dt, box, options, var_name,
var_dict, res) )
var_dict, res, weights_file, target_grid) )
# compute weights file if necessary
if ti == 0:
if not os.path.exists(weights_file):
sellatlon_ARPEGE(inp_file, out_file, dt, box, options,
var_name, var_dict, res, weights_file,
target_grid)
# run function serial or parallel
if n_tasks > 1:
......
......@@ -2,11 +2,11 @@
export GRIB_DEFINITION_PATH=/mnt/lustre01/sw/rhel6-x64/eccodes/definitions
dirin=/work/ka1081/DYAMOND/FV3-3.25km/2016081100/
for d in 11; do
#for d in 11; do
echo $d
#cdo -P 8 -f nc -O -fldmean -sellonlatbox,-50,-30,0,20 -seltimestep,1/8 ${dirin}ql_plev_C3072_360x180.fre.nc FV3_qc_dc_mean_201608${d}.nc
#cdo -P 8 -f nc -O -fldmean -sellonlatbox,-50,-30,0,20 -seltimestep,1/8 ${dirin}h_plev_C3072_360x180.fre.nc FV3_h_dc_mean_201608${d}.nc
#cdo timavg FV3_h_dc_mean_20160811.nc FV3_h_dc_mean_201608${d}_mean.nc
#cdo timavg FV3_qc_dc_mean_20160811.nc FV3_qc_dc_mean_201608${d}_mean.nc
done
#done
gridtype = lonlat
xsize = 445
ysize = 278
xname = lon
xlongname = "longitude"
xunits = "degrees_east"
yname = lat
ylongname = "latitude"
yunits = "degrees_north"
xfirst = -95
xinc = 0.036
yfirst = -24
yinc = 0.036
#!/bin/bash
n_jobs=12
grid_size_keyword=global_0.022
target_grid_nc=0.022_grid.nc
#grid_size_keyword=global_0.02
#target_grid_nc=0.02_grid.nc
#grid_size_keyword=global_0.05
#target_grid_nc=0.05_grid.nc
griddes=griddes.arpege1
weights_file=weights_ARPEGE_test
split_file=/work/ka1081/2019_06_Hackathon_Mainz/christoph_heim/newdata/ARPEGE-NH_2.5/tmp/dirtmp_201608100000/split.t119.l5400.grb.0.3.4.gp
out_file=H_test.nc
#cdo -O -P $n_jobs --cellsearchmethod spherepart genycon,$grid_size_keyword -setgrid,$griddes \
# -setgridtype,regular $split_file $weights_file
#cdo -O -P $n_jobs -f nc -topo,$grid_size_keyword $target_grid_nc
cdo -O -f nc4 -sellonlatbox,-95,-79,-24,-14 -remap,$target_grid_nc,$weights_file \
-setgrid,$griddes -setgridtype,regular $split_file $out_file
#!/bin/bash
mkdir -p tmp_fv3
mkdir -p fv3
fv3_dir=/work/ka1081/DYAMOND/FV3-3.25km
out_dir=fv3
tmp_dir=tmp_fv3
#cdo -P 12 -genycon,${out_dir}/0.10_grid.nc \
# -setgrid,${tmp_dir}/gridspec.nc:2 \
# -collgrid,gridtype=unstructured ${fv3_dir}/2016080100/lhflx_15min.tile?.nc \
# ${tmp_dir}/fv3_weights_to_0.10deg.nc
cdo -O -f nc -topo,global_0.10 ${tmp_dir}/0.10_grid.nc
cdo -P 12 -O remapnn,${tmp_dir}/0.10_grid.nc \
-setgrid,${tmp_dir}/grid_spec.tile5.nc \
${fv3_dir}/2016080100/t2m_15min.tile5.nc \
${out_dir}/test.nc
#cdo -P 12 -O setgrid,${tmp_dir}/grid_spec.tile5.nc \
# -collgrid,gridtype=unstructured ${fv3_dir}/2016080100/t2m_15min.tile5.nc \
# ${out_dir}/test2.nc
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment