My “~/.gmt/sessions” is filling up substantially and then causing my PBS job to crash. I have a lot files that I’m creating. Please let me know if you have any ideas on how I might optmise my pygmt calls in pygmt_circumpolar_map (below). I’m out of ideas.
Using the following:
1.) pretty irrelevant to GMT, but just so their is awareness the pygmt script is being called by a PBS job
#!/bin/bash
#PBS -N plot-afim-maps2
#PBS -P jk72
#PBS -l walltime=36:00:00
#PBS -q normalbw
#PBS -l mem=256GB
#PBS -l ncpus=28
#PBS -l storage=gdata/cj50+gdata/jk72+gdata/ik11+gdata/hh5+gdata/rt52+gdata/gb6
#PBS -M daniel.atwater@utas.edu.au
module purge
module use /g/data/hh5/public/modules
module load conda/analysis3
python3 ~/AFIM/src/python/PBS_jobs/plot_afim_maps2.py
2.) python script:
import sys, os
sys.path.insert(0, '/home/581/da1339/AFIM/src/python')
import AFIM
import pandas as pd
import xarray as xr
from datetime import timedelta, date, datetime
afim_anal = AFIM.analysis()
dt0_str = "1993-09-30"
dtN_str = "1994-03-01"
var_names = ["aice","hs","uvel","strintx","strairx","strength"]
vmins = [-1,-3,-2,-2,-2,-1e6]
vmaxs = [-x for x in vmins]
cmap = "cmocean/balance"
run_names = ["mocn-d100-cice-d100","wocn-d100-cice-d100","docn-d100-cice-d100"]#"docn-mod-bath-JFprb"
dt_rng = pd.date_range(start=dt0_str, end=dtN_str)
for run_name in run_names:
for var_name in var_names:
for dt_obj in dt_rng:
dt_str = dt_obj.strftime("%Y-%m-%d")
P_ERA5 = afim_anal.ERA5_dict['P_frcg'].format(yr=dt_obj.strftime('%Y'))
dt0_ERA5 = f"{dt_str} 00:00"
dtN_ERA5 = f"{dt_str} 23:00"
cice = afim_anal.load_data(run_name, dt0_str=dt_str, single_file=True, concat_vars=False)
era5 = xr.open_dataset(P_ERA5).sel(time=slice(dt0_ERA5, dtN_ERA5)).mean(dim="time")
afim_anal.pygmt_circumpolar_map(run_name, cice, dt_str, var_name,
show_fig=False,
ow_fig=False,
sic_mask=True,
show_sst=True,
sst_vmin=-2,
sst_vmax=4,
plt_fill_style="c0.15c")
afim_anal.pygmt_regional_map(run_name, cice, dt_str, var_name,
show_fig=False,
ow_fig=False,
sic_mask=True,
show_sst=True,
sst_vmin=-2,
sst_vmax=4,
show_bath=False,
show_coast_shape=False,
da_uatm=era5.wndewd,
da_vatm=era5.wndnwd,
plt_fill_style="c0.15c",
vector_string="v.1c+ea",
vector_scale=50,
vector_skip=75)
3.) uses the following method – I’ve only included pygmt_circumpolar_map as pygmt_regional_map is very similar:
def pygmt_circumpolar_map(self, run_name, ds, dt_str, plot_variable,
sic_mask=False,
plot_fi=False,
da_fi=None,
fia_val=None,
da_uatm=None,
da_vatm=None,
hemisphere=None,
P_save=None,
show_fig=False,
show_sia=False,
sia_val=None,
show_sst=False,
sst_vmin=-2,
sst_vmax=2,
show_bath=False,
bath_cntrs=None,
show_coast_shape=False,
ow_fig=False,
plt_fill_style="c0.2c",
dif_map=False,
cmap=None,
vmin=None,
vmax=None,
var_units_str=None,
vector_string="v1c+bc+ea",
vector_scale=10,
vector_skip=30,
**kwargs):
txt = None
pygmt.config(FONT_HEADING="20p,Courier-Bold")
pygmt.config(PROJ_LENGTH_UNIT="cm")
self.method_init['pygmt_regional_panel_map'] = True
hemisphere = hemisphere if hemisphere is not None else self.hemisphere
hemisphere = self.define_hemisphere(hemisphere)
var_info = self.plot_var_dict[plot_variable]
self.define_run_and_model_name(run_name)
if not var_info["plot"]:
print(f"Variable {plot_variable} set in JSON file to *NOT* plot. Continuing to next variable.")
return
if plot_variable not in ds.variables:
print(f"Variable {plot_variable} does not exist in the dataset")
return
self.define_cice_coords(ds)
G_lon, G_lat = ds[self.CICE_lon_coord_name].values, ds[self.CICE_lat_coord_name].values
G_lon = self.convert_longitudes_to_360(G_lon)
G_lon_flat = G_lon.flatten()
G_lat_flat = G_lat.flatten()
if show_bath:
P_bath = self.anal_dict[run_name]['P_bath']
print(f"loading bathymetry in file: {P_bath}")
bath_flat = xr.open_dataset(P_bath).Bathymetry.values.flatten()
if bath_cntrs is None:
bath_cntrs = np.concatenate([np.arange(0, 500, 100)])
if show_coast_shape:
P_cst = self.P_coast
print(f"loading coastline in file: {P_cst}")
gdf = gpd.read_file(P_cst)
gdf = gdf.to_crs(epsg=4326)
self.pygmt_figure_essentials(var_info, plot_fi=plot_fi, **kwargs)
cmap = cmap if cmap is not None else self.pygmt_cmap
vmin = vmin if vmin is not None else self.pygmt_vmin
vmax = vmax if vmax is not None else self.pygmt_vmax
var_units_str = var_units_str if var_units_str is not None else self.pygmt_units
self.pygmt_dataset_preparation(run_name, ds, var_name=plot_variable, mask=sic_mask, da_fi=da_fi)
if P_save is None:
self.construct_graphical_output_path(dt_str=dt_str, var_name=plot_variable, run_name=run_name, circumpolar=True, fi=plot_fi, dif_map=dif_map)
else:
self.P_fig = os.path.join(P_save, f"{dt_str}_{reg_name}_{plot_variable}.png")
if os.path.exists(self.P_fig) and not ow_fig:
print(f"{self.P_fig} exists 'ow_fig' is set to FALSE")
print(f" WE ARE NOT OVER-WRITING THIS FIGURE")
return
print(f"Creating figure: {self.P_fig}")
fig = pygmt.Figure()
map_scale = f"g180/-85+w250k"
if self.hemisphere=='south':
ext = [-180,180,-90,-55]
else:
ext = [-180,180,55,90]
fig.basemap(region=ext, projection=f"A0/-90/15c", map_scale=map_scale, frame=["x30g10", "y15g5", f"wsne+t{run_name}-{dt_str}"],)
if show_coast_shape:
print(f"plotting coast shape file")
for surface_type, color in zip(['land', 'ice shelf', 'ice tongue', 'rumple'],
['dimgray', 'lightblue', 'powderblue', 'lightsteelblue']):
subset_gdf = gdf[gdf['surface'] == surface_type]
if not subset_gdf.empty:
fig.plot(data=subset_gdf, pen=f"1p,{color}")
else:
fig.coast(shorelines="0.25p,black,solid", resolution='f', water="black", land="grey", region=ext)
if show_sst:
dt_obj = datetime.strptime(dt_str, "%Y-%m-%d").date()
if run_name=='aom2-era5':
P_sst = os.path.join(self.D_dict['AOM2_jk72'],'ERA5','ocean','0p25','daily',f"{dt_obj.strftime('%Y')}_ocean_daily.nc")
sst_flat = xr.open_dataset(P_sst).surface_pot_temp.sel(time=dt_str,method='nearest').values.flatten()-273.15
elif run_name=='aom2-jra55':
P_sst = os.path.join(self.D_dict['AOM2_jk72'],'JRA55','ocean','0p25','daily',f"{dt_obj.strftime('%Y')}_ocean_daily.nc")
sst_flat = xr.open_dataset(P_sst).sst.sel(time=dt_str,method='nearest').values.flatten()-273.15
else:
sst_flat = ds['sst'].values.flatten()
nan_mask = ~np.isnan(sst_flat)
sst_plt = sst_flat[nan_mask]
lon_plt = G_lon_flat[nan_mask]
lat_plt = G_lat_flat[nan_mask]
pygmt.makecpt(cmap="cmocean/thermal", series=[sst_vmin,sst_vmax])
fig.plot(x=lon_plt, y=lat_plt, style=plt_fill_style, cmap=True, fill=sst_plt)
fig.colorbar(position="JMR+o1c/0c+w7c/0.5c", frame=[f"x+lSST", f"y+lC"]) #+n+mc
# MAIN VARIABLE
if dif_map:
pygmt.makecpt(cmap=cmap, series=[vmin, vmax])
else:
pygmt.makecpt(cmap=cmap, series=[vmin, vmax], transparency=vmin)
#transparency = (self.pygmt_da - self.pygmt_da.min()) / (self.pygmt_da.max() - self.pygmt_da.min()) * 100 # Scale to 0-100%
fig.plot(x=self.pygmt_lon, y=self.pygmt_lat, style=plt_fill_style, cmap=True, fill=self.pygmt_da)#, transparency=transparency)
fig.colorbar(position="JBC+w15c/0.5c+h", frame=[f"x+l{plot_variable}", f"y+l{var_units_str}"])
if show_bath:
print(f"plotting bathymetry")
fig.contour(x=G_lon_flat, y=G_lat_flat, z=bath_flat, levels=bath_cntrs, pen="0.25p,black")#, region=reg_cfg['ext'], projection=f"A{reg_cfg['MC']}/-90/15c")
if plot_fi:
txt = f"FIA: {fia_val} ({self.fic_scale/1e6:.1E})km²"
elif show_sia:
txt = f"SIA: {sia_val:2.2f} ({self.sic_scale/1e6:.1E})km²"
if fia_val is not None:
txt = f"{txt} and FIA: {fia_val} ({self.fic_scale/1e6:.1E})km²"
if txt is not None:
fig.text(x=ext[0], y=ext[2], text=txt, font="12p,Helvetica-Bold,black", justify="CB")
if da_uatm is not None and da_vatm is not None:
print(f"plotting winds")
uatm = da_uatm.values*units('m/s')
vatm = da_vatm.values*units('m/s')
wspd = wind_speed(uatm,vatm).magnitude.flatten()[::vector_skip]/vector_scale
wdir = wind_direction(uatm,vatm).magnitude.flatten()[::vector_skip]
wlon = G_lon_flat[::vector_skip]
wlat = G_lat_flat[::vector_skip]
fig.plot(x=wlon, y=wlat, direction=[wdir, wspd], style=vector_string, pen='0.25,black') #pen='0.25,black') #f'v1c+bc+ea+a30'#{vector_scale}c+e+a40+gblack+h0.3+p1p,black
fig.savefig(self.P_fig)
print(f"finished saving figure")
if show_fig:
fig.show()