# Import the tools we are going to need today:
import matplotlib.pyplot as plt  # plotting library
import numpy as np  # numerical library
import xarray as xr  # netCDF library
import pandas as pd  # tabular library
import cartopy  # Map projections libary
import cartopy.crs as ccrs  # Projections list
import cartopy.feature as cfeature
# Some defaults:
plt.rcParams['figure.figsize'] = (12, 5)  # Default plot size

ds = xr.open_dataset('ERA5_LowRes_Monthly_t2m.nc') - 273.15 # convert to °C

ds_sst = xr.open_dataset('ERA5_LowRes_Monthly_sst.nc')

ds_prec = xr.open_dataset('ERA5_LowRes_Monthly_tp.nc') * 1000 #convert to mm

ds_wind = xr.open_dataset('ERA5_LowRes_Monthly_uvslp.nc')

ds = xr.merge([ds, ds_sst, ds_wind, ds_prec])

# Open csv as dataframe with matching length
df = pd.read_csv('co2_mm_gl.csv', skiprows=38, skipfooter=57, parse_dates={'date' : [0, 1]}, index_col='date')

# Add needed column to dataset
ds['co2'] = ('time', df['average'])
ds

C:\Users\michl\AppData\Local\Temp\ipykernel_12320\2067514920.py:12: ParserWarning: Falling back to the 'python' engine because the 'c' engine does not support skipfooter; you can avoid this warning by specifying engine='python'.
  df = pd.read_csv('co2_mm_gl.csv', skiprows=38, skipfooter=57, parse_dates={'date' : [0, 1]}, index_col='date')
C:\Users\michl\AppData\Local\Temp\ipykernel_12320\2067514920.py:12: UserWarning: Could not infer format, so each element will be parsed individually, falling back to `dateutil`. To ensure parsing is consistent and as-expected, please specify a format.
  df = pd.read_csv('co2_mm_gl.csv', skiprows=38, skipfooter=57, parse_dates={'date' : [0, 1]}, index_col='date')

<xarray.Dataset>
Dimensions:    (longitude: 480, latitude: 241, time: 480)
Coordinates:
  * longitude  (longitude) float32 -179.6 -178.9 -178.1 ... 178.1 178.9 179.6
  * latitude   (latitude) float32 90.0 89.25 88.5 87.75 ... -88.5 -89.25 -90.0
  * time       (time) datetime64[ns] 1979-01-01 1979-02-01 ... 2018-12-01
Data variables:
    t2m        (time, latitude, longitude) float32 -28.44 -28.44 ... -23.82
    sst        (time, latitude, longitude) float32 ...
    u10        (time, latitude, longitude) float32 ...
    v10        (time, latitude, longitude) float32 ...
    msl        (time, latitude, longitude) float32 ...
    tp         (time, latitude, longitude) float32 0.2937 0.2937 ... 0.3377
    co2        (time) float64 336.6 337.3 337.9 338.3 ... 406.6 408.1 409.2
Attributes:
    Conventions:  CF-1.6
    history:      2019-11-18 09:36:58 GMT by grib_to_netcdf-2.14.0: /opt/ecmw...

array([-179.625, -178.875, -178.125, ...,  178.125,  178.875,  179.625],
      dtype=float32)

array([ 90.  ,  89.25,  88.5 , ..., -88.5 , -89.25, -90.  ], dtype=float32)

array(['1979-01-01T00:00:00.000000000', '1979-02-01T00:00:00.000000000',
       '1979-03-01T00:00:00.000000000', ..., '2018-10-01T00:00:00.000000000',
       '2018-11-01T00:00:00.000000000', '2018-12-01T00:00:00.000000000'],
      dtype='datetime64[ns]')

