Global Climate Change

Complex surface temperature models have been created to model changes over many years. Here we will read in such a file that models global surface temperature for years 1920 through 2100. We will use temperatures from some years to train a neural network model that accepts temperatures as input and predicts the year. This will be tested on other years.

In [1]:
import numpy as np
import matplotlib.pyplot as plt
import mpl_toolkits.basemap as bm
import neuralnetworksA4 as nn
%matplotlib inline
import glob
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-1-efb2388c0068> in <module>
      1 import numpy as np
      2 import matplotlib.pyplot as plt
----> 3 import mpl_toolkits.basemap as bm
      4 import neuralnetworksA4 as nn
      5 get_ipython().run_line_magic('matplotlib', 'inline')

ModuleNotFoundError: No module named 'mpl_toolkits.basemap'
In [2]:
import netCDF4 as nc
In [69]:
dataset = nc.Dataset('hist_rcp85.TREFHT.001.192001-210012.nc')
dataset.variables.keys()
Out[69]:
odict_keys(['time', 'time_bnds', 'lon', 'lat', 'TREFHT'])
In [70]:
dataset.variables['TREFHT']
Out[70]:
<class 'netCDF4._netCDF4.Variable'>
float32 TREFHT(time, lat, lon)
    long_name: Reference height temperature
    units: K
    cell_methods: time: mean
unlimited dimensions: time
current shape = (181, 192, 288)
filling off
In [71]:
temps = dataset.variables['TREFHT'][:]
temps.shape
Out[71]:
(181, 192, 288)
In [72]:
lons = dataset.variables['lon'][:]
lats = dataset.variables['lat'][:]
lons.shape, lats.shape
Out[72]:
((288,), (192,))
In [73]:
plt.imshow(temps[0, :, :])
Out[73]:
<matplotlib.image.AxesImage at 0x7fb454a6fe10>
In [74]:
def drawOnGlobe(basemap, lats, lons, data, cmap='coolwarm', vmin=None, vmax=None):
    basemap.drawcoastlines()
    # basemap.drawparallels(np.arange(-90.,120.,30.))
    dshifted, newlons = bm.shiftgrid(180., data, lons, start=False)
    lonsmg, latsmg = np.meshgrid(newlons, lats)
    x, y = basemap(lonsmg, latsmg)
    # cs = basemap.contourf(x,y,dshifted)
    # maxmag = np.max(np.abs(data))
    if vmin:
        basemap.pcolormesh(x, y, dshifted, cmap=cmap, vmin=vmin, vmax=vmax)
    else:
        basemap.pcolormesh(x, y, dshifted, cmap=cmap)
    plt.colorbar(fraction=0.03, pad=0.02)
In [75]:
basemap = bm.Basemap(projection='robin', lat_0=-90, lon_0=0.0)
In [76]:
vmin = np.min(temps)
vmax = np.max(temps)
In [77]:
plt.figure(figsize=(20, 10))
plt.subplot(1, 2, 1)
drawOnGlobe(basemap, lats, lons, temps[0, :, :], vmin=vmin, vmax=vmax)
plt.subplot(1, 2, 2)
drawOnGlobe(basemap, lats, lons, temps[180, :, :], vmin=vmin, vmax=vmax)
/s/parsons/e/fac/anderson/anaconda3/lib/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:1623: MatplotlibDeprecationWarning: The get_axis_bgcolor function was deprecated in version 2.0. Use get_facecolor instead.
  fill_color = ax.get_axis_bgcolor()
/s/parsons/e/fac/anderson/anaconda3/lib/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:3413: MatplotlibDeprecationWarning: The ishold function was deprecated in version 2.0.
  b = ax.ishold()
