{ "cells": [ { "cell_type": "markdown", "id": "385d21fd-6f63-41d8-8633-79721f90fe8d", "metadata": {}, "source": [ "# LiDAR point cloud of the Trail Islands in British Columbia, Canada\n", "\n", "This is a LiDAR point cloud from the [LidarBC](https://www2.gov.bc.ca/gov/content/data/geographic-data-services/lidarbc) project, sliced to the small [Trail Islands](https://apps.gov.bc.ca/pub/bcgnws/names/21973.html) North of Vancouver to reduce the data size. The islands have some nice looking topography and their isolated nature creates problems for some interpolation methods.\n", "\n", "Coordinates of the original data are in NAD83 UTM zone 10 and the vertical datum for elevation is CGVD2013 (based on the Canadian geoid).\n", "\n", "The horizontal datum was changed to WGS84 and coordinates un-projected to longitude, latitude.\n", "\n", "License: [Open Government Licence - British Columbia](https://www2.gov.bc.ca/gov/content/data/open-data/open-government-licence-bc)" ] }, { "cell_type": "code", "execution_count": 1, "id": "34fa5814-2463-44eb-97a1-9aa57670e3df", "metadata": {}, "outputs": [], "source": [ "import os\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pandas as pd\n", "import verde as vd\n", "import pooch\n", "import pyproj\n", "import pygmt\n", "import laspy " ] }, { "cell_type": "markdown", "id": "15c63575-49cf-4503-86fe-abd6f455f552", "metadata": {}, "source": [ "## Download the tile we need\n", "\n", "The islands are covered by the tile `bc_092g041` of the dataset. Download it with Pooch and cache it." ] }, { "cell_type": "code", "execution_count": 2, "id": "50ee303e-927d-4ba8-8d3e-e3c8849505ec", "metadata": {}, "outputs": [], "source": [ "fname = pooch.retrieve(\n", " url=\"https://nrs.objectstore.gov.bc.ca/gdwuts/092/092g/2019/pointcloud/bc_092g041_4_2_2_xyes_8_utm10_2019.laz\",\n", " known_hash=\"md5:5acad29ff5ba1f2a09453ed16c3004da\",\n", ")" ] }, { "cell_type": "code", "execution_count": 3, "id": "fbc2e32b-9b7c-4404-925c-9e00a7c12ff6", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "size: 83.871764 Mb\n" ] } ], "source": [ "print(f\"size: {os.path.getsize(fname) / 1e6} Mb\")" ] }, { "cell_type": "markdown", "id": "86ca37f7-c83c-4741-941b-27c29f6d70f8", "metadata": {}, "source": [ "## Read the data\n", "\n", "Use laspy to read in the data from the file. We'll read it in chunks and only take the ground reflectance points to avoid overloading our memory." ] }, { "cell_type": "code", "execution_count": 4, "id": "6bbd81fa-0f94-46e0-b458-a324295b6e2b", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", " | easting | \n", "northing | \n", "height | \n", "
---|---|---|---|
0 | \n", "441500.56 | \n", "5479185.08 | \n", "-1.14 | \n", "
1 | \n", "441500.54 | \n", "5479184.50 | \n", "-1.15 | \n", "
2 | \n", "441500.43 | \n", "5479183.88 | \n", "-1.10 | \n", "
3 | \n", "441500.84 | \n", "5479184.13 | \n", "-1.14 | \n", "
4 | \n", "441501.26 | \n", "5479184.38 | \n", "-1.13 | \n", "
... | \n", "... | \n", "... | \n", "... | \n", "
2489382 | \n", "442012.88 | \n", "5477949.65 | \n", "-0.37 | \n", "
2489383 | \n", "442012.58 | \n", "5477949.15 | \n", "-0.47 | \n", "
2489384 | \n", "442014.16 | \n", "5477950.85 | \n", "-0.06 | \n", "
2489385 | \n", "442013.24 | \n", "5477949.33 | \n", "-0.40 | \n", "
2489386 | \n", "442015.80 | \n", "5477949.79 | \n", "-0.69 | \n", "
2489387 rows × 3 columns
\n", "\n", " | easting | \n", "northing | \n", "height | \n", "
---|---|---|---|
0 | \n", "441030.75 | \n", "5478941.38 | \n", "-0.67 | \n", "
1 | \n", "441032.74 | \n", "5478942.59 | \n", "-0.96 | \n", "
2 | \n", "441029.94 | \n", "5478940.32 | \n", "-0.78 | \n", "
3 | \n", "441031.37 | \n", "5478941.18 | \n", "-0.61 | \n", "
4 | \n", "441031.89 | \n", "5478941.49 | \n", "-0.62 | \n", "
... | \n", "... | \n", "... | \n", "... | \n", "
829728 | \n", "441484.31 | \n", "5478352.06 | \n", "-1.24 | \n", "
829729 | \n", "441484.07 | \n", "5478350.77 | \n", "-1.25 | \n", "
829730 | \n", "441483.66 | \n", "5478349.07 | \n", "-1.19 | \n", "
829731 | \n", "441482.76 | \n", "5478347.58 | \n", "-1.21 | \n", "
829732 | \n", "441482.49 | \n", "5478347.16 | \n", "-1.20 | \n", "
829733 rows × 3 columns
\n", "\n", " | longitude | \n", "latitude | \n", "elevation_m | \n", "
---|---|---|---|
0 | \n", "-123.8137526 | \n", "49.4602634 | \n", "-0.67 | \n", "
1 | \n", "-123.8137253 | \n", "49.4602745 | \n", "-0.96 | \n", "
2 | \n", "-123.8137636 | \n", "49.4602538 | \n", "-0.78 | \n", "
3 | \n", "-123.8137440 | \n", "49.4602617 | \n", "-0.61 | \n", "
4 | \n", "-123.8137369 | \n", "49.4602645 | \n", "-0.62 | \n", "
... | \n", "... | \n", "... | \n", "... | \n", "
829728 | \n", "-123.8074072 | \n", "49.4550069 | \n", "-1.24 | \n", "
829729 | \n", "-123.8074104 | \n", "49.4549953 | \n", "-1.25 | \n", "
829730 | \n", "-123.8074158 | \n", "49.4549800 | \n", "-1.19 | \n", "
829731 | \n", "-123.8074280 | \n", "49.4549665 | \n", "-1.21 | \n", "
829732 | \n", "-123.8074316 | \n", "49.4549627 | \n", "-1.20 | \n", "
829733 rows × 3 columns
\n", "\n", " | longitude | \n", "latitude | \n", "elevation_m | \n", "
---|---|---|---|
0 | \n", "-123.813753 | \n", "49.460263 | \n", "-0.67 | \n", "
1 | \n", "-123.813725 | \n", "49.460274 | \n", "-0.96 | \n", "
2 | \n", "-123.813764 | \n", "49.460254 | \n", "-0.78 | \n", "
3 | \n", "-123.813744 | \n", "49.460262 | \n", "-0.61 | \n", "
4 | \n", "-123.813737 | \n", "49.460265 | \n", "-0.62 | \n", "
... | \n", "... | \n", "... | \n", "... | \n", "
829728 | \n", "-123.807407 | \n", "49.455007 | \n", "-1.24 | \n", "
829729 | \n", "-123.807410 | \n", "49.454995 | \n", "-1.25 | \n", "
829730 | \n", "-123.807416 | \n", "49.454980 | \n", "-1.19 | \n", "
829731 | \n", "-123.807428 | \n", "49.454966 | \n", "-1.21 | \n", "
829732 | \n", "-123.807432 | \n", "49.454963 | \n", "-1.20 | \n", "
829733 rows × 3 columns
\n", "