array([[[-28.443253 , -28.443253 , -28.443253 , ..., -28.443253 ,
         -28.443253 , -28.443253 ],
        [-28.795761 , -28.80458  , -28.813385 , ..., -28.778137 ,
         -28.785187 , -28.790466 ],
        [-28.4926   , -28.508469 , -28.526093 , ..., -28.446777 ,
         -28.460876 , -28.476746 ],
        ...,
        [-31.897858 , -31.885529 , -31.876709 , ..., -31.933105 ,
         -31.926056 , -31.911957 ],
        [-31.37262  , -31.35675  , -31.339127 , ..., -31.43431  ,
         -31.411392 , -31.390244 ],
        [-30.223434 , -30.223434 , -30.223434 , ..., -30.223434 ,
         -30.223434 , -30.223434 ]],

       [[-31.70398  , -31.70398  , -31.70398  , ..., -31.70398  ,
         -31.70398  , -31.70398  ],
        [-32.24156  , -32.24861  , -32.259186 , ..., -32.229218 ,
         -32.232742 , -32.236267 ],
        [-32.521805 , -32.53415  , -32.551773 , ..., -32.486557 ,
         -32.500656 , -32.51123  ],
...
        [-37.497498 , -37.479874 , -37.4693   , ..., -37.541565 ,
         -37.53099  , -37.513367 ],
        [-37.6015   , -37.583862 , -37.564484 , ..., -37.661423 ,
         -37.64203  , -37.62088  ],
        [-36.91057  , -36.91057  , -36.91057  , ..., -36.91057  ,
         -36.91057  , -36.91057  ]],

       [[-21.588669 , -21.588669 , -21.588669 , ..., -21.588669 ,
         -21.588669 , -21.588669 ],
        [-21.766693 , -21.771973 , -21.779022 , ..., -21.763153 ,
         -21.763153 , -21.764923 ],
        [-22.13153  , -22.142105 , -22.15445  , ..., -22.108627 ,
         -22.113907 , -22.124481 ],
        ...,
        [-26.181885 , -26.160736 , -26.151932 , ..., -26.225952 ,
         -26.218903 , -26.201279 ],
        [-25.0327   , -25.015076 , -24.997452 , ..., -25.083817 ,
         -25.066193 , -25.050323 ],
        [-23.823593 , -23.823593 , -23.823593 , ..., -23.823593 ,
         -23.823593 , -23.823593 ]]], dtype=float32)

[55526400 values with dtype=float32]

[55526400 values with dtype=float32]

[55526400 values with dtype=float32]

[55526400 values with dtype=float32]

# Compute the temporal mean temperature
T_mean = ds['t2m'].mean(dim='time')

# Compute the zonal anomaly map of average temperature
T_star = ds['t2m'].mean(dim='time') - ds['t2m'].mean(dim=['time', 'longitude'])

# Compute the monthly average temperature for each month
T_M = ds['t2m'].groupby('time.month').mean(dim='time')

# Compute the average monthly temperature range
T_M_range = T_M.max(dim='month') - T_M.min(dim='month')

# Plot the temporal mean temperature
T_mean.plot()
plt.title('Temporal Mean Temperature')
plt.show()

# Plot the zonal anomaly map of average temperature
T_star.plot()
plt.title('Zonal Anomaly Map of Average Temperature')
plt.show()

# Plot montly average temperature (annual cycle)
T_M.mean(dim=['longitude', 'latitude']).plot(figsize=[12,5])
plt.title('Annual cycle of mean temperature')
plt.show()

levels = np.linspace(0,60,11)
cmap = 'YlGnBu'

# Plot the average monthly temperature range
T_M_range.plot(cmap='Reds')
plt.title('Average Monthly Temperature Range')
plt.show()

month_dict = {
    1: 'January',
    2: 'February',
    3: 'March',
    4: 'April',
    5: 'May',
    6: 'June',
    7: 'July',
    8: 'August',
    9: 'September',
    10: 'October',
    11: 'November',
    12: 'December'
}


# Compute the average daily precipitation for each month
monthly_precip = ds['tp'].groupby('time.month').mean('time')

# Define the levels and colormap for the plot
levels = [0.5, 1, 2, 3, 4, 5, 7, 10, 15, 20, 40]
cmap = 'YlGnBu'

# Create a single figure with two subplots
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12,10), subplot_kw={'projection':ccrs.PlateCarree()})

# Plot average daily precipitation for January
month = 1
monthly_precip.sel(month=month).plot(ax=ax1, levels=levels, cmap=cmap, transform=ccrs.PlateCarree())
ax1.add_feature(cfeature.COASTLINE)
gridlines = ax1.gridlines(draw_labels=True)
gridlines.bottom_labels  = False
gridlines.right_labels = False
ax1.set_title(f'Average Daily Precipitation in {month_dict[month]}', rotation=90, x=-0.1, y=0)

