{ "cells": [ { "cell_type": "markdown", "metadata": { "id": "SZeETtIsOEYf" }, "source": [ "## Processing COPC Data with PDAL, STAC, and Planetary Computer\n", "\n", "This scenario demonstrates how to use the PDAL Python bindings to produce some simple raster products from a COPC file. \n", "\n", "### COPC\n", "\n", "COPC is [Cloud Optimized Point Cloud](https://copc.io). COPC files are LASzip LAZ 1.4 data stored as a clustered octree. COPC files allow applications to select data for a window or a resolution, and allow them to limit how much data they must fetch, decompress, and process. You can read more about COPC at https://lidarmag.com/2021/12/27/cloud-native-geospatial-lidar-with-the-cloud-optimized-point-cloud/\n", "\n", "### PDAL\n", "\n", "PDAL supports reading and writing COPC with [readers.copc](https://pdal.io/stages/readers.copc.html) and [writers.copc](https://pdal.io/stages/writers.copc.html). \n", "\n", "### PDAL-Python\n", "\n", "The [PDAL Python bindings](https://https://github.com/PDAL/python/) allow programmatic composition of [PDAL pipelines](https://https://pdal.io/pipeline.html). \n", "\n", "### What are we doing?\n", "\n", "* Select some data from Chicago using Planetary Computer's STAC API\n", "* Write a relative intensity image (RII) \n", "* Extract a HeightAboveGround surface and estimate the tallest point from the ground in the scene" ] }, { "cell_type": "markdown", "metadata": { "id": "kY0ikB6JQ2G7", "tags": [] }, "source": [ "### Setup\n", "\n", "We have some imports and helper code to define first.\n", "\n", "**Imports**" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "id": "kydbokRLpzsm" }, "outputs": [], "source": [ "import pdal\n", "import pystac_client\n", "import planetary_computer\n", "import PIL\n", "import pyproj" ] }, { "cell_type": "markdown", "metadata": { "id": "_GOto7pvWhVy", "tags": [] }, "source": [ "**Useful code for estimating the UTM zone of a point**\n", "\n", "Estimate the UTM zone of the point using [PyPROJ's method](https://gis.stackexchange.com/a/423614)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Estimate our UTM zone\n", "from pyproj import CRS\n", "from pyproj.aoi import AreaOfInterest\n", "from pyproj.database import query_utm_crs_info\n", "\n", "\n", "def get_utm(point):\n", " longitude, latitude = point.x, point.y\n", " buffer = 0.001\n", " utm_crs_list = query_utm_crs_info(\n", " datum_name=\"WGS 84\",\n", " area_of_interest=AreaOfInterest(\n", " west_lon_degree=longitude - buffer,\n", " south_lat_degree=latitude - buffer,\n", " east_lon_degree=longitude + buffer,\n", " north_lat_degree=latitude + buffer,\n", " ),\n", " )\n", " utm_crs = CRS.from_epsg(utm_crs_list[0].code)\n", " return utm_crs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Figure out where to query\n", "\n", "* Define a GeoJSON point geometry for [The Bean](https://en.wikipedia.org/wiki/Cloud_Gate)\n", "* Buffer that in UTM by 400m\n", "* Reproject it back to EPSG:4326 so we can use that to query the STAC API\n", "* Plot it so it looks ok\n" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "tags": [] }, "outputs": [ { "data": { "text/html": [ "