---

Python Code Version 3.0 for

Climate Mathematics: Theory and Applications - Chapter 9

---

Samuel S.P. Shen and Richard C.J. Somerville

---

This Version 3.0 is authored by Briana Ramirez, edited by Samuel Shen. Liu Yang, Sandra Villamar, and Joaquin Stawsky contributed codes to this version.

Video tutorial for the python code can be found at the following URL: https://www.youtube.com/channel/UC7D9i0kBMzPTHyEhU0h6W9g

This version is based upon the previous version described in the following box.

######################################################################################################################
#This Python Code for Climate Science is written for the book entitled "Climate Mathematics: Theory and Applications"#
#A Cambridge University Press book authored by SSP Shen and RCJ Somerville in July 2019                              #
#The Python codes were based on the R codes written by Samuel Shen Distinguished Professor,                          #
#San Diego State University, USA and were translated from R by Louis Selstad, Stephen Shen,                          #
#Gregori Clarke, and Dakota Newmann and edited by Samuel Shen.                                                       #
######################################################################################################################

#FIRST TIME Python users*****
#These package need to be installed (on the terminal or anaconda interface) before importing them below. 

#Follow this tutorial for package installation before
# https://towardsdatascience.com/importerror-no-module-named-xyz-45e4a5339e1b

#Change your file path to the folder where your downloaded data is stored
#MAC HELP: https://support.apple.com/guide/mac-help/go-directly-to-a-specific-folder-on-mac-mchlp1236/mac
#PC HELP: https://www.sony.com/electronics/support/articles/00015251
import os
# os.chdir("/Users/sshen/climmath/data")
os.chdir('/Users/HP/Documents/sshen/climmath/data')

%%javascript
IPython.OutputArea.prototype._should_scroll = function(lines) {
    return false;
}

#Style Dictionary to standardize plotting scheme between different python scripts 
import matplotlib.pyplot as plt

styledict = {'xtick.labelsize':20,
             'xtick.major.size':9,
             'xtick.major.width':1,
             'ytick.labelsize':20,
             'ytick.major.size':9,
             'ytick.major.width':1,
             'legend.framealpha':0.0,
             'legend.fontsize':15,
             'axes.labelsize':20,
             'axes.titlesize':25,
             'axes.linewidth':2,
             'figure.figsize':(12,8),
             'savefig.format':'jpg'}
plt.rcParams.update(**styledict)

#Function that creates personalized discrete Colormap
import numpy as np 
from matplotlib import cm as cm1
from matplotlib.colors import ListedColormap, to_rgba

def newColMap(colors):
    """
    This function creates a new color map from a list of colors given
    as a parameter. Recommended length of list of colors is at least 6.
    """
    first = np.repeat([to_rgba(colors[0])], 2, axis = 0)
    last = np.repeat([to_rgba(colors[-1])], 2, axis = 0)
    v = cm1.get_cmap('viridis', 16*(len(colors)-2))
    newcolors = v(np.linspace(0, 1, 16*(len(colors)-2)))
    for (i, col) in enumerate(colors[1:-1]):
        newcolors[16*i : 16*(i+1), :] = to_rgba(col)
    return ListedColormap(np.append(np.append(first,newcolors, axis=0), last, axis=0))

---

Chapter 9: Python Graphics for Climate Science

---

import numpy as np 
import pandas as pd
import warnings
import os
import plotly.express as px 
import sympy as sm
import math as m
import netCDF4 as nc

import cartopy 
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter

import matplotlib
from matplotlib import cm as cm1
from matplotlib import pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
from matplotlib.patches import Polygon
from matplotlib import animation as animation
from matplotlib.animation import ArtistAnimation
from mpl_toolkits.basemap import Basemap, cm, shiftgrid, addcyclic

from datetime import datetime

from scipy import optimize as opt
from scipy.ndimage import gaussian_filter
import time
warnings.filterwarnings("ignore")

---

9.1 Two-Dimensional Line Plots and Setups of Margins and Labels

---

Two Different Time Series on the Same Plot

