Commit 9ff75c2b authored by Modellers Operational's avatar Modellers Operational

Initial commit, only surface plots working

parents
import matplotlib as plt
plt.use('Agg')
import os
import socket
import sys
import numpy as np
import socket
import multiprocessing
from pathlib import Path
from cmocean import cm
import PyFVCOM as pf
def plot_var(idx):
fig, ax = plt.pyplot.subplots(1,2, figsize=(20,15))
DATA = np.squeeze(var_surf_data[varlist[0]][idx, :])
# surface plot
pc = ax[0].tripcolor(X, Y, trinodes, DATA, shading='flat', cmap=cmap_to_use)
for this_cross_section_ind in np.arange(0, len(cross_section_points)-1):
this_cross_section = cross_section_points[this_cross_section_ind:this_cross_section_ind + 2]
ax[0].plot([this_cross_section[0][0], this_cross_section[1][0]], [this_cross_section[0][1], this_cross_section[1][1]], linewidth=1, color='green')
ax[0].set_xlim(plot_lims[0:2])
ax[0].set_ylim(plot_lims[2:])
ax[0].set_title(''.join(Times[idx, :-7].astype(str)).replace('T', ' '))
pc.set_clim(clim[0], clim[1])
cb = fig.colorbar(pc, extend='both')
#cb.set_label(label)
# cross section plot
this_c_plot = pf.plot.CrossPlotter(filereader, figure=fig, axes=ax[1], vmin=clim[0], vmax=clim[1], cmap=cmap_to_use)
this_c_plot.cross_plot_x_pcolor = c_plot.cross_plot_x_pcolor
this_c_plot.cross_plot_y_pcolor = c_plot.cross_plot_y_pcolor
this_c_plot.sel_points = c_plot.sel_points
this_c_plot.sample_points = c_plot.sample_points
this_c_plot.chan_x = c_plot.chan_x
this_c_plot.chan_y = c_plot.chan_y
this_c_plot.xlim_vals = c_plot.xlim_vals
this_c_plot.ylim_vals = c_plot.ylim_vals
this_c_plot.wet_points_data = c_plot.wet_points_data
this_c_plot.sample_points_ind = c_plot.sample_points_ind
if var=='tchl':
var_to_plot=this_c_plot._var_prep('P1_Chl',idx)
var_to_plot+=this_c_plot._var_prep('P2_Chl',idx)
var_to_plot+=this_c_plot._var_prep('P3_Chl',idx)
var_to_plot+=this_c_plot._var_prep('P4_Chl',idx)
this_c_plot.plot_pcolor_field(var_to_plot.T,idx)
else:
this_c_plot.plot_pcolor_field(var, idx)
fig.savefig('{}_{:04d}.png'.format(baseName, idx + plot_no_start), bbox_inches='tight', pad_inches=0.2, dpi=360)
plt.pyplot.close() # to free memory
def nan_extend(in_array):
in_shape = in_array.shape
if len(in_shape) == 3:
nan_ext = np.empty([in_shape[0],in_shape[1],1])
elif len(in_shape)==2:
nan_ext = np.empty([in_shape[0],1])
nan_ext[:] = np.NAN
return np.append(in_array, nan_ext, axis=len(in_shape)-1)
if __name__ == '__main__':
## setup params
# model run and variable settings
fvcom_fname = sys.argv[1]
var = sys.argv[2]
cross_section_file = sys.argv[3]
clim_min = sys.argv[4]
clim_max = sys.argv[5]
varlist = [var]
layer_depth = {'siglay': '[' + str(layer_depth_no) + ']'}
exec('cmap_to_use = cm.' + cmap_str)
# plot area and cross section settings
plot_lims_scratch = np.load('plot_lims.npy')
plot_lims = [np.min(plot_lims_scratch[:,0]), np.max(plot_lims_scratch[:,0]), np.min(plot_lims_scratch[:,1]), np.max(plot_lims_scratch[:,1])]
cross_section_points = np.load('cross_section.npy')
noisy = False
show_range = False # plot the location of extreme values
# retrieve mesh and define the points for cross section
mesh_data = pf.read.FileReader(fvcom_fname)
cross_section_list = []
for i in np.arange(0, len(cross_section_points)-1):
cross_section_list.append(cross_section_points[i:i+2,:])
filereader = pf.read.FileReader(fvcom_fname)
c_plot = pf.plot.CrossPlotter(filereader)
c_plot.cross_section_init(cross_section_list, dist_res=0.1,wetting_and_drying=False)
cross_plot_x_pcolor = c_plot.cross_plot_x_pcolor
cross_plot_y_pcolor = c_plot.cross_plot_y_pcolor
sel_node = c_plot.sel_points
sample_nodes = c_plot.sample_points
sample_nodes_ind = c_plot.sample_points_ind
chan_x = c_plot.chan_x
chan_y = c_plot.chan_y
wet_points_data = c_plot.wet_points_data
# get the variables
if var=='tchl':
out_data = pf.read.FileReader(fvcom_fname, ['P1_Chl','P2_Chl','P3_Chl','P4_Chl'], dims=layer_depth)
var_surf_data = {}
var_surf_data['tchl'] = out_data['P1_Chl'] + out_data['P2_Chl'] + out_data['P3_Chl'] + out_data['P4_Chl']
del out_data
else:
var_surf_data = pf.read.FileReader(fvcom_fname, varlist, dims=layer_depth, noisy=noisy)
# decide plot limits
if clim_min == 'auto':
clim = [np.floor(np.nanmin(var_surf_data)), np.ceil(np.nanmax(var_surf_data))]
else:
clim = [float(clim_min), float(clim_max)]
## run the loop
time_inds = np.arange(len(mesh_data.times.time_dt))
# Launch the parallel plotting and then close the pool ready for the
# next variable.
pool = multiprocessing.Pool(2)
pool.map(plot_var, time_inds)
pool.close()
[command]
default = today_output=/${ARCHIVE_DIR}/${FORECAST_DAY}/${GRID_NAME}_0001.nc; echo ${today_output}
OIFS=$IFS
IFS=","
cross_section_files_arr=(${CROSS_SECTION_FILESTR})
var_arr=(${CROSS_SECTION_VARS})
for this_cross_section_params in "${cross_section_files_arr[@]}"; do
this_cross_section_file=${this_cross_section_params[0]}
this_cross_section_name=${this_cross_section_params[1]}
for this_var_params in "${var_arr[@]}"; do
IFS=";"
this_var_details=(${this_var_params})
this_var=${this_var_params[0]}
this_var_lower_lim=${this_var_params[1]}
this_var_upper_lim=${this_var_params[2]}
echo ${this_var}
echo ${this_var_lower_lim}
echo ${this_var_upper_lim}
python3 plot_cross_section.py ${today_output} ${this_var} ${this_cross_section_file} ${this_cross_section_name}
this_var_plot_dir=${OUTPUT_DIR}/${FORECAST_DAY}/cross_sections/${this_var}
echo ${this_var_plot_dir}
mkdir -p /${this_var_plot_dir}
mv *.png /${this_var_plot_dir}
done
done
IFS=${OIFS}
import matplotlib as mpl
mpl.use('Agg')
import sys
import multiprocessing
import numpy as np
from cmocean import cm
import PyFVCOM as pf
import pmltools as pt
labels = {'q2': 'Turbulent kinetic energy $(m^{2}s^{-2})$',
'l': 'Turbulent macroscale $(m^{3}s^{-2})$',
'q2l': 'Turbulent kinetic\nenergy x turblent\nmacroscale ($cm^{3}s^{-2}$)',
'tke': 'Turbulent kinetic energy $(m^{2}s^{-2})$',
'viscofh': 'Horizontal Turbulent Eddy Viscosity $(m^{2}s^{-1})$',
'teps': 'Turbulent kinetic\nenergy x turblent\nmacroscale ($cm^{3}s^{-2}$)',
'tauc': 'Bed shear stress $(m^{2}s^{-2})$',
'temp': 'Temperature ($\degree C$)',
'salinity': 'Salinity (PSU)',
'zeta': 'Surface elevation (m)',
'uv': 'Speed $(ms^{-1})$',
'uava': 'Depth averaged speed $(ms^{-1})$',
'uvanomaly': 'Speed anomaly $(ms^{-1})$',
'direction': 'Direction $(\degree)$',
'O3_c': 'Carbonate total dissolved\ninorganic carbon $(mmol C/m^3)$',
'O3_pH': 'Carbonate pH',
'O3_TA': 'Total alkalinity $(umol/kg)$',
'O3_fair': 'Carbonate air-sea flux of $CO_{2} (mmol C/m^{2}/d)$',
'volume': 'Node-based control water column volume $(m^{3})$'}
def plot_var(idx):
plot = pf.plot.Plotter(fvcom, figsize=(23, 18), cmap=cmap, cb_label=label, extend=extension, res=None)
plot.plot_field(np.squeeze(getattr(fvcom.data, var))[idx, level, :])
plot.tripcolor_plot.set_clim(clim[0], clim[1])
plot.axes.set_title(fvcom.time.Times[idx][:-7].replace('T', ' '))
suffix = ''
plot.figure.savefig('{}_{:04d}.png'.format(var, idx + 1),
bbox_inches='tight',
pad_inches=0.2,
dpi=120)
plot.close()
fname = sys.argv[1]
var = sys.argv[2]
clim = [float(sys.argv[3]), float(sys.argv[4])]
sub_lowerleft = [float(sys.argv[5]), float(sys.argv[6])]
sub_upperright = [float(sys.argv[7]), float(sys.argv[8])]
print(fname)
cmap = pt.plotting.pmlcmaps(var)
pool_size = 4
fvcom = pf.read.FileReader(fname, [var], subset=np.asarray([sub_lowerleft, sub_upperright]))
label = labels[var]
extension = pt.plotting.colourbar_extension(*clim, getattr(fvcom.data, var).min(), getattr(fvcom.data, var).max())
level = 0
time_indices = range(fvcom.dims.time)
# Launch the parallel plotting and then close the pool ready for the
# next variable.
pool = multiprocessing.Pool(pool_size)
pool.map(plot_var, time_indices)
pool.close()
[command]
default = today_output=/${ARCHIVE_DIR}/${FORECAST_DAY}/${GRID_NAME}_0001.nc; echo ${today_output}
OIFS=$IFS
IFS=","
subset_arr=(${SUBSET_LLS})
var_arr=(${SUBSET_VARS})
IFS=";"
for this_subset_params in "${subset_arr[@]}"; do
subset_name=${this_subset_params[0]}
subset_lon_lower_left=${this_subset_params[1]}
subset_lat_lower_left=${this_subset_params[2]}
subset_lon_upper_right=${this_subset_params[3]}
subset_lat_upper_right=${this_subset_params[4]}
for this_var_params in "${var_arr[@]}"; do
this_var=${this_var_params[0]}
this_var_lower_lim=${this_var_params[1]}
this_var_upper_lim=${this_var_params[2]}
echo ${this_var}
echo ${this_var_lower_lim}
echo ${this_var_upper_lim}
python3 plot_var.py ${today_output} ${this_var} ${this_var_lower_lim} ${this_var_upper_lim} ${subset_lon_lower_left} ${subset_lat_lower_left} ${subset_lon_upper_right} ${subset_lat_upper_right}
this_var_plot_dir=${OUTPUT_DIR}/${FORECAST_DAY}/subset_plots/${subset_name}/${this_var}
echo ${this_var_plot_dir}
mkdir -p /${this_var_plot_dir}
mv *.png /${this_var_plot_dir}
done
done
IFS=${OIFS}
import matplotlib as mpl
mpl.use('Agg')
import sys
import multiprocessing
import numpy as np
from cmocean import cm
import PyFVCOM as pf
import pmltools as pt
labels = {'q2': 'Turbulent kinetic energy $(m^{2}s^{-2})$',
'l': 'Turbulent macroscale $(m^{3}s^{-2})$',
'q2l': 'Turbulent kinetic\nenergy x turblent\nmacroscale ($cm^{3}s^{-2}$)',
'tke': 'Turbulent kinetic energy $(m^{2}s^{-2})$',
'viscofh': 'Horizontal Turbulent Eddy Viscosity $(m^{2}s^{-1})$',
'teps': 'Turbulent kinetic\nenergy x turblent\nmacroscale ($cm^{3}s^{-2}$)',
'tauc': 'Bed shear stress $(m^{2}s^{-2})$',
'temp': 'Temperature ($\degree C$)',
'salinity': 'Salinity (PSU)',
'zeta': 'Surface elevation (m)',
'uv': 'Speed $(ms^{-1})$',
'uava': 'Depth averaged speed $(ms^{-1})$',
'uvanomaly': 'Speed anomaly $(ms^{-1})$',
'direction': 'Direction $(\degree)$',
'O3_c': 'Carbonate total dissolved\ninorganic carbon $(mmol C/m^3)$',
'O3_pH': 'Carbonate pH',
'O3_TA': 'Total alkalinity $(umol/kg)$',
'O3_fair': 'Carbonate air-sea flux of $CO_{2} (mmol C/m^{2}/d)$',
'volume': 'Node-based control water column volume $(m^{3})$'}
def plot_var(idx):
plot = pf.plot.Plotter(fvcom, figsize=(23, 18), cmap=cmap, cb_label=label, extend=extension, res=None)
plot.plot_field(np.squeeze(getattr(fvcom.data, var))[idx, level, :])
plot.tripcolor_plot.set_clim(clim[0], clim[1])
plot.axes.set_title(fvcom.time.Times[idx][:-7].replace('T', ' '))
suffix = ''
plot.figure.savefig('{}_{:04d}.png'.format(var, idx + 1),
bbox_inches='tight',
pad_inches=0.2,
dpi=120)
plot.close()
fname = sys.argv[1]
var = sys.argv[2]
clim = [float(sys.argv[3]), float(sys.argv[4])]
print(fname)
cmap = pt.plotting.pmlcmaps(var)
pool_size = 4
fvcom = pf.read.FileReader(fname, [var])
label = labels[var]
extension = pt.plotting.colourbar_extension(*clim, getattr(fvcom.data, var).min(), getattr(fvcom.data, var).max())
level = 0
time_indices = range(fvcom.dims.time)
# Launch the parallel plotting and then close the pool ready for the
# next variable.
pool = multiprocessing.Pool(pool_size)
pool.map(plot_var, time_indices)
pool.close()
[command]
default = today_output=${PATH_PREFIX}${FVCOM_OUTPUT_DIR}/${FORECAST_DAY}/${FVCOM_GRID_NAME}_0001.nc; echo ${today_output}
OIFS=$IFS
IFS=","
var_arr=(${FULL_MODEL_SURFACE_VARS})
for this_var_params in "${var_arr[@]}"; do
IFS=";"
this_var_details=(${this_var_params})
this_var=${this_var_details[0]}
this_var_lower_lim=${this_var_details[1]}
this_var_upper_lim=${this_var_details[2]}
echo ${this_var}
echo ${this_var_lower_lim}
echo ${this_var_upper_lim}
python3 plot_var.py ${today_output} ${this_var} ${this_var_lower_lim} ${this_var_upper_lim}
this_var_plot_dir=${OUTPUT_DIR}/${FORECAST_DAY}/surface_${this_var}
echo ${this_var_plot_dir}
mkdir -p /${this_var_plot_dir}
mv *.png /${this_var_plot_dir}
done
IFS=${OIFS}
[jinja2:suite.rc]
## Run properties
INITIAL_START_DATE='2019-04-21T00:00:00Z'
FINAL_CYCLE_POINT='NONE'
MAIL_TO='mbe@pml.ac.uk'
NODES=1
USE_CETO=False
REMOTE_USER='modop'
TRIGGER_SUITE='fvcom_tamar'
TRIGGER_TASK='transfer_data'
# MODEL SETUP
FVCOM_GRID_NAME='tamar_v2'
FVCOM_OUTPUT_DIR='data/sthenno1/scratch/modop/Model/FVCOM_tamar/output'
FULL_MODEL_SURFACE_VARS='temp;8;18,salinity;31;36,uv;0;1.5,zeta;-4;4'
SUBSET_VARS=''
SUBSET_LLS=''
CROSS_SECTION_VARS=''
CROSS_SECTION_FILESTR=''
OUTPUT_DIR='data/sthenno1/scratch/modop/Model/FVCOM_tamar/plots'
access-list=mbe
owner=mbe
project=test_suite
sub-project=A
title=hmm_title
#!jinja2
[cylc]
UTC mode = True # Ignore DST
abort if any task fails = False
[scheduling]
initial cycle point = {{INITIAL_START_DATE}}
{%- if FINAL_CYCLE_POINT not in ['NONE','None'] %}
final cycle point = {{FINAL_CYCLE_POINT}}
{%- endif %}
[[special tasks]]
clock-trigger = start_cycle(PT0M)
[[dependencies]]
[[[P1D]]]
graph = """
finish_plot[-P1D]:finish => start_cycle => suite_trigger <{{TRIGGER_SUITE}}::{{TRIGGER_TASK}}> => start_plots
start_plots => surface_plots & subset_surface_plots & cross_section_plots => finish_plot
"""
[runtime]
[[root]]
env-script = eval $(rose task-env --cycle-offset=P1D)
script = rose task-run --verbose
[[[job]]]
execution time limit = PT3H
[[[events]]]
mail events = submission timeout, execution timeout, failed
mail to = {{MAIL_TO}}
submission timeout = P1D
[[[environment]]]
{% if USE_CETO %}
PATH_PREFIX='/pml'
{% else %}
PATH_PREFIX='/'
{% endif %}
FVCOM_GRID_NAME={{FVCOM_GRID_NAME}}
FVCOM_OUTPUT_DIR={{FVCOM_OUTPUT_DIR}}
FULL_MODEL_SURFACE_VARS={{FULL_MODEL_SURFACE_VARS}}
SUBSET_VARS={{SUBSET_VARS}}
SUBSET_LLS={{SUBSET_LLS}}
CROSS_SECTION_VARS={{CROSS_SECTION_VARS}}
CROSS_SECTION_FILESTR={{CROSS_SECTION_FILESTR}}
OUTPUT_DIR={{OUTPUT_DIR}}
FORECAST_DAY=$(rose date --print-format='%Y-%m-%d' $CYLC_TASK_CYCLE_POINT)
FORECAST_YEAR=$(rose date --print-format='%Y' $CYLC_TASK_CYCLE_POINT)
FORECAST_MONTH=$(rose date --print-format='%m' $CYLC_TASK_CYCLE_POINT)
FVCOM_FILE={{FVCOM_OUTPUT_DIR}}/${FORECAST_DAY}/${FVCOM_GRID_NAME}_0001.nc
[[start_cycle]]
script = """
echo 'Starting plotting cycle'
"""
[[start_plots]]
script = """
echo 'Starting plots'
"""
[[suite_trigger]]
script =""
[[[suite state polling]]]
interval = PT10M
max-polls = 1440
[[[job]]]
execution retry delays = 3*PT15M
{% if USE_CETO %}
[[slurm_job]]
[[[job]]]
batch system = slurm
submission polling intervals = PT10S
execution polling intervals = PT10S, PT1M
[[[directives]]]
--nodes = {{NODES}}
--ntasks-per-node=20
--threads-per-core=1
--time=24:00:00
[[[remote]]]
host = login.ceto.npm.ac.uk
owner = {{REMOTE_USER}}
[[remote_job]]
[[[remote]]]
host = login.ceto.npm.ac.uk
owner = {{REMOTE_USER}}
{% else %}
[[slurm_job]]
[[remote_job]]
{% endif %}
[[surface_plots]]
inherit = slurm_job
[[subset_surface_plots]]
inherit = slurm_job
[[cross_section_plots]]
inherit = slurm_job
[[finish_plot]]
script = """
echo 'Finish plotting cycle'
"""
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