notebooks/1. AHN preprocessing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c3fd9730-0dc5-4520-a485-b2445c984820",
   "metadata": {},
   "source": [
    "# Pre-processing of LAS tiles\n",
    "\n",
    "Pre-processing consists of two parts:\n",
    "1. Clipping of AHN point clouds to match target point clouds;\n",
    "2. Ground and building surface extraction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57f542f8-21bf-4382-bfb9-11b11b05802b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Uncomment to load the local package rather than the pip-installed version.\n",
    "# Add project src to path.\n",
    "import set_path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8adb4173-257d-49f2-9b20-c5af54f642b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import modules.\n",
    "import numpy as np\n",
    "import time\n",
    "\n",
    "import upcp.preprocessing.ahn_preprocessing as ahn_preprocessing\n",
    "import upcp.utils.ahn_utils as ahn_utils\n",
    "import upcp.utils.las_utils as las_utils\n",
    "import upcp.utils.plot_utils as plot_utils"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f62d5f5-6bd7-4ec0-9232-2b2f79fe907d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set-up.\n",
    "\n",
    "# The point cloud tile for which a corresponding AHN tile should be created.\n",
    "tile_code = '2397_9705'\n",
    "filename = '../datasets/pointcloud/filtered_' + tile_code + '.laz'\n",
    "\n",
    "# The folder where AHN data will be saved.\n",
    "ahn_data_folder = '../datasets/ahn/'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de6d646d-51ab-49ac-9555-9a74138ce08d",
   "metadata": {},
   "source": [
    "## 1a. Clip a single file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4deb32af-cb1d-4856-857a-d16ef6d2215e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download the corresponding AHN cloud from \n",
    "# https://www.arcgis.com/apps/Embed/index.html?appid=a3dfa5a818174aa787392e461c80f781\n",
    "!wget https://download.pdok.nl/rws/ahn3/v1_0/laz/C_25DN2.LAZ -P ../datasets/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05d4ee80-dd40-4796-a650-9745959412da",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the AHN point cloud.\n",
    "# WARNING: this may require 16+ GB of RAM, depending on the size of the cloud.\n",
    "ahn_cloud_file = '../datasets/C_25DN2.LAZ'\n",
    "ahn_cloud = las_utils.read_las(ahn_cloud_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a115c01-bd26-4648-932d-8ca4f5289aa5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Clip the corresponding AHN point cloud tile.\n",
    "start = time.time()\n",
    "ahn_preprocessing.clip_ahn_las_tile(ahn_cloud, filename, out_folder=ahn_data_folder, buffer=1)\n",
    "end = time.time()\n",
    "print(f'Tile clipped in {end-start:.2f} seconds.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56a36c6c-fd39-478c-8b08-dfcb7622aad9",
   "metadata": {},
   "source": [
    "## 1b. Clip all files in a folder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90e42e65-7096-439f-94c6-2dbad3fa6743",
   "metadata": {},
   "outputs": [],
   "source": [
    "las_folder = '../datasets/pointcloud/'\n",
    "\n",
    "# Set resume=True to skip existing files in the out_folder.\n",
    "ahn_preprocessing.clip_ahn_las_folder(ahn_cloud, las_folder, out_folder=ahn_data_folder, buffer=1, resume=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "369980ec-43b0-4817-b80f-f005c20ba978",
   "metadata": {},
   "source": [
    "## 2a. Surface extraction for a single file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63abda59-af1e-43ce-9414-f1df601ff089",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the clipped tile that we have just created.\n",
    "ahn_las_file = ahn_data_folder + '/ahn_' + tile_code + '.laz'\n",
    "\n",
    "start = time.time()\n",
    "# This will create the AHN ground and building surfaces, with a 0.1m resolution, and save them as an .npz file.\n",
    "_file = ahn_preprocessing.process_ahn_las_tile(ahn_las_file, out_folder=ahn_data_folder)\n",
    "end = time.time()\n",
    "print(f'Processed LAS tile in {end-start:.2f}s.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19462d26-e59c-40b8-ab1f-2f93d5f35390",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the result.\n",
    "ahn_reader = ahn_utils.NPZReader(ahn_data_folder)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86a5b5c4-4412-4e6c-ba9d-800610d14043",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Optional: plot the resulting surfaces.\n",
    "%matplotlib widget\n",
    "plot_utils.plot_ahn_sidebyside(tile_code, ahn_reader)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75c2ae56-5929-49bd-a394-e73392734e36",
   "metadata": {},
   "source": [
    "This will look like:\n",
    "\n",
    "![Surface Plot](https://github.com/Amsterdam-AI-Team/Urban_PointCloud_Processing/raw/main/media/examples/surface_plot.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce1793e2-878e-4224-becb-37476e3e25a1",
   "metadata": {},
   "source": [
    "## 2b. Surface extraction for a folder of ahn cloud tiles"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "889ab9b4-a304-4afa-a715-7b783aa9ad2b",
   "metadata": {},
   "source": [
    "We have a batch script to do this, it allows using multiple CPU threads to speed things up.\n",
    "\n",
    "```sh\n",
    "$ cd ../scripts\n",
    "$ python ahn_batch_processor.py --in_folder ../datasets/ahn/ --out_folder ../datasets/ahn/ --resume --workers 4\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "927c31d9-08a6-4b7a-8bb3-aefc7ccfd2fc",
   "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.10.8"
  },
  "vscode": {
   "interpreter": {
    "hash": "e6888fd1df2ad0cf70295ff83cef12f2beafc15e4cff2fac0fbf5a36213f1b1c"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}