#Set up Variables 
time = np.arange(2001,2011)
Tmean = np.array([12.06, 11.78, 11.81, 11.72, 12.02,
                  12.36, 12.03, 11.27, 11.33, 11.66])
Prec = np.array([737.11, 737.87, 774.95, 844.55, 764.03,
                 757.43, 741.17, 793.50, 820.42, 796.80])

#Fig 9.1
# set up figure
fig, ax = plt.subplots(1,1,figsize=(12, 8))

# plot data 
ax.plot(time, Tmean, 'o-r', label="$T_{mean}$")

# set labels 
ax.set_ylabel("$T_{mean}$ [$\degree$C]", size=18)
ax.set_xlabel("Year", size=18)
ax.set_title("Contiguous U.S. Annual Mean Temperature and Total Precipitation", size=18, fontweight="bold")
ax.grid()
plt.yticks(rotation=90, size=18)
plt.xticks(size=18)

# creates twin axis 
ax1 = ax.twinx()

# plot of twin axis 
ax1.plot(time, Prec, 'o-b', label="Precipitation")
ax1.set_ylabel("Precipitation [mm]")

plt.yticks(rotation=90, size=18)

# create legend 
hand1, lab1 = ax.get_legend_handles_labels()
hand2, lab2 = ax1.get_legend_handles_labels()

ax.legend(handles=hand1+hand2, labels=lab1+lab2,
          loc='upper left', prop={'size': 20});

ax1.tick_params(direction='out')
ax.tick_params(direction='out')

# show plot
plt.show()

Figure Setups: Margins, Fonts, Mathematical Symbols, and More

#Fig 9.2
x = 0.25*np.arange(-30, 31)

# find constraint on variables
yPositive = np.where(np.sin(x) > 0)[0]
yNegative = np.where(np.sin(x) < 0)[0]

# set up figure
fig, ax = plt.subplots(figsize=(14,8))
plt.rcParams['hatch.color'] = 'blue'

# create bar plot 
ax.bar(x[yPositive], np.sin(x[yPositive]), width=.1, color='red')
ax.bar(x[yNegative], np.sin(x[yNegative]), hatch='.', fc='white', color="blue")

# add annotations
ax.text(-5, .7, "$Sine\ Waves$", size=75, color='green')
ax.set_xlabel(r"Angle in Radians: $\theta - \phi_0$", color='red',
              labelpad=-75)
ax.set_xticks([2*i for i in range(-3, 4)])
ax.set_xlim(x.min(), x.max())
ax.set_ylabel(r"$y = \sin(\theta - \hat{\phi})$", color='blue')
ax.set_yticks([-1., 0., .5, 1.])
ax.spines['bottom'].set_position(('data',-1.4))
ax.spines['bottom'].set_bounds(-6,6)

# twin axis 
ax1 = ax.twinx()

# plot on twin axis 
ax1.plot(x, np.sin(x), alpha=0)

# set labels
ax1.set_yticks([-1., -.5, 0., .5, 1.])
ax1.set_yticklabels(['A', 'B', 'C', 'D', 'E'])
ax1.tick_params(width=1, length=9)

fig.text(0.28,.95,"Text outside of the figure on side 3", size=25, fontweight="bold")

# show plot 
fig.show()

#Fig 9.3
#Illustrating how to adjust font size, axis labels space, and margins.

# create Random Variable
randVals = np.random.standard_normal(200)
x = np.linspace(0, 10, randVals.size)

# set up Figure 
fig, ax = plt.subplots(1, figsize=(12,8))
fig.suptitle("Normal Random Values", fontsize=25, y=.93)

# plot data 
ax.scatter(x, randVals, alpha=1, facecolors='none', edgecolors='k')

# add plot labels 
ax.set_title("Subtitle: 200 Random Values", y=-.2, fontsize=20)
ax.set_ylabel("Random Values")
plt.yticks(rotation=90)

ax.set_yticks([i for i in range(-2, 3)])
ax.set_xlabel("Time")

# show plot and grid
plt.grid()
fig.show()

