import pandas as pd
import numpy as np
import math
import statistics as stats
import scipy.stats as scistats
import matplotlib.pyplot as plt
import matplotlib.colors as color
import os
import sys
import statsmodels.api as sm
from pandas import read_table
from netCDF4 import Dataset as ds
import shapely
import warnings
warnings.filterwarnings("ignore")
from matplotlib import cm as cm1
from numpy.linalg import multi_dot
import netCDF4 as nc
from matplotlib.colors import LinearSegmentedColormap
from mpl_toolkits.basemap import Basemap, cm, shiftgrid, addcyclic
from matplotlib.patches import Polygon
from scipy.ndimage import gaussian_filter
# Style Dictionary to standarize plotting
# scheme between different python scripts
styledict = {'xtick.labelsize':21,
'xtick.major.size':9,
'xtick.major.width':1,
'ytick.labelsize':21,
'ytick.major.size':9,
'ytick.major.width':1,
'legend.framealpha':0.0,
'legend.fontsize':15,
'axes.labelsize':22,
'axes.titlesize':22,
'axes.linewidth':2,
'figure.figsize':(12,8),
'savefig.format':'jpg'}
plt.rcParams.update(**styledict)
# define font style for text
font = {'family': 'serif',
'color': 'blue',
'weight': 'normal',
'size': 12,
}
#Go to the working directory
import os
os.chdir("/Users/sshen/climstats")
import netCDF4 as nc
#import .nc data
nc = nc.Dataset("data/air.mon.anom.nc")
# extrat data for variables
lon = nc.variables['lon'][:]
lat = nc.variables['lat'][:]
time = nc.variables['time'][:]
sat = nc.variables['air'][:]
#sat = surface air temp anomalies: mon, 5deg grid
# 1674 months from Jan 1880-June 2019: 139yrs 6mons
print(np.shape(sat)) #(1674, 36, 72): mon, lat, lon
#print(nc)
Dec = list(np.arange(11, 1674, 12))#Dec index
Decsat = sat[Dec, :, :] #extract Dec data
# space-time data matrix STsat
N = 72 * 36
P = len(Dec)
STsat = np.zeros((N,P))
for k in range(P):
STsat[:,k] = Decsat[k,:,:].flatten()
print(STsat.shape)#(2592, 139) 139 Decembers
#print(STsat[0:4, 0:4])
#Space-time dataframe with lat, lon, time
# repeat each value in lat 72 times
LAT = np.repeat(lat,72).reshape(2592,1)
# repeat lon 36 times
LON = np.tile(lon,36).reshape(2592,1)
STanom0 = np.hstack((LAT, LON, STsat))
yrlist = list(np.arange(1880, 2019))
cnames = ["LAT", "LON"] + yrlist
STanom = pd.DataFrame(STanom0, columns=cnames)
#Select the zonal band
n1 = np.where((STanom["LAT"] > -4) & (STanom["LAT"] < 0)\
& (STanom["LON"] > 0) & (STanom["LON"] < 360))
P1 = 81 #Dec 1961 in STsat space-time data matrix
P2 = 130 #Dec 2010
P = P2 - P1 + 1 #50 years
print(STsat.shape) #The space-time data
da0 = STsat[n1,P1:(P2 + 1)]
da1 = da0[0,:,:]
(1674, 36, 72) (2592, 139) (2592, 139)
STanom #Display the space-time data frame
LAT | LON | 1880 | 1881 | 1882 | 1883 | 1884 | 1885 | 1886 | 1887 | ... | 2009 | 2010 | 2011 | 2012 | 2013 | 2014 | 2015 | 2016 | 2017 | 2018 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | -87.5 | 2.5 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | ... | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 |
1 | -87.5 | 7.5 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | ... | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 |
2 | -87.5 | 12.5 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | ... | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 |
3 | -87.5 | 17.5 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | ... | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 |
4 | -87.5 | 22.5 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | ... | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 |
... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
2587 | 87.5 | 337.5 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | ... | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 |
2588 | 87.5 | 342.5 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | ... | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 |
2589 | 87.5 | 347.5 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | ... | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 |
2590 | 87.5 | 352.5 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | ... | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 |
2591 | 87.5 | 357.5 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | ... | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 | -9.969210e+36 |
2592 rows × 141 columns
# define AreaFrac as a float for future manipulation
AreaFrac = np.sqrt(np.cos(LAT[n1]*np.pi/180))
#AreaFrac = float(np.sqrt(np.cos(Lat1[0]*np.pi/180)))
# subtract row means from dat1 to get math anomalies
dat2 = da1 - da1.mean(axis=1, keepdims=True)
# multiply each value by AreaFrac
Banddat = AreaFrac * dat2
#print(Banddat.shape)
covBand = np.dot(Banddat, Banddat.transpose())
# check max and min values
print(np.max(covBand))
print(np.min(covBand))
covBand62 = covBand #To be used for Fig 6.2 later
90.67207862062392 -6.18168468483208
#plot the figure: Fig. 6.1a
int = np.linspace(-7,92,81) # define levels for plot
fig, ax = plt.subplots()
# define color bar
myColMap = LinearSegmentedColormap.from_list(name='my_list',
colors=['navy','darkgreen','lime','yellow',
'pink','red','maroon'], N=100)
# plot covariance plot
Lon1 = np.arange(2.5, 360, 5)
colormap = plt.contourf(Lon1,Lon1,covBand,
levels = int, cmap = myColMap)
plt.title("Covariance of SAT Anomalies on \
Zonal Band 2.5$\degree$S",
size = 21, pad = 15)
plt.xlabel("Longitude", labelpad = 15)
plt.ylabel("Longitude", labelpad = 15)
# plot color bar
cbar = plt.colorbar(colormap, shrink=0.9, pad = 0.05)
plt.text(372, 348, '[$\degree$C]$^2$', size = 20)
cbar.set_ticks([0,20,40,60,80])
plt.savefig("fig0601a.eps")# save figure
#Covert missing values into NaN
STanom[STanom < -999] = pd.NA
STanom.iloc[0:2, 0:5]
LAT | LON | 1880 | 1881 | 1882 | |
---|---|---|---|---|---|
0 | -87.5 | 2.5 | NaN | NaN | NaN |
1 | -87.5 | 7.5 | NaN | NaN | NaN |
plt.figure(figsize=(9, 7))
n1 = np.where((STanom.iloc[:, 0] > -90) & (STanom.iloc[:, 0] < 90)
& (STanom.iloc[:, 1] > 0) & (STanom.iloc[:, 1] < 4))
n1 = np.concatenate(n1)
P1 = 83 # Time index for Dec 1961
P2 = 133 # Dec 2010
P = P2 - P1 # =50 Decembers
dat1 = STanom.iloc[n1, P1:P2] # 1961-2010
Lat1 = STanom.iloc[n1, 0]
Lon1 = STanom.iloc[n1, 1]
AreaFac = np.sqrt(np.cos(Lat1 * np.pi / 180))
dat2 = dat1.transpose() - np.array(np.mean(dat1, axis=1))
Banddat = AreaFac[:, np.newaxis] * dat2.transpose()
covBand = np.dot(Banddat, Banddat.T)
#print(dat2)
# Display max and min of covariance matrix
print("Max:", np.nanmax(covBand))
print("Min:", np.nanmin(covBand))
# Creating the filled contour plot
intervals = np.linspace(-35, 120, 81)
myColMap = LinearSegmentedColormap.from_list(name='my_list',
colors=['navy','darkgreen','lime','yellow',
'pink','red','maroon'], N=100)
ticks = np.arange(-90, 91, 30)
plt.contourf(Lat1, Lat1, covBand, levels=intervals, cmap=myColMap)
clb = plt.colorbar(label='', ticks=np.arange(-40, 120, 20),
pad = 0.05, shrink=0.9)
clb.ax.set_title('[°C]²')
plt.title("SAT Meridional Covariance: 2.5°E")
plt.xlabel("Latitude")
plt.ylabel("Latitude")
plt.xticks(ticks)
plt.yticks(ticks)
plt.grid()
plt.show()
Max: 116.8833093892827 Min: -34.57828535281515
sat.shape
mapmat = sat[1416, :, :]
sat[1410:1419, 0, 0]
masked_array(data=[--, --, --, --, --, --, --, --, --], mask=[ True, True, True, True, True, True, True, True, True], fill_value=-9.96921e+36, dtype=float32)
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import matplotlib.colors as mcolors
import netCDF4 as nc
plt.figure(figsize=(12, 6))
#import .nc data
nc = nc.Dataset("data/air.mon.anom.nc")
# extrat data for variables
lon = nc.variables['lon'][:]
lat = nc.variables['lat'][:]
time = nc.variables['time'][:]
sat = nc.variables['air'][:] #The temp anomalies data
# sat is a 3D array, and lon, lat are longitude and latitude
mapmat = sat[1415, :, :] #Dec 1997
mapmat = np.clip(mapmat, -5, 5)
# Create the color palette
rgb_palette = mcolors.LinearSegmentedColormap.from_list(
"custom", ['black', 'blue', 'darkgreen', 'green', 'white',
'yellow', 'pink', 'red', 'maroon'])
# Initialize the Basemap
m = Basemap(projection='cyl', llcrnrlat=-90,
urcrnrlat=90, llcrnrlon=0, urcrnrlon=360)
# Create meshgrid for plotting
lon, lat = np.meshgrid(lon, lat)
x, y = m(lon, lat)
# Plot the data
cf = m.contourf(x, y, mapmat, levels=np.linspace(-5, 5, 31),
cmap=myColMap, extend='both')
# Add map features
m.drawcoastlines()
m.drawparallels(np.arange(-90., 91., 30.),
labels=[1,0,0,0], fontsize = 15)
m.drawmeridians(np.arange(0., 361., 60.),
labels=[0,0,0,1], fontsize = 15)
# Add the two thick dashed lines
plt.plot([0,360], [-2.5, 2.5],
linestyle='dashed', color = 'k', linewidth = 5)
plt.plot([2.5, 2.5], [-90, 90],
linestyle='dashed', color = 'k', linewidth = 5)
# Title and labels
plt.title("December 1997 SAT Anomalies", fontsize=20)
clb = plt.colorbar(cf, label='', ticks=np.arange(-5, 5.1, 1),
pad = 0.03, shrink=0.8)
clb.ax.set_title('°C')
# Save the figure
plt.savefig("fig0601c.eps", format='eps')
plt.show()