{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preprocessing of data for SfM/MVS\n",
    "In this notebook, we import the data collected in the fied with the workflow described in Casella et al. (Year) and do some preprocessing steps to prepare them for analysis in Agisoft Metashape.\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 \"Agisoft_Input\" where all the results of preprocessing are stored."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Libraries\n",
    "import os\n",
    "import glob\n",
    "import pandas as pd\n",
    "import os\n",
    "from PIL import Image as PILImage\n",
    "import exifread\n",
    "from IPython.display import display, Image\n",
    "from datetime import datetime\n",
    "import os\n",
    "from PIL import Image, UnidentifiedImageError\n",
    "from PIL.ExifTags import TAGS\n",
    "import pandas as pd\n",
    "from datetime import datetime, timedelta\n",
    "import xml.etree.ElementTree as ET\n",
    "import shutil \n",
    "import geopandas as gpd\n",
    "from shapely.geometry import Point\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib_scalebar.scalebar import ScaleBar\n",
    "\n",
    "# Folder where processed data is stored\n",
    "in_folder = 'Data/07_08_2020/Agisoft_Input' \n",
    "os.makedirs(in_folder, exist_ok=True)\n",
    "\n",
    "# Folder where CSV files from the data are stored\n",
    "csv_folder = 'Data/07_08_2020/Echosounder' # data from portable echosounder\n",
    "tide_folder = 'Data/07_08_2020/Tide' # data from nearest tide gauge\n",
    "\n",
    "# Folder where the photo to sync with gpx is stored\n",
    "photo_sync='Data/07_08_2020/Camera/Camera_sync/Sync.JPG' # Photo showing date and time of the GPS device\n",
    "\n",
    "# Directory containing photos\n",
    "photo_dir = 'Data/07_08_2020/Camera/all_photos' # Directory where all pictures that will be used as inputs to metashape are stored\n",
    "\n",
    "# GPX file path\n",
    "gpx_file = 'Data/07_08_2020/GPS/GPX/20200807-063735-0004665-003482.gpx' # Path and filename of the GPX file exported from the GPS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Import and process tidal data\n",
    "After some format adjustments, the timestamps of the water level data from the nearest tide gauge are compared with the timestamps of the echosounder data, extracting the water levels at the time of survey. Then, the average of water levels during the time of survey is calculated and a \"Z_correction.txt\" file is written that reports the camera elevation that will be inserted in Agisoft Metashape. The average water level is used to correct the echosounder bathymetry, and a \"tide_corrected_bathymetry.csv\" file is saved."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Define column names for bathymetric data\n",
    "bathy_columns = ['Latitude (dec degrees)', 'Longitude (dec degrees)', 'Depth (m)', 'UNIX Time']\n",
    "\n",
    "# Define column names for tide data\n",
    "tide_data_col = 'DATA'\n",
    "tide_time_col = 'ORA'\n",
    "tide_value_col = 'VALORE'\n",
    "\n",
    "# Process bathymetric data\n",
    "bathy_dirname = os.path.join(os.path.realpath(''), csv_folder)\n",
    "bathy_files = glob.glob(os.path.join(bathy_dirname, \"*.csv\"))\n",
    "\n",
    "if not bathy_files:\n",
    "    raise ValueError(\"No CSV files found in the specified bathymetric data folder.\")\n",
    "\n",
    "bathy_li = []\n",
    "\n",
    "for filename in bathy_files:\n",
    "    bathy_raw = pd.read_csv(filename, index_col=None, header=None)\n",
    "    bathy_li.append(bathy_raw)\n",
    "\n",
    "if len(bathy_li) > 1:\n",
    "    bathy_raw = pd.concat(bathy_li, axis=0, ignore_index=True)\n",
    "else:\n",
    "    bathy_raw = bathy_li[0]\n",
    "\n",
    "bathy_raw.columns = bathy_columns\n",
    "bathy_raw['Time (UTC)'] = pd.to_datetime(bathy_raw['UNIX Time'], unit='ms')\n",
    "bathy_raw = bathy_raw.drop(bathy_raw[(bathy_raw[bathy_columns[0]] == 0) & \n",
    "                                     (bathy_raw[bathy_columns[1]] == 0)].index)\n",
    "\n",
    "# Process tide data\n",
    "tide_dirname = os.path.join(os.path.realpath(''), tide_folder)\n",
    "tide_files = glob.glob(os.path.join(tide_dirname, \"*.csv\"))\n",
    "\n",
    "if not tide_files:\n",
    "    raise ValueError(\"No CSV files found in the specified tide data folder.\")\n",
    "\n",
    "tide_li = []\n",
    "\n",
    "for filename in tide_files:\n",
    "    tide = pd.read_csv(filename, index_col=None, header=4, sep=';')\n",
    "    tide_li.append(tide)\n",
    "\n",
    "if len(tide_li) > 1:\n",
    "    tide = pd.concat(tide_li, axis=0, ignore_index=True)\n",
    "else:\n",
    "    tide = tide_li[0]\n",
    "\n",
    "tide[tide_data_col] = tide[tide_data_col].astype(str)\n",
    "tide[tide_time_col] = tide[tide_time_col].astype(str)\n",
    "tide['TideGauge_Time_UTC'] = tide[tide_data_col].str.cat(tide[tide_time_col], sep=' ')\n",
    "tide['TideGauge_Time_UTC'] = pd.to_datetime(tide['TideGauge_Time_UTC'])\n",
    "\n",
    "start = str(min(bathy_raw['Time (UTC)']))\n",
    "end = str(max(bathy_raw['Time (UTC)']))\n",
    "\n",
    "tide = tide[tide['TideGauge_Time_UTC'].between(start, end)]\n",
    "tide[tide_value_col] = tide[tide_value_col].str.replace(',', '.').astype(float)\n",
    "\n",
    "# Calculate tide-corrected depth\n",
    "bathy_raw['Tide-corrected depth (m)'] = (bathy_raw[bathy_columns[2]] - tide[tide_value_col].mean()) * -1\n",
    "\n",
    "# Reset index\n",
    "bathy_raw.reset_index(drop=True, inplace=True)\n",
    "\n",
    "# Construct the full file path for saving\n",
    "filename = os.path.join(os.path.realpath(''), in_folder, 'tide_corrected_bathy.csv')\n",
    "\n",
    "# Save the processed dataframe to a CSV file\n",
    "bathy_raw.to_csv(filename, index=False)\n",
    "\n",
    "tide_correction=round(tide[tide_value_col].mean(),3)\n",
    "\n",
    "# Define the file path\n",
    "file_path = os.path.join(in_folder, 'Z_correction.txt')\n",
    "\n",
    "# Write the text to the file\n",
    "with open(file_path, 'w') as file:\n",
    "    file.write(f'The Z value to insert into camera coordinates in Agisoft Metashape is {tide_correction}')\n",
    "\n",
    "tide\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Syncing GPS time and camera time\n",
    "The code reads the image that is provided for camera and GPS syncing, and displays the image. The user is requested to insert the date and time as shown in the photo. Then, the code uses the information to build a \"camera_coordinates.txt\" file that associates via timestamps a position from the GPX file to each camera. This file can be used in the \"Reference\" panel of Agisoft Metashape to assign camera positions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to read EXIF data and extract the timestamp\n",
    "def get_exif_timestamp(photo_sync):\n",
    "    with open(photo_sync, 'rb') as f:\n",
    "        tags = exifread.process_file(f)\n",
    "        timestamp = tags.get('EXIF DateTimeOriginal', None)\n",
    "        if timestamp:\n",
    "            return str(timestamp)\n",
    "        else:\n",
    "            return \"No EXIF timestamp found\"\n",
    "\n",
    "# Display the photo\n",
    "img = PILImage.open(photo_sync)\n",
    "display(img)\n",
    "\n",
    "# Extract and display the timestamp from EXIF data\n",
    "exif_timestamp = get_exif_timestamp(photo_sync)\n",
    "print(f\"Original EXIF Timestamp: {exif_timestamp}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Insert the time you see in the photo\n",
    "year=2020\n",
    "month=8\n",
    "day=7\n",
    "hour=6\n",
    "minutes=38\n",
    "seconds=0\n",
    "\n",
    "# Convert EXIF timestamp to datetime object\n",
    "if exif_timestamp != \"No EXIF timestamp found\":\n",
    "    exif_timestamp_dt = datetime.strptime(exif_timestamp, \"%Y:%m:%d %H:%M:%S\")\n",
    "    # Create datetime object for user-provided timestamp\n",
    "    user_timestamp_dt = datetime(year, month, day, hour, minutes, seconds)\n",
    "    \n",
    "    # Calculate the time difference\n",
    "    time_difference = exif_timestamp_dt - user_timestamp_dt\n",
    "    \n",
    "    # Display the time difference\n",
    "    print(f\"Time difference: {time_difference}\")\n",
    "else:\n",
    "    print(\"No EXIF timestamp found in the photo.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to extract EXIF data from images\n",
    "def get_exif_data(image_path):\n",
    "    try:\n",
    "        image = Image.open(image_path)\n",
    "        exif_data = {}\n",
    "        info = image._getexif()\n",
    "        if info:\n",
    "            for tag, value in info.items():\n",
    "                decoded = TAGS.get(tag, tag)\n",
    "                exif_data[decoded] = value\n",
    "        return exif_data\n",
    "    except UnidentifiedImageError:\n",
    "        return None\n",
    "\n",
    "# Function to extract datetime from EXIF data\n",
    "def get_datetime(exif_data):\n",
    "    datetime_str = exif_data.get('DateTimeOriginal') or exif_data.get('DateTime')\n",
    "    if datetime_str:\n",
    "        return datetime.strptime(datetime_str, '%Y:%m:%d %H:%M:%S')\n",
    "    return None\n",
    "\n",
    "\n",
    "photo_timestamps = []\n",
    "\n",
    "# Iterate over photos and extract timestamps\n",
    "for photo in os.listdir(photo_dir):\n",
    "    photo_path = os.path.join(photo_dir, photo)\n",
    "    if not photo.lower().endswith(('.png', '.jpg', '.jpeg', '.tiff', '.bmp', '.gif')):\n",
    "        continue\n",
    "    exif_data = get_exif_data(photo_path)\n",
    "    if exif_data:\n",
    "        photo_datetime = get_datetime(exif_data)\n",
    "        if photo_datetime:\n",
    "            photo_timestamps.append((photo, photo_datetime))\n",
    "\n",
    "# Convert to DataFrame\n",
    "photo_df = pd.DataFrame(photo_timestamps, columns=['Label', 'Timestamp'])\n",
    "\n",
    "# Adjust photo timestamps\n",
    "photo_df['Adjusted Timestamp'] = photo_df['Timestamp'] - time_difference\n",
    "\n",
    "# Function to parse GPX file\n",
    "def parse_gpx(gpx_file):\n",
    "    tree = ET.parse(gpx_file)\n",
    "    root = tree.getroot()\n",
    "    ns = {'default': 'http://www.topografix.com/GPX/1/1'}\n",
    "    gpx_data = []\n",
    "    for trkpt in root.findall('.//default:trkpt', ns):\n",
    "        lat = float(trkpt.get('lat'))\n",
    "        lon = float(trkpt.get('lon'))\n",
    "        time_str = trkpt.find('default:time', ns).text\n",
    "        try:\n",
    "            time = datetime.strptime(time_str, '%Y-%m-%dT%H:%M:%S.%fZ')\n",
    "        except ValueError:\n",
    "            time = datetime.strptime(time_str, '%Y-%m-%dT%H:%M:%SZ')\n",
    "        gpx_data.append((lat, lon, time))\n",
    "    return pd.DataFrame(gpx_data, columns=['Latitude', 'Longitude', 'Timestamp'])\n",
    "\n",
    "gpx_df = parse_gpx(gpx_file)\n",
    "\n",
    "# Function to find nearest GNSS point\n",
    "def find_nearest_gnss(photo_time):\n",
    "    time_diffs = gpx_df['Timestamp'] - photo_time\n",
    "    min_diff = min(time_diffs, key=abs)\n",
    "    nearest_gnss = gpx_df.loc[time_diffs.abs().idxmin()]\n",
    "    return nearest_gnss, min_diff\n",
    "\n",
    "# Add GNSS data to photos\n",
    "output_data = []\n",
    "for index, row in photo_df.iterrows():\n",
    "    nearest_gnss, time_diff = find_nearest_gnss(row['Adjusted Timestamp'])\n",
    "    output_data.append([\n",
    "        row['Label'],\n",
    "        nearest_gnss['Latitude'],\n",
    "        nearest_gnss['Longitude'],\n",
    "        tide_correction,\n",
    "        1,\n",
    "        0.2,\n",
    "        abs(time_diff.total_seconds())\n",
    "    ])\n",
    "\n",
    "# Convert to DataFrame\n",
    "output_df = pd.DataFrame(output_data, columns=['Label', 'Latitude', 'Longitude', 'Altitude', 'Hrz Accuracy','Vrt Accuracy', 'Time Difference'])\n",
    "\n",
    "\n",
    "# Output file path\n",
    "output_file = os.path.join(in_folder, 'camera_coordinates.txt')\n",
    "output_df.to_csv(output_file, index=False, sep=',')\n",
    "\n",
    "# Print summary\n",
    "print(f\"Total images processed: {len(photo_df)}\")\n",
    "print(f\"Total images geotagged: {len(output_df[output_df['Time Difference'] < 1])}\")\n",
    "print(f\"Total images not geotagged: {len(output_df[output_df['Time Difference'] >= 1])}\")\n",
    "print(f\"Time and date of first image: {photo_df['Timestamp'].min()}\")\n",
    "print(f\"Time and date of last image: {photo_df['Timestamp'].max()}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Display data\n",
    "The data is then displayed in a map to ensure that the processing was successful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert DataFrame to GeoDataFrame\n",
    "geometry = [Point(xy) for xy in zip(output_df['Longitude'], output_df['Latitude'])]\n",
    "geo_df = gpd.GeoDataFrame(output_df, geometry=geometry)\n",
    "\n",
    "# Convert bathy_raw to a GeoDataFrame\n",
    "geometry = [Point(xy) for xy in zip(bathy_raw['Longitude (dec degrees)'], bathy_raw['Latitude (dec degrees)'])]\n",
    "bathy_geo_df = gpd.GeoDataFrame(bathy_raw, geometry=geometry)\n",
    "\n",
    "# Plotting the bathymetric points, the output_df points, and the minimum bounding rectangle\n",
    "fig, ax = plt.subplots(figsize=(10, 10))\n",
    "\n",
    "# Plot the bathymetric points\n",
    "bathy_geo_df.plot(ax=ax, marker='o', color='gray', markersize=20, label='Echosounder Points')\n",
    "\n",
    "# Plot the points from output_df\n",
    "geo_df.plot(ax=ax, marker='o', color='blue', markersize=5, label='Camera locations')\n",
    "\n",
    "title = f\"Survey time: {day} {month} {year}, {hour:02d}:{minutes:02d}\"\n",
    "plt.title(title)\n",
    "\n",
    "# Adding title and legend\n",
    "plt.legend()\n",
    "plt.xlabel('Longitude')\n",
    "plt.ylabel('Latitude')\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  }
 ],
 "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": 4
}