# read in data 
NOAATemp = np.loadtxt("aravg.ann.land_ocean.90S.90N.v4.0.1.2016.txt",skiprows=0)

# set up variables
x=NOAATemp[:,0]
y=NOAATemp[:,1]
z = [-99] * x.size

def forz(z,y):
    for i in range(2,134):
        
        rslt = [(y[i-2],y[i-1],y[i],y[i+1],y[i+2])]
        z[i] = np.mean(rslt)
    return z

zz = forz(z,y)

def n1func():
    n1 = []
    for i in range(0,137):
        if y[i] >= 0:
            n1.append(i)
    return n1

def n2func():
    n2 = []
    for i in range(0,137):
        if y[i] < 0:
            n2.append(i)
    return n2

n1 = n1func()
n2 = n2func()

x1= x[n1]
y1= y[n1]
x2 = x[n2]
y2 = y[n2]
x3 = x[2:134]
y3 = z[2:134]

#Fig 9.4
# fancy plot of the NOAAGlobalTemp time series
plt.bar(x1,y1,color='r', width=0.5)
plt.bar(x2,y2,color='b', width=0.5)
plt.plot(x3,y3,color='k')
plt.xlim(1877,2020)

plt.xlabel('Time')
plt.ylabel('Temperature [$\degree$C]')
plt.title('NOAA Global Average Annual Mean Temperature Anomalies', fontweight="bold", size=20)
plt.grid()
plt.show()

#Fig 9.5
# Contiguous United States (a) annual mean temperature and (b) annual total precipitation 

# plot the US temp and prec time series on the same figure 
fig, ax = plt.subplots(2,1, figsize=(10,8))

# add additional axis 
ax[0].axes.get_xaxis().set_visible(False)

# plot first subplot data
ax[0].plot(time, Tmean, 'ro-', alpha=.6, label="(a)")

ax[0].set_yticks([11.4, np.mean([11.4, 11.8]), 11.8
                  , np.mean([11.8, 12.2]), 12.2])
ax[0].set_yticklabels([11.4, None, 11.8, None, 12.2])
ax[0].set_ylabel("$T_{mean}$ [$\degree$C]")
ax[0].text(2004,12.0,"US Annual Mean Temperature"
           , fontsize = "xx-large", fontweight="bold")

ax[0].text(2001,12.3,"(a)", size=18)

# plot second subplot data 
ax[1].plot(time, Prec, 'bo-', alpha=.6, label="(b)")

ax[1].text(2004,810,"US Annual Total Precipitation"
           , fontsize = "xx-large", fontweight="bold")
plt.ylabel("Precipitation [mm]")
ax[1].set_yticks([740, 760, 780, 800, 820, 840])
ax[1].set_yticklabels([740, None, 780, None, 820, None])

#set labels 
plt.text(2001,840,"(b)", size=18)

# ensures no overlap in plots 
fig.tight_layout(pad=-1.5)

# show plot
fig.show()

---

9.2 Color Contour Maps

---

Basic Principals for a Python Contour Plot

#9.2.1
# set up variables
x, y = np.linspace(-1, 1, 25), np.linspace(-1, 1, 25)
z = np.random.standard_normal(size=(25,25))

# set up figure 
fig, ax = plt.subplots(1,figsize=(12,12))

# plot contour plot
mymap = ax.contourf(x,y,z, cmap=cm1.get_cmap('jet'))

# create and set labels
ticklabs = [-1., None, -.5, None, 0., None, .5, None, 1.]
ax.set_yticklabels(ticklabs)
ax.set_xticklabels(ticklabs)

ax.set_ylabel("$y$")
ax.set_xlabel("$x$")

ax.set_title("Filled Contour Plot of Normal Random Values", size = 22)
fig.colorbar(mymap)

# show plot
fig.show()

Plot Contour Color Maps for Random Values on a Map

#Fig 9.6
# plot a 5-by-5 grid of global map of standard normal random values
# set up variables
lat = np.linspace(-90, 90, 37)
lon = np.linspace(0, 360, 72)