# Plot average daily precipitation for August
month = 8
monthly_precip.sel(month=month).plot(ax=ax2, levels=levels, cmap=cmap, transform=ccrs.PlateCarree())
ax2.add_feature(cfeature.COASTLINE)
gridlines = ax2.gridlines(draw_labels=True)
gridlines.top_labels  = False
gridlines.right_labels = False
ax2.set_title(f'Average Daily Precipitation in {month_dict[month]}', rotation=90, x=-0.1, y=0)

plt.savefig('rain.png')
plt.show()

ds = ds.rename({'u10':'Zonal Winds', 'v10':'Meridional Winds'})
# Compute the temporal and zonal average of sea-level pressure
avg_slp = ds['msl'].mean(dim=['time', 'longitude'])

# Convert it to hPa
avg_slp_hpa = avg_slp / 100

fig, ax1 = plt.subplots(figsize=(10, 6))

color = 'tab:red'
ax1.axhline(y=1013.25, color='r', linestyle='--', label='standard atmosphere pressure')  # standard atmosphere pressure
ax1.set_xlabel('Latitude')
ax1.set_ylabel('Pressure (hPa)', color=color)
ax1.plot(avg_slp_hpa['latitude'], avg_slp_hpa, color=color, label='sea-lvl-pressure')  # specify latitude as x-axis
ax1.tick_params(axis='y', labelcolor=color)


ax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis

colors = ['blue', 'lightblue']
variables = ['Zonal Winds', 'Meridional Winds']
ax2.set_ylabel('Speed (m/s)')  # we already handled the x-label with ax1
ax2.axhline(y=0, color='b', linestyle='--', linewidth=1, label='0-line')
#ax1.axvline(x=)
for var, color in zip(variables, colors):
    avg_var = ds[var].mean(dim=['time', 'longitude'])
    ax2.plot(avg_var['latitude'], avg_var, color=color, label=var)  # specify latitude as x-axis
    ax2.tick_params(axis='y', labelcolor=color)

# Adjust the y-axis limits so that s_msl corresponds to 0
ax1.set_ylim(bottom=983, top=1043.5)
ax2.set_ylim(bottom=-8, top=8)

fig.tight_layout()  # otherwise the right y-label is slightly clipped
plt.title('Temporal and Zonal Average of Sea-Level Pressure and Wind Speeds')

# Add legends
ax1.legend(loc='upper left')
ax2.legend(loc='upper right')
plt.savefig('msl.png')
plt.show()

def slope(ds):
    y = ds.values
    x = np.arange(len(y))
    return np.polyfit(x, y, 1)[0]

annual_avg = ds['co2'].resample(time='Y').mean()
pre_idust = annual_avg.sel(time=slice('1980', '1985'))
idust = annual_avg.sel(time=slice('2013', '2018'))
co2_pi = pre_idust.mean()
co2_i = idust.mean()
print(f'pre industrial incrase (1980-1985) = {slope(pre_idust):.2f} ppm/year')
print(f'industrial incrase in (2013-2018) = {slope(idust):.2f} ppm/year')
print(f'pre industrial co2 concentratin (1980-1985): {co2_pi.values:.2f} ppm')
print(f'industrial co2 concentratin (2013-2018): {co2_i.values:.2f} ppm')
# Compute the global temp for every timestep resample by years and get means
avg_temp = ds['t2m'].mean(dim=['longitude', 'latitude']).resample(time='Y').mean()

# Initialize figure
fig, ax1 = plt.subplots(figsize=(10, 6))

# Plot monthly and yearly average of CO2 concentration
color = 'tab:blue'
ax1.set_xlabel('Year')
ax1.set_ylabel('CO2 Concentration (ppm)', color=color)
ln1 = ax1.plot(ds['time'], ds['co2'], linewidth=1, color='black', label='Monthly CO2 Concentration')
ln2 = ax1.plot(annual_avg['time'], annual_avg, linewidth=3, color=color, label='Annual Mean CO2')
ax1.tick_params(axis='y', labelcolor=color)

# Create second y-axis
ax2 = ax1.twinx()