/s/parsons/e/fac/anderson/anaconda3/lib/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:3422: MatplotlibDeprecationWarning: axes.hold is deprecated.
    See the API Changes document (http://matplotlib.org/api/api_changes.html)
    for more details.
  ax.hold(b)
In [104]:
basemap = bm.Basemap(projection='cyl', lat_0=-90, lon_0=0.0)
In [105]:
plt.figure(figsize=(20, 10))
plt.subplot(1, 2, 1)
drawOnGlobe(basemap, lats, lons, temps[0, :, :], vmin=vmin, vmax=vmax)
plt.subplot(1, 2, 2)
drawOnGlobe(basemap, lats, lons, temps[180, :, :], vmin=vmin, vmax=vmax)
/s/parsons/e/fac/anderson/anaconda3/lib/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:3413: MatplotlibDeprecationWarning: The ishold function was deprecated in version 2.0.
  b = ax.ishold()
/s/parsons/e/fac/anderson/anaconda3/lib/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:3422: MatplotlibDeprecationWarning: axes.hold is deprecated.
    See the API Changes document (http://matplotlib.org/api/api_changes.html)
    for more details.
  ax.hold(b)
In [78]:
temps.shape
Out[78]:
(181, 192, 288)
In [79]:
rows = np.arange(180).reshape((-1,10))
rows
Out[79]:
array([[  0,   1,   2,   3,   4,   5,   6,   7,   8,   9],
       [ 10,  11,  12,  13,  14,  15,  16,  17,  18,  19],
       [ 20,  21,  22,  23,  24,  25,  26,  27,  28,  29],
       [ 30,  31,  32,  33,  34,  35,  36,  37,  38,  39],
       [ 40,  41,  42,  43,  44,  45,  46,  47,  48,  49],
       [ 50,  51,  52,  53,  54,  55,  56,  57,  58,  59],
       [ 60,  61,  62,  63,  64,  65,  66,  67,  68,  69],
       [ 70,  71,  72,  73,  74,  75,  76,  77,  78,  79],
       [ 80,  81,  82,  83,  84,  85,  86,  87,  88,  89],
       [ 90,  91,  92,  93,  94,  95,  96,  97,  98,  99],
       [100, 101, 102, 103, 104, 105, 106, 107, 108, 109],
       [110, 111, 112, 113, 114, 115, 116, 117, 118, 119],
       [120, 121, 122, 123, 124, 125, 126, 127, 128, 129],
       [130, 131, 132, 133, 134, 135, 136, 137, 138, 139],
       [140, 141, 142, 143, 144, 145, 146, 147, 148, 149],
       [150, 151, 152, 153, 154, 155, 156, 157, 158, 159],
       [160, 161, 162, 163, 164, 165, 166, 167, 168, 169],
       [170, 171, 172, 173, 174, 175, 176, 177, 178, 179]])
In [80]:
trainRows = rows[::2,:]
testRows = rows[1::2,:]
trainRows, testRows            
Out[80]:
(array([[  0,   1,   2,   3,   4,   5,   6,   7,   8,   9],
        [ 20,  21,  22,  23,  24,  25,  26,  27,  28,  29],
        [ 40,  41,  42,  43,  44,  45,  46,  47,  48,  49],
        [ 60,  61,  62,  63,  64,  65,  66,  67,  68,  69],
        [ 80,  81,  82,  83,  84,  85,  86,  87,  88,  89],
        [100, 101, 102, 103, 104, 105, 106, 107, 108, 109],
        [120, 121, 122, 123, 124, 125, 126, 127, 128, 129],
        [140, 141, 142, 143, 144, 145, 146, 147, 148, 149],
        [160, 161, 162, 163, 164, 165, 166, 167, 168, 169]]),
 array([[ 10,  11,  12,  13,  14,  15,  16,  17,  18,  19],
        [ 30,  31,  32,  33,  34,  35,  36,  37,  38,  39],
        [ 50,  51,  52,  53,  54,  55,  56,  57,  58,  59],
        [ 70,  71,  72,  73,  74,  75,  76,  77,  78,  79],
        [ 90,  91,  92,  93,  94,  95,  96,  97,  98,  99],
        [110, 111, 112, 113, 114, 115, 116, 117, 118, 119],
        [130, 131, 132, 133, 134, 135, 136, 137, 138, 139],
        [150, 151, 152, 153, 154, 155, 156, 157, 158, 159],
        [170, 171, 172, 173, 174, 175, 176, 177, 178, 179]]))
In [81]:
trainRows = trainRows.reshape(-1)
testRows = testRows.reshape(-1)
In [82]:
testRows = np.array(testRows.tolist() + [180])
trainRows, testRows
Out[82]:
(array([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  20,  21,  22,
         23,  24,  25,  26,  27,  28,  29,  40,  41,  42,  43,  44,  45,
         46,  47,  48,  49,  60,  61,  62,  63,  64,  65,  66,  67,  68,
         69,  80,  81,  82,  83,  84,  85,  86,  87,  88,  89, 100, 101,
        102, 103, 104, 105, 106, 107, 108, 109, 120, 121, 122, 123, 124,
        125, 126, 127, 128, 129, 140, 141, 142, 143, 144, 145, 146, 147,
        148, 149, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169]),
 array([ 10,  11,  12,  13,  14,  15,  16,  17,  18,  19,  30,  31,  32,
         33,  34,  35,  36,  37,  38,  39,  50,  51,  52,  53,  54,  55,
         56,  57,  58,  59,  70,  71,  72,  73,  74,  75,  76,  77,  78,
         79,  90,  91,  92,  93,  94,  95,  96,  97,  98,  99, 110, 111,
        112, 113, 114, 115, 116, 117, 118, 119, 130, 131, 132, 133, 134,
        135, 136, 137, 138, 139, 150, 151, 152, 153, 154, 155, 156, 157,
        158, 159, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180]))
In [83]:
Xtrain = temps[trainRows, :, :].reshape((len(trainRows), -1))
Ttrain = (trainRows + 1920).reshape((-1, 1))
Xtest = temps[testRows, :, :].reshape((len(testRows), -1))
Ttest = (testRows + 1920).reshape((-1, 1))
Xtrain.shape, Ttrain.shape, Xtest.shape, Ttest.shape
Out[83]:
((90, 55296), (90, 1), (91, 55296), (91, 1))
In [84]:
def rmse(A, B):
    return np.sqrt(np.mean((A - B)**2))
In [91]:
nnet = nn.NeuralNetwork(Xtrain.shape[1], [5], 1)
In [92]:
nnet.train(Xtrain, Ttrain, 200, verbose=True)
nnet
SCG: Iteration 20 fValue Eval 0.0418089956278 Scale 134.217728
SCG: Iteration 40 fValue Eval 0.0134210862166 Scale 0.131072
SCG: Iteration 60 fValue Eval 0.00495722179261 Scale 1e-06
SCG: Iteration 80 fValue Eval 0.00173022090185 Scale 9.5367431640625e-13
SCG: Iteration 100 fValue Eval 0.00108589636888 Scale 1e-15
SCG: Iteration 120 fValue Eval 0.000887136742876 Scale 1e-15
SCG: Iteration 140 fValue Eval 0.000289956375343 Scale 1e-15
SCG: Iteration 160 fValue Eval 0.000185848275202 Scale 1e-15
SCG: Iteration 180 fValue Eval 0.00014925721939 Scale 1e-15
SCG: Iteration 200 fValue Eval 0.000125005590987 Scale 1e-15
Out[92]:
NeuralNetwork(55296, [5], 1)
   Network was trained for 201 iterations that took 23.9856 seconds. Final error is 0.011180589921250452.
In [93]:
plt.plot(nnet.getErrors())
Out[93]:
[<matplotlib.lines.Line2D at 0x7fb454784198>]
In [94]:
Ytrain = nnet.use(Xtrain)
Ytest = nnet.use(Xtest)
rmse(Ytrain, Ttrain), rmse(Ytest, Ttest)
Out[94]:
(0.81777691182631362, 16.26077973593193)
In [95]:
plt.subplot(2, 1, 1)
plt.plot(Ttrain)
plt.plot(Ytrain)
plt.subplot(2, 1, 2)
plt.plot(Ttest)
plt.plot(Ytest)
Out[95]:
[<matplotlib.lines.Line2D at 0x7fb454727208>]
In [96]:
for hiddens in [[nu]*nl for nu in [1, 2, 5, 10, 20] for nl in [1, 2, 3, 4, 5, 10]]:
    nnet = nn.NeuralNetwork(Xtrain.shape[1], hiddens, 1)
    nnet.train(Xtrain, Ttrain, 100)
    print(rmse(Ttrain, nnet.use(Xtrain)), rmse(Ttest, nnet.use(Xtest)), hiddens)
14.3665567449 18.0403176703 [1]
22.5339019716 26.3277565172 [1, 1]
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-96-9795b0538098> in <module>()
      1 for hiddens in [[nu]*nl for nu in [1, 2, 5, 10, 20] for nl in [1, 2, 3, 4, 5, 10]]:
      2     nnet = nn.NeuralNetwork(Xtrain.shape[1], hiddens, 1)
----> 3     nnet.train(Xtrain, Ttrain, 100)
      4     print(rmse(Ttrain, nnet.use(Xtrain)), rmse(Ttest, nnet.use(Xtest)), hiddens)

~/res/atmos/libby/neuralnetworksA4.py in train(self, X, T, nIterations, verbose, weightPrecision, errorPrecision, saveWeightsHistory)
    143                             verbose=verbose,
    144                             ftracep=True,
--> 145                             xtracep=saveWeightsHistory)
    146 
    147         self._unpack(scgresult['x'])

~/res/atmos/libby/mlutilities.py in scg(x, f, gradf, *fargs, **params)
    488                 fold = fnew
    489                 gradold = gradnew
--> 490                 gradnew = gradf(x, *fargs)
    491                 ## If the gradient is zero then we are done.
    492                 if np.dot(gradnew, gradnew) == 0:

~/res/atmos/libby/neuralnetworksA4.py in _gradientF(self, w, X, T)
    108                             Z[Zi-1].T @ delta))
    109             dVs.insert(0, dV)