# create data to be plotted 
mapmat = np.random.standard_normal(size=(lat.size, lon.size))

levels = np.linspace(-3, 3, 81)

# set up figure 
fig, ax = plt.subplots(figsize=(14, 7))

ax = plt.subplot(1,1,1, projection=cartopy.crs.PlateCarree(central_longitude=180))
ax.coastlines()
ax.add_feature(cartopy.feature.BORDERS)
ax.set_global()
ax.set_title("Standard Normal Random Values")

# create desired Colormap 
colors = ["darkred","firebrick","firebrick","indianred","lightcoral","lightcoral","lightpink"
          ,"bisque","moccasin","yellow","yellow","greenyellow","yellowgreen","limegreen","limegreen"
          ,"mediumspringgreen","palegreen","honeydew","lavender","mediumpurple",
         "blueviolet","blueviolet","mediumorchid","mediumvioletred","darkmagenta","indigo","black"]
colors.reverse()
myColMap = newColMap(colors)

# plot data
contf = ax.contourf(lon-180, lat, mapmat, levels, cmap=myColMap)

# add labels and colorbar
colbar = plt.colorbar(contf, drawedges=True, aspect=12)
colbar.set_ticks([i for i in range(-3, 4)])
ax.set_aspect('auto')

ax.set_xticks([0, 50, 100, 150, 200, 250, 300, 350]
              ,crs=cartopy.crs.PlateCarree())
ax.set_xticklabels([0, 50, 100, 150, 200, 250, 300, 350])
ax.set_xlabel("Longitude")

ax.set_yticks([-50, 0, 50], crs= cartopy.crs.PlateCarree())
ax.set_ylabel("Latitude")

# show plot 
fig.show()

#Fig 9.7
# plot a 5-by-5 grid regional map to cover USA and Canada 
# set up varibles 
lat3 = np.linspace(10, 70, 13)
lon3 = np.linspace(230, 295, 14)
mapdat = np.random.standard_normal(size=(lat3.size,lon3.size))

# set up figure
fig, ax = plt.subplots(figsize=(14,10))

ax = plt.subplot(111, projection=cartopy.crs.PlateCarree())
ax.set_extent([230, 295, 10, 65])
ax.coastlines()
ax.add_feature(cartopy.feature.BORDERS)
ax.set_title("Standard Normal Random Values")

# plot data 
contf = ax.contourf(lon3, lat3, mapdat, levels, cmap=myColMap)

# add colorbar and labels 
colbar = plt.colorbar(contf, drawedges=True, aspect=12)
colbar.set_ticks([i for i in range(-3, 4)])
ax.set_aspect('auto')

ax.set_yticks([10*i for i in range(1, 8)],crs=cartopy.crs.PlateCarree())
ax.set_ylabel("Latitude")

ax.set_xticks([10*i for i in range(-13, -6)], crs=cartopy.crs.PlateCarree())
ax.set_xticklabels([10*i for i in range(23,30)])
ax.set_xlabel("Longitude")

ax.tick_params(length=9)

# show plot 
fig.show()

9.2.3 Plot Contour Maps from Climate Model Data in NetCDF Files

# Import Data
datamat = nc.Dataset("air.mon.mean.nc")

# Preparing data for plotting 
Lon = datamat.variables['lon'][:]
Lat = datamat.variables['lat'][:]
Time = datamat.variables['time']
precnc = datamat.variables['air']

#Fig 9.8
# set up figure
fig, ax = plt.subplots(1, figsize=(12,8))

# plot data 
ax.plot(np.linspace(90, -90, precnc.shape[1]), precnc[0, :, 64], '-k',
       linewidth=4)

# set labels and ticks
ax.set_ylabel("Temperature [$\degree$C]")
ax.set_yticks([10*i for i in range(-3, 3)])
ax.set_xlabel("Latitude")
ax.set_xticks([-50, 0, 50])
ax.set_title("90S-90N Temperature [$\degree$C] \nalong a Meridional Line at 160$\degree$E: January 1948",
             size = 22, pad = 10)
