Home
If you're just getting started, be sure to check README.md.
How do we run this thing? Lets fire up a Python3 script and give it a whirl.
We begin by creating a embm.Model instance and let it spin up by letting the model run for about 208 days ("model time" not 208 days IRL, each time step represents about 1/2 hour in the model, by default):
import numpy as np
import embm
m = embm.Model()
m.step(10000)
This takes within a minute to run on my old laptop.
That's fantastic, but lets get some data and see what's actually being simulated when the model spins up. First, we need to reset the model to its starting conditions:
m.reset()
t_hist, q_hist, p_hist = m.step(10000, trace = True)
Ha! With the trace argument set to True, we're asking the model to give us the global mean of its air temperature, specific humidity and precipitation for each of its 10000 time steps. If you have matplotlib and Basemap installed, we can get some basic plots. If you just have matplotlib and not Basemap, that's okay, you can do much of the same plotting with pylab.pcolor(), pylab.pcolormesh(), etc. You'll just be missing the nice map projections.
Anyways, lets import our modules and plot the data.
import pylab as plt
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
plt.figure()
plt.subplot(2, 1, 1)
plt.plot(p_hist)
plt.grid()
plt.ylabel("Global mean precipitation (m/yr)")
plt.xlabel("Model step")
plt.subplot(2, 1, 2)
plt.plot(t_hist)
plt.grid()
plt.ylabel("Global mean temperature (K)")
plt.xlabel("Model step")
plt.tight_layout()
plt.show()
This will give you a graphic with two subplots that look something like this:
You can see how the model becomes stable after a few thousand steps.
We can also plot the global distribution of air temperature and precipitation across the globe using a Robinson projection, this time pulling data directly from the model instance m.
# Air temperature.
plt.figure()
m = Basemap(projection='robin',lon_0=0,resolution='c')
lonmesh, latmesh = np.meshgrid(m.lon_range, m.lat_range)
x, y = m(lonmesh, latmesh)
m.drawcoastlines()
# draw parallels and meridians.
m.drawparallels(np.arange(-90.,120.,30.))
m.drawmeridians(np.arange(0.,360.,60.))
m.drawmapboundary()
m.pcolormesh(x, y, m.t[1], vmin = 240, vmax = 300)
cb = m.colorbar()
cb.set_label("Air temperature (K)")
plt.title("EMBM modeled")
plt.show()
# Precipitation.
plt.figure()
m = Basemap(projection='robin',lon_0=0,resolution='c')
lonmesh, latmesh = np.meshgrid(m.lon_range, m.lat_range)
x, y = m(lonmesh, latmesh)
m.drawcoastlines()
# draw parallels and meridians.
m.drawparallels(np.arange(-90.,120.,30.))
m.drawmeridians(np.arange(0.,360.,60.))
m.drawmapboundary()
m.pcolormesh(x, y, m.p, cmap = plt.cm.Blues, vmin = 0, vmax = 5)
cb = m.colorbar()
cb.set_label("Precipitation (m/yr)")
plt.title("EMBM modeled")
plt.show()
The air temperature map looks something like this:
while our precipitation map should be along the lines of this:
Interesting, I wonder how well this compares with an "observed" climatology, or a dataset of observations averaged over a 30 year period. I'm going to pull this reanalysis data (NCEP/NCAR R1) in from a netCDF file, this step isn't going to be described here (maybe I'll write that in later). The important thing is that the reanalysis data arrays have the same shape and global grid as our m.p and m.t[1]. obs_precip_lat, obs_precip_lon, and obs_temp_lat, obs_temp_lon are the latitude and longitude for the reanalysis precipitation and temperature data, much like m.lat_range and m.lon_range. The reanalysis data arrays, themselves, are obs_precip and obs_temp.
Lets see what the observation data look like:
plt.figure()
m = Basemap(projection='robin',lon_0=0,resolution='c')
lonmesh, latmesh = np.meshgrid(obs_temp_lon, obs_temp_lat)
x, y = m(lonmesh, latmesh)
m.drawcoastlines()
m.drawparallels(np.arange(-90.,120.,30.))
m.drawmeridians(np.arange(0.,360.,60.))
m.drawmapboundary()
m.pcolormesh(x, y, obs_temp)
cb = m.colorbar()
cb.set_label("Air temperature (K)")
plt.title("Observed climatology")
plt.show()
plt.figure()
m = Basemap(projection='robin',lon_0=0,resolution='c')
lonmesh, latmesh = np.meshgrid(obs_precip_lon, obs_precip_lat)
x, y = m(lonmesh, latmesh)
m.drawcoastlines()
m.drawparallels(np.arange(-90.,120.,30.))
m.drawmeridians(np.arange(0.,360.,60.))
m.drawmapboundary()
m.pcolormesh(x, y, obs_precip)
cb = m.colorbar()
cb.set_label("Precipitation (m/yr)")
plt.title("Observed climatology")
plt.show()
The resulting reanalysis temperature and precipitation maps are below.
Interesting. The temperature looks similar but the precipitation is quite different. Lets get some anomaly maps for each of these variables:
plt.figure()
m = Basemap(projection='robin',lon_0=0,resolution='c')
lonmesh, latmesh = np.meshgrid(obs_temp_lon, obs_temp_lat)
x, y = m(lonmesh, latmesh)
m.drawcoastlines()
m.drawparallels(np.arange(-90.,120.,30.))
m.drawmeridians(np.arange(0.,360.,60.))
m.drawmapboundary()
m.pcolormesh(x, y, m.t[1] - obs_temp, cmap = plt.cm.RdBu_r)
cb = m.colorbar()
cb.set_label("Air temperature anomaly (K)")
plt.title("EMBM modeled - Climatology")
plt.show()
plt.figure()
m = Basemap(projection='robin',lon_0=0,resolution='c')
lonmesh, latmesh = np.meshgrid(obs_precip_lon, obs_precip_lat)
x, y = m(lonmesh, latmesh)
m.drawcoastlines()
m.drawparallels(np.arange(-90.,120.,30.))
m.drawmeridians(np.arange(0.,360.,60.))
m.drawmapboundary()
m.pcolormesh(x, y, m.t[1] - obs_precip, cmap = plt.cm.RdBu_r)
cb = m.colorbar()
cb.set_label("Precipitation anomaly (m/yr)")
plt.title("EMBM modeled - Climatology")
plt.show()
The two plotted maps are here:
Here we're really starting to get an idea of the EMBM's skill and why it might be hitting or missing in some areas. The temperature is off near the Andes, Himalayas and Antarctica. The EMBM is also colder than the reanalysis record over Eurasia. That being said the general pattern for air temperature isn't terrible considering how simple the EMBM is. Precipitation is another story. Precipitation has a reputation for being very difficult to accurately reproduce in models. It's clear from the plot above that the EMBM isn't even close. The simple precipitation relationship in the model does not reliably reproduce the precipitation patterns in the reanalysis data.
Lets get another summary by plotting a zonal profile (i.e. average across the longitudes) of the observation and model data. For temperature:
plt.figure()
plt.subplot(3, 1, 1)
plt.plot(m.lat_range, m.t[1].mean(axis = 1))
plt.ylim([230, 300])
plt.xlabel("Latitude")
plt.ylabel("Air temperature (K)")
plt.grid()
plt.title("EMBM modeled profile")
plt.subplot(3, 1, 2)
plt.plot(obs_temp_lat, obs_temp.mean(axis = 1))
plt.ylim([230, 300])
plt.xlabel("Latitude")
plt.ylabel("Air temperature (K)")
plt.grid()
plt.title("Observed climatology profile")
plt.subplot(3, 1, 3)
plt.plot(m.lat_range, m.t[1].mean(axis = 1) - (obs_temp.mean(axis = 1)))
plt.xlabel("Latitude")
plt.ylabel("Temperature anomaly (K)")
plt.title("EMBM modeled - Climatology anomaly")
plt.grid()
plt.tight_layout()
plt.show()
Which gives us:
This shows that our model has a warm bias in the southern latitudes and a cold bias in the northern latitudes.
The same profile plot but comparing EMBM and reanalysis precipitation:
plt.figure()
plt.subplot(3, 1, 1)
plt.plot(m.lat_range, m.p.mean(axis = 1))
plt.xlabel("Latitude")
plt.ylabel("Precipitation (m/yr)")
plt.ylim((0, 2))
plt.grid()
plt.title("EMBM modeled profile")
plt.subplot(3, 1, 2)
plt.plot(obs_precip_lat, obs_precip)
plt.xlabel("Latitude")
plt.ylabel("Precipitation (m/yr)")
plt.ylim((0, 2))
plt.grid()
plt.title("Observed climatology profile")
plt.subplot(3, 1, 3)
plt.plot(m.lat_range, m.p.mean(axis = 1) - obs_precip)
plt.xlabel("Latitude")
plt.ylabel("Precipitation\nanomaly (m/yr)")
plt.title("EMBM modeled - Climatology anomaly")
plt.grid()
plt.tight_layout()
plt.show()
Which produces this plot:
Again, we see that the EMBM has a precipitation problem. In general, the EMBM under predicts precipitation. The model clearly struggles most near the equator and tropics where we see the most extreme positive and negative anomalies, relative to the reanalysis data.