lidar-processing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2e1b91fc-82c5-488a-bd91-b1105ddf6b90",
   "metadata": {},
   "source": [
    "# Access to Swisstopo Lidar Data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "195b43d2-b46c-4d6f-a6e9-e6fd151e8160",
   "metadata": {},
   "source": [
    "This notebook shows and explain how to get the low resolution Lidar data from Swisstopo and how to process them.\n",
    "\n",
    "Most example are taken or inspired from [laspy documentation](https://laspy.readthedocs.io/en/latest/index.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# These library are needed for the notebook\n",
    "!pip install numpy          # must have\n",
    "!pip install pandas         # must have\n",
    "!pip install laspy          # must have\n",
    "\n",
    "!pip install matplotlib     # nice to have"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Get the data"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Go to [this](https://www.swisstopo.admin.ch/fr/geodata/height/surface3d.html#technische_details) website, to the 'swissSURFACE3D - Accès aux géodonnées' part and select the area of data that you want. \n",
    "\n",
    "For this example we select by commune and selct the commune of Fribourg (FR). As option, we keep the format LIDAR (ZIP), coordonates system MN95 (default) and the actual state. \n",
    "4 links are made available to download. For the sake of the example, we only download the first one here: 'swisssurface3d_2019_2576-1183_2056_5728.las.zip'. Click on 'Download'.\n",
    "\n",
    "You now have a file 'swisssurface3d_2019_2576-1183_2056_5728.las.zip' locally of around 130MB (or any other las.zip file). Make it accessible from your notebook."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Process the data"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# TODO: adapt this path to access the downloaded data from your notebook\n",
    "SOURCE_LIDAR_ZIP = 'swisssurface3d_2019_2576-1183_2056_5728.las.zip'\n",
    "TARGET_LIDAR_FOLDER = 'fribourg_part_1' "
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Unzip the file and get its path"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "import zipfile\n",
    "\n",
    "with zipfile.ZipFile(SOURCE_LIDAR_ZIP, 'r') as zip_ref:\n",
    "    zip_ref.extractall(TARGET_LIDAR_FOLDER)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "import os\n",
    "las_path = os.path.join(TARGET_LIDAR_FOLDER, os.listdir(TARGET_LIDAR_FOLDER)[0])\n",
    "las_path"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Open Lidar data\n",
    "The file is stored in [las format](https://laspy.readthedocs.io/en/latest/intro.html) which is a public standard format for exchanging 3D point data. There are different python librairies to process this type of data. Here we will use [laspy](https://pypi.org/project/laspy/)."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "import laspy as lp\n",
    "import numpy as np"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Get metadata information"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "#  Read and show metadate about point cloud\n",
    "with lp.open(las_path) as f:\n",
    "    # Get metadata information\n",
    "    print(f\"Point format:             {f.header.point_format}\")\n",
    "    print(f\"Number of points:         {f.header.point_count}\")\n",
    "    print(f\"Variable length record:   {f.header.vlrs}\\n\")\n",
    "    \n",
    "    # Read the point cloud data\n",
    "    las = f.read()\n",
    "    print(f\"Las file:                 {las}\")\n",
    "    point_format = las.point_format\n",
    "    print(f\"Dimensions:               {list(point_format.dimension_names)}\\n\")\n",
    "    \n",
    "    # Offset and scaling for projections (discussed later)\n",
    "    print(f\"Points offset:            {f.header.offset}\")\n",
    "    print(f\"Points scaling:           {f.header.scale}\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Process Lidar Data"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "import pandas as pd"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "las = lp.read(las_path)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Get pandas dataset of the point cloud coordinates and class\n",
    "df = pd.DataFrame({'X': las.X, 'Y': las.Y, 'Z':las.Z, 'Class': np.array(las.classification)})\n",
    "df.head()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "From there, you can do anything that you can usually do in pandas. You can modify the dataset to add other columns from the point cloud."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Plot Lidar data"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Subsample (to have less points)\n",
    "print(f'Previous size: {df.shape}')\n",
    "df = df.sample(frac=0.50, replace=False, random_state=1)\n",
    "\n",
    "# Keep a smaller region\n",
    "df = df[df['X'] > np.mean(df['X'])]\n",
    "df = df[df['Y'] > np.mean(df['Y'])]\n",
    "print(f'New size: {df.shape}')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(15, 15))\n",
    "ax = plt.axes(projection='3d')\n",
    "ax.scatter3D(df['X'], df['Y'], df['Z'], c =  df['Class'], cmap='viridis', s=1);"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Reproject to another Reference Coordinate System"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Unfortunatelly with these data, the variable length record which should contain the reference coordinate system (see above in metadata part) is empty. \n",
    "\n",
    "Based on the additonal documentation of swissSURFACE3D, we know that every downloaded tiles names contain an id which corresponds to the kilometric coordinates in the reference system [MN95] (https://www.swisstopo.admin.ch/fr/connaissances-faits/mensuration-geodesie/cadres-de-reference/local/mn95.html#234_1607609199106) in the south-west angle.\n",
    "\n",
    "Hence our `swisssurface3d_2019_2576-1183_2056_5728.las` would be the tile recorded in 2019 at the coordinate (2576, 1183) as can be checked below:"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "#  Read and show metadate about point cloud\n",
    "with lp.open(las_path) as f:\n",
    "    # Offset and scaling for projections ()\n",
    "    print(f\"Points offset:            {f.header.offset}\")\n",
    "    print(f\"Points scaling:           {f.header.scale}\")\n",
    "    x_offset, y_offset, z_offset = f.header.offset\n",
    "    x_scale, y_scale, z_scale = f.header.scale"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "To reconstruct the original point cloud in MN95 coordinates in the system CH1903+, the coordinates are calculated using the following formula(s): \\\n",
    "x = (x_int * x_scale) + x_offset \\\n",
    "y = (y_int * y_scale) + y_offset \\\n",
    "z = (z_int * z_scale) + z_offset"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "df['X'] = df['X'] * x_scale + x_offset\n",
    "df['Y'] = df['Y'] * y_scale + y_offset\n",
    "df['Z'] = df['Z'] * z_scale + z_offset"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "df.rename(columns = {'X': 'Easting', 'Y': 'Northing', 'Z': 'Elevation'}, inplace = True)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "df.head() # in CH1903+ coordinates"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "To reproject in any other coordinate reference system, more complicated formulas are required. They can be found online and applied to the points."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d2e542e-47fb-4f90-b3a4-1e57f0627730",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}