ax.tick_params(length=9, width=1)
ax.grid()

# show plot 
fig.show()

#Fig 9.9(a)
#Compute and plot climatology and standard deviation from Jan 1984 to Dec 2015

# set up variables for map
# number of colors that correspond to data values
contour_levels = np.linspace(-50,50,81)

# create color palette
myColMap = LinearSegmentedColormap.from_list(name='my_list', colors=['black','blue','darkgreen',
                                            'green','white','yellow','pink','red','maroon'], N=256)

JMon = precnc[12*np.arange(67)]
sdmat = np.std(precnc, axis=0)
climmat = np.mean(JMon, axis=0)
levels1 = np.linspace(-50, 50, 81)
levels2 = np.linspace(0, 20, 81)

# set up figure
fig, ax = plt.subplots(figsize=(14,7))

ax = plt.subplot(111, projection=cartopy.crs.PlateCarree(central_longitude=180))
ax.coastlines()
ax.add_feature(cartopy.feature.BORDERS, edgecolor='black')
ax.set_title("NCEP RA 1948-2015 January Climatology [$\degree$C]", size = 22, pad = 10)

# plot data
contf = ax.contourf(Lon-180, Lat, climmat, levels1, cmap=myColMap)

# add colorbar and labels
colbar = plt.colorbar(contf, drawedges=True, aspect=12, ticks=[])
colbar.set_ticks([20*i for i in range(-2, 3)])
ax.set_aspect('auto')

ax.set_yticks([50*i for i in range(-1, 2)],crs=cartopy.crs.PlateCarree())
ax.set_ylabel("Latitude")

ax.set_xticks([50*i for i in range(8)], crs=cartopy.crs.PlateCarree())
ax.set_xticklabels([50*i for i in range(8)])
ax.set_xlabel("Longitude")

ax.tick_params(length=9, width=1)
plt.text(198, 93, "[$\degree$C]", size=20)

# show plot 
fig.show()

#Fig 9.9(b)
# plot the standard deviation
# set up figure 
fig, ax = plt.subplots(figsize=(14,7))

ax = plt.subplot(111, projection=cartopy.crs.PlateCarree(central_longitude=180))
ax.coastlines()
ax.add_feature(cartopy.feature.BORDERS, edgecolor='black')
ax.set_title("NCEP 1948-2015 Jan. SAT RA Standard Deviation [$\degree$C]",size = 22, pad = 10)

# plot data
contf = ax.contourf(Lon-180, Lat, sdmat, levels2, cmap=myColMap)

# add colorbar and labels
colbar = plt.colorbar(contf, drawedges=True, aspect=12)
colbar.set_ticks([5*i for i in range(5)])
ax.set_aspect('auto')

ax.set_yticks([50*i for i in range(-1, 2)],crs=cartopy.crs.PlateCarree())
ax.set_ylabel("Latitude")

ax.set_xticks([50*i for i in range(8)], crs=cartopy.crs.PlateCarree())
ax.set_xticklabels([50*i for i in range(8)])
ax.set_xlabel("Longitude")

ax.tick_params(length=9, width=1)
plt.text(198, 93, "[$\degree$C]", size=20)

# show and save figure 
fig.savefig("CH9;JanuarySTD.jpg", bbox_inches='tight');

Plot Data for Displaying Climate Features

#Fit 9.10
# plot the January 1983 temperature using the above setup
mapmat = climmat

# set up figure 
fig, ax = plt.subplots(figsize=(14,7))

ax = plt.subplot(111, projection=cartopy.crs.PlateCarree(central_longitude=180))
ax.coastlines()
ax.add_feature(cartopy.feature.BORDERS, edgecolor='black')
ax.set_title("January 1983 Surface Air Temperature [$\degree$C]", size = 22, pad = 10)

# plot data 
contf = ax.contourf(Lon-180, Lat, mapmat, levels1, cmap=myColMap)

