{ "cells": [ { "cell_type": "markdown", "id": "9a18f94b", "metadata": {}, "source": [ "# Postprocessing of SfM/MVS data\n", "In this notebook, we postprocess the outputs of Agisoft Metashape and compared them with the data collected in the field with the workflow described in Casella et al. (Year). \n", "\n", "> CITATION (WHEN AVAILABLE)\n", "\n", "## 1. Import libraries and define source folders\n", "First, we need to import the libraries needed and indicate in which folders the data are stored. The notebook creates a folder called \"Postprocessing_output\" where all the results of postprocessing are stored." ] }, { "cell_type": "code", "execution_count": 1, "id": "62cde758", "metadata": {}, "outputs": [ { "ename": "ImportError", "evalue": "dlopen(/Users/alessiorovere/opt/anaconda3/envs/Shallow_water_bathy/lib/python3.8/site-packages/rasterio/_version.cpython-38-darwin.so, 0x0002): Library not loaded: @rpath/libpoppler.126.dylib\n Referenced from: /Users/alessiorovere/opt/anaconda3/envs/Shallow_water_bathy/lib/libgdal.32.3.6.2.dylib\n Reason: tried: '/Users/alessiorovere/opt/anaconda3/envs/Shallow_water_bathy/lib/libpoppler.126.dylib' (no such file), '/Users/alessiorovere/opt/anaconda3/envs/Shallow_water_bathy/lib/python3.8/site-packages/rasterio/../../../libpoppler.126.dylib' (no such file), '/Users/alessiorovere/opt/anaconda3/envs/Shallow_water_bathy/lib/python3.8/site-packages/rasterio/../../../libpoppler.126.dylib' (no such file), '/Users/alessiorovere/opt/anaconda3/envs/Shallow_water_bathy/bin/../lib/libpoppler.126.dylib' (no such file), '/Users/alessiorovere/opt/anaconda3/envs/Shallow_water_bathy/bin/../lib/libpoppler.126.dylib' (no such file), '/usr/local/lib/libpoppler.126.dylib' (no such file), '/usr/lib/libpoppler.126.dylib' (no such file, not in dyld cache)", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mImportError\u001b[0m Traceback (most recent call last)", "Cell \u001b[0;32mIn[1], line 7\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mshapely\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mgeometry\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m Point\n\u001b[1;32m 6\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mmatplotlib\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mpyplot\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mplt\u001b[39;00m\n\u001b[0;32m----> 7\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mrasterio\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mplot\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m show\n\u001b[1;32m 8\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mnumpy\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mnp\u001b[39;00m\n\u001b[1;32m 9\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mrasterio\u001b[39;00m\n", "File \u001b[0;32m~/opt/anaconda3/envs/Shallow_water_bathy/lib/python3.8/site-packages/rasterio/__init__.py:28\u001b[0m\n\u001b[1;32m 24\u001b[0m os\u001b[38;5;241m.\u001b[39madd_dll_directory(os\u001b[38;5;241m.\u001b[39mpath\u001b[38;5;241m.\u001b[39mabspath(p))\n\u001b[1;32m 27\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mrasterio\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01m_show_versions\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m show_versions\n\u001b[0;32m---> 28\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mrasterio\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01m_version\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m gdal_version, get_geos_version, get_proj_version\n\u001b[1;32m 29\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mrasterio\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mcrs\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m CRS\n\u001b[1;32m 30\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mrasterio\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mdrivers\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m driver_from_extension, is_blacklisted\n", "\u001b[0;31mImportError\u001b[0m: dlopen(/Users/alessiorovere/opt/anaconda3/envs/Shallow_water_bathy/lib/python3.8/site-packages/rasterio/_version.cpython-38-darwin.so, 0x0002): Library not loaded: @rpath/libpoppler.126.dylib\n Referenced from: /Users/alessiorovere/opt/anaconda3/envs/Shallow_water_bathy/lib/libgdal.32.3.6.2.dylib\n Reason: tried: '/Users/alessiorovere/opt/anaconda3/envs/Shallow_water_bathy/lib/libpoppler.126.dylib' (no such file), '/Users/alessiorovere/opt/anaconda3/envs/Shallow_water_bathy/lib/python3.8/site-packages/rasterio/../../../libpoppler.126.dylib' (no such file), '/Users/alessiorovere/opt/anaconda3/envs/Shallow_water_bathy/lib/python3.8/site-packages/rasterio/../../../libpoppler.126.dylib' (no such file), '/Users/alessiorovere/opt/anaconda3/envs/Shallow_water_bathy/bin/../lib/libpoppler.126.dylib' (no such file), '/Users/alessiorovere/opt/anaconda3/envs/Shallow_water_bathy/bin/../lib/libpoppler.126.dylib' (no such file), '/usr/local/lib/libpoppler.126.dylib' (no such file), '/usr/lib/libpoppler.126.dylib' (no such file, not in dyld cache)" ] } ], "source": [ "# Libraries\n", "import os\n", "import geopandas as gpd\n", "import pandas as pd\n", "from shapely.geometry import Point\n", "import matplotlib.pyplot as plt\n", "from rasterio.plot import show\n", "import numpy as np\n", "import rasterio\n", "from sklearn.model_selection import train_test_split\n", "from scipy.optimize import minimize\n", "\n", "# Folder where processed data is stored\n", "out_folder = 'Data/07_08_2020/Postprocessing_output' \n", "os.makedirs(out_folder, exist_ok=True)\n", "\n", "# Link to files exported from Agisoft Metashape\n", "bathy = 'Data/07_08_2020/Agisoft_Input/tide_corrected_bathy.csv' # data from portable echosounder\n", "dbm = 'Data/07_08_2020/Agisoft_processing/DBM07082020.tif' # Digital bathymetric model from Agisoft Metashape\n", "cloudpoint = 'Data/07_08_2020/Agisoft_processing/PointCloud07082020.laz' # \n", "\n", "#Define the EPSG of your products:\n", "EPSG=32632" ] }, { "cell_type": "markdown", "id": "78111e5d", "metadata": {}, "source": [ "## 2. Check bathymetric data\n", "The bathymetric data and the echosounder points are plotted to check for consistency." ] }, { "cell_type": "code", "execution_count": null, "id": "6bbdbdc7", "metadata": {}, "outputs": [], "source": [ "# Make a geodataframe with the echosounder data\n", "bathy_df = pd.read_csv(bathy)\n", "bathy_df['geometry'] = bathy_df.apply(lambda row: Point(row['Longitude (dec degrees)'], row['Latitude (dec degrees)']), axis=1)\n", "bathy_gdf = gpd.GeoDataFrame(bathy_df, geometry='geometry')\n", "bathy_gdf.set_crs(epsg=4326, inplace=True)\n", "# Transform the GeoDataFrame to EPSG 32623\n", "bathy_gdf = bathy_gdf.to_crs(epsg=EPSG)\n", "\n", "\n", "# Open the DBM using rasterio\n", "with rasterio.open(dbm) as dbm_src:\n", " dbm_data = dbm_src.read(1) # Read the first band\n", " dbm_transform = dbm_src.transform\n", " dbm_crs = dbm_src.crs\n", " dbm_bounds = dbm_src.bounds\n", "\n", " # Replace NoData values with NaN for better handling\n", " dbm_data = np.where(dbm_data == -9999, np.nan, dbm_data)\n", "\n", " # Sample DBM values at the bathymetric points' locations\n", " bathy_gdf['DBM_Value'] = [\n", " value[0] for value in dbm_src.sample(bathy_gdf['geometry'].apply(lambda geom: (geom.x, geom.y)))\n", " ]\n", "\n", "# Replace -9999 with NaN in the sampled DBM values\n", "bathy_gdf['DBM_Value'] = bathy_gdf['DBM_Value'].replace(-9999, np.nan)\n", "\n", "# Remove points where DBM value is NaN\n", "bathy_gdf = bathy_gdf[~np.isnan(bathy_gdf['DBM_Value'])]\n", "\n", "# Plot the DBM and the filtered bathymetric points\n", "fig, ax = plt.subplots(figsize=(10, 10))\n", "show(dbm_data, transform=dbm_transform, ax=ax, cmap='viridis')\n", "bathy_gdf.plot(ax=ax, color='k', markersize=10, label='Filtered Bathymetric Points')\n", "ax.set_title('DBM with Filtered Bathymetric Points')\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "daeeccc2", "metadata": {}, "source": [ "## 3. Calibration of DBM\n", "The echosounder points are divided into two subsets (of equal value). One is used to find the optimal shift that minimizes the RMSE. The second subset is used as validation. The difference between validation echosounder points and DBM depths at the same locations is then plotted as histogram." ] }, { "cell_type": "code", "execution_count": null, "id": "4083f007", "metadata": {}, "outputs": [], "source": [ "# Replace -9999 with NaN in the sampled DBM values\n", "bathy_gdf['DBM_Value'] = bathy_gdf['DBM_Value'].replace(-9999, np.nan)\n", "\n", "# Remove points where DBM value is NaN\n", "bathy_gdf = bathy_gdf[~np.isnan(bathy_gdf['DBM_Value'])]\n", "\n", "# Step 1: Split the data into 50% for training and 50% for validation\n", "train_gdf, val_gdf = train_test_split(bathy_gdf, test_size=0.5, random_state=42)\n", "\n", "# Step 2: Define the function to minimize (RMSE)\n", "def rmse_shift(shift, dbm_values, depth_values):\n", " # Apply the shift to the DBM values\n", " shifted_dbm = dbm_values + shift\n", " # Calculate RMSE between shifted DBM and actual depth values\n", " rmse = np.sqrt(np.mean((shifted_dbm - depth_values) ** 2))\n", " return rmse\n", "\n", "# Step 3: Find the optimal shift\n", "result = minimize(rmse_shift, x0=[0], args=(train_gdf['DBM_Value'].values, train_gdf['Tide-corrected depth (m)'].values))\n", "optimal_shift = result.x[0]\n", "print(f\"Optimal Shift: {optimal_shift:.4f} meters\")\n", "\n", "# Step 4: Apply the optimal shift to the entire DBM\n", "corrected_dbm = dbm_data + optimal_shift\n", "\n", "# Step 5: Sample the corrected DBM values at the validation points\n", "def get_row_col(transform, x, y):\n", " col, row = ~transform * (x, y)\n", " return int(row), int(col)\n", "\n", "corrected_values = []\n", "for geom in val_gdf['geometry']:\n", " row, col = get_row_col(dbm_transform, geom.x, geom.y)\n", " if 0 <= row < corrected_dbm.shape[0] and 0 <= col < corrected_dbm.shape[1]:\n", " corrected_values.append(corrected_dbm[row, col])\n", " else:\n", " corrected_values.append(np.nan)\n", "\n", "val_gdf['Corrected_DBM_Value'] = corrected_values\n", "\n", "# Remove any points where sampling failed (i.e., outside the bounds)\n", "val_gdf = val_gdf.dropna(subset=['Corrected_DBM_Value'])\n", "\n", "val_gdf['Corrected_Difference'] = val_gdf['Tide-corrected depth (m)'] - val_gdf['Corrected_DBM_Value']\n", "val_gdf['UnCorrected_Difference'] = val_gdf['Tide-corrected depth (m)'] - val_gdf['DBM_Value']\n", "\n", "corrected_rmse = np.sqrt(np.mean(val_gdf['Corrected_Difference'] ** 2))\n", "print(f\"Corrected RMSE: {corrected_rmse:.2f} meters\")\n", "\n", "# Step 6: Create a histogram of differences for the corrected validation set\n", "plt.figure(figsize=(10, 6))\n", "plt.hist(val_gdf['Corrected_Difference'], bins=30, color='skyblue', edgecolor='black')\n", "plt.title(f'Histogram of Corrected Depth Differences (Validation Set)\\nRMSE: {corrected_rmse:.2f} meters')\n", "plt.xlabel('Corrected Difference (Tide-corrected depth - Corrected DBM_Value) [m]')\n", "plt.ylabel('Frequency')\n", "plt.grid(True)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "04c5c11a", "metadata": {}, "source": [ "## 4. RMSE by depth bin\n", "The validation data is divided into depth bins and the corresponding RMSE for each bin is calculated." ] }, { "cell_type": "code", "execution_count": null, "id": "8ce65675", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "from sklearn.metrics import mean_squared_error\n", "\n", "# Create a DataFrame\n", "data = pd.DataFrame({\n", " 'Observed': val_gdf['Tide-corrected depth (m)'],\n", " 'Predicted': val_gdf['Corrected_DBM_Value']\n", "})\n", "\n", "# Define bin limits dynamically based on the min and max of the observed data, rounding to integers\n", "min_depth = np.floor(data['Observed'].min()).astype(int) # Round down to nearest integer\n", "max_depth = np.ceil(data['Observed'].max()).astype(int) # Round up to nearest integer\n", "\n", "# Create integer bins using np.arange\n", "bins = np.arange(min_depth, max_depth + 1) # Ensure max depth is included\n", "\n", "# Apply the binning\n", "data['Depth Bin'] = pd.cut(data['Observed'], bins, include_lowest=True)\n", "\n", "# Calculate RMSE for each bin, skipping empty bins\n", "def calculate_rmse(group):\n", " if len(group) > 0:\n", " return np.sqrt(mean_squared_error(group['Observed'], group['Predicted']))\n", " else:\n", " return np.nan # Return NaN for empty bins\n", "\n", "# Apply the function to each bin\n", "bin_rmse = data.groupby('Depth Bin').apply(calculate_rmse)\n", "\n", "# Show RMSE by depth bin, excluding empty bins\n", "print(bin_rmse.dropna())" ] }, { "cell_type": "markdown", "id": "c97e1724", "metadata": {}, "source": [ "## 5. Export raster\n", "The depth-calibrated raster is exported for further analysis." ] }, { "cell_type": "code", "execution_count": null, "id": "1a91e6e7", "metadata": {}, "outputs": [], "source": [ "import rasterio\n", "from rasterio.transform import from_origin\n", "\n", "# Step 1: Export the corrected DBM as a GeoTIFF\n", "corrected_dbm_path = os.path.join(out_folder, 'corrected_dbm.tif')\n", "\n", "# Create the corrected DBM GeoTIFF\n", "with rasterio.open(\n", " corrected_dbm_path,\n", " 'w',\n", " driver='GTiff',\n", " height=corrected_dbm.shape[0],\n", " width=corrected_dbm.shape[1],\n", " count=1,\n", " dtype=corrected_dbm.dtype,\n", " crs=dbm_crs,\n", " transform=dbm_transform,\n", ") as dst:\n", " dst.write(corrected_dbm, 1)\n", "\n", "print(f\"Corrected DBM exported to {corrected_dbm_path}\")\n", "\n", "# Step 2: Export the validation GeoDataFrame (val_gdf) to CSV\n", "val_gdf_path = os.path.join(out_folder, 'validation_points.csv')\n", "\n", "# Export to CSV\n", "val_gdf.to_csv(val_gdf_path, index=False)\n", "\n", "print(f\"Validation GeoDataFrame exported to {val_gdf_path}\")" ] }, { "cell_type": "markdown", "id": "15ae36ed", "metadata": {}, "source": [ "## 6. Cloud point confidence\n", "The cloud point is processed and a text file containing the confidence bins is exported for further inspection." ] }, { "cell_type": "code", "execution_count": null, "id": "be436307", "metadata": {}, "outputs": [], "source": [ "import laspy\n", "import numpy as np\n", "import os\n", "\n", "# Specify the output file for the histogram data\n", "output_txt_file = os.path.join(out_folder, 'confidence_histogram.txt')\n", "\n", "# Open the LAZ file\n", "with laspy.open(cloudpoint) as laz_file:\n", " # Read the point data\n", " points = laz_file.read()\n", "\n", " # Extract the 'confidence' field if it exists\n", " if 'confidence' in points.point_format.dimension_names:\n", " confidence_values = points['confidence']\n", "\n", " # Get unique confidence values and their counts\n", " unique, counts = np.unique(confidence_values, return_counts=True)\n", "\n", " # Calculate the percentage for each unique confidence value\n", " total_points = len(confidence_values)\n", " percentages = (counts / total_points) * 100\n", "\n", " # Write the histogram data (confidence value and percentage) to the output file\n", " with open(output_txt_file, 'w') as txt_file:\n", " txt_file.write(\"Confidence_Value,Percentage\\n\")\n", " for value, percentage in zip(unique, percentages):\n", " txt_file.write(f\"{value},{percentage:.2f}\\n\")\n", "\n", " print(f\"Histogram data has been exported to {output_txt_file}\")\n", " else:\n", " print(\"The 'confidence' field was not found in the LAZ file.\")" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.8.19" } }, "nbformat": 4, "nbformat_minor": 5 }