--> 110             delta = (delta @ self.Vs[Vi][1:, :].T) * (1 - Z[Zi-1]**2)
    111         return self._pack(dVs, dW)
    112 

KeyboardInterrupt: 
In [100]:
hiddens = [20]
nnet = nn.NeuralNetwork(Xtrain.shape[1], hiddens, 1)
nnet.train(Xtrain, Ttrain, 200, verbose=True)
Ytrain = nnet.use(Xtrain)
Ytest = nnet.use(Xtest)
print(rmse(Ttrain, Ytrain), rmse(Ttest, Ytest), hiddens)
SCG: Iteration 20 fValue Eval 0.0291903344948 Scale 16.777216
SCG: Iteration 40 fValue Eval 0.00491355480406 Scale 0.262144
SCG: Iteration 60 fValue Eval 0.00331268684174 Scale 5e-07
SCG: Iteration 80 fValue Eval 0.000253606479732 Scale 4.76837158203125e-13
SCG: Iteration 100 fValue Eval 5.64014842614e-05 Scale 1e-15
SCG: Iteration 120 fValue Eval 1.54152564995e-05 Scale 1e-15
SCG: Iteration 140 fValue Eval 2.76506270074e-06 Scale 1e-15
SCG: Iteration 160 fValue Eval 1.69329507172e-06 Scale 1e-15
SCG: Iteration 180 fValue Eval 6.43660439427e-07 Scale 1e-15
SCG: Iteration 200 fValue Eval 4.0427110885e-07 Scale 1e-15
0.0465057314084 10.1746564552 [20]
In [102]:
plt.plot(nnet.getErrors())
Out[102]:
[<matplotlib.lines.Line2D at 0x7fb4545d5860>]
In [103]:
plt.subplot(2, 1, 1)
plt.plot(Ttrain)
plt.plot(Ytrain)
plt.subplot(2, 1, 2)
plt.plot(Ttest)
plt.plot(Ytest)
Out[103]:
[<matplotlib.lines.Line2D at 0x7fb455044b00>]

