Climate Modeling và Big Data

# Climate Model Ensemble Analysis
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter

# Load CMIP6 model ensemble data
# This example uses pre-processed ensemble data
def load_cmip6_ensemble():
    # In practice, you would load from multiple model files
    # Here we simulate ensemble data
    models = ["CESM2", "GFDL-ESM4", "MPI-ESM1-2-HR", "UKESM1-0-LL"]
    
    # Simulate temperature change data (2081-2100 vs 1995-2014)
    lat = np.linspace(-90, 90, 180)
    lon = np.linspace(-180, 180, 360)
    
    ensemble_data = []
    for i, model in enumerate(models):
        # Simulate model-specific warming patterns
        temp_change = np.random.normal(2.5 + i*0.3, 0.5, (len(lat), len(lon)))
        # Add realistic spatial patterns
        lat_grid, lon_grid = np.meshgrid(lat, lon, indexing="ij")
        temp_change += 1.5 * np.cos(np.radians(lat_grid))**2  # Arctic amplification
        
        da = xr.DataArray(
            temp_change,
            coords={"lat": lat, "lon": lon},
            dims=["lat", "lon"],
            name=f"tas_change_{model}"
        )
        ensemble_data.append(da)
    
    return xr.concat(ensemble_data, dim="model")

# Load ensemble data
ensemble = load_cmip6_ensemble()

# Calculate ensemble statistics
ensemble_mean = ensemble.mean(dim="model")
ensemble_std = ensemble.std(dim="model")
model_agreement = (ensemble > 0).sum(dim="model") / len(ensemble.model)

# Create visualization
fig = plt.figure(figsize=(15, 12))

# Ensemble mean temperature change
ax1 = plt.subplot(2, 2, 1, projection=ccrs.PlateCarree())
im1 = ax1.contourf(ensemble_mean.lon, ensemble_mean.lat, ensemble_mean, 
                   levels=np.linspace(0, 6, 13), cmap="RdYlBu_r", 
                   transform=ccrs.PlateCarree())
ax1.add_feature(cfeature.COASTLINE)
ax1.add_feature(cfeature.BORDERS)
ax1.set_global()
ax1.set_title("Ensemble Mean Temperature Change (°C)")
cbar1 = plt.colorbar(im1, ax=ax1, shrink=0.6)

# Ensemble standard deviation (uncertainty)
ax2 = plt.subplot(2, 2, 2, projection=ccrs.PlateCarree())
im2 = ax2.contourf(ensemble_std.lon, ensemble_std.lat, ensemble_std,
                   levels=np.linspace(0, 2, 11), cmap="viridis",
                   transform=ccrs.PlateCarree())
ax2.add_feature(cfeature.COASTLINE)
ax2.add_feature(cfeature.BORDERS)
ax2.set_global()
ax2.set_title("Ensemble Standard Deviation (°C)")
cbar2 = plt.colorbar(im2, ax=ax2, shrink=0.6)

# Model agreement
ax3 = plt.subplot(2, 2, 3, projection=ccrs.PlateCarree())
im3 = ax3.contourf(model_agreement.lon, model_agreement.lat, model_agreement,
                   levels=np.linspace(0, 1, 11), cmap="RdYlGn",
                   transform=ccrs.PlateCarree())
ax3.add_feature(cfeature.COASTLINE)
ax3.add_feature(cfeature.BORDERS)
ax3.set_global()
ax3.set_title("Model Agreement (fraction agreeing on warming)")
cbar3 = plt.colorbar(im3, ax=ax3, shrink=0.6)

# Regional analysis
ax4 = plt.subplot(2, 2, 4)
# Arctic vs Tropical warming
arctic_warming = ensemble.sel(lat=slice(70, 90)).mean(dim=["lat", "lon"])
tropical_warming = ensemble.sel(lat=slice(-20, 20)).mean(dim=["lat", "lon"])

ax4.boxplot([arctic_warming.values, tropical_warming.values], 
           labels=["Arctic (70-90°N)", "Tropics (20°S-20°N)"])
ax4.set_ylabel("Temperature Change (°C)")
ax4.set_title("Regional Warming Comparison")
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Print summary statistics
print("Climate Model Ensemble Analysis Summary:")
print(f"Global mean warming: {ensemble_mean.mean().values:.2f} ± {ensemble_std.mean().values:.2f}°C")
print(f"Arctic amplification factor: {arctic_warming.mean().values / tropical_warming.mean().values:.2f}")
print(f"High agreement regions (>80%): {(model_agreement > 0.8).sum().values / model_agreement.size * 100:.1f}% of globe")