Skip to content

GitLab

Commit bad25369 authored by v99 Workshop's avatar v99 Workshop
Browse files

Initial commit

parents
Hide whitespace changes
Inline Side-by-side
00_plot.py 0 → 100755
View file @ bad25369
import os
from datetime import datetime, timedelta
import matplotlib.pyplot as plt
import xarray as xr
import numpy as np
from cdo import Cdo
if __name__ == '__main__':
base_path_sh = '/work/ka1081/2019_06_Hackathon_Mainz/shweta'
base_path_ch = '/work/ka1081/2019_06_Hackathon_Mainz/christoph_heim'
models = {
'ICON':{
'bp':base_path_sh,
'res':[80,40,20,10,5,2.5],
#'res':[10,5],
'var_key':{'QV':'q', 'QC':'param212.1.0', 'T':'t'},
'var_file':{'QV':'qv', 'QC':'tot_qc_dia', 'T':'t'},
'vcoord':'height',
},
'NICAM':{
'bp':base_path_ch,
'res':[7,3.5],
'var_key':{'QV':'q', 'QC':'ms_qc', 'T':'ms_tem'},
'var_file':{'QV':'qv', 'QC':'ms_qc', 'T':'ms_tem'},
'vcoord':'lev',
},
'SAM':{
'bp':base_path_ch,
'res':[4],
'var_key':{'QV':'QV', 'QC':'QC', 'T':'TABS'},
'var_file':{'QV':'QV', 'QC':'QC', 'T':'TABS'},
'vcoord':'z',
},
}
cdo = Cdo()
var_name = 'T'
#var_name = 'QC'
#dt0 = datetime(2016,8,11)
#dt1 = datetime(2016,8,12)
dt_range = [datetime(2016,8,11), datetime(2016,8,12),
datetime(2016,8,13), datetime(2016,8,14)]
#dt_range = [datetime(2016,8,13)]
fig,axes = plt.subplots(2,2, figsize=(12,12))
for ti,dt in enumerate(dt_range):
ri = ti % 2
ci = (ti - ri) / 2
handles = []
labels = []
for mi,mk in enumerate(models.keys()):
#mk = 'ICON'
print(mk)
for res in models[mk]['res']:
#res = models[mk]['res'][0]
#res = 80
model_name = mk + '-' + str(res) + 'km'
dt_str = str(dt.year) + str(dt.month) + str(dt.day)
dt_str = '{:%Y%m%d}'.format(dt)
file_name = os.path.join(models[mk]['bp'], 'SCu', 'data',
model_name, dt_str,
models[mk]['var_file'][var_name] + '.nc')
print(file_name)
#field = cdo.fldmean(input=file_name, returnXArray='field')
#dataset = xr.open_dataset(cdo.fldmean(input=file_name,
# options='-f nc'))
dataset = xr.open_dataset(cdo.fldmean(input=file_name,
options='-f nc'))
#dataset[models[mk]['var_key'][var_name]].plot()
# TODO: show to Ralf
#field = dataset[models[mk]['var_key'][var_name]]
#print(field)
#try:
if mk == 'ICON':
#kind = 18
height = dataset.coords[models[mk]['vcoord']].values[-18:]
time = dataset.coords['time'].values[-18:]
values = dataset[models[mk]['var_key'][var_name]].\
mean('time').values.squeeze()[-18:]
else:
height = dataset.coords[models[mk]['vcoord']].values
time = dataset.coords['time'].values
values = dataset[models[mk]['var_key'][var_name]].\
mean('time').values.squeeze()
if mk == 'SAM':
values *= 0.001
dataset.close()
#except:
# pass
## fixes that are not nice..
#print(height)
#print(len(height))
#print(np.arange(-len(height),-1,0))
#quit()
if mk == 'ICON':
ICON_hgt = np.loadtxt('ICON_hgt.txt')
ICON_hgt = ICON_hgt[:-1] + np.diff(ICON_hgt)/2
height = ICON_hgt[-len(height):]
print(len(height))
print(values.shape, height.shape)
label = mk + '_' + str(res)
labels.append(label)
ax = axes[ri,ci]
line, = ax.plot(values,height, label=label)
#line, = ax.plot(1000*values,height, label=label)
handles.append(line)
ax.set_title('{:%Y-%m-%d}'.format(dt))
ax.set_ylim(0,3000)
ax.set_xlim(280,305)
ax.set_xlabel('Temperature [K]')
#ax.set_xlabel('QC [g kg-1]')
ax.set_ylabel('Height [m]')
#fig.legend(handles=handles, labels=labels, bbox_to_achnor=[0.9,0.9])
plt.legend(handles=handles, fontsize=9)
#fig.subplots_adjust(right=0.9)
plt.show()
#plt.savefig('qc_prof.pdf', format='pdf')
#plt.savefig('temp_prof.pdf', format='pdf')
#plt.close(fig.number)
00_plot_t_z.py 0 → 100644
View file @ bad25369
import os
from datetime import datetime, timedelta
import matplotlib.pyplot as plt
import xarray as xr
import numpy as np
from cdo import Cdo
if __name__ == '__main__':
base_path_sh = '/work/ka1081/2019_06_Hackathon_Mainz/shweta'
base_path_ch = '/work/ka1081/2019_06_Hackathon_Mainz/christoph_heim'
models = {
'ICON':{
'bp':base_path_sh,
'res':[80,40,20,10,5],#,2.5],
'res':[5],
'var_key':{'QV':'q', 'QC':'param212.1.0', 'T':'t'},
'var_file':{'QV':'qv', 'QC':'tot_qc_dia', 'T':'t'},
'vcoord':'height',
},
'NICAM':{
'bp':base_path_ch,
'res':[7,3.5],
'var_key':{'QV':'q', 'QC':'ms_qc', 'T':'ms_tem'},
'var_file':{'QV':'qv', 'QC':'ms_qc', 'T':'ms_tem'},
'vcoord':'lev',
},
}
cdo = Cdo()
var_name = 'QC'
#dt0 = datetime(2016,8,11)
#dt1 = datetime(2016,8,12)
dt_range = [datetime(2016,8,11), datetime(2016,8,12),
datetime(2016,8,13)]
dt_range = [datetime(2016,8,11)]
for dt in dt_range:
#fig,ax = plt.subplots(1,2, figsize=(10,10))
fig,ax = plt.subplots(1,1, figsize=(10,10))
handles = []
for mi,mk in enumerate(models.keys()):
mk = 'ICON'
print(mk)
for res in models[mk]['res']:
#res = models[mk]['res'][0]
#res = 80
model_name = mk + '-' + str(res) + 'km'
dt_str = str(dt.year) + str(dt.month) + str(dt.day)
dt_str = '{:%Y%m%d}'.format(dt)
file_name = os.path.join(models[mk]['bp'], 'SCu', 'data',
model_name, dt_str,
models[mk]['var_file'][var_name] + '.nc')
print(file_name)
#field = cdo.fldmean(input=file_name, returnXArray='field')
#dataset = xr.open_dataset(cdo.fldmean(input=file_name,
# options='-f nc'))
dataset = xr.open_dataset(cdo.fldmean(input=file_name,
options='-f nc'))
#dataset[models[mk]['var_key'][var_name]].plot()
# TODO: show to Ralf
#field = dataset[models[mk]['var_key'][var_name]]
#print(field)
#try:
valuesc=[]
if mk == 'ICON':
#kind = 18
height = dataset.coords[models[mk]['vcoord']].values[-18:]
time = dataset.coords['time'].values[-18:]
values = dataset[models[mk]['var_key'][var_name]].\
mean('time').values.squeeze()[-18:]
valuesc.append(values)
else:
height = dataset.coords[models[mk]['vcoord']].values
time = dataset.coords['time'].values
values = dataset[models[mk]['var_key'][var_name]].\
mean('time').values.squeeze()
dataset.close()
#except:
# pass
## fixes that are not nice..
print(height)
print(len(height))
#print(np.arange(-len(height),-1,0))
#quit()
if mk == 'ICON':
ICON_hgt = np.loadtxt('ICON_hgt.txt')
ICON_hgt = ICON_hgt[:-1] + np.diff(ICON_hgt)/2
height = ICON_hgt[-len(height):]
print(len(height))
#print(values)
print(height)
#print(time)
#quit()
print(values.shape, height.shape)
#ax[mi].plot(values,height)
label = mk + '_' + str(res)
line, = ax.contourf(time,height,valuesc, label=label)
# line, = ax.plot(values,height, label=label)
handles.append(line)
#ax.set_ylim(0,3000)
plt.legend(handles=handles)
plt.show()
plt.close(fig.number)
01_nicam.sh 0 → 100755
View file @ bad25369
#!/bin/bash
out_base='/work/ka1081/2019_06_Hackathon_Mainz/christoph_heim/data'
#lon0=-90
#lon1=-80
#lat0=-24
#lat1=-14
lon0=-90
lon1=-89
lat0=-24
lat1=-23
var_names=(ms_qc ms_qv ms_tem)
var_names=(ms_qc)
ress=(3.5 7)
ress=(7)
source ~/config
month=08
days=(11 12 13 14 15 16 17 18 19 20)
#inds=(0 1 2 3 4 5 6 7 8 9)
model_name=NICAM-${res}km
model_path=$store/${box_type}/data/$model_name
mkdir -p $model_path
cd $datad/$model_name
for ind in ${inds[@]} ; do
dt0=2016${months[$ind]}${days[$ind]}
day1=${days[$ind]}
dt1=2016${months[$ind]}$(( $day1 + 1))
echo $dt0
mkdir -p $model_path/${dt0}
cd ${dt0}.000000-${dt1}.000000
cdo -P $njobs -sellonlatbox,-90,-80,-24,-14 \
-sellevidx,1/18 \
${var_name}.nc $model_path/$dt0/${var_name}.nc
cd ../
done
02_sam.py 0 → 100755
View file @ bad25369
"""
author: Christoph Heim
date: 20.06.2019
"""
import os, glob, subprocess
import numpy as np
from datetime import datetime, timedelta
from multiprocessing import Pool
from cdo import Cdo
def run_cdo(inp_file, out_file, dt, var_name):
#print('-sellevidx,1/28 -setdate,'+'{:%Y-%m-%d}'.format(dt)+
# ' -settime,'+'{:%H:%M:%S}'.format(dt)+' '+inp_file)
#print('-sellevidx,1/28 -settaxis,'+'{:%Y-%m-%d,%H:%M:%S,3hour}'.format(dt)+
# ' '+inp_file)
tmp_file = '../SCu/data/SAM_RAW/'+var_name+'_tmp_{:%d%H}'.format(dt)+'.nc'
bash_command = ('ncatted -O -a units,lon,o,c,degrees_east ' +
'-a units,lat,o,c,degrees_north '+inp_file+
' '+tmp_file)
process = subprocess.Popen(bash_command.split(), stdout=subprocess.PIPE)
output, error = process.communicate()
ofile = cdo.sellonlatbox(270,280,-24,-14,
#input=('-sellevidx,1/28 '+inp_file),
#input=('-sellevidx,1/28 -setdate,'+'{:%Y-%m-%d}'.format(dt)+
# ' -settime,'+'{:%H:%M:%S}'.format(dt)+' '+inp_file),
input=('-sellevidx,1/28 -settaxis,'+
'{:%Y-%m-%d,%H:%M:%S,3hour}'.format(dt)+
' '+tmp_file),
output=out_file)
if __name__ == '__main__':
#################
var_names = ['QC', 'QV', 'TABS']
spl_inds = {
'QC':{'start':-16,'end':-6},
'TABS':{'start':-18,'end':-8},
}
model = 'SAM'
res = 4
out_base_dir = os.path.join('/work','ka1081','2019_06_Hackathon_Mainz',
'christoph_heim','SCu','data')
data_dir = os.path.join('/work','ka1081','DYAMOND')
dt = 7.5
base_time = datetime(2016,8,1)
#################
cdo = Cdo()
var_name = 'TABS'
#var_name = 'QC'
model_name = model + '-' + str(res) + 'km'
inp_dir = os.path.join(data_dir, model_name, 'OUT_3D')
model_path = os.path.join(out_base_dir, model + '-' + str(res)+'km')
print(os.path.join(inp_dir, '*_'+var_name))
inp_files_glob = glob.glob(os.path.join(inp_dir, '*_'+var_name+'.nc'))
times = [base_time + timedelta(seconds=dt*int(
f[spl_inds[var_name]['start']:spl_inds[var_name]['end']])) \
for f in inp_files_glob]
#times = np.sort(np.asarray(times))
for ind,inp_file in enumerate(inp_files_glob):
dt = times[ind]
#print('{:%Y%m%d}'.format(dt))
if dt >= datetime(2016,8,11,0) and dt < datetime(2016,8,15,0):
print(inp_file)
print(dt)
print()
out_file = os.path.join(model_path,'{:%Y%m%d}'.format(dt),
var_name+'_{:%H}'.format(dt)+'.nc')
print(out_file)
if not os.path.exists(out_file):
run_cdo(inp_file, out_file, dt, var_name)
#quit()
03_fv3.sh 0 → 100755
View file @ bad25369
#!/bin/bash
export GRIB_DEFINITION_PATH=/mnt/lustre01/sw/rhel6-x64/eccodes/definitions
dirin=/work/ka1081/DYAMOND/FV3-3.25km/2016081100/
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
ICON_hgt.txt 0 → 100755
View file @ bad25369
75000.00000
72363.54688
69842.38281
67357.79688
64946.44531
62606.29688
60335.46484
58132.16797
55976.21484
53877.92969
51824.68359
49826.94922
47890.75000
46014.77734
44197.79297
42438.62500
40736.15234
39089.29688
37497.04688
35958.42578
34472.50781
33038.39844
31655.24805
30322.25000
29038.62109
27803.62305
26617.34961
25488.96289
24416.90820
23408.79688
22460.81445
21569.37500
20731.10742
19942.83789
19201.58594
18504.54492
17849.08203
17232.71289
16653.10742
16108.07422
15595.54883
15113.59473
14660.38672
14234.20898
13821.52832
13421.51855
13021.52441
12621.51465
12221.51660
11821.54004
11421.52832
11021.51660
10621.50391
10221.51758
9821.51562
9421.51562
9021.52734
8621.55371
8221.53125
7821.50977
7421.53271
7021.50977
6621.52734
6221.53662
5821.54248
5421.53809
5033.73242
4659.99219
4300.15576
3954.21582
3622.06934
3303.77051
2999.27954
2708.58911
2431.72412
2168.61621
1919.35303
1683.92688
1462.60999
1255.32886
1062.16406
883.60358
719.54578
570.43378
436.51855
318.25900
216.48891
131.87152
65.69055
20.03601
0.00000
Session.vim 0 → 100755
View file @ bad25369
let SessionLoad = 1
if &cp | set nocp | endif
vnoremap  :norm
map ## 24i#
map co 0i#$
let s:cpo_save=&cpo
set cpo&vim
nmap gx <Plug>NetrwBrowseX
map oo ok0
vmap r "_dP
map un 0x$
nnoremap <silent> <Plug>NetrwBrowseX :call netrw#NetrwBrowseX(expand("<cWORD>"),0)
map <F9> gt
map <F8> gT
imap jk 
let &cpo=s:cpo_save
unlet s:cpo_save
set autowrite
set background=dark
set backspace=indent,eol,start
set directory=~/tmp,/var/tmp,/tmp
set expandtab
set fileencodings=ucs-bom,utf-8,latin1
set guicursor=n-v-c:block,o:hor50,i-ci:hor15,r-cr:hor30,sm:block,a:blinkon0
set helplang=en
set hidden
set hlsearch
set ignorecase
set incsearch
set laststatus=2
set mouse=a
set ruler
set shiftwidth=4
set showcmd
set showmatch
set smartcase
set softtabstop=4
set statusline=%f
set tabstop=4
set viminfo='20,\"50
set wildignore=*.pyc
let s:so_save = &so | let s:siso_save = &siso | set so=0 siso=0
let v:this_session=expand("<sfile>:p")
silent only
cd /mnt/lustre02/work/ka1081/2019_06_Hackathon_Mainz/christoph_heim/scripts
if expand('%') == '' && !&modified && line('$') <= 1 && getline(1) == ''
let s:wipebuf = bufnr('%')
endif
set shortmess=aoO
badd +0 01_nicam.sh
badd +0 02_sam.py
badd +0 03_fv3.sh
argglobal
silent! argdel *
argadd 01_nicam.sh
set stal=2
edit 01_nicam.sh
set splitbelow splitright
set nosplitbelow
set nosplitright
wincmd t
set winheight=1 winwidth=1
argglobal
setlocal keymap=
setlocal noarabic
setlocal noautoindent
setlocal backupcopy=
setlocal nobinary
setlocal nobreakindent
setlocal breakindentopt=
setlocal bufhidden=
setlocal buflisted
setlocal buftype=
setlocal nocindent
setlocal cinkeys=0{,0},0),:,0#,!^F,o,O,e
setlocal cinoptions=
setlocal cinwords=if,else,while,do,for,switch
setlocal colorcolumn=
setlocal comments=s1:/*,mb:*,ex:*/,://,b:#,:%,:XCOMM,n:>,fb:-
setlocal commentstring=#%s
setlocal complete=.,w,b,u,t,i
setlocal concealcursor=
setlocal conceallevel=0
setlocal completefunc=
setlocal nocopyindent
setlocal cryptmethod=
setlocal nocursorbind
setlocal nocursorcolumn
setlocal nocursorline
setlocal define=
setlocal dictionary=
setlocal nodiff
setlocal equalprg=
setlocal errorformat=
setlocal expandtab
if &filetype != 'sh'
setlocal filetype=sh
endif
setlocal foldcolumn=0
setlocal foldenable
setlocal foldexpr=0
setlocal foldignore=#
setlocal foldlevel=0
setlocal foldmarker={{{,}}}
setlocal foldmethod=manual
setlocal foldminlines=1
setlocal foldnestmax=20
setlocal foldtext=foldtext()
setlocal formatexpr=
setlocal formatoptions=tcq
setlocal formatlistpat=^\\s*\\d\\+[\\]:.)}\\t\ ]\\s*
setlocal grepprg=
setlocal iminsert=0
setlocal imsearch=0
setlocal include=
setlocal includeexpr=
setlocal indentexpr=
setlocal indentkeys=0{,0},:,0#,!^F,o,O,e
setlocal noinfercase
setlocal iskeyword=@,48-57,_,192-255,.
setlocal keywordprg=
setlocal nolinebreak
setlocal nolisp
setlocal lispwords=
setlocal nolist
setlocal makeprg=
setlocal matchpairs=(:),{:},[:]
setlocal modeline
setlocal modifiable
setlocal nrformats=octal,hex
set number
setlocal number
setlocal numberwidth=4
setlocal omnifunc=
setlocal path=
setlocal nopreserveindent
setlocal nopreviewwindow