notebooks/2. Create GeoTiff raster.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "49ddb8fa-69f7-4b82-8cb3-7e26e4ebf626",
   "metadata": {},
   "source": [
    "# Convert annotated point cloud to a GeoTiff raster\n",
    "\n",
    "The input for this notebok is the pre-processed and filtered cloud that results from notebook `1. Tree filter.ipynb`.\n",
    "\n",
    "We create a raster with three bands:\n",
    "1. `n_points`: the number of tree-points per cell\n",
    "2. `hag`: the 80th percentile height above ground (HAG) for that cell\n",
    "3. `nap`: the 80th percentile NAP height for that cell\n",
    "\n",
    "This raster can be used for further analysis in standard GIS tools."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3a935d0-dfe7-40bb-a2ae-abb7440134ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pathlib\n",
    "\n",
    "import gvl.geotiff_utils as gt_utils\n",
    "import gvl.helper_functions as utils\n",
    "import laspy\n",
    "import numpy as np\n",
    "from tqdm.notebook import tqdm\n",
    "from upcp.utils import clip_utils\n",
    "\n",
    "import set_path"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "3aab804a-13d7-4990-891d-3a186fb7a3da",
   "metadata": {},
   "source": [
    "## Settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b20f31c0-ac8a-486d-975e-98f4353482f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "DATA_FOLDER = pathlib.Path(\"../data\")\n",
    "\n",
    "input_dir = DATA_FOLDER / \"ahn4_trees\"\n",
    "output_dir = DATA_FOLDER / \"ahn4_tree_grid\"\n",
    "N = 4  # Number of digits in tilecode format\n",
    "TILE_WIDTH = 50  # Tile width in meters\n",
    "\n",
    "EPSG = 28992  # RD coordinates\n",
    "grid_res = 0.5  # 50cm resolution grid\n",
    "\n",
    "include_overlap = False  # Add the existing overlap buffer around the tile\n",
    "BUFFER = 5\n",
    "\n",
    "# Our classification\n",
    "UNKNOWN = 0\n",
    "TREE = 1\n",
    "NOISE = 2\n",
    "OTHER = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18c23393-55f3-4680-8046-92ba68ee5016",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create output folder\n",
    "output_dir.mkdir(parents=True, exist_ok=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "af74836f-a1d3-43e9-ba29-36ada57273b0",
   "metadata": {},
   "source": [
    "## Main loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b695c0f9-06a5-4bdc-b3ad-4eaa346f3110",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_files = list(input_dir.glob(\"trees*.laz\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4e2a3f9-06e4-4779-9d6f-67fd995f5238",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check existing output files (ignore this cell to re-run for all tiles)\n",
    "existing_files = list(output_dir.glob(\"raster*.tif\"))\n",
    "existing_codes = {\n",
    "    utils.get_tilecode_from_filename(file.name, n_digits=N) for file in existing_files\n",
    "}\n",
    "\n",
    "input_files = [\n",
    "    file\n",
    "    for file in input_files\n",
    "    if utils.get_tilecode_from_filename(file.name, n_digits=N) not in existing_codes\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25ab0627-87cd-4bea-a4b0-7691b299b9d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "pbar = tqdm(input_files, unit=\"file\", smoothing=0)\n",
    "\n",
    "for file in pbar:\n",
    "    tilecode = utils.get_tilecode_from_filename(file.name, n_digits=N)\n",
    "    pbar.set_postfix_str(tilecode)\n",
    "\n",
    "    # Load LAS data\n",
    "    las = laspy.read(file)\n",
    "    mask = las.label == TREE\n",
    "    if np.count_nonzero(mask) == 0:\n",
    "        continue\n",
    "    points_xyz = np.vstack((las.x[mask], las.y[mask], las.z[mask])).T\n",
    "    hag = las.hag[mask]\n",
    "\n",
    "    if include_overlap:\n",
    "        ((x_min, y_max), (x_max, y_min)) = utils.get_bbox_from_tile_code(\n",
    "            tilecode, padding=BUFFER, width=TILE_WIDTH, height=TILE_WIDTH\n",
    "        )\n",
    "        tilewidth = x_max - x_min  # Assume square tile\n",
    "    else:\n",
    "        ((x_min, y_max), (x_max, y_min)) = utils.get_bbox_from_tile_code(\n",
    "            tilecode, padding=0, width=TILE_WIDTH, height=TILE_WIDTH\n",
    "        )\n",
    "        tilewidth = x_max - x_min  # Assume square tile\n",
    "        mask = clip_utils.rectangle_clip(points_xyz, (x_min, y_min, x_max, y_max))\n",
    "        if np.count_nonzero(mask) == 0:\n",
    "            continue\n",
    "        points_xyz = points_xyz[mask, :]\n",
    "        hag = hag[mask]\n",
    "\n",
    "    n_points, hag, nap = gt_utils.get_tile_grid(\n",
    "        points_xyz, hag, x_min, y_min, tilewidth, grid_res\n",
    "    )\n",
    "\n",
    "    output_file = output_dir / f\"raster_{tilecode}.tif\"\n",
    "    gt_utils.save_geotiff(\n",
    "        output_file.as_posix(),\n",
    "        n_points.T,\n",
    "        hag.T,\n",
    "        nap.T,\n",
    "        x_min,\n",
    "        y_min,\n",
    "        grid_res,\n",
    "        epsg=EPSG,\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94f596de",
   "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.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}