Project 1 – Oxygen Budget

!pwd

!ls

!ls erda_mount/

!ls /home/jovyan/erda_mount/theocean.nc # Verify the file exists first

!pip install xarray

import xarray as xr

# Open and inspect the NetCDF file
ds = xr.open_dataset('/home/jovyan/erda_mount/theocean.nc')

# Show dataset details
print(ds)

# Display variables specifically
print("\nData Variables:")
print(ds.data_vars)

import xarray as xr

# Load the dataset
ds = xr.open_dataset('/home/jovyan/erda_mount/theocean.nc')

# Check if 'DIC' exists in the dataset
if 'DIC' in ds.variables:
    print("'DIC' variable found!\n")
    
    # Print detailed info about the 'DIC' variable
    print(ds['DIC'])
else:
    print("'DIC' variable not found in the dataset.") 

import xarray as xr

# Open the NetCDF file
ds = xr.open_dataset('/home/jovyan/erda_mount/theocean.nc')

# Print all variables and their long_name 
print("=== Variables and Their Descriptions ===\n")
for var in ds.variables:
    long_name = ds[var].attrs.get('long_name', 'No long_name')
    print(f"{var:20}: {long_name}")

import xarray as xr

# Open dataset
ds = xr.open_dataset('/home/jovyan/erda_mount/theocean.nc')

# Loop through all variables
for var in ds.variables:
    da = ds[var]
    
    print(f"{da.dtype} {var}({', '.join(f'{dim}={ds.dims[dim]}' for dim in da.dims)});")
    
    for attr, val in da.attrs.items():
        print(f"  :{attr} = {val!r};")
    
    print()  # Blank line between variables

import xarray as xr
import matplotlib.pyplot as plt

# Load the dataset
ds = xr.open_dataset('/home/jovyan/erda_mount/theocean.nc')

# Select surface DIC (assuming variable DIC is 4D: time, z, y, x)
surface_DIC = ds['DIC'].isel(time=-1, z_t=0)

# Plotting the map of sea surface DIC
fig, ax = plt.subplots(figsize=(6, 5), dpi = 300)
surface_DIC.plot.pcolormesh(
    ax=ax,
    cmap='coolwarm',  # choose colormap here
    shading='auto',
    cbar_kwargs={'label': 'DIC (mol m$^{-3}$)'}
)

ax.set_title('Sea Surface DIC')
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')

plt.tight_layout()

# Save the figure at the save folder
plt.savefig('./sea_surface_DIC.png', transparent=True)

plt.show()

import numpy as np

ds = xr.open_dataset('/home/jovyan/erda_mount/theocean.nc')

# Print coordinates specifically
print("\n--- Coordinate Variables ---")
for coord in ds.coords:
    print(f"{coord}: {ds.coords[coord].dims} {ds.coords[coord].shape}")

# Print variable names and their dimensions
print("\n--- Data Variables ---")
for var in ds.data_vars:
    print(f"{var}: {ds[var].dims} {ds[var].shape}")

# Print attributes of TLAT and TLONG if they exist
if 'TLAT' in ds:
    print("\nTLAT attributes:", ds['TLAT'].attrs)
if 'TLONG' in ds:
    print("\nTLONG attributes:", ds['TLONG'].attrs)

# Print random sample lat/lon points to verify
print("\n--- Sample Latitude and Longitude Points ---")
TLAT = ds['TLAT'].values
TLONG = ds['TLONG'].values

nlat, nlon = TLAT.shape
for _ in range(5):
    i = np.random.randint(0, nlat)
    j = np.random.randint(0, nlon)
    print(f"Index (i={i}, j={j}) -> Lat: {TLAT[i, j]:.2f}°, Lon: {TLONG[i, j]:.2f}°")

# Surface map (latitude-longitude at surface level z_t=0)
dic_surface = ds['DIC'].isel(time=time_index, z_t=0)
lat_2d = ds['TLAT'].values
lon_2d = ds['TLONG'].values
lon_2d = np.where(lon_2d > 180, lon_2d - 360, lon_2d)  # Convert to [-180, 180]

import matplotlib.pyplot as plt
import numpy as np
from matplotlib import colors

# Create sample data
N = 100
X, Y = np.mgrid[-3:3:complex(0, N), -2:2:complex(0, N)]
Z1 = np.exp(-X**2 - Y**2)
Z2 = np.exp(-(X - 1)**2 - (Y - 1)**2)
Z = (Z1 - Z2) * 2

# Normalize options
vmin, vmax = -1.0, 1.0

