Merge pull request #159 from megandevlan/updateLandAtmTCI

Update land Terrestrial Coupling Notebook to add scatter plot.

nblibrary/lnd/Global_TerrestrialCouplingIndex_VisualCompareObs.ipynb

@@ -17,7 +17,7 @@
     "- Note: Built to use monthly output; ideally, CI should be based on daily data. \n",
     "- Optional: Comparison against FLUXNET obs\n",
     "<br><br>\n",
-    "Notebook created by [email protected]; Last update: 2 Aug 2024 "
+    "Notebook created by [email protected]; Last update: 11 Dec 2024 "
    ]
   },
   {
@@ -923,6 +923,87 @@
     "    return"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6003435f-f16d-4c8b-80da-40f9cb571560",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def plotScatter(seasonstr, caseSel=None):\n",
+    "    node_lats = uxgrid.face_lat.values\n",
+    "    node_lons = uxgrid.face_lon.values\n",
+    "\n",
+    "    predictions = []\n",
+    "\n",
+    "    fig, axs = plt.subplots(1, 1, figsize=(8, 8))\n",
+    "\n",
+    "    CI_model = couplingIndex_DS[\"CouplingIndex\"].sel(season=seasonstr)\n",
+    "\n",
+    "    iSeason = np.where(seasons == seasonstr)[0]\n",
+    "    iStations = np.where(np.isfinite(terraCI_fluxnetConverted[:, iSeason]) == True)[0]\n",
+    "\n",
+    "    for iPoint in range(len(iStations)):\n",
+    "        this_lon = lon_fluxnet[iStations[iPoint]]  # lon1\n",
+    "        this_lat = lat_fluxnet[iStations[iPoint]]  # lat1\n",
+    "        obs_point = np.array((this_lon, this_lat))\n",
+    "\n",
+    "        # Get subset of relevant points\n",
+    "        i = np.where(\n",
+    "            (node_lats >= (this_lat - 2))\n",
+    "            & (node_lats <= (this_lat + 2))\n",
+    "            & (node_lons >= (this_lon - 2))\n",
+    "            & (node_lons <= (this_lon + 2))\n",
+    "        )[0]\n",
+    "\n",
+    "        minDistance = 100\n",
+    "        for iSelClose in range(len(i)):\n",
+    "            # Find point in uxarray? Use euclidian distance\n",
+    "            distance = np.linalg.norm(\n",
+    "                obs_point - np.array((node_lons[i[iSelClose]], node_lats[i[iSelClose]]))\n",
+    "            )\n",
+    "\n",
+    "            if (distance < minDistance) & (\n",
+    "                np.isfinite(\n",
+    "                    couplingIndex_DS[\"CouplingIndex\"]\n",
+    "                    .sel(season=\"JJA\")\n",
+    "                    .values[i[iSelClose]]\n",
+    "                )\n",
+    "                == True\n",
+    "            ):\n",
+    "                minDistance = distance\n",
+    "                selLon = node_lons[i[iSelClose]]\n",
+    "                selLat = node_lats[i[iSelClose]]\n",
+    "                selLoc = i[iSelClose]\n",
+    "\n",
+    "        predictions = np.append(\n",
+    "            predictions,\n",
+    "            couplingIndex_DS[\"CouplingIndex\"].sel(season=seasonstr).values[selLoc],\n",
+    "        )\n",
+    "        axs.plot(\n",
+    "            terraCI_fluxnetConverted[iPoint, iSeason],\n",
+    "            couplingIndex_DS[\"CouplingIndex\"].sel(season=seasonstr).values[selLoc],\n",
+    "            \"bo\",\n",
+    "            alpha=0.5,\n",
+    "        )\n",
+    "\n",
+    "    axs.set_xlabel(\"FLUXNET CI Value\", fontsize=12)\n",
+    "    axs.set_ylabel(\"CESM CI Value\", fontsize=12)\n",
+    "    axs.set_title(\n",
+    "        \"Individual station CI vs. nearest gridcell CI: \" + seasonstr, fontsize=14\n",
+    "    )\n",
+    "    axs.set_xlim([-25, 25])\n",
+    "    axs.set_ylim([-25, 25])\n",
+    "    axs.plot(np.arange(-25, 26), np.arange(-25, 26), \"k--\")\n",
+    "\n",
+    "    rmse = np.sqrt(\n",
+    "        ((predictions - terraCI_fluxnetConverted[iStations, iSeason]) ** 2).mean()\n",
+    "    )\n",
+    "    axs.text(0.05, 0.95, \"RMSE: \" + str(rmse), transform=axs.transAxes)\n",
+    "\n",
+    "    return axs"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -941,10 +1022,16 @@
     "if len(caseNames) == 1:\n",
     "    plotTCI_case(\"JJA\", None)\n",
     "    plotTCI_case(\"DJF\", None)\n",
+    "\n",
+    "    plotScatter(\"JJA\", None)\n",
+    "    plotScatter(\"DJF\", None)\n",
     "else:\n",
     "    for iCase in range(len(caseNames)):\n",
     "        plotTCI_case(\"JJA\", iCase)\n",
-    "        plotTCI_case(\"DJF\", iCase)"
+    "        plotTCI_case(\"DJF\", iCase)\n",
+    "\n",
+    "        plotScatter(\"JJA\", iCase)\n",
+    "        plotScatter(\"DJF\", iCase)"
    ]
   }
  ],