# add colorbar and labels
colbar = plt.colorbar(contf, drawedges=True, aspect=12)
colbar.set_ticks([20*i for i in range(-2, 3)])
ax.set_aspect('auto')

ax.set_yticks([50*i for i in range(-1, 2)],crs=cartopy.crs.PlateCarree())
ax.set_ylabel("Latitude")

ax.set_xticks([50*i for i in range(8)], crs=cartopy.crs.PlateCarree())
ax.set_xticklabels([50*i for i in range(8)])
ax.set_xlabel("Longitude")

plt.text(198,93,"[$\degree$C]", size=20)

ax.tick_params(length=9, width=1)

# show figure 
fig.show()

#Fig 9.11
# NOTE: the color distribution compared to the figures in the textbook look slightly different
# due to differences between R and Python ColorMapping characteristics

# plot the January 1983 anomaly from NCEP data
# create variables
levels1 = np.linspace(-6, 6, 81)
anomat = precnc[420,:,:] - climmat

## suppress values into [-6,6] range
# for i in range(73):
#     for j in range(144):
#         if anomat[i,j] > 6:
#             anomat[i,j] = 6
#         if anomat[i,j] < -6:
#             anomat[i,j] = -6

## another method to suppress anomat values into [-6,6]
anomat = np.maximum(np.minimum(anomat, 6), -6) 

# set up figure            
fig, ax = plt.subplots(figsize=(14,7))

ax = plt.subplot(111, projection=cartopy.crs.PlateCarree(central_longitude=180))
ax.coastlines()
ax.add_feature(cartopy.feature.BORDERS, edgecolor='black')
ax.set_title("January 1983 Surface Air Temperature Anomaly [$\degree$C]", size = 22, pad = 10)

# plot data
contf = ax.contourf(Lon-180, Lat, anomat, levels1, cmap=myColMap)

# add colorbar and labels 
colbar = plt.colorbar(contf, drawedges=True, aspect=12)
colbar.set_ticks([2*i for i in range(-4, 4)])
ax.set_aspect('auto')

ax.set_yticks([50*i for i in range(-1, 2)],crs=cartopy.crs.PlateCarree())
ax.set_ylabel("Latitude")

ax.set_xticks([50*i for i in range(8)], crs=cartopy.crs.PlateCarree())
ax.set_xticklabels([50*i for i in range(8)])
ax.set_xlabel("Longitude")

plt.text(198,93,"[$\degree$C]", size=20)
ax.tick_params(length=9, width=1)

# show figure 
fig.show()

Plot the January 1983 Temperature Anomaly from NCEP data

#Fig 9.12(a)
# create and prepare variables 
myslice = range(144)
anomat = anomat[:,myslice]
int1 = np.linspace(-6,6,81)

# set up figure
fig, ax = plt.subplots(figsize=(14,6))

ax = plt.subplot(111, projection=cartopy.crs.PlateCarree(central_longitude=200))
ax.coastlines()
ax.add_feature(cartopy.feature.BORDERS, edgecolor='black')
ax.set_title("January 1983 Surface Air Temperature Anomaly [$\degree$C]", size = 22, pad = 10)

# plot contour 
contf = ax.contourf(Lon-200, Lat, anomat, int1, cmap=myColMap) 
# subtract longitude by 200 to adjust for perspective change when zooming in

# add colorbar and labels
colbar = plt.colorbar(contf, drawedges=True, aspect=12)
colbar.set_ticks([2*i for i in range(-3, 4)])
ax.set_aspect('auto')

ax.set_yticks([20*i for i in range(-3, 3)],crs=cartopy.crs.PlateCarree())
ax.set_ylabel("Latitude")

ax.set_xticks([50*i for i in range(8)], crs=cartopy.crs.PlateCarree())
ax.set_xticklabels([50*i for i in range(8)])
ax.set_xlabel("Longitude")

ax.set_extent([100,300,50,-50],crs=cartopy.crs.PlateCarree())
ax.tick_params(length=9, width=1)

plt.text(111.5,52,"[$\degree$C]", size=20)

# show figure
fig.show()

