Skip to content

GitLab

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

NICAM, SAM, ICON, MPAS, IFS and UM done.

parent 48720275
Hide whitespace changes
Inline Side-by-side
03_icon.py
View file @ b0eda3bf
......@@ -4,7 +4,7 @@
description: Extract lat-lon box of data from model ICON.
author: Christoph Heim
date created: 27.06.2019
date changed: 05.07.2019
date changed: 07.07.2019
usage: arguments:
1st: n jobs for multiprocessing pool
python: 3.5.2
......@@ -16,13 +16,17 @@ from datetime import datetime, timedelta
from multiprocessing import Pool
from pathlib import Path
from cdo import Cdo
from utilities import Timer
from utilities import Timer, write_grid_file
###############################################################################
def comp_weights_file(target_grid, weights_file, inp_file, grid_def_file):
def comp_weights_file(target_grid, weights_file, inp_file, grid_def_file,
res, box):
"""
"""
if not os.path.exists(target_grid):
write_grid_file(box, target_grid, res)
print('Compute weights file')
ofile = cdo.gennn(target_grid, input=(' -setgrid,'+grid_def_file+
' '+inp_file), output=weights_file,
......@@ -72,7 +76,8 @@ if __name__ == '__main__':
'christoph_heim','newdata')
# vert box to subselect
box = {'vert0':73-14,'vert1':91-14}
box = {'lon0': 265, 'lon1': 281, 'lat0': -24, 'lat1': -14,
'vert0':73-14,'vert1':91-14}
# name of model
model_name = 'ICON'
......@@ -83,6 +88,7 @@ if __name__ == '__main__':
# model resolutions [km] of simulations
ress = [10,5,2.5]
ress = [10,2.5]
#ress = [5]
#ress = [10]
......@@ -106,17 +112,17 @@ if __name__ == '__main__':
10: {'grid_def_file':os.path.join(grid_def_base_dir,
'10km/icon_grid_0025_R02B08_G.nc'),
'weights_file':'weights_10km',
'target_grid':'grids/latlon_0.10_deg',
'target_grid':'grids/latlon_dx_10.00km',
},
5: {'grid_def_file':os.path.join(grid_def_base_dir,
'5km_2/icon_grid_0015_R02B09_G.nc'),
'weights_file':'weights_5km',
'target_grid':'grids/latlon_0.05_deg',
'target_grid':'grids/latlon_dx_05.00km',
},
2.5:{'grid_def_file':os.path.join(grid_def_base_dir,
'2.5km/icon_grid_0017_R02B10_G.nc'),
'weights_file':'weights_2.5km',
'target_grid':'grids/latlon_0.025_deg',
'target_grid':'grids/latlon_dx_02.50km',
},
}
###########################################################################
......@@ -155,7 +161,6 @@ if __name__ == '__main__':
'_' + str(res),'tmp')
Path(out_tmp_dir).mkdir(parents=True, exist_ok=True)
# prepare grid interpolation
# icon grid definition file
grid_def_file = grid_dict[res]['grid_def_file']
......@@ -182,7 +187,8 @@ if __name__ == '__main__':
if not os.path.exists(weights_file):
comp_weights_file(target_grid, weights_file,
inp_file, grid_def_file)
inp_file, grid_def_file,
res, box)
# run function serial or parallel
......
05_mpas.py
View file @ b0eda3bf
......@@ -4,7 +4,7 @@
description: Extract lat-lon box of data from model MPAS.
author: Christoph Heim
date created: 05.07.2019
date changed: 05.07.2019
date changed: 07.07.2019
usage: arguments:
1st: n jobs for multiprocessing pool
python: 3.5.2
......@@ -16,12 +16,27 @@ from datetime import datetime, timedelta
from multiprocessing import Pool
from pathlib import Path
from cdo import Cdo
from utilities import Timer
from utilities import Timer, write_grid_file
###############################################################################
#from netCDF4 import Dataset
def sellatlon_MPAS(inp_file, out_file, dt, box, i_recompute):
def comp_weights_file(target_grid, weights_file, inp_file, grid_def_file,
res, box):
"""
"""
if not os.path.exists(target_grid):
write_grid_file(box, target_grid, res)
print('Compute weights file')
ofile = cdo.gennn(target_grid,
input=(' -setgrid,mpas:'+grid_def_file+
' '+inp_file),
output=weights_file,
options='-P 48')
def sellatlon_MPAS(inp_file, out_file, dt, box, i_recompute, res, var_dict,
weights_file, target_grid):
TM = Timer()
TM.start('total')
......@@ -33,40 +48,39 @@ def sellatlon_MPAS(inp_file, out_file, dt, box, i_recompute):
TM.start('cdo')
TM.stop('cdo')
else:
print('\t\t{:%Y%m%d%H}'.format(dt))
split = os.path.split(out_file)
tmp_file = os.path.join(split[0],'tmp_'+split[1])
print(inp_file)
print(tmp_file)
# first step
TM.start('prepr')
#ofile = cdo.setattribute('*@axis="txz"',
# input=(
# '-selname,qc '+inp_file
# ),
# output=tmp_file)
if not os.path.exists(tmp_file):
ofile = cdo.setattribute('*@axis="txz"',
input=(
'-selname,'+var_dict['file']+' '+inp_file
),
output=tmp_file, options='-f nc4')
if not os.path.exists(weights_file):
comp_weights_file(target_grid, weights_file,
tmp_file, grid_def_file,
res, box)
TM.stop('prepr')
#TODO
grid_def = '/work/ka1081/2019_06_Hackathon_Mainz/falko/MPAS_7.5km_grid.nc'
# cdo
TM.start('cdo')
print('\t{:%Y%m%d%H}'.format(dt))
ofile = cdo.sellonlatbox(
box['lon0'],box['lon1'],
box['lat0'],box['lat1'],
ofile = cdo.remap(target_grid, weights_file,
input=(
'-setgrid,mpas:'+grid_def+
' ' + tmp_file,
),
#input=('-sellevidx,'+str(box['vert0'])+'/'+
# str(box['vert1'])+' '+inp_file),
output=out_file, options='-f nc4')
#print('\t\t{:%Y%m%d%H} completed'.format(dt))
' -sellevidx,'+
str(box['vert0'])+'/'+str(box['vert1'])+
' -setgrid,mpas:'+grid_def_file+
' '+tmp_file),
output=out_file, options='-f nc')
TM.stop('cdo')
TM.stop('total')
......@@ -85,21 +99,22 @@ if __name__ == '__main__':
# lat lon vert box to subselect
box = {'lon0': 265, 'lon1': 281, 'lat0': -24, 'lat1': -14,
'vert0':1,'vert1':20}
'vert0':1,'vert1':22}
# name of model
model_name = 'MPAS'
# variables to extract
var_names = ['QC', 'T']
var_names = ['QC']
#var_names = ['QC']
# model resolutions [km] of simulations
ress = [7.5]
ress = [7.5, 3.75]
#ress = [3.75]
# date range
first_date = datetime(2016,8,10)
last_date = datetime(2016,8,10)
last_date = datetime(2016,8,19)
# recompute cdo
i_recompute = 0
......@@ -108,9 +123,23 @@ if __name__ == '__main__':
# MPAS SPECIFIC SETTINGS
###########################################################################
grid_def_base_dir = os.path.join('/work','ka1081',
'2019_06_Hackathon_Mainz', 'falko')
var_dict = {
'QC':{'file':'qc','type':'history'},
'T':{'file':'ta',},
'T':{'file':'temperature','type':'history'},
}
grid_dict = {
3.75:{'grid_def_file':os.path.join(grid_def_base_dir,
'MPAS_3.75km_grid.nc'),
'weights_file':'weights_3.75km',
'target_grid':'grids/latlon_dx_03.75km',
},
7.5: {'grid_def_file':os.path.join(grid_def_base_dir,
'MPAS_7.5km_grid.nc'),
'weights_file':'weights_7.5km',
'target_grid':'grids/latlon_dx_07.50km',
},
}
inc_min = {'history':180}
###########################################################################
......@@ -152,15 +181,34 @@ if __name__ == '__main__':
'_' + str(res),'tmp')
Path(out_tmp_dir).mkdir(parents=True, exist_ok=True)
# prepare grid interpolation
# mpas grid definition file
grid_def_file = grid_dict[res]['grid_def_file']
# weights file to compute once for a grid
weights_file = os.path.join(out_tmp_dir,
grid_dict[res]['weights_file'])
# target grid on which to interpolate the model output
target_grid = grid_dict[res]['target_grid']
# find times and files that should be extracted
# and prepare arguments for function
args = []
for dt in dt_range:
inp_file = glob.glob(os.path.join(inp_dir,var_dict[var_name]['type']+
'.{:%Y-%m-%d_%H.%M.%S}.nc'.format(dt)))[0]
inp_file = glob.glob(os.path.join(
inp_dir,var_dict[var_name]['type']+
'.{:%Y-%m-%d_%H.%M.%S}.nc'.format(dt)))[0]
out_file = os.path.join(out_tmp_dir,
var_name+'_{:%Y%m%d}'.format(dt)+'.nc')
args.append( (inp_file, out_file, dt, box, i_recompute) )
args.append( (inp_file, out_file, dt, box, i_recompute,
res, var_dict[var_name], weights_file,
target_grid) )
#if not os.path.exists(weights_file):
# comp_weights_file(target_grid, weights_file,
# inp_file, grid_def_file,
# res, box)
# run function serial or parallel
if n_tasks > 1:
......
06_ifs.py
View file @ b0eda3bf
......@@ -4,7 +4,7 @@
description: Extract lat-lon box of data from model IFS.
author: Christoph Heim
date created: 05.07.2019
date changed: 05.07.2019
date changed: 07.07.2019
usage: arguments:
1st: n jobs for multiprocessing pool
python: 3.5.2
......@@ -19,7 +19,7 @@ from cdo import Cdo
from utilities import Timer
###############################################################################
def sellatlon_IFS(inp_file, out_file, dt, box, i_recompute):
def sellatlon_IFS(inp_file, out_file, dt, box, i_recompute, var_dict):
TM = Timer()
TM.start('total')
......@@ -30,25 +30,32 @@ def sellatlon_IFS(inp_file, out_file, dt, box, i_recompute):
TM.stop('cdo')
else:
print('\t\t{:%Y%m%d%H}'.format(dt))
split = os.path.split(out_file)
tmp_file = os.path.join(split[0],'tmp_'+split[1])
# cdo
TM.start('cdo')
#print('\t{:%Y%m%d%H} cdo'.format(dt))
if not os.path.exists(tmp_file):
cdo.select('name='+var_dict['file'],
input=(inp_file),
options='--eccodes',
output=tmp_file)
ofile = cdo.sellonlatbox(
box['lon0'],box['lon1'],
box['lat0'],box['lat1'],
input=('-sellevidx,'+str(box['vert0'])+'/'+
str(box['vert1'])+' -settaxis,'+
'{:%Y-%m-%d,%H:%M:%S,3hour}'.format(dt)+
' '+nco_file),
output=out_file)
os.remove(tmp_file)
input=(
'-sellevel,'+str(box['vert0'])+'/'+
str(box['vert1'])+
' -setgridtype,regular'+
' '+tmp_file),
output=out_file, options='-f nc')
print('\t\t{:%Y%m%d%H} completed'.format(dt))
TM.stop('cdo')
TM.stop('total')
......@@ -67,22 +74,22 @@ if __name__ == '__main__':
# lat lon vert box to subselect
box = {'lon0': 265, 'lon1': 281, 'lat0': -24, 'lat1': -14,
'vert0':1,'vert1':28}
'vert0':105,'vert1':137}
# name of model
model_name = 'IFS'
# variables to extract
var_names = ['QC', 'T']
#var_names = ['QC']
var_names = ['T']
# model resolutions [km] of simulations
ress = [4,9]
ress = [9]
#ress = [9]
# date range
first_date = datetime(2016,8,10)
last_date = datetime(2016,8,10)
last_date = datetime(2016,8,19)
# recompute cdo
i_recompute = 0
......@@ -92,10 +99,12 @@ if __name__ == '__main__':
# IFS SPECIFIC SETTINGS
###########################################################################
var_dict = {
'QC':{'file':'QC',
'loc':'OUT_3D',},
'T':{'file':'TABS',
'loc':'OUT_3D',},
'QC':{'file':'clwc',
'group':'mars_out_ml_moist',},
'T':{'file':'t',
'group':'gg_mars_out_ml_upper_sh',},
#'PHI':{'file':'z',
# 'group':'mars_out',},
}
base_time = datetime(2016,8,1)
......@@ -121,8 +130,7 @@ if __name__ == '__main__':
#res = 4
sim_name = model_name + '-' + str(res) + 'km'
inp_dir = os.path.join(raw_data_dir, sim_name,
var_dict[var_name]['loc'])
inp_dir = os.path.join(raw_data_dir, sim_name)
# directory for final model output (after mergetime)
out_dir = os.path.join(out_base_dir, model_name + '_' + str(res))
Path(out_dir).mkdir(parents=True, exist_ok=True)
......@@ -134,10 +142,11 @@ if __name__ == '__main__':
# find times and files that should be extracted
inp_files_glob = glob.glob(os.path.join(inp_dir,
'*_'+var_dict[var_name]['file']+'.nc'))
times = [base_time + timedelta(seconds=dt*int(
f[var_dict[var_name]['fntime'][0]:
var_dict[var_name]['fntime'][1]]))
var_dict[var_name]['group']+'.*'))
inp_files_glob = [f for f in inp_files_glob if not '.tar.gz' in f]
times = [base_time + timedelta(hours=int(
f.split('.')[1]))
for f in inp_files_glob]
use_times = [dt for dt in times if dt >= first_date and
dt < last_date+timedelta(days=1)]
......@@ -152,7 +161,7 @@ if __name__ == '__main__':
out_file = os.path.join(out_tmp_dir,
var_name+'_{:%Y%m%d%H}'.format(use_times[i])+'.nc')
args.append( (inp_file, out_file, use_times[i], box,
i_recompute) )
i_recompute, var_dict[var_name]) )
# run function serial or parallel
if n_tasks > 1:
......
grids/latlon_0.025_deg → grids/latlon_0.075_deg 100755 → 100644
View file @ b0eda3bf
gridtype = lonlat
xsize = 400
ysize = 400
xsize = 238
ysize = 149
xname = lon
xlongname = "longitude"
xunits = "degrees_east"
......@@ -8,6 +8,6 @@ yname = lat
ylongname = "latitude"
yunits = "degrees_north"
xfirst = 265
xinc = 0.025
xinc = 0.067
yfirst = -24
yinc = 0.025
yinc = 0.067
grids/latlon_0.10_deg → grids/latlon_dx_02.500km 100755 → 100644
View file @ b0eda3bf
gridtype = lonlat
xsize = 100
ysize = 100
xsize = 712
ysize = 445
xname = lon
xlongname = "longitude"
xunits = "degrees_east"
......@@ -8,6 +8,6 @@ yname = lat
ylongname = "latitude"
yunits = "degrees_north"
xfirst = 265
xinc = 0.10
xinc = 0.022
yfirst = -24
yinc = 0.10
yinc = 0.022
grids/latlon_dx_02.50km 0 → 100644
View file @ b0eda3bf
gridtype = lonlat
xsize = 712
ysize = 445
xname = lon
xlongname = "longitude"
xunits = "degrees_east"
yname = lat
ylongname = "latitude"
yunits = "degrees_north"
xfirst = 265
xinc = 0.022
yfirst = -24
yinc = 0.022
grids/latlon_dx_07.50km 0 → 100644
View file @ b0eda3bf
gridtype = lonlat
xsize = 238
ysize = 149
xname = lon
xlongname = "longitude"
xunits = "degrees_east"
yname = lat
ylongname = "latitude"
yunits = "degrees_north"
xfirst = 265
xinc = 0.067
yfirst = -24
yinc = 0.067
grids/latlon_0.05_deg → grids/latlon_dx_10.000km 100755 → 100644
View file @ b0eda3bf
gridtype = lonlat
xsize = 200
ysize = 200
xsize = 178
ysize = 112
xname = lon
xlongname = "longitude"
xunits = "degrees_east"
......@@ -8,6 +8,6 @@ yname = lat
ylongname = "latitude"
yunits = "degrees_north"
xfirst = 265
xinc = 0.05
xinc = 0.09
yfirst = -24
yinc = 0.05
yinc = 0.09
ifs_script_0
View file @ b0eda3bf
......@@ -10,6 +10,6 @@ for d in 11; do
#cdo remapbil,grids/latlon_0.10_deg -setgridtype,regular test.nc out.nc
cdo --eccodes select,name=clwc ${dirin}mars_out_ml_moist.240 tmp.grb
cdo -f nc setgridtype,regular tmp.grb test.nc
cdo -f nc -sellonlatbox,-90,-80,-34,-24 -setgridtype,regular tmp.grb test.nc
done
mpas_script
View file @ b0eda3bf
......@@ -3,7 +3,7 @@
dirin=/work/ka1081/DYAMOND/MPAS-7.5km/
for d in 11; do
#cdo setattribute,*@axis="txz" -selname,qc ${dirin}history.2016-08-11_00.00.00.nc temp.nc
cdo setattribute,*@axis="txz" -selname,qc ${dirin}history.2016-08-11_00.00.00.nc temp.nc
#cdo -P 16 -f nc4 -O -fldmean -sellonlatbox,-50,-30,0,20 \
# -setgrid,mpas:/work/ka1081/2019_06_Hackathon_Mainz/falko/MPAS_7.5km_grid.nc temp.nc test.nc
......
utilities.py
View file @ b0eda3bf
......@@ -4,7 +4,7 @@
description: Useful stuff
author: Christoph Heim
date created: 27.06.2019
date changed: 27.06.2019
date changed: 07.07.2019
usage: arguments:
python: 3.5.2
"""
......@@ -56,3 +56,36 @@ class Timer:
str(np.round(100*value/comp_time,n_decimal_perc)) +
'\t%\t' + str(np.round(value/60,n_decimal_min)) + ' \tmin')
def write_grid_file(box, file, dx_km):
rE = 6371
dx_deg = dx_km / rE / np.pi * 180
lat0 = box['lat0']
lat1 = box['lat1']
lon0 = box['lon0']
lon1 = box['lon1']
xsize = (lon1 - lon0)/dx_deg
ysize = (lat1 - lat0)/dx_deg
with open(file, 'w') as f:
f.write('gridtype = lonlat\n')
f.write('xsize = '+str(int(np.ceil(xsize)))+'\n')
f.write('ysize = '+str(int(np.ceil(ysize)))+'\n')
f.write('xname = lon'+'\n')
f.write('xlongname = "longitude"'+'\n')
f.write('xunits = "degrees_east"'+'\n')
f.write('yname = lat'+'\n')
f.write('ylongname = "latitude"'+'\n')
f.write('yunits = "degrees_north"'+'\n')
f.write('xfirst = '+str(lon0)+'\n')
f.write('xinc = '+str(np.round(dx_deg,3))+'\n')
f.write('yfirst = '+str(lat0)+'\n')
f.write('yinc = '+str(np.round(dx_deg,3))+'\n')
if __name__ == '__main__':
box = {'lon0': 265, 'lon1': 281, 'lat0': -24, 'lat1': -14,
'vert0':1,'vert1':22}
file = 'test_grid'
write_grid_file(box, file, dx_km=1)
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