Skip to content

Instantly share code, notes, and snippets.

@aradhakrishnanGFDL

aradhakrishnanGFDL/tw_atmos_ts_monthly_sfc_ocean_updated.ipynb

Created April 28, 2026 19:21
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save aradhakrishnanGFDL/391b23f4a6e4cd0f34cf001d27d075b3 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save aradhakrishnanGFDL/391b23f4a6e4cd0f34cf001d27d075b3 to your computer and use it in GitHub Desktop.
Download ZIP
tw_atmos_ts_monthly_sfc_ocean_updated.ipynb
Display the source blob
Display the rendered blob
Raw
tw_atmos_ts_monthly_sfc_ocean_updated.ipynb
{
"cells": [
{
"cell_type": "markdown",
"id": "72eedf38",
"metadata": {},
"source": [
"## About the analysis\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "cd28b307",
"metadata": {},
"outputs": [],
"source": [
"diag_name = \"atw_atmos_ts_monthly_sfc_ocean\"\n",
"diag_desc = \"Annual-mean climatological diagnostics\"\n",
"diag_lead_author = \"Andrew Wittenberg\"\n",
"contributors = [\"Soelem Bhuiyan\", \"Aparna Radhakrishnan\"]\n",
"\n",
"more_info = \"\"\"\n",
"================================================================================\n",
"INPUT REQUIREMENTS & DOCUMENTATION\n",
"================================================================================\n",
"This script evaluates climatological biases by comparing high-resolution climate \n",
"model outputs against observational datasets. It regrids the model to the obs \n",
"grid using bilinear interpolation after applying a native land mask.\n",
"\n",
"Required Input Data:\n",
" 1. Model Data: Monthly time-series NetCDF files (e.g., atmos.*.t_surf.nc).\n",
" 2. Obs Data: Monthly observational NetCDF file (e.g., OISST).\n",
" 3. Static File: Model grid specification file (atmos.static.nc) containing \n",
" the 'land_mask' variable.\n",
"================================================================================\n",
"\"\"\""
]
},
{
"cell_type": "markdown",
"id": "904133dc",
"metadata": {},
"source": [
"## Jupyter Notebook Workflow Outline\n",
"\n",
"The analysis is built upon the following procedures:\n",
"\n",
"---\n",
"\n",
"### <a id=\"section1\"></a>Section 1: Prepare Diagnostics Settings\n",
"- Set the parameters \n",
"- Set the runtime configurations \n",
"\n",
"---\n",
"\n",
"### <a id=\"section2\"></a>Section 2: Load Data from Catalog Files\n",
"Two tested options available (users can override with documentation):\n",
"- Intake-ESM option \n",
"This section also handles retrieving files from TAPE file using dmget utility\n",
"\n",
"---\n",
"\n",
"### <a id=\"section3\"></a>Section 3: Diagnostics Computation\n",
"- Step 1. Compute climatology \n",
"- Step 2. Apply land masking \n",
"- Step 3. Regrid \n",
"- Step 4. Spatial slicing \n",
"- Step 5. Statistics computation \n",
"\n",
"---\n",
"\n",
"### <a id=\"section4\"></a>Section 4: Results Visualization\n",
"- Plot specifications \n",
"- Surface temperature bias maps (SPEAR-HI vs OISST) \n",
"\n",
"---\n",
"\n",
"### <a id=\"section5\"></a>Section 5: Citations and More About this Diagnostics\n",
"- TBA\n"
]
},
{
"cell_type": "markdown",
"id": "101c1837-1c75-469d-87af-30a0d0edde15",
"metadata": {},
"source": [
"## Section 1: Prepare Diagnostics Setup "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6a5582c2-e975-46c9-8bec-b0a1e72bd39f",
"metadata": {},
"outputs": [],
"source": [
"import xarray as xr\n",
"import numpy as np\n",
"import xesmf as xe\n",
"import matplotlib.pyplot as plt\n",
"import cartopy.crs as ccrs\n",
"import cartopy.feature as cfeature\n",
"from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter\n",
"import time\n",
"import os\n",
"import sys\n",
"import logging\n",
"\n",
"# 1. Standardize the Path \n",
"# If we can't find 'diagnostics', we add the parent directory to sys.path\n",
"try:\n",
" import diagnostics\n",
"except ImportError:\n",
" parent_dir = os.path.abspath(os.path.join(os.getcwd(), '..'))\n",
" if parent_dir not in sys.path:\n",
" sys.path.append(parent_dir)\n",
"\n",
"# 2. Unified Imports\n",
"try:\n",
" # Attempt absolute imports (Standard for Project Root/Pipelines)\n",
" from diagnostics.utils import dhelper, dmget\n",
" from diagnostics.utils.logger import init_logger # Assuming location based on your snippet\n",
" from diagnostics.utils.dhelper import load_settings, dict_to_settings\n",
"except ModuleNotFoundError:\n",
" # Fallback to relative/local imports (Standard for Interactive/Subfolder)\n",
" from utils import dhelper, dmget, logger\n",
" from utils.dhelper import load_settings, dict_to_settings\n",
" from utils.logger import init_logger\n",
"\n",
"init_logger(f\"{diag_name}_analysis.log\",level=logging.DEBUG) #TODO use analysis script name as prefix for log file name.\n",
"logger = logging.getLogger(__name__) # Things we could do if we know the debug level INFO DEBUG WARNING ERROR CRITICAL\n",
"logging.info(f\"Analysis log is written to {diag_name}_analysis.log\")\n",
"t0 = time.perf_counter() if logger.isEnabledFor(logging.INFO) else None #track time if DEBUG, otherwise set to None to avoid unnecessary perf_counter calls\n",
"\n",
"# Silence the xarray warnings about future updates\n",
"xr.set_options(use_new_combine_kwarg_defaults=True)\n"
]
},
{
"cell_type": "markdown",
"id": "20829a13",
"metadata": {},
"source": [
"## General settings"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "73c3efbb",
"metadata": {},
"outputs": [],
"source": [
"#If the mode is `interactive`, the user will manually input the needed info of the POD options, otherwise the info will be read from the default `settings.jsonc` file and the [analysis_script]_config.json\n",
"# Define a mode \n",
"mode = \"interactive\" #TODO this is not really needed for us when we use snakemake, but we can keep it for now to allow for interactive use of the notebook. We can remove it later if we decide it's not needed.\n",
"\n",
"if(mode != \"interactive\"):\n",
" case_config_path = \"atw_atmos_ts_monthly_sfc_ocean_updated_config.json\"\n",
" settings_path = \"settings.json\"\n",
" "
]
},
{
"cell_type": "markdown",
"id": "9afc358c",
"metadata": {},
"source": [
"## Diagnostics settings"
]
},
{
"cell_type": "markdown",
"id": "17fab358",
"metadata": {},
"source": [
"##### The following simply helps us load the settings and config files and return the variables in an intuitive batch-friendly way\n",
"##### Refer to settings.json and atw_atmos_ts_monthly_sfc_ocean_updated_config.json in diagnostics directory"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3b358602",
"metadata": {},
"outputs": [],
"source": [
"if( mode != \"interactive\"):\n",
" # 1. Load and convert in one line\n",
" case_info = dhelper.load_settings(case_config_path) #case_info\n",
" sd = dhelper.load_settings(settings_path) #sd\n",
"\n",
" # 2. Extract your favorite handles\n",
" dimensions = sd.dimensions\n",
" settings = sd.settings\n",
" varlist = sd.varlist"
]
},
{
"cell_type": "markdown",
"id": "5ab6d16c",
"metadata": {},
"source": [
"### **Quick Reference: Configuration Variables**\n",
"\n",
"| Prefix | Source File | Description | Code Example |\n",
"| :--- | :--- | :--- | :--- |\n",
"| **`case_info`** | `case_info.json` | **Paths/System:** File locations and directories. | `case_info.DATA_CATALOG` `case_info.OBS_FILE_PATH`|\n",
"| **`case_info.case_list`** | `case_info.json` | **Model Info:** Case metadata (e.g., `SPEAR_HI_8`). | `case_info.case_list.SPEAR_HI_8.enddate` |\n",
"| **`dimensions`** | `settings.json` | **Slicing:** Coordinate bounds (lat, lon, time). | `slice(*dimensions.lat.LAT_BOUNDS)` |\n",
"| **`varlist`** | `settings.json` | **Variables:** Variable names and unit offsets. | `varlist.t_surf.MODEL_OFFSET` |\n",
"| **`settings`** | `settings.json` | **Plotting:** Titles, labels, and contour levels. | `levels=settings.SST_SHADE_LEVELS` |"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "330ba338",
"metadata": {},
"outputs": [],
"source": [
"#These are the runtime configuration values.\n",
"#For e.g. which experiments to run diagnostis on, where the data is located, etc. \n",
"#If you're not running this interactively, these values will be overridden by the atw_atmos_ts_monthly_sfc_ocean_updated_config.json file. \n",
"# If you are running this interactively, these values will be used instead of the json file.\n",
"if mode == \"interactive\":\n",
" print(\"Applying Interactive case settings...\")\n",
" \n",
" # 1. Define the dictionary to match the JSON structure exactly\n",
" case_info_dict = {\n",
" \"pod_list\": [\n",
" \"atw_atmos_ts_monthly_sfc_ocean_updated\"\n",
" ],\n",
" \"case_list\": {\n",
" \"SPEAR-HI_8\": {\n",
" \"convention\": \"GFDL\",\n",
" \"enddate\": \"03981231000000\",\n",
" \"model\": \"SPEAR_HI_8\",\n",
" \"startdate\": \"00010101000000\",\n",
" \"frequency\": \"mon\",\n",
" \"static_vars\": [\"land_mask\", \"frac_ocean\"]\n",
" }\n",
" },\n",
" \"DATA_CATALOG\": \"/home/a1r/git/atw/atw_diags/cats/ciheim/48-SPEAR-cat.json\",\n",
" \"OBS_DATA_ROOT\": \"/home/atw/data/reynolds/v2/noaa.oisst.v2/\",\n",
" \"OBS_FILE_PATH\": \"/home/atw/data/reynolds/v2/noaa.oisst.v2/sst_janstart.nc\",\n",
" \"WORK_DIR\": \"/home/a1r/git/atw/atw_diags/outputs/\",\n",
" \"OUTPUT_DIR\": \"/home/a1r/git/atw/atw_diags/outputs/\"\n",
" }\n",
"\n",
" # 2. Convert to Namespace for dot notation\n",
" case_info = dict_to_settings(case_info_dict)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3861f668-7ce6-4b20-8e46-581ab74f8a75",
"metadata": {},
"outputs": [],
"source": [
"# ==============================================================================\n",
"# SETTINGS These normally do not get changed by the user. \n",
"#Plot settings, variable settings, dimension settings, etc. that are specific for the diagnostic script are mentioned here.\n",
"# These are settings that are used by the diagnostic script. They can be overridden by the user when running in interactive mode, \n",
"# but they are not expected to be changed by the user. Use the settings.json file in the diagnostics directory to set these values for the interactive mode. \n",
"# The settings in the settings.json file will override the default settings defined here if you don't run in interactive mode. \n",
"# If you run in interactive mode, the settings in the settings.json file will be ignored and the settings defined here will be used instead.\n",
"# ==============================================================================\n",
"if (mode) == \"interactive\":\n",
" print(\"Interactive settings\")\n",
" settings = {\n",
" \"settings\": {\n",
" \"description\": \"SPEAR Model vs OISST Surface Temperature Analysis\",\n",
" \"driver\": \"atw_atmos_ts_monthly_sfc_ocean.py\",\n",
" \"long_name\": \"Monthly Surface Temperature Diagnostics\",\n",
" \"convention\": \"GFDL\",\n",
" \"PLOT_TITLE_MODEL\": \"SPEAR Model (Years 0001-0100)\",\n",
" \"PLOT_TITLE_OBS\": \"OISST Observations (1982-2016)\",\n",
" \"CBAR_LABEL\": \"SST (°C)\",\n",
" \"SST_LINE_LEVELS\": [\n",
" 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33\n",
" ],\n",
" \"SST_SHADE_LEVELS\": [\n",
" 16.0, 16.5, 17.0, 17.5, 18.0, 18.5, 19.0, 19.5, 20.0, 20.5, 21.0, \n",
" 21.5, 22.0, 22.5, 23.0, 23.5, 24.0, 24.5, 25.0, 25.5, 26.0, 26.5, \n",
" 27.0, 27.5, 28.0, 28.5, 29.0, 29.5, 30.0, 30.5, 31.0, 31.5, 32.0, \n",
" 32.5, 33.0\n",
" ],\n",
" \"BIAS_LINE_LEVELS\": [-5, 5.5, 0.5],\n",
" \"BIAS_TICK_INTERVAL\": [-5, 6, 1],\n",
" \"BIAS_SHADE_LEVELS\": [\n",
" -5.0, -4.75, -4.5, -4.25, -4.0, -3.75, -3.5, -3.25, -3.0, -2.75, \n",
" -2.5, -2.25, -2.0, -1.75, -1.5, -1.25, -1.0, -0.75, -0.5, -0.25, \n",
" 0.0, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, \n",
" 3.0, 3.25, 3.5, 3.75, 4.0, 4.25, 4.5, 4.75, 5.0\n",
" ]\n",
" },\n",
" \"dimensions\": {\n",
" \"lat\": {\n",
" \"standard_name\": \"latitude\",\n",
" \"axis\": \"Y\",\n",
" \"LAT_BOUNDS\": [-20, 20]\n",
" },\n",
" \"lon\": {\n",
" \"standard_name\": \"longitude\",\n",
" \"axis\": \"X\",\n",
" \"LON_BOUNDS\": [25, 385]\n",
" },\n",
" \"time\": {\n",
" \"standard_name\": \"time\",\n",
" \"MODEL_TIME_SLICE\": [0, 1200],\n",
" \"OBS_TIME_SLICE\": [\"1982\", \"2016\"]\n",
" }\n",
" },\n",
" \"varlist\": {\n",
" \"vars\": {\n",
" \"freq\": \"mon\",\n",
" \"realm\": \"atmos\",\n",
" \"chunk_freq\": \"10yr\",\n",
" \"MODEL_VAR\": \"t_surf\",\n",
" \"OBS_VAR\": \"SST\",\n",
" \"MODEL_OFFSET\": -273.15,\n",
" \"dimensions\": [\"time\", \"lat\", \"lon\"]\n",
" },\n",
" \"STATIC_VARS\": {\n",
" \"freq\": \"fx\",\n",
" \"realm\": \"atmos\",\n",
" \"STATIC_VARS\": [\"land_mask\", \"frac_ocean\"],\n",
" \"dimensions\": [\"lat\", \"lon\"]\n",
" }\n",
" }\n",
" }\n",
"\n",
" # Convert the settings dictionary into an object so values can be accessed\n",
" # using dot notation instead of dictionary keys (e.g., settings.OBS_PATH)\n",
" settings_ns = dhelper.dict_to_settings(settings)\n",
" case_info_ns = dhelper.dict_to_settings(case_info)\n",
" # Extract your favorite handles\n",
" dimensions = settings_ns.dimensions\n",
" settings = settings_ns.settings\n",
" varlist = settings_ns.varlist\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1ed205af",
"metadata": {},
"outputs": [],
"source": [
"#Time section\n",
"if t0 is not None:\n",
" logger.debug(f\"Section 1 Process completed in {time.perf_counter() - t0:.4f} seconds\")"
]
},
{
"cell_type": "markdown",
"id": "8ec3f9b3-d93a-4448-9389-bbb6e044a7a8",
"metadata": {},
"source": [
"## Section 2 Loading Data from Catalogs\n",
"lets figure out ds_model with catalogs"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "765e5ad9",
"metadata": {},
"outputs": [],
"source": [
"#TODO handle time and profiling better\n",
"\n",
"t02 = time.perf_counter()\n",
"import intake, intake_esm"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8d19c8ca-883c-4c27-baab-5947dd32e9d0",
"metadata": {},
"outputs": [],
"source": [
"cat = intake.open_esm_datastore(case_info.DATA_CATALOG)\n",
"\n",
"#Exception - TODO to be addressed at the catalog level\n",
"df = cat.df \n",
"df['table_id'] = df['table_id'].fillna('NA')\n",
"df['table_id'].value_counts(dropna=False)\n",
"\n",
"cat_subset = cat.search(\n",
" variable_id=varlist.vars.MODEL_VAR,\n",
" frequency=varlist.vars.freq,\n",
" realm=varlist.vars.realm,\n",
" chunk_freq=varlist.vars.chunk_freq,\n",
" cell_methods=\"ts\", \n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0a03dafb-4be5-4c62-8d93-5ea949fb2a6c",
"metadata": {},
"outputs": [],
"source": [
"cat_subset"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "52bb1403",
"metadata": {},
"outputs": [],
"source": [
"#dmgets the files before we proceed\n",
"#TODO check if the filesystem is TAPE or not\n",
"\n",
"try:\n",
" # Works when running from the Project Root (e.g., Snakemake/Papermill)\n",
" from diagnostics.utils import dmget\n",
"except ModuleNotFoundError:\n",
" try:\n",
" # Works when running interactively from inside the 'diagnostics' folder\n",
" from utils import dmget\n",
" except ModuleNotFoundError:\n",
" # Fallback: add the path manually if both fail\n",
" import sys\n",
" import os\n",
" sys.path.append(os.path.abspath(os.path.join(os.getcwd(), '..')))\n",
" from diagnostics.utils import dmget\n",
"\n",
"#from diagnostics.utils import dmget\n",
"logging.info(\"Dmgetting files...\")\n",
"dmstatus = cat_subset.df[\"path\"].apply(dmget.dmgetmagic)\n",
"logging.info(\"done\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e29d108c-e2d3-4765-b7b6-135692f6d6d4",
"metadata": {},
"outputs": [],
"source": [
"dset_dict = cat_subset.to_dataset_dict(\n",
" xarray_open_kwargs={\"decode_times\": True,\"chunks\": {}}\n",
")\n",
"dset_dict.keys()\n",
"logging.info(\"to_dataset_dict complete\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3ab2b993-ca86-4b1f-804a-8f04f2e109ab",
"metadata": {},
"outputs": [],
"source": [
"#TODO don't hardcode this, sanity check dset_dict and retrieve from that with index\n",
"ds_model = dset_dict['SPEAR_HI_8.SPEAR_c384_OM4p08_Control_1990_A13.mon.NA.atmos.10yr']"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "cfd0f58b-974b-4da6-9b71-3b5498cadf5f",
"metadata": {},
"outputs": [],
"source": [
"ds_model"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f430ebb8",
"metadata": {},
"outputs": [],
"source": [
"# ==========================================\n",
"# 1A. LOAD DATA\n",
"# ==========================================\n",
"logging.info(\"[1/7] Loading data...\")\n",
"\n",
"ds_model = ds_model.isel(time=dimensions.time.MODEL_TIME_SLICE)\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "26ba509e",
"metadata": {},
"outputs": [],
"source": [
"\n",
"# ==========================================\n",
"# 1B. LOAD OBS DATA\n",
"# ==========================================\n",
"\n",
"ds_obs = xr.open_dataset(case_info.OBS_FILE_PATH) #orig OBS_PATH"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5accaf25",
"metadata": {},
"outputs": [],
"source": [
"\n",
"# ==========================================\n",
"# 1C. LOAD STATIC DATA\n",
"# ==========================================\n",
"\n",
"#\"STATIC_VARS\": [\"land_mask\",\"frac_ocean\"]\n",
"\n",
"cat_subset = cat.search(\n",
" variable_id=\"fixed\",\n",
" frequency=varlist.STATIC_VARS.freq,\n",
" realm=varlist.STATIC_VARS.realm,\n",
" table_id = varlist.STATIC_VARS.freq, #fx stands for fixed variables, which includes static variables like land_mask and frac_ocean\n",
")\n",
"dset_dict_static = cat_subset.to_dataset_dict(\n",
" xarray_open_kwargs={\"decode_times\": True,\"chunks\": {}}\n",
")\n",
"dset_dict_static.keys()\n",
"ds_static = dset_dict_static['SPEAR_HI_8.SPEAR_c384_OM4p08_Control_1990_A13.fx.fx.atmos']"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "cf55c075",
"metadata": {},
"outputs": [],
"source": [
"if t02 is not None: logger.debug(f\"Step 2 : {time.perf_counter() - t02:.4f}s\"); t02 = time.perf_counter()"
]
},
{
"cell_type": "markdown",
"id": "08c035df",
"metadata": {},
"source": [
"## Section 3 Diagnostics Computation"
]
},
{
"cell_type": "markdown",
"id": "79e08452",
"metadata": {},
"source": [
"###1. Computing Climatology"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8a140ba9",
"metadata": {},
"outputs": [],
"source": [
"# ==========================================\n",
"# 2. COMPUTING CLIMATOLOGY\n",
"# ==========================================\n",
"logging.info(\"[2/7] Computing climatologies...\") #TODO redo the numbering of the steps in the logging messages\n",
"t0 = time.perf_counter()\n",
"\n",
"model_clim = ds_model[varlist.vars.MODEL_VAR].mean(dim='time', keep_attrs=True)\n",
"model_clim = model_clim + varlist.vars.MODEL_OFFSET \n",
"obs_clim = ds_obs[varlist.vars.OBS_VAR].sel(time=slice(*dimensions.time.OBS_TIME_SLICE)).mean(dim='time', keep_attrs=True)\n",
"\n",
"model_clim = model_clim.compute()\n",
"obs_clim = obs_clim.compute()\n"
]
},
{
"cell_type": "markdown",
"id": "4f0e7272",
"metadata": {},
"source": [
"####2.Apply Land Masking"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2a65a700",
"metadata": {},
"outputs": [],
"source": [
"# ==========================================\n",
"# 3. NATIVE LAND MASK\n",
"# ==========================================\n",
"logging.info(\"[3/7] Applying native model land mask...\")\n",
"\n",
"try:\n",
" ds_static # = xr.open_dataset(settings.STATIC_PATH)\n",
" if 'land_mask' in ds_static:\n",
" model_clim = model_clim.where(ds_static['land_mask'] < 0.5)\n",
" elif 'frac_ocean' in ds_static:\n",
" model_clim = model_clim.where(ds_static['frac_ocean'] > 0.5)\n",
" else:\n",
" print(\"WARNING: 'land_mask' or 'frac_ocean' not found in static file. Land bleed may occur.\")\n",
"except FileNotFoundError:\n",
" print(f\"WARNING: Could not find {settings.STATIC_PATH}. Land bleed will occur!\")\n"
]
},
{
"cell_type": "markdown",
"id": "95fdde25",
"metadata": {},
"source": [
"####3. Regridding"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "90bc8a24",
"metadata": {},
"outputs": [],
"source": [
"# ==========================================\n",
"# 4. REGRIDDING\n",
"# ==========================================\n",
"\n",
"logging.info(\"[4/7] Evaluating grids and regridding Model to Obs...\")\n",
"t0 = time.perf_counter()\n",
"\n",
"# 1. Calculate the resolution (total pixels) of each grid\n",
"source_pixels = model_clim.lat.size * model_clim.lon.size\n",
"target_pixels = obs_clim.lat.size * obs_clim.lon.size\n",
"\n",
"# 2. The Switch Logic (Replicating Andrew's def_regrid_ax_str.jnl)\n",
"if source_pixels > target_pixels:\n",
" regrid_method = 'conservative' # Equivalent to Ferret's @ave\n",
" logger.info(\"-> High-to-Low resolution detected. Switch = 1 (CONSERVATIVE)\")\n",
"else:\n",
" regrid_method = 'bilinear' # Equivalent to Ferret's @lin\n",
" logger.info(\"-> Low-to-High resolution detected. Switch = 0 (BILINEAR)\")\n",
"\n",
"# 3. Execute the regridding using the auto-selected method\n",
"regridder = xe.Regridder(model_clim, obs_clim, regrid_method, periodic=True)\n",
"model_regridded = regridder(model_clim)\n"
]
},
{
"cell_type": "markdown",
"id": "79267229",
"metadata": {},
"source": [
"####4. Spatial Slicing"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4a9c297a",
"metadata": {},
"outputs": [],
"source": [
"# ==========================================\n",
"# 5. SPATIAL SLICING & MASKING\n",
"# ==========================================\n",
"logging.info(\"[5/7] Isolating spatial domain and applying common mask...\")\n",
"\n",
"model_domain = model_regridded.sel(lat=slice(dimensions.lat.LAT_BOUNDS[0], dimensions.lat.LAT_BOUNDS[1]))\n",
"obs_domain = obs_clim.sel(lat=slice(dimensions.lat.LAT_BOUNDS[0], dimensions.lat.LAT_BOUNDS[1]))\n",
"\n",
"common_mask = model_domain.notnull() & obs_domain.notnull()\n",
"\n",
"a_var = obs_domain.where(common_mask) # Ferret Variable A (Obs)\n",
"b_var = model_domain.where(common_mask) # Ferret Variable B (Model)\n",
"bma_var = b_var - a_var # Ferret Bias (Model minus Obs)\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "719e0e0b",
"metadata": {},
"source": [
"####5. Statistics Computation"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a90444e3",
"metadata": {},
"outputs": [],
"source": [
"# ==========================================\n",
"# 6. SPHERICAL STATISTICS (Cosine Weighting)\n",
"# ==========================================\n",
"logging.info(\"[6/7] Calculating spherically-weighted statistics...\")\n",
"\n",
"weights = np.cos(np.deg2rad(a_var.lat))\n",
"weights.name = \"weights\"\n",
"\n",
"def get_weighted_stats(da, w):\n",
" daw = da.weighted(w)\n",
" ave = daw.mean(skipna=True).values\n",
" std = daw.std(skipna=True).values\n",
" vmax = da.max(skipna=True).values\n",
" vmin = da.min(skipna=True).values\n",
" return ave, std, vmax, vmin\n",
"\n",
"a_ave, a_std, a_max, a_min = get_weighted_stats(a_var, weights)\n",
"b_ave, b_std, b_max, b_min = get_weighted_stats(b_var, weights)\n",
"bma_ave, bma_std, bma_max, bma_min = get_weighted_stats(bma_var, weights)\n",
"\n",
"ab_rmsd = np.sqrt(((bma_var)**2).weighted(weights).mean(skipna=True)).values\n",
"covar = ((b_var - b_ave) * (a_var - a_ave)).weighted(weights).mean(skipna=True)\n",
"ab_corr = (covar / (b_var.weighted(weights).std(skipna=True) * a_var.weighted(weights).std(skipna=True))).values\n",
"\n",
"logger.info(f\"Stats Check -> Bias Ave: {bma_ave:.3f}, RMSD: {ab_rmsd:.3f}, Corr: {ab_corr:.2f}\")\n",
"\n",
"if t02 is not None: logger.debug(f\"Section 3 : {time.perf_counter() - t02:.4f}s\"); t02 = time.perf_counter()"
]
},
{
"cell_type": "markdown",
"id": "7b324190",
"metadata": {},
"source": [
"## Section 4 Results Visualization"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4e86ee36",
"metadata": {},
"outputs": [],
"source": [
"# ==========================================\n",
"# 7. EXECUTION & PLOTTING \n",
"# ==========================================\n",
"logging.info(\"[7/7] Plotting data...\")\n",
"\n",
"fig, axes = plt.subplots(3, 1, figsize=(10, 12), subplot_kw={'projection': ccrs.PlateCarree(central_longitude=205)})\n",
"\n",
"plt.subplots_adjust(right=0.75, bottom=0.08, hspace=0.3)\n",
"\n",
"for ax in axes:\n",
" ax.set_extent([dimensions.lon.LON_BOUNDS[0], dimensions.lon.LON_BOUNDS[1], dimensions.lat.LAT_BOUNDS[0], dimensions.lat.LAT_BOUNDS[1]], crs=ccrs.PlateCarree())\n",
" ax.set_aspect('auto', adjustable='box')\n",
" ax.add_feature(cfeature.LAND, facecolor='black', zorder=2)\n",
" ax.coastlines(color='white', linewidth=0.5, zorder=3)\n",
" ax.set_xticks(np.arange(dimensions.lon.LON_BOUNDS[0], dimensions.lon.LON_BOUNDS[1]+1, 60), crs=ccrs.PlateCarree())\n",
" ax.set_yticks(np.arange(dimensions.lat.LAT_BOUNDS[0], dimensions.lat.LAT_BOUNDS[1]+1, 10), crs=ccrs.PlateCarree())\n",
" ax.xaxis.set_major_formatter(LongitudeFormatter())\n",
" ax.yaxis.set_major_formatter(LatitudeFormatter())\n",
" ax.gridlines(color='black', linestyle='-', alpha=0.3, zorder=4)\n",
"\n",
"# --- Panel 1: Obs ---\n",
"cs1 = a_var.plot.contourf(ax=axes[0], transform=ccrs.PlateCarree(), levels=settings.SST_SHADE_LEVELS, cmap='turbo', add_colorbar=False, extend='both')\n",
"line1 = a_var.plot.contour(ax=axes[0], transform=ccrs.PlateCarree(), levels=settings.SST_LINE_LEVELS, colors='black', linewidths=0.8)\n",
"axes[0].clabel(line1, inline=True, fontsize=8, fmt='%1.0f')\n",
"axes[0].set_title(settings.PLOT_TITLE_OBS)\n",
"fig.colorbar(cs1, ax=axes[0], label=settings.CBAR_LABEL, pad=0.02, shrink=0.85)\n",
"\n",
"axes[0].text(1.18, 0.5, f\"ave\\n{a_ave:.1f}\\n\\nstd\\n{a_std:.2f}\\n\\nmax\\n{a_max:.1f}\\n\\nmin\\n{a_min:.1f}\", \n",
" transform=axes[0].transAxes, va='center', ha='center', fontsize=10, family='monospace', clip_on=False)\n",
"\n",
"# --- Panel 2: Model ---\n",
"cs2 = b_var.plot.contourf(ax=axes[1], transform=ccrs.PlateCarree(), levels=settings.SST_SHADE_LEVELS, cmap='turbo', add_colorbar=False, extend='both')\n",
"line2 = b_var.plot.contour(ax=axes[1], transform=ccrs.PlateCarree(), levels=settings.SST_LINE_LEVELS, colors='black', linewidths=0.8)\n",
"axes[1].clabel(line2, inline=True, fontsize=8, fmt='%1.0f')\n",
"axes[1].set_title(settings.PLOT_TITLE_MODEL)\n",
"fig.colorbar(cs2, ax=axes[1], label=settings.CBAR_LABEL, pad=0.02, shrink=0.85)\n",
"\n",
"axes[1].text(1.18, 0.5, f\"ave\\n{b_ave:.1f}\\n\\nstd\\n{b_std:.2f}\\n\\nmax\\n{b_max:.1f}\\n\\nmin\\n{b_min:.1f}\", \n",
" transform=axes[1].transAxes, va='center', ha='center', fontsize=10, family='monospace', clip_on=False)\n",
"\n",
"# --- Panel 3: Bias ---\n",
"cs3 = bma_var.plot.contourf(ax=axes[2], transform=ccrs.PlateCarree(), levels=settings.BIAS_SHADE_LEVELS, cmap='RdBu_r', add_colorbar=False, extend='both')\n",
"line_widths = [1.5 if lev == 0 else 0.8 for lev in settings.BIAS_LINE_LEVELS]\n",
"line3 = bma_var.plot.contour(ax=axes[2], transform=ccrs.PlateCarree(), levels=settings.BIAS_LINE_LEVELS, colors='black', linewidths=line_widths)\n",
"axes[2].clabel(line3, inline=True, fontsize=8, fmt='%1.1f')\n",
"axes[2].set_title(\"Difference (Model - Obs)\")\n",
"fig.colorbar(cs3, ax=axes[2], label='Bias (°C)', pad=0.02, shrink=0.85, ticks=settings.BIAS_TICK_INTERVAL)\n",
"\n",
"axes[2].text(1.18, 0.5, f\"ave\\n{bma_ave:.3f}\\n\\nstd\\n{bma_std:.3f}\\n\\nmax\\n{bma_max:.2f}\\n\\nmin\\n{bma_min:.2f}\", \n",
" transform=axes[2].transAxes, va='center', ha='center', fontsize=10, family='monospace', clip_on=False)\n",
"\n",
"axes[2].text(0.25, -0.25, f\"corr(a,b) = {ab_corr:.2f}\", transform=axes[2].transAxes, ha='center', fontsize=12, family='monospace', clip_on=False)\n",
"axes[2].text(0.75, -0.25, f\"RMSD(a,b) = {ab_rmsd:.3f}\", transform=axes[2].transAxes, ha='center', fontsize=12, family='monospace', clip_on=False)\n",
"\n",
"# Save Figure\n",
"plt.savefig(case_info.OUTPUT_DIR + \"t_surf_final_replica.png\", dpi=600, bbox_inches='tight')\n",
"plt.show()\n",
"plt.close()\n",
"\n",
"if t02 is not None: logger.debug(f\"Section 4 : {time.perf_counter() - t02:.4f}s\"); t02 = time.perf_counter()\n",
"if t0 is not None: logger.debug(f\"TOTAL TIME: {time.perf_counter() - t0:.4f}s\"); t0 = time.perf_counter()\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.14.4"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment