diffplot_final.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8d6aa4a-f0aa-4e98-8268-427f1dd2b86d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "import gc\n",
    "import os\n",
    "import netCDF4 as nc4\n",
    "import numpy as np\n",
    "import numpy.ma as ma\n",
    "import pandas as pd\n",
    "import geopy as gpy\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.tri as tri\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n",
    "from shapely import Polygon, box"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0256a32-5b7e-4487-8193-6188210d3bf2",
   "metadata": {},
   "outputs": [],
   "source": [
    "withoffset = nc4.Dataset('/scratch/07174/soelem/test_jason/fort.63.nc')\n",
    "withoutoffset = nc4.Dataset('/work2/07174/soelem/data/fort.63.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "faa47871-cdae-4f8d-bf6f-fee3d5db49d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "lat = withoffset.variables['y'][:].data\n",
    "lon = withoffset.variables['x'][:].data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf1fd7c6-de80-41d2-96b4-9d2e07618f4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "start_date = pd.Timestamp('2023-07-16 00:00:00')\n",
    "end_date = pd.Timestamp('2023-10-31 23:00:00')\n",
    "time = pd.date_range(start_date, end_date, freq='H')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72d6cfa1-49fa-4e83-bd5c-fb313e65ccc9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the directory where you want to save the files\n",
    "output_dir = '/scratch/07174/soelem/results'\n",
    "\n",
    "#Map Extent Bounds\n",
    "lat1,lat2 = 7, 50\n",
    "lon1,lon2 = -100, -58\n",
    "cmap =  mpl.cm.get_cmap('seismic')\n",
    "for i in range(len(time)):\n",
    "    #Basic Plot parameters\n",
    "    ssh_with = withoffset.variables['zeta'][i,:].data\n",
    "    ssh_without = withoutoffset.variables['zeta'][i,:].data\n",
    "    diff = ssh_with - ssh_without\n",
    "    ts = pd.to_datetime(str(time[i]))\n",
    "    fig  = plt.figure(figsize=(10,6))\n",
    "    proj = ccrs.PlateCarree()\n",
    "    ax = plt.axes(projection = proj)\n",
    "    ax.set_extent([lon1, lon2, lat1, lat2], crs=proj)\n",
    "    cmap =  mpl.cm.get_cmap('seismic')\n",
    "    levels = np.arange(-0.9, 0.9, 0.01)\n",
    "    triang = tri.Triangulation(lon, lat, triangles=withoffset.variables['element'][:,:]-1)\n",
    "    if ma.isMaskedArray(diff):  # Check if diff is a masked array\n",
    "        if diff.mask.any():  # Then check if any element is masked\n",
    "            point_mask_indices = np.where(diff.mask)\n",
    "            tri_mask = np.any(np.in1d(withoffset.variables['element'][:,:]-1, point_mask_indices).reshape(-1, 3), axis=1)\n",
    "            triang.set_mask(tri_mask)\n",
    "    diff_array = np.array(diff.data)  # Convert memoryview to numpy array\n",
    "    diff_array[diff_array<-999] = np.nan  # Now you can do the comparison and assignment\n",
    "    if diff.size == 0 or np.all(np.isnan(diff)):  # Check if data is empty or all NaNs\n",
    "        print(f\"Skipping iteration {i} due to invalid data.\")\n",
    "        continue  # Skip this iteration if data is not valid\n",
    "    cf = ax.tricontourf(triang, diff.data, transform = proj, cmap = cmap, levels=levels,extend='both')\n",
    "    ax.coastlines(resolution='10m', linewidth=0.5)\n",
    "    ax.add_feature(cfeature.LAND, facecolor='lightgrey')\n",
    "    cbar = fig.colorbar(cf)\n",
    "    cbar.set_label('WL Difference (m)', size=25)\n",
    "    cbar.ax.tick_params(labelsize=15)\n",
    "    gl = ax.gridlines(draw_labels=True, linewidth=1, color='gray', alpha=0.5, linestyle='--')\n",
    "    \n",
    "    # Customize the gridlines\n",
    "    gl.xlabels_top = False  # Hide the top x-axis labels\n",
    "    gl.ylabels_right = False  # Hide the right y-axis labels\n",
    "    gl.xlabel_style = {'size': 10}  # Set the x-axis label font size\n",
    "    gl.ylabel_style = {'size': 10, 'rotation': 90}  # Set the y-axis label font size and rotation\n",
    "    #norm = plt.Normalize(min, max)\n",
    "    #cbar = plt.colorbar(norm=norm, cmap=cmap)\n",
    "    # plt.title ('2014-07-01 10:00', size=15)\n",
    "    savename = os.path.join(output_dir, ts.strftime('%Y-%m-%d %H.%M'))\n",
    "    plt.title(savename, size=15)\n",
    "    title = str(f'{savename}.png')\n",
    "    try:\n",
    "        fig.savefig(title, dpi=100, bbox_inches='tight')\n",
    "    except ValueError as e:\n",
    "        print(f\"An error occurred while saving the figure in iteration {i}: {e}\")\n",
    "        continue  # Skip the iteration if an error occurred during saving\n",
    "    \n",
    "    plt.close('all')  # Close all figures to ensure memory is freed\n",
    "    del ssh_with, ssh_without, diff, cf, fig, ax, triang  # Delete large variables explicitly\n",
    "    gc.collect()  # Suggest to the Python garbage collector to free up memory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import cartopy.crs as ccrs\n",
    "import numpy as np\n",
    "import matplotlib.tri as tri\n",
    "import cartopy.feature as cfeature\n",
    "\n",
    "# Assuming you have already loaded your datasets and extracted necessary variables\n",
    "# Specify the timestep you're interested in\n",
    "timestep = 400  # Replace with your specific timestep\n",
    "\n",
    "# Define the levels for contouring\n",
    "levels = np.arange(-0.9, 0.9, 0.01)  # Adjust as needed\n",
    "\n",
    "# Create the figure and subplots\n",
    "fig, axs = plt.subplots(1, 2, figsize=(10, 10), subplot_kw={'projection': ccrs.PlateCarree()}, gridspec_kw={'wspace': 0.05})\n",
    "\n",
    "# Plot for 'withoutoffset' dataset on the left\n",
    "ssh_without = withoutoffset.variables['zeta'][timestep, :].data\n",
    "triang_without = tri.Triangulation(lon, lat, triangles=withoutoffset.variables['element'][:, :] - 1)\n",
    "cf_without = axs[0].tricontourf(triang_without, ssh_without, levels=levels, cmap='viridis', extend='both')\n",
    "axs[0].coastlines()\n",
    "axs[0].add_feature(cfeature.LAND, facecolor='darkgrey')\n",
    "gl_without = axs[0].gridlines(draw_labels=True, xlocs=np.arange(-180, 181, 10), ylocs=np.arange(-90, 91, 10))\n",
    "gl_without.ylabels_right = False\n",
    "gl_without.xlabels_top = False\n",
    "gl_without.xlabel_style = {'size': 12}\n",
    "gl_without.ylabel_style = {'rotation': 90, 'size': 12}\n",
    "\n",
    "# Plot for 'withoffset' dataset on the right\n",
    "ssh_with = withoffset.variables['zeta'][timestep, :].data\n",
    "triang_with = tri.Triangulation(lon, lat, triangles=withoffset.variables['element'][:, :] - 1)\n",
    "cf_with = axs[1].tricontourf(triang_with, ssh_with, levels=levels, cmap='viridis', extend='both')\n",
    "axs[1].coastlines()\n",
    "axs[1].add_feature(cfeature.LAND, facecolor='darkgrey')\n",
    "gl_with = axs[1].gridlines(draw_labels=True, xlocs=np.arange(-180, 181, 10), ylocs=np.arange(-90, 91, 10))\n",
    "gl_with.ylabels_left = False\n",
    "gl_with.ylabels_right = False\n",
    "gl_with.xlabels_top = False\n",
    "gl_without.xlabel_style = {'size': 12}\n",
    "\n",
    "# Position the colorbar outside and below the plot boundaries\n",
    "cbar = fig.colorbar(cf_with, ax=axs.ravel().tolist(), orientation='horizontal', pad=0.05, fraction=0.05)\n",
    "cbar.set_label('Water Level (m)', size=20)\n",
    "\n",
    "# Adjust layout\n",
    "plt.tight_layout()\n",
    "\n",
    "# Save the figure\n",
    "plt.savefig(\"water_levels_comparison.png\", dpi=600)\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "adcsat",
   "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.9.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}