#Zoom in to a specific lat-lon region: Pacific

# prepare data 
JMon = list(12*np.arange(68))
climmat = np.zeros((144,73))
for i in range(144):
    for j in range(73):
        climmat[i,j] = np.mean(precnc[JMon,j,i])

anomat = precnc[420,:,:] - climmat.T

# suppress values into the (-5,5) range 
for i in range(73):
    for j in range(144):
        if anomat[i,j] > 5:
            anomat[i,j] = 5
        if anomat[i,j] < -5:
            anomat[i,j] = -5
            
matdiff = anomat
int1 = np.linspace(-5,5,81)

#Fig 9.12(b)
# set up figure 
fig, ax = plt.subplots(figsize=(14,8))

ax = plt.subplot(111, projection=cartopy.crs.PlateCarree(central_longitude=250))
ax.coastlines()
ax.set_title("January 1983 Surface Air Temperature Anomaly [$\degree$C]", size = 22, pad = 10)

# plot data 
contf = ax.contourf(Lon-250, Lat, matdiff, int1, cmap=myColMap)
# subtract longitude by 250 to adjust for perspective change when zooming in

# add colorbar and labels 
colbar = plt.colorbar(contf, drawedges=True, aspect=12)
colbar.set_ticks([2*i for i in range(-3, 4)])
ax.set_aspect('auto')

ax.set_yticks([10*i for i in range(2, 7)],crs=cartopy.crs.PlateCarree())
ax.set_ylabel("Latitude")

ax.set_xticks([200+i*(10) for i in range(3,11)], crs=cartopy.crs.PlateCarree())
ax.set_xticklabels([200+i*(10) for i in range(3,11)])
ax.set_xlabel("Longitude")

ax.set_extent([230,300,20,60],crs=cartopy.crs.PlateCarree())
ax.add_feature(cartopy.feature.BORDERS, edgecolor='black')
ax.tick_params(length=9, width=1)

plt.text(54,61,"[$\degree$C]", size=20)

# show figure 
fig.show()

datamat.close()

---

9.3 Plot Wind Velocity Field on a Map

---

#Fig 9.13
# prepare variables 
lat = np.repeat(np.linspace(-75,75,num = 6),12)
lon = np.repeat(np.linspace(-165,165,num=12),6)

# create wind vector variables
u = np.repeat([-1,1,-1,-1,1,-1],864)
v = np.repeat([1,-1,1,-1,1,-1],864)

X, Y = np.meshgrid(lon, lat)

# set up figure
fig, ax = plt.subplots(figsize=(14,8))

ax = plt.axes(projection=cartopy.crs.PlateCarree())
ax.coastlines()
ax.add_feature(cartopy.feature.BORDERS, edgecolor='black')

# plot every wind vector in a for loop
for i in lon:
    for j in lat:
        plt.plot(i,j, markersize=3, marker='o', color='black')

# create the quiver plot
ax.quiver(X,Y,u,v, scale=65, color="blue")

# add labels and text
# plt.title("Polar High", color = 'green', size = 25)
ax.text(-165,80,"Polar High", color="green", size = 20)
ax.text(-165,30,"Subtropical High", color="green", size = 20)
ax.text(-165,0,"Intertropical Convergence Zone (ITCZ)", color="green", size = 20)
ax.text(-165,-30,"Subtropical High", color="green", size = 20)

ax.axis([-180,180,-90,90])
ax.set_xticks([-165,-105,-45,15,75,135], crs=cartopy.crs.PlateCarree())
ax.set_yticks([-75,-45,-15,15,45,75], crs=cartopy.crs.PlateCarree())

# format lon and lat 
lon_formatter = LongitudeFormatter(degree_symbol='', dateline_direction_label=True)
lat_formatter = LatitudeFormatter(degree_symbol='')
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)

# show figure 
plt.show()

#-- open netcdf file
ncd = nc.Dataset('uvclm95to05.nc')

#-- read variable
u10 = ncd.variables['u'][0,:,:]
v10 = ncd.variables['v'][0,:,:]
lat = ncd.variables['lat'][::-1]
lon = ncd.variables['lon'][:]

#Fig 9.14 
#The NOAA sea wind field of January 1, 1995. 

## THIS FIGURE IS BEING ADJUSTED SEPARATELY AND WILL BE UPDATED SHORTLY ##

ncd.close()

---

9.4 GG-Plot for Data

# create dict of state names and labels
state_names = {"AL":"Alabama","AK":"Alaska","AZ":"Arizona","AR":"Arkansas","CA":"California","CO":"Colorado"
               ,"CT":"Connecticut","DE":"Delaware","FL":"Florida","GA":"Georgia","HI":"Hawaii","ID":"Idaho"
               ,"IL":"Illinois","IN":"Indiana","IA":"Iowa","KS":"Kansas","KY":"Kentucky","LA":"Louisiana"
               ,"ME":"Maine","MD":"Maryland","MA":"Massachusetts","MI":"Michigan","MN":"Minnesota","MS":"Mississippi"
               ,"MO":"Missouri","MT":"Montana","NE":"Nebraska","NV":"Nevada","NH":"New Hampshire"
               ,"NJ":"New Jersey","NM":"New Mexico","NY":"New York","NC":"North Carolina","ND":"North Dakota"
               ,"OH":"Ohio","OK":"Oklahoma","OR":"Oregon","PA":"Pennsylvania","RI":"Rhode Island","SC":"South Carolina"
               ,"SD":"South Dakota","TN":"Tennessee","TX":"Texas","UT":"Utah","VT":"Vermont","VA":"Virginia"
               ,"WA":"Washington","WV":"West Virginia","WI":"Wisconsin","WY":"Wyoming"}

# create state Abbreviations
state_abr = ["AL", "AK", "AZ", "AR", "CA", "CO", "CT", "DC", "DE", "FL", "GA", 
          "HI", "ID", "IL", "IN", "IA", "KS", "KY", "LA", "ME", "MD", 
          "MA", "MI", "MN", "MS", "MO", "MT", "NE", "NV", "NH", "NJ", 
          "NM", "NY", "NC", "ND", "OH", "OK", "OR", "PA", "RI", "SC", 
          "SD", "TN", "TX", "UT", "VT", "VA", "WA", "WV", "WI", "WY"]
state_df = pd.DataFrame({"State":list(state_names.values()),"Abbreviation":list(state_names.keys())})

#Fig 9.15
# set up and create figure 
fig = px.choropleth(state_df,  # Input Pandas DataFrame
                    locations="Abbreviation",  # DataFrame column with locations
                    color="State",  # DataFrame column with color values
                    locationmode = 'USA-states',
                   color_continuous_scale = 'rainbow') # Set to plot as US States

# update layout of interactive plot 
fig.update_layout(
    title_text = 'Color Map of the United States', # Create a Title
    geo_scope='usa',  # Plot only the USA instead of globe
)

# show figure
fig.show()

---

9.5 Animation

---

import turtle
import numpy as np
# animation will open in another window

# window setup
width = 400
height = 500
window = turtle.Screen()
window.setup(width, height)
# window.setup()
window.title('Free Fall Ball')

# ball characteristics
ball = turtle.Turtle()
ball.penup()
ball.color("green")
ball.shape("circle")
# ball.speed('slowest')

g = -9.81 # meters/second^2 = gravitational acceleration
v = 0 # meters/second = initial velocity
y = 490 # meters = initial height
time = np.linspace(0,20.5,100) # 0 to 10 seconds

try:
    while y > 0.0:
        for t in time:
    #         print(t,y)
            v0 = v
            y0 = y
            y = y0 + (v0 * t) + (0.5 * g * t**2)
            v = np.sqrt(v0**2 + 2*g*(y-y0))
            ball.sety(y)
            ball.speed(1 - y/490)
            window.update()

    window.bye()
    # exit window
    
    turtle.done()
except Exception: 
    # so the animation does not end in an error
    pass

# turtle.done()
# uncomment above if program does not stop alone