For Rotation Students#

Depending on the project, you may want to make maps, scatter plots, or animations using Python or MatLab. Below are some example scripts. You can run your scripts on your local machine or on Compute1 . For help, contact any current doctoral students.

Example Python scripts#

You are encouraged to read the Python Data Science Handbook if you are new to python. Commonly used python libraries inlcude xarray, pandas, numpy, matplotlib and cartopy.

Script to make maps from GCHP output:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
Created on Thu Feb 16 13:02:46 2023
@author: yanyu
edit: Haihui

import matplotlib.pyplot as plt
import as ccrs  # cartopy version must be >=0.19
import xarray as xr
import cartopy.feature as cfeature
import calendar

InDir = 'your/input/file/name'
OutDir = 'your/output/file/name'

for mon in range(1,13):
    ds = xr.open_dataset(InDir + 'your_file_name_{}.nc'.format(mon))'default')
    plt.figure(figsize=(5, 3))
    left = 0.1   # Adjust the left position as needed
    bottom = 0.2   # Adjust the bottom position as needed
    width = 0.8   # Adjust the width as needed
    height = 0.8   # Adjust the height as needed
    ax = plt.axes([left, bottom, width, height], projection=ccrs.Miller())

    ax.set_extent([-140, 160, -60, 60], crs=ccrs.PlateCarree())# World map without Arctic and Antarctic region

    for face in range(6):
        x = ds.corner_lons.isel(nf=face)
        y = ds.corner_lats.isel(nf=face)
        v = ds.PM25.isel(nf=face) #Change species as needed
        im = ax.pcolormesh(x, y, v, transform=ccrs.PlateCarree(), vmin=0, vmax=100)

    month_str = calendar.month_name[mon]
    ax.text(0.45, 0.1, '{}'.format(month_str), fontsize=10, transform=ax.transAxes)

    plt.title('figure title')
    plt.colorbar(im,label="$PM_{2.5}$ concentrations (\u03bcg/m$\mathregular{^3}$)", orientation="horizontal", pad=0.01, fraction=0.040)
    plt.savefig(OutDir + 'Your_figure_name_{}.png'.format(mon), dpi=500)

Script to make scatter plots comparing simulations to observations:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
Created on Thu Feb 16 16:41:56 2023

@author: yanyu
edit: Haihui
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
from scipy import stats, odr

InDir = 'Your/input/path'
OutDir = 'Yout/output/path'

species = pd.Series(['PM25', 'SO4', 'NIT', 'NH4'])
compr = pd.read_csv(InDir + 'Your_file_name.csv')

site =

fig, ax = plt.subplots(figsize=(5, 5))

unique_species = species.unique()
palette = sns.color_palette('Paired', len(unique_species))
color_dict = dict(zip(unique_species, palette))
colors = [color_dict[s] for s in species]

plt.scatter(PM25,PM25_sim,s=10, label=r'PM$_{2.5}$', c=colors[0])
plt.scatter(SO4,SO4_sim,s=10, label='Sulfate',  c=colors[1])
plt.scatter(NIT,NIT_sim,s=10, label='Nitrate',  c=colors[2])
plt.scatter(NH4,NH4_sim,s=10, label='Ammonium',  c=colors[3])

plt.grid(True, alpha=0.2)

# Add the legend
font_legend = {'family': 'Arial', 'size': 10}
plt.legend(loc='lower right', prop=font_legend, frameon=False)

x = PM25 #Observation
y = PM25_sim #Simulation

data = odr.Data(x, y)
odr_obj = odr.ODR(data, odr.unilinear)
output =
slope = output.beta[0]
offset = output.beta[1]
rsq = (stats.linregress(x, y)[2]) ** 2
plt.plot(np.linspace(0, max(y)+10), slope*np.linspace(0, max(y)+10)+offset,
        linestyle='-', color='black')
plt.plot([0, max(y)+13], [0, max(y)+13], color='black', linestyle='--')

if offset > 0:
    t = plt.text(0.02, 0.98, r'$y = {:.2f} x + {:.2f}$''\n''$R^2 = {:.2f}$''\n'.format(slope, offset, rsq)
                    + r'$N$' + f' = {len(x)}', transform=ax.transAxes, va='top', fontsize=10)
    t = plt.text(0.02, 0.98, r'$y = {:.2f} x {:.2f}$''\n''$R^2 = {:.2f}$''\n'.format(slope, offset, rsq)
                    + r'$N$' + f' = {len(x)}', transform=ax.transAxes, va='top', fontsize=10)

font = {'fontname': 'Arial'}
plt.xlabel('Observed  (\u03bcg/m$\mathregular{^3}$)', fontsize=10, fontdict=font)
plt.ylabel('Simulated  (\u03bcg/m$\mathregular{^3}$)', fontsize=10, fontdict=font)

ax.set_xlim(left=0, right=max(y)+13)
ax.set_ylim(bottom=0, top=max(y)+13)
ax.set_xticks(np.arange(0, max(y)+13, 10))
ax.set_yticks(np.arange(0, max(y)+13, 10))
plt.savefig(OutDir + 'Your_figure_name.png', dpi = 500)

Example MatLab scripts#

Matlab provides many powerful toolboxes and built-in functions, which are well documented on its website. You can also use the help function to look up the definition and examples of a built-in function. For example, type help interp2 in your command window to learn about the interp2 funciton, which is commonly used for mapping.

Example script to make maps:

clear % clear variables in workspace, comment out if you don't want to do so
close all % close all figure windows, comment out if you don't want to do so

SimYear = 2019;
InDir = 'Your/Input/Dir';

for Mon = 1:12
fname = sprintf('%s/your_file_name_%d%.2d.mat',InDir,SimYear,Mon);
load (fname,'mapdata_PM25','mapdata_AOD','lat','lon')


% sub figure 1
spec = 'PM2.5'; units = '\mug/m^3'; rng = [0 80];
Label = 'PM_{2.5}';

% sub figure 2
spec = 'AOD'; units = 'unitless'; rng = [0 0.8];
Label = 'AOD';

% sub figure 3
spec = 'ETA'; units = '\mug/m^3'; rng = [0 300];
Label = '\eta';
mapdata = mapdata_PM25./mapdata_AOD;

% Save figure

% below is a user-built function called in the previous codes:

function map = make_submap(Position,mapdata,lat,lon,spec,rng,units,Label)

    fz = 14; % define font size

    map = worldmap([-50 60],[-150 160]);
    surfm(lat, lon, mapdata);
    load coastlines

    % define a color scheme if you don't want to use the default one
    cm = flip(cbrewer('div','RdYlBu',12,'spline'));

    % adding color bar
    cb1=colorbar('vertical', 'fontsize',fz, 'fontweight', 'bold','Ticks',rng(1):(rng(2)-rng(1))/4:rng(2));
    colorbar_label = sprintf('%s (%s)',spec,units);

    % adding label to the figure
    text(max(xlim)-0.1*(max(xlim)-min(xlim)), min(ylim)+0.1*(max(ylim)-min(ylim)),...
            Region,'fontsize',fz,'fontweight','bold','Horiz','right', 'Vert','bottom');

Running your scripts on Compute1#

You can run your python scripts on compute1 interactively or on the background by submitting an interactive or batch job using the geoschem/gcpy:latest docker. Refer to Learn to use Compute1 for more instructions.

To run MatLab on Compute1, refer to Using MatLab.