uupdate_ground_truth_prepare_simulation_study

Publication/src/AUVSimulator/CTDSimulator.py

@@ -43,17 +43,80 @@ def __init__(self, random_seed: int = 0,
         self.grid_sinmod_tree = KDTree(self.grid_sinmod)
 
         # Set up essential parameters
+        self.ar1_corr = .965
+        self.sigma = sigma
         l_range = 700
         eta = 4.5 / l_range
         t0 = time()
         dm = cdist(self.grid_sinmod, self.grid_sinmod)
-        cov = sigma ** 2 * ((1 + eta * dm) * np.exp(-eta * dm))
-        L = np.linalg.cholesky(cov)
+        cov = self.sigma ** 2 * ((1 + eta * dm) * np.exp(-eta * dm))
+        self.L = np.linalg.cholesky(cov)
         print("Cholesky decomposition takes: ", time() - t0)
 
+        self.construct_ground_truth_field()
+
     def construct_ground_truth_field(self) -> None:
+        """
+        This function constructs the ground truth field given all the timestamps and locations.
+        """
+        print("Start constructing ground truth field!")
+        x0 = (self.salinity_sinmod[0, :, :].flatten() +
+             (self.L @ np.random.randn(len(self.L)).reshape(-1, 1)).flatten())
+        xt_1 = x0
+        xt = x0
+        self.mu_truth = np.empty([0, len(xt)])
+        self.mu_truth = np.vstack((self.mu_truth, xt))
+
+        import matplotlib.pyplot as plt
+        from matplotlib.pyplot import get_cmap
+        from matplotlib.gridspec import GridSpec
+        figpath = os.getcwd() + "/../../../../OneDrive - NTNU/MASCOT_PhD/Projects" \
+                                "/GOOGLE/Docs/fig/Sim_2DNidelva/Simulator/groundtruth/"
+        t0 = time()
+        for i in range(1, len(self.timestamp_sinmod)):
+            print("Time string: ", datetime.fromtimestamp(self.timestamp_sinmod[i]).strftime("%Y-%m-%d %H:%M:%S"))
+            fig = plt.figure(figsize=(35, 10))
+            gs = GridSpec(1, 3, figure=fig)
+            ax = fig.add_subplot(gs[0])
+            im = ax.scatter(self.grid_sinmod[:, 1], self.grid_sinmod[:, 0], c=xt, cmap=get_cmap("BrBG", 10),
+                        vmin=10, vmax=30)
+            plt.colorbar(im)
+            ax.set_title("Ground truth at time: " +
+                         datetime.fromtimestamp(self.timestamp_sinmod[i]).strftime("%Y-%m-%d %H:%M:%S"))
+            ax.set_xlabel("East (m)")
+            ax.set_ylabel("North (m)")
+
+            ax = fig.add_subplot(gs[1])
+            im = ax.scatter(self.grid_sinmod[:, 1], self.grid_sinmod[:, 0], c=self.salinity_sinmod[i, :, :].flatten(),
+                        cmap=get_cmap("BrBG", 10), vmin=10, vmax=30)
+            plt.colorbar(im)
+            ax.set_title("SINMOD at time: " +
+                         datetime.fromtimestamp(self.timestamp_sinmod[i]).strftime("%Y-%m-%d %H:%M:%S"))
+            ax.set_xlabel("East (m)")
+            ax.set_ylabel("North (m)")
+
+            ax = fig.add_subplot(gs[2])
+            im = ax.scatter(self.grid_sinmod[:, 1], self.grid_sinmod[:, 0],
+                            c=xt - self.salinity_sinmod[i, :, :].flatten(),
+                            cmap=get_cmap("RdBu", 10), vmin=-5, vmax=5)
+            plt.colorbar(im)
+            ax.set_title("Difference at time: " +
+                            datetime.fromtimestamp(self.timestamp_sinmod[i]).strftime("%Y-%m-%d %H:%M:%S"))
+            ax.set_xlabel("East (m)")
+            ax.set_ylabel("North (m)")
+
+            plt.savefig(figpath + "P_{:03d}.png".format(i))
+            plt.close("all")
+
+            xt = self.salinity_sinmod[i, :, :].flatten() + \
+                 self.ar1_corr * (xt_1 - self.salinity_sinmod[i-1, :, :].flatten()) + \
+                 (np.sqrt(1 - self.ar1_corr ** 2) * (self.L @ np.random.randn(len(self.L)).reshape(-1, 1))).flatten()
+            xt_1 = xt
+            self.mu_truth = np.vstack((self.mu_truth, xt))
 
-        pass
+        self.mu_truth
+        print("Time consumed: ", time() - t0)
+        xt_1
 
     def get_salinity_at_dt_loc(self, dt: float, loc: np.ndarray) -> Union[np.ndarray, None]:
         """

Publication/truth.ipynb

@@ -0,0 +1,59 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# This script is used to pre-generate the ground truth field for a specific day\n",
+    "\n",
+    "Author: Yaolin Ge\n",
+    "Email: [email protected]\n",
+    "Date: 2023-08-28\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from src.WGS import WGS\n",
+    "import os\n",
+    "import netCDF4\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "from shapely.geometry import Polygon, Point\n",
+    "\n",
+    "filepath = \"/Users/yaolin/OneDrive - NTNU/MASCOT_PhD/Data/Nidelva/SINMOD_DATA/samples/samples_2022.05.11.nc\"\n",
+    "path_border = \"src/csv/polygon_border.csv\"\n",
+    "path_obstacle = \"src/csv/polygon_obstacle.csv\"\n",
+    "\n",
+    "sinmod = netCDF4.Dataset(filepath)\n",
+    "plg_border = pd.read_csv(path_border).to_numpy()\n",
+    "plg_obstacle = pd.read_csv(path_obstacle).to_numpy()\n",
+    "plg_b_sh = Polygon(plg_border)\n",
+    "plg_o_sh = Polygon(plg_obstacle)\n",
+    "\n",
+    "timestamp = sinmod['time']\n",
+    "lat = np.array(sinmod['gridLats'])\n",
+    "lon = np.array(sinmod['gridLons'])\n",
+    "depth = np.array(sinmod['depth'])\n",
+    "salinity = np.array(sinmod['salinity'])\n",
+    "\n",
+    "sal = np.nanmean(salinity[:, 0, :, :], axis=0)\n",
+    "lat = lat.flatten()\n",
+    "lon = lon.flatten()\n",
+    "sal = sal.flatten()"
+   ]
+  }
+ ],
+ "metadata": {
+  "language_info": {
+   "name": "python"
+  },
+  "orig_nbformat": 4
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}