norms = {
    'Linear': None,
    'LogNorm': colors.LogNorm(vmin=Z[Z > 0].min(), vmax=Z.max()),
    'PowerNorm (quartic γ=0.25)': colors.PowerNorm(gamma=0.25, vmin=vmin, vmax=vmax),
    'PowerNorm (γ=0.5)': colors.PowerNorm(gamma=0.5, vmin=vmin, vmax=vmax),
    'TwoSlopeNorm (center=0)': colors.TwoSlopeNorm(vmin=vmin, vcenter=0, vmax=vmax),
    'Quantile Bins': colors.BoundaryNorm(boundaries=np.quantile(Z, np.linspace(0, 1, 10)), ncolors=256)
}

# Plotting
plt.rcParams['text.color'] = '#d6ebff'
plt.rcParams['axes.labelcolor'] = '#d6ebff'
plt.rcParams['xtick.color'] = '#d6ebff'
plt.rcParams['ytick.color'] = '#d6ebff'
plt.rcParams['axes.edgecolor'] = '#d6ebff'

fig, axs = plt.subplots(nrows=3, ncols=2, figsize=(14, 12), dpi = 300)
axs = axs.ravel()
fig.patch.set_alpha(0)
for ax in axs:
    ax.set_facecolor('none')

for ax, (name, norm) in zip(axs, norms.items()):
    # LogNorm only works on positive values
    if name == 'LogNorm':
        pcm = ax.pcolormesh(X, Y, np.abs(Z), norm=norm, cmap='coolwarm', shading='auto')
    else:
        pcm = ax.pcolormesh(X, Y, Z, norm=norm, cmap='coolwarm', shading='auto')

    fig.colorbar(pcm, ax=ax)
    ax.set_title(name)
    ax.set_aspect('equal')

plt.suptitle("Colorbar Normalization Examples", fontsize=16)
plt.tight_layout(rect=[0, 0, 1, 0.97])
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle

# --- Generate clean sine signal ---
N = 20
dt = 1.0
A = 1.0
f = 0.1
omega = 2 * np.pi * f
timesteps = np.linspace(1, N, N)
signal = A * np.sin(omega * timesteps)

# --- Add noise ---
np.random.seed(1)
noise = 0.3 * (np.random.rand(N) - 0.5)
noisy_signal = signal + noise

# --- Boxcar smoothing ---
def boxcar_smooth(x, w):
    smoothed = np.zeros_like(x)
    half_w = w // 2
    for i in range(len(x)):
        left = max(0, i - half_w)
        right = min(len(x), i + half_w + 1)
        smoothed[i] = np.mean(x[left:right])
    return smoothed

w = 5
smoothed_signal = boxcar_smooth(noisy_signal, w=w)

# --- Plot with styled settings ---
axis_color = '#d6ebff'
noise_color = 'red'
smooth_color = '#66c2ff'
box_color = '#00ffff'

fig, ax = plt.subplots(figsize=(10, 5), dpi=300)
fig.patch.set_alpha(0.0)
ax.set_facecolor('none')

# Scatter noisy and smoothed data
ax.scatter(timesteps, noisy_signal, color=noise_color, s=30, label='Noisy Data', marker='X')
ax.scatter(timesteps, smoothed_signal, color=smooth_color, s=50, label='Smoothed')

