03_sightlines_dtm.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "99bcb5f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import rioxarray\n",
    "import math\n",
    "import shapely\n",
    "import xvec\n",
    "\n",
    "from shapely.geometry import Point, MultiPoint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "0a21db39",
   "metadata": {},
   "outputs": [],
   "source": [
    "sightlines = pd.read_pickle(\"my_results/method_1/d06_sightlines_geometries.pickle\")\n",
    "sightlines = sightlines[[\"start_point\", \"end_point\", \"sight_line_points\"]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0523e664",
   "metadata": {},
   "source": [
    "# Convert tuples with coordinates to shapely geometry\n",
    "Should eventually be deprecated by already using shapely geometries\n",
    "in the first place."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "fa705303",
   "metadata": {},
   "outputs": [],
   "source": [
    "def convert_to_point(row):\n",
    "    return Point(row)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d467e829",
   "metadata": {},
   "outputs": [],
   "source": [
    "def convert_to_multipoint(row):\n",
    "    if row:\n",
    "        return MultiPoint(row)\n",
    "    else:\n",
    "        return None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e28ee0a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert start and end point to geometry\n",
    "tuple_columns = [\"start_point\", \"end_point\"]\n",
    "\n",
    "for col in tuple_columns:\n",
    "    sightlines[col+\"_geom\"] = sightlines[col].apply(convert_to_point)\n",
    "\n",
    "\n",
    "# Convert sight line points to geometry\n",
    "sightlines[\"sl_points\"] = sightlines[\"sight_line_points\"].apply(convert_to_multipoint)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc0ba283",
   "metadata": {},
   "source": [
    "# Load raster dtm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f6e743be",
   "metadata": {},
   "outputs": [],
   "source": [
    "raster = rioxarray.open_rasterio(\"data/dtm/dtm_raw/RGEALTI_FXX_1025_6285_MNT_LAMB93_IGN69.asc\")\n",
    "raster = raster.drop_vars(\"band\").squeeze()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51d28053",
   "metadata": {},
   "source": [
    "# Extract z coordinates from raster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c203ee34",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract z coords from raster\n",
    "z_start = raster.drop_vars(\"spatial_ref\").xvec.extract_points(points=sightlines[\"start_point_geom\"], x_coords=\"x\",y_coords=\"y\")\n",
    "\n",
    "# Convert do gpd\n",
    "z_start = z_start.xvec.to_geopandas().rename(columns={0: \"z\"})\n",
    "\n",
    "# Append z values to points\n",
    "z_start[\"start_point_3d\"] = z_start.apply(lambda row: Point(row[\"geometry\"].x, row[\"geometry\"].y, row[\"z\"]), axis=1)\n",
    "z_start = z_start.drop(columns=[\"geometry\", \"z\"])\n",
    "\n",
    "\n",
    "\n",
    "# Extract z coords from raster\n",
    "z_end = raster.drop_vars(\"spatial_ref\").xvec.extract_points(points=sightlines[\"end_point_geom\"], x_coords=\"x\",y_coords=\"y\")\n",
    "\n",
    "# Convert do gpd\n",
    "z_end = z_end.xvec.to_geopandas().rename(columns={0: \"z\"})\n",
    "\n",
    "# Append z values to points\n",
    "z_end[\"end_point_3d\"] = z_end.apply(lambda row: Point(row[\"geometry\"].x, row[\"geometry\"].y, row[\"z\"]), axis=1)\n",
    "z_end = z_end.drop(columns=[\"geometry\", \"z\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "8f119e3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "z_points_list = []\n",
    "\n",
    "for row in sightlines[\"sl_points\"]:\n",
    "    if row is not None:\n",
    "        points = row.geoms\n",
    "\n",
    "        z_points = raster.drop_vars(\"spatial_ref\").xvec.extract_points(points=points, x_coords=\"x\",y_coords=\"y\")\n",
    "        z_points = z_points.xvec.to_geopandas().rename(columns={0: \"z\"})\n",
    "\n",
    "        z_points[\"geometry\"] = z_points.apply(lambda row: Point(row[\"geometry\"].x, row[\"geometry\"].y, row[\"z\"]), axis=1)\n",
    "        z_points = z_points.drop(columns=\"z\")\n",
    "        \n",
    "        multipoint = MultiPoint(z_points[\"geometry\"].tolist())\n",
    "\n",
    "    else:\n",
    "        multipoint = None\n",
    "\n",
    "    z_points_list.append(multipoint)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f3a110e",
   "metadata": {},
   "source": [
    "## Put results together"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "2f8f8ff1",
   "metadata": {},
   "outputs": [],
   "source": [
    "sightlines = pd.concat([z_start, z_end], axis=1)\n",
    "\n",
    "sightlines = sightlines.rename(columns={\"start_point_3d\": \"sl_start\",\n",
    "                                          \"end_point_3d\": \"sl_end\"})\n",
    "\n",
    "sightlines[\"sl_points\"] = z_points_list"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf47ff54",
   "metadata": {},
   "source": [
    "# Compute slope function\n",
    "Finally, yaaay."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "de597a50",
   "metadata": {},
   "outputs": [],
   "source": [
    "NODATA_RASTER = raster.rio.nodata\n",
    "\n",
    "def compute_slope(road_row):\n",
    "    start = road_row.sl_start   # Point z\n",
    "    end = road_row.sl_end       # Point z\n",
    "    slp = road_row.sl_points    # Multipoint z\n",
    "\n",
    "    if slp == None:\n",
    "        # Case when there is no sight line point (e.g. when the road is too short)\n",
    "        # just computes slope between start and end\n",
    "        if start.z == NODATA_RASTER or end.z == NODATA_RASTER:\n",
    "            # Case when there is at least one invalid z coord\n",
    "            return 0, 0, 0, False\n",
    "        slope_percent = abs(start.z - end.z) / shapely.distance(start, end)\n",
    "        slope_degree = math.degrees(math.atan(slope_percent))\n",
    "\n",
    "        return slope_percent, slope_degree, 1, True\n",
    "    \n",
    "    # From Multipoint z to Point z list\n",
    "    slp_list = [p for p in slp.geoms]\n",
    "\n",
    "    points = []\n",
    "\n",
    "    points.append(start)\n",
    "    # From Point z list to all points list\n",
    "    for p in slp_list:\n",
    "        points.append(p)\n",
    "    points.append(end)\n",
    "\n",
    "    # number of points\n",
    "    nb_points = len([start]) + len([end]) + len(slp_list)\n",
    "\n",
    "    # temporary variables to store inter slope values\n",
    "    sum_slope_percent = 0\n",
    "    sum_slope_radian = 0\n",
    "    sum_nb_points = 0\n",
    "\n",
    "    # if there is one or more sight line points\n",
    "    for i in range (1, nb_points-1):\n",
    "        a = points[i-1]\n",
    "        b = points[i+1]\n",
    "\n",
    "        if a.z == NODATA_RASTER or b.z == NODATA_RASTER:\n",
    "            # Case when there is no valid z coord in slpoint\n",
    "            continue\n",
    "            \n",
    "        sum_nb_points += 1\n",
    "        inter_slope_percent = abs(a.z - b.z) / shapely.distance(a, b)\n",
    "\n",
    "        sum_slope_percent += inter_slope_percent\n",
    "        sum_slope_radian += math.atan(inter_slope_percent)\n",
    "\n",
    "    if sum_nb_points == 0:\n",
    "        # Case when no slpoint has a valid z coord\n",
    "        # Unable to compute slope\n",
    "        return 0, 0, 0, False\n",
    "        \n",
    "    # compute mean of inter slopes\n",
    "    slope_percent = sum_slope_percent/sum_nb_points\n",
    "    slope_degree = math.degrees(sum_slope_radian/sum_nb_points)\n",
    "\n",
    "    return slope_degree, slope_percent, sum_nb_points, True"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "246d682e",
   "metadata": {},
   "source": [
    "# Compute results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "563729ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "slope_values = []\n",
    "\n",
    "for _, road_row in sightlines.iterrows():\n",
    "\n",
    "    # display([p for p in road_row.sl_points.geoms])\n",
    "    # break\n",
    "\n",
    "    slope_degree, slope_percent, n_slopes, slope_valid = compute_slope(road_row)\n",
    "\n",
    "    slope_values.append([slope_degree, slope_percent, n_slopes, slope_valid])\n",
    "\n",
    "    df_slopes = pd.DataFrame(slope_values,\n",
    "                             columns=[\"slope_degree\",\n",
    "                                      \"slope_percent\",\n",
    "                                      \"n_slopes\",\n",
    "                                      \"slope_valid\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "3fb51854",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>slope_degree</th>\n",
       "      <th>slope_percent</th>\n",
       "      <th>n_slopes</th>\n",
       "      <th>slope_valid</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2.305050</td>\n",
       "      <td>0.040300</td>\n",
       "      <td>28</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.094379</td>\n",
       "      <td>0.019111</td>\n",
       "      <td>33</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3.159987</td>\n",
       "      <td>0.055233</td>\n",
       "      <td>15</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2.404723</td>\n",
       "      <td>0.042030</td>\n",
       "      <td>34</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.878624</td>\n",
       "      <td>0.015341</td>\n",
       "      <td>13</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2162</th>\n",
       "      <td>0.284191</td>\n",
       "      <td>0.004960</td>\n",
       "      <td>7</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2163</th>\n",
       "      <td>4.664951</td>\n",
       "      <td>0.081627</td>\n",
       "      <td>5</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2164</th>\n",
       "      <td>9.312918</td>\n",
       "      <td>0.166963</td>\n",
       "      <td>9</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2165</th>\n",
       "      <td>6.194412</td>\n",
       "      <td>0.109414</td>\n",
       "      <td>31</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2166</th>\n",
       "      <td>7.294105</td>\n",
       "      <td>0.132081</td>\n",
       "      <td>26</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>2167 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "      slope_degree  slope_percent  n_slopes  slope_valid\n",
       "0         2.305050       0.040300        28         True\n",
       "1         1.094379       0.019111        33         True\n",
       "2         3.159987       0.055233        15         True\n",
       "3         2.404723       0.042030        34         True\n",
       "4         0.878624       0.015341        13         True\n",
       "...            ...            ...       ...          ...\n",
       "2162      0.284191       0.004960         7         True\n",
       "2163      4.664951       0.081627         5         True\n",
       "2164      9.312918       0.166963         9         True\n",
       "2165      6.194412       0.109414        31         True\n",
       "2166      7.294105       0.132081        26         True\n",
       "\n",
       "[2167 rows x 4 columns]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_slopes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "5077645c",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_slopes.to_parquet(\"my_results/method_3/d06_dtm_result.parquet\")"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "encoding": "# -*- coding: utf-8 -*-",
   "formats": "ipynb,py:light"
  },
  "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}