# Plot yearly average temperature and convert it from Kelvin to Celsius
color = 'tab:red'
ax2.set_ylabel('Temperature (°C)', color=color)
ln3 = ax2.plot(avg_temp['time'], avg_temp, linewidth=3, color=color, label='Annual Mean Temperature')
ax2.tick_params(axis='y', labelcolor=color)

# Plot title
plt.title('Annual Average Global CO2 Concentration and 2m Temperature')

# Add a grid
plt.grid(True)

# Add a common legend
lns = ln1 + ln2 + ln3
labs = [l.get_label() for l in lns]
ax1.legend(lns, labs, loc='upper left')

plt.show()

pre industrial incrase (1980-1985) = 1.33 ppm/year
industrial incrase in (2013-2018) = 2.52 ppm/year
pre industrial co2 concentratin (1980-1985): 342.00 ppm
industrial co2 concentratin (2013-2018): 401.38 ppm

def compute_anomalies(ds, variable, start_year, end_year):
    """Compute the zonal anomaly map of average for a given variable and time period.

    Parameters:
    ds (xarray.Dataset): The dataset containing the variable data.
    variable (str): The variable to compute anomalies for.
    start_year (str): The start year of the period.
    end_year (str): The end year of the period.

    Returns:
    xarray.DataArray: The anomalies for the specified variable and period.
    """
    # Compute the long-term mean
    long_term_mean = ds[variable].mean(dim='time')

    # Compute the mean for the specified period
    period_mean = ds.sel(time=slice(start_year, end_year))[variable].mean(dim='time')

    # Compute the anomalies
    anomalies = period_mean - long_term_mean

    return anomalies

variables = ['t2m', 'sst', 'tp']
periods = {'El Nino' : ('1997-06-01', '1998-05-01'),
           'El Nino2' : ('1982-07-01', '1983-07-01'),
           'La Nina' : ('1999-06-01', '2000-05-01'),
           'La Nina2' : ('2010-07-01', '2011-04-01')
           }
col = len(periods)
row = len(variables)

fig, axs = plt.subplots(row, col, figsize=(col*12, row*5))

# Add a title for the whole figure
fig.suptitle('Strong El Nino / La Nina seasons anomalies against longterm mean', fontsize=32)

for i, variable in enumerate(variables):
    for j, (period_name, (start_year, end_year)) in enumerate(periods.items()):
        # Compute anomalies
        anomalies = compute_anomalies(ds, variable, start_year, end_year)
        
        # trash outliers to get a better readable plot
        mask = anomalies > 10
        count = mask.sum().values
        print(f'{period_name}:{count}')
        
        #more deatail
        anomalies = anomalies.where(anomalies < 10, np.nan)

        # Choose a colormap
        cmap = 'RdBu_r' if variable != 'tp' else 'RdBu'

        # Plot the anomalies
        anomalies.plot(ax=axs[i, j], cmap=cmap)
        
        # Add a title for each row
        if j == 0:
            axs[i, j].set_ylabel(f'{variable.capitalize()}', rotation=90, fontsize=16, )#labelpad=30)
        
        # Add a title for each column
        if i == 0:
            axs[i, j].set_title(f'{period_name} ({start_year},{end_year})', fontsize=16, pad=20)

plt.tight_layout(rect=[0, 0, 1, 0.96])  # Adjust for the suptitle
plt.savefig('anomalies.png')
plt.show()

El Nino:0
El Nino2:0
La Nina:0
La Nina2:0
El Nino:0
El Nino2:0
La Nina:0
La Nina2:0
El Nino:12
El Nino2:6
La Nina:0
La Nina2:0

Group project: the Climate System¶

Instructions¶

Part 1 - temperature climatology¶

Part 2 - Precipitation climatology¶

Part 3: sea-level pressure and surface winds¶

Part 4: temperature change and CO $_2$ concentrations¶

Part 5: variability and ENSO (open research question)¶

Group project: the Climate System¶

Instructions¶

Part 1 - temperature climatology¶

Part 2 - Precipitation climatology¶

Part 3: sea-level pressure and surface winds¶

Part 4: temperature change and CO2_2 concentrations¶

Part 5: variability and ENSO (open research question)¶

Part 4: temperature change and CO $_2$ concentrations¶