# Draw example boxcar windows
for center in [3, 10, 17]:
    left = max(0, center - w // 2)
    right = min(N, center + w // 2 + 1)
    x_min = timesteps[left]
    x_max = timesteps[right - 1]
    y_min = min(noisy_signal[left:right]) - 0.1
    y_max = max(noisy_signal[left:right]) + 0.1
    rect = Rectangle((x_min, y_min), x_max - x_min, y_max - y_min,
                     edgecolor=box_color, facecolor='none', lw=2)
    ax.add_patch(rect)

ax.set_title("Boxcar Smoothing", color=axis_color)
ax.set_xlabel("Time", color=axis_color)
ax.set_ylabel("Signal", color=axis_color)
ax.set_ylim(-1.5, 1.5)
ax.grid(True, linestyle='--', color='#3c5c78', alpha=0.6)

# Legend
legend = ax.legend()
legend.get_frame().set_alpha(0.3)
for text in legend.get_texts():
    text.set_color(axis_color)


# Axis ticks and spines
ax.tick_params(colors=axis_color)
for spine in ax.spines.values():
    spine.set_color(axis_color)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# --- Generate clean sine signal ---
N = 1000           # Number of samples
dt = 1.0 / 100.0    # 100 samples/sec
A = 5.0             # Amplitude
f = 0.5             # Frequency in Hz
omega = 2 * np.pi * f
timesteps = np.linspace(0.0, N*dt, N)
signal = A * np.sin(omega * timesteps)

# --- Add noise ---
np.random.seed(0)  # For reproducibility
noise = 3 * (np.random.rand(N) - 0.5)  # Random noise centered around 0
noisy_signal = signal + noise

# --- Boxcar smoothing ---
def boxcar_smooth(x, w):
    """Smooth 1D array x using a boxcar (moving average) of width w."""
    smoothed = np.zeros_like(x)
    half_w = w // 2
    for i in range(len(x)):
        left = max(0, i - half_w)
        right = min(len(x), i + half_w + 1)
        smoothed[i] = np.mean(x[left:right])
    return smoothed

smoothed_signal = boxcar_smooth(noisy_signal, w=13)

# --- Plot noisy vs smoothed with transparent background and styled text ---
axis_color = '#d6ebff'  # light blue for axis
noise_color = '#aeeafc' # lighter color for noisy points

fig, axs = plt.subplots(1, 2, figsize=(14, 5), dpi=300)
fig.patch.set_alpha(0.0)

# Plot 1: Noisy Signal
axs[0].set_facecolor('none')
axs[0].plot(timesteps, noisy_signal, '.', color=noise_color, alpha=0.7)  # Change color here
axs[0].set_xlim(0, 4)
axs[0].set_ylim(-8, 8)
axs[0].set_title("Noisy Signal", color=axis_color)
axs[0].tick_params(colors=axis_color)
axs[0].set_xlabel("Time", color=axis_color)
axs[0].set_ylabel("Amplitude", color=axis_color)
for spine in axs[0].spines.values():
    spine.set_color(axis_color)

# Plot 2: Smoothed Signal
axs[1].set_facecolor('none')
axs[1].plot(timesteps, smoothed_signal, 'r')
axs[1].set_xlim(0, 4)
axs[1].set_ylim(-8, 8)
axs[1].set_title("Smoothed Signal (Boxcar w=13)", color=axis_color)
axs[1].tick_params(colors=axis_color)
axs[1].set_xlabel("Time", color=axis_color)
axs[1].set_ylabel("Amplitude", color=axis_color)
for spine in axs[1].spines.values():
    spine.set_color(axis_color)

plt.tight_layout()
plt.show()

import matplotlib.pyplot as plt
import numpy as np
from matplotlib import colors

# --- Generate Sample 2D Data ---
N = 100
X, Y = np.mgrid[-3:3:complex(0, N), -2:2:complex(0, N)]
Z1 = np.exp(-X**2 - Y**2)
Z2 = np.exp(-(X - 1)**2 - (Y - 1)**2)
Z = (Z1 - Z2) * 2

# --- Boxcar Smoothing 2D ---
def boxcar_smooth_2d(data, w):
    smoothed = np.zeros_like(data)
    half_w = w // 2
    for i in range(data.shape[0]):
        for j in range(data.shape[1]):
            left_i = max(0, i - half_w)
            right_i = min(data.shape[0], i + half_w + 1)
            left_j = max(0, j - half_w)
            right_j = min(data.shape[1], j + half_w + 1)
            smoothed[i, j] = np.mean(data[left_i:right_i, left_j:right_j])
    return smoothed

Z_smoothed = boxcar_smooth_2d(Z, w=50)

# --- Quantile Bin Normalization ---
boundaries = np.quantile(Z, np.linspace(0, 1, 10))
quantile_norm = colors.BoundaryNorm(boundaries=boundaries, ncolors=256)

# --- Plotting Style Settings ---
plt.rcParams['text.color'] = '#d6ebff'
plt.rcParams['axes.labelcolor'] = '#d6ebff'
plt.rcParams['xtick.color'] = '#d6ebff'
plt.rcParams['ytick.color'] = '#d6ebff'
plt.rcParams['axes.edgecolor'] = '#d6ebff'

fig, axs = plt.subplots(1, 2, figsize=(14, 6), dpi=300)
fig.patch.set_alpha(0)
for ax in axs:
    ax.set_facecolor('none')

# --- Plot Original ---
pcm1 = axs[0].pcolormesh(X, Y, Z, norm=quantile_norm, cmap='coolwarm', shading='auto')
fig.colorbar(pcm1, ax=axs[0])
axs[0].set_title('Original Data', fontsize=14)
axs[0].set_aspect('equal')

# --- Plot Smoothed ---
pcm2 = axs[1].pcolormesh(X, Y, Z_smoothed, norm=quantile_norm, cmap='coolwarm', shading='auto')
fig.colorbar(pcm2, ax=axs[1])
axs[1].set_title('Smoothed Data (Boxcar w=50)', fontsize=14)
axs[1].set_aspect('equal')

plt.suptitle('Boxcar Smoothing on 2D Heatmap (Quantile Normalization)', fontsize=16)
plt.tight_layout(rect=[0, 0, 1, 0.96])
plt.show()