How is this network able to estimate the year so well? What parts of the globe are most relevant?

Let's look at the weights in the first layer units, displayed on the globe like we have displayed the temperature data.

In [108]:
plt.figure(figsize=(15, 15))
V = nnet.Vs[0]
nPlots = V.shape[1]
nCols = int(np.sqrt(nPlots))
if nCols * nCols != nPlots:
    nCols += 1
basemap = bm.Basemap(projection='cyl', lat_0=-90, lon_0=0.0)
vmin, vmax = np.min(V[1:, :]), np.max(V[1:, :])
for i in range(nPlots):
    plt.subplot(nCols, nCols, i + 1)
    d = V[1:, i].reshape(len(lats), len(lons))
    drawOnGlobe(basemap, lats, lons, d, cmap='RdYlGn', vmin=vmin, vmax=vmax)
    plt.axis('auto')
    plt.axis('off')
    plt.title('Output W = {:.3f}'.format(nnet.W[1 + i, 0]))
/s/parsons/e/fac/anderson/anaconda3/lib/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:3413: MatplotlibDeprecationWarning: The ishold function was deprecated in version 2.0.
  b = ax.ishold()
/s/parsons/e/fac/anderson/anaconda3/lib/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:3422: MatplotlibDeprecationWarning: axes.hold is deprecated.
    See the API Changes document (http://matplotlib.org/api/api_changes.html)
    for more details.
  ax.hold(b)