<jupyter_start><jupyter_text>Read SST data, Create and Save nino3 time seriesIn this notebook, we will finish these tasks* subsample SST data over the nino3 area* calcualte climatology* calculate anomalies* calculate regional mean* plot regional mean SST time series* save data  1. Load basic libraries<jupyter_code>%matplotlib inline
import numpy as np
from numpy import nonzero

import matplotlib.pyplot as plt         # to generate plots
from mpl_toolkits.basemap import Basemap # plot on map projections
import matplotlib.dates as mdates

import datetime             
from netCDF4 import Dataset # http://unidata.github.io/netcdf4-python/
from netCDF4 import netcdftime
from netcdftime import utime<jupyter_output><empty_output><jupyter_text>2. Set and read input NetCDF file info  2.1 Read data<jupyter_code>ncfile = 'data\skt.mon.mean.nc'

fh     = Dataset(ncfile, mode='r') # file handle, open in read only mode
lon    = fh.variables['lon'][:]
lat    = fh.variables['lat'][:]
nctime = fh.variables['time'][:]
t_unit = fh.variables['time'].units
skt    = fh.variables['skt'][:]

try :
    t_cal = fh.variables['time'].calendar
except AttributeError : # Attribute doesn't exist
    t_cal = u"gregorian" # or standard

fh.close() # close the file<jupyter_output><empty_output><jupyter_text>2.2 Parse time<jupyter_code>utime = netcdftime.utime(t_unit, calendar = t_cal)
datevar = utime.num2date(nctime)
print(datevar.shape)

datevar[0:5]<jupyter_output>(687L,)<jupyter_text>3. Subregion for nino3 area* Lat: -5 ~ 5* Lon: 210 ~ 270  3.1 Get indices of time, lat and lon over the nino3 area<jupyter_code>idx_lat_n3  = (lat>=-5.0) * (lat<=5.0)
idx_lon_n3  = (lon>=210.0) * (lon<=270.0)<jupyter_output><empty_output><jupyter_text>* time: 1970-1999<jupyter_code>years = np.array([idx.year for idx in datevar])
idx_tim_n3 = (years>=1970) * (years<=1999)<jupyter_output><empty_output><jupyter_text>Get Index using np.nonzero<jupyter_code>idxtim = nonzero(idx_tim_n3)[0]
#idxlat = nonzero(idx_lat_n3)[0]
idxlon = nonzero(idx_lon_n3)[0]
idxlon<jupyter_output><empty_output><jupyter_text>3.2 Extract data over nino3 areaUse logical indexing is fine for 1D array; however, a little funny for multiple dimension array.<jupyter_code>lat_n3 = lat[idx_lat_n3]
lon_n3 = lon[idx_lon_n3]
dates_n3  = datevar[idx_tim_n3]

skt_n3 = skt[idx_tim_n3, :, :][:,idx_lat_n3,:][:,:,idx_lon_n3]
print(skt_n3.shape)
print(dates_n3.shape)<jupyter_output>(360L, 6L, 33L)
(360L,)<jupyter_text>4. Calculate region means  4.1 Transform skt_n3 from months|lat|lon => 12|year|lat|lon => year|12|lat|lon<jupyter_code>skt_n3 = np.reshape(skt_n3, (12,30,6,33), order='F')
skt_n3 = np.transpose(skt_n3, (1, 0, 2, 3))
skt_n3.shape<jupyter_output><empty_output><jupyter_text>4.2 Calculate monthly climatology<jupyter_code>clima_skt_n3 = np.mean(skt_n3, axis=0)
clima_skt_n3.shape<jupyter_output><empty_output><jupyter_text>4.3 Calculate anomaly of SST over nino3 area<jupyter_code>num_repeats = 30 # 30 years
clima_skt_n3 = np.vstack([clima_skt_n3]*num_repeats)
clima_skt_n3.shape
clima_skt_n3 = np.reshape(clima_skt_n3, (12,30,6,33),order='F')
clima_skt_n3 = np.transpose(clima_skt_n3, (1, 0, 2, 3))
clima_skt_n3.shape
ssta  = skt_n3-clima_skt_n3
ssta2 = np.reshape(ssta,(30,12,6*33), order='F') # 30x12x198
ssta3 = np.mean(ssta2, axis=2); # 30x12
ssta3.shape
ssta_series = np.reshape(ssta3.T,(12*30,1), order='F'); # 1x360
ssta_series.shape<jupyter_output><empty_output><jupyter_text>5. Have a beautiful look<jupyter_code>import matplotlib.dates as mdates
from matplotlib.dates import MonthLocator, WeekdayLocator, DateFormatter
import matplotlib.ticker as ticker

fig, ax = plt.subplots(1, 1 , figsize=(15,5))
ax.plot(dates_n3, ssta_series)
ax.set_ylim((-4,4))

#horiz_line_data = np.array([0 for i in np.arange(len(dates_n3))])
#ax.plot(dates_n3, horiz_line_data, 'r--') 
ax.axhline(0, color='r')

ax.set_title('NINO3 SSTA 1970-1999')
ax.set_ylabel(['$^oC$'])
ax.set_xlabel('Date')

# rotate and align the tick labels so they look better
fig.autofmt_xdate()

# use a more precise date string for the x axis locations in the toolbar
ax.fmt_xdata = mdates.DateFormatter('%Y')<jupyter_output><empty_output><jupyter_text>6. Save data<jupyter_code>np.savez('data/ssta.nino3.30y.npz', ssta_series=ssta_series)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13125436350942518, 'fraction_numerical': 0.04863858505934373, 'mean_word_length': 6.487804878048781, 'pattern_counts': {'":': 0, '<': 50, '<?xml version=': 0, '>': 52, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '417837', 'license': 'MIT', 'n_tokens_mistral': 1885, 'n_tokens_neox': 1769, 'n_words': 356, 'path': 'ex02-Read SST data, create and save nino3 time series.ipynb', 'repo_name': 'monocilindro/Python-Practical-Application-on-Climate-Variability-Studies', 'size': '4297'}