{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Visualizing MAAP STAC Dataset with MosaicJSON\n", "\n", "Authors: Samuel Ayers (UAH), Alex Mandel (Development Seed), Aimee Barciauskas (Development Seed)\n", "\n", "Date: July 23, 2021\n", "\n", "Description: In this notebook, we visualize SRTM Cloud-Optimized GeoTIFFs (COGs) from MAAP's STAC using a generated mosaic from MAAP's TiTiler." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Run This Notebook\n", "To access and run this tutorial within MAAP's Algorithm Development Environment (ADE), please refer to the [\"Getting started with the MAAP\"](https://docs.maap-project.org/en/latest/getting_started/getting_started.html) section of our documentation.\n", "\n", "Disclaimer: it is highly recommended to run a tutorial within MAAP's ADE, which already includes packages specific to MAAP, such as maap-py. Running the tutorial outside of the MAAP ADE may lead to errors." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Note About the Data\n", "\n", "The NASA Shuttle Radar Topographic Mission (SRTM) has provided digital elevation data (DEMs) for over 80% of the globe. This data is currently distributed free of charge by USGS and is available for download from the National Map Seamless Data Distribution System, or the USGS FTP site.\n", "\n", "At MAAP, we've converted this elevation data into Cloud-Optimized GeoTIFFs (COGs) so they can be efficiently queried and visualized. These COGs are available in the [MAAP STAC](https://stac.maap-project.org/collections/SRTMGL1_COD/items)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Additional Resources\n", "- [USGS EROS Archive - Shuttle Radar Topography Mission (SRTM)](https://www.usgs.gov/centers/eros/science/usgs-eros-archive-digital-elevation-shuttle-radar-topography-mission-srtm)\n", "- [TiTiler API Documentation - MosaicJSON](https://titiler.maap-project.org/docs#/MosaicJSON)\n", "- [Github - MosaicJSON](https://github.com/developmentseed/mosaicjson-spec)\n", "- [Working with MosaicJSON - TiTiler (Development Seed Example)](https://developmentseed.org/titiler/examples/notebooks/Working_with_MosaicJSON/)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Importing and Installing Packages\n", "\n", "To be able to run this notebook you'll need the following requirements:\n", "\n", "- rasterio\n", "- folium\n", "- requests\n", "- pystac-client\n", "- cogeo-mosaic (Optional)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If the packages below are not installed already, uncomment the following cell:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "tags": [] }, "outputs": [], "source": [ "# %pip install rasterio\n", "# %pip install folium\n", "# %pip install requests\n", "# %pip install pystac-client\n", "# %pip install cogeo-mosaic --pre" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "tags": [] }, "outputs": [], "source": [ "import requests\n", "\n", "from pystac_client import Client\n", "\n", "from pprint import pprint\n", "\n", "from rasterio.features import bounds as featureBounds\n", "\n", "from folium import Map, TileLayer, GeoJson" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fetch SRTM COG STAC Items" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "tags": [] }, "outputs": [], "source": [ "cat = Client.open('https://stac.maap-project.org/')\n", "\n", "items = list(\n", " cat.search(\n", " collections=\"SRTMGL1_COD\",\n", " bbox=[4, 42, 16, 48],\n", " max_items=120\n", " ).items_as_dicts(),\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Map the Data Bounds" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "tags": [] }, "outputs": [ { "data": { "text/html": [ "
Make this Notebook Trusted to load map: File -> Trust Notebook
" ], "text/plain": [ "" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "geojson = {\"type\": \"FeatureCollection\", \"features\": items}\n", "\n", "bounds = featureBounds(geojson)\n", "zoom_start = 5\n", "\n", "m = Map(\n", " tiles=\"OpenStreetMap\",\n", " location=((bounds[1] + bounds[3]) / 2, (bounds[0] + bounds[2]) / 2),\n", " zoom_start=zoom_start,\n", ")\n", "\n", "geo_json = GeoJson(\n", " data=geojson,\n", " style_function=lambda x: {\n", " \"opacity\": 1,\n", " \"dashArray\": \"1\",\n", " \"fillOpacity\": 0,\n", " \"weight\": 1,\n", " },\n", ")\n", "geo_json.add_to(m)\n", "m" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create Mosaic" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We're using the TiTiler deployed by MAAP" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "tags": [] }, "outputs": [], "source": [ "titiler_endpoint = \"https://titiler.maap-project.org\" # MAAP titiler endpoint" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To begin, we'll pull our COG:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'band_descriptions': [['b1', '']],\n", " 'band_metadata': [['b1', {}]],\n", " 'bounds': [15.99986111111111,\n", " 47.999861111111116,\n", " 17.000138888888888,\n", " 49.00013888888889],\n", " 'colorinterp': ['gray'],\n", " 'count': 1,\n", " 'driver': 'GTiff',\n", " 'dtype': 'int16',\n", " 'height': 3601,\n", " 'maxzoom': 12,\n", " 'minzoom': 8,\n", " 'nodata_type': 'Nodata',\n", " 'nodata_value': -32768.0,\n", " 'overviews': [2, 4, 8, 16],\n", " 'width': 3601}\n" ] } ], "source": [ "cog_info = requests.get(\n", " f\"{titiler_endpoint}/cog/info\",\n", " params={\n", " \"url\": geojson[\"features\"][0][\"assets\"][\"cog_default\"][\"href\"],\n", " },\n", ").json()\n", "pprint(cog_info)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we will create the mosaic using the geojson from above. SRTMGL1 COGs have a \"cog default\" asset type, so we can create a mosaic of these type=\"image/tiff\" assets. We can get access to these assets by using the accessor function." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "tags": [] }, "outputs": [], "source": [ "from cogeo_mosaic.mosaic import MosaicJSON\n", "\n", "mosaicdata = MosaicJSON.from_features(\n", " geojson.get(\"features\"),\n", " minzoom=cog_info[\"minzoom\"],\n", " maxzoom=cog_info[\"maxzoom\"],\n", " accessor=lambda feature: feature[\"assets\"][\"cog_default\"][\"href\"],\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Alternatively, we can build a MosaicJSON using a list of files within MAAP's ADE by converting the local \"my-private\", \"my-public\", and \"shared-buckets\" paths to their respective AWS S3 prefixes like so:\n", "\n", "```python\n", "from maap.maap import MAAP\n", "maap = MAAP(maap_host='api.maap-project.org')\n", "username = maap.profile.account_info()[\"username\"]\n", "\n", "local_path = \"/local/path/to/files/\"\n", "\n", "files = glob.glob(os.path.join(dps_output, \"*.tif\"), recursive=False)\n", "\n", "def local_to_s3(url):\n", " ''' A Function to convert local paths to s3 urls'''\n", " if url.startswith(\"my-private-bucket\"):\n", " return url.replace(\"my-private-bucket\", f\"s3://maap-ops-workspace/{username}\")\n", "\n", " if url.startswith(\"my-public-bucket\"):\n", " return url.replace(\"my-public-bucket\", f\"s3://maap-ops-workspace/shared/{username}\")\n", "\n", " if url.startswith(\"shared-buckets\"):\n", " return url.replace(\"shared-buckets\", \"s3://maap-ops-workspace/shared\")\n", "\n", "tiles = [local_to_s3(file) for file in files]\n", "\n", "mosaicdata = MosaicJSON.from_urls(tiles, minzoom=9, maxzoom=16)\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we will upload the mosaicjson to the TiTiler:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'id': 'a6739144-b9ab-4ef6-82c3-f092100858cb',\n", " 'links': [{'href': 'https://titiler.maap-project.org/mosaics/a6739144-b9ab-4ef6-82c3-f092100858cb',\n", " 'rel': 'self',\n", " 'title': 'Self',\n", " 'type': 'application/json'},\n", " {'href': 'https://titiler.maap-project.org/mosaics/a6739144-b9ab-4ef6-82c3-f092100858cb/mosaicjson',\n", " 'rel': 'mosaicjson',\n", " 'title': 'MosaicJSON',\n", " 'type': 'application/json'},\n", " {'href': 'https://titiler.maap-project.org/mosaics/a6739144-b9ab-4ef6-82c3-f092100858cb/tilejson.json',\n", " 'rel': 'tilejson',\n", " 'title': 'TileJSON',\n", " 'type': 'application/json'},\n", " {'href': 'https://titiler.maap-project.org/mosaics/a6739144-b9ab-4ef6-82c3-f092100858cb/tiles/{z}/{x}/{y}',\n", " 'rel': 'tiles',\n", " 'title': 'Tiles',\n", " 'type': 'application/json'},\n", " {'href': 'https://titiler.maap-project.org/mosaics/a6739144-b9ab-4ef6-82c3-f092100858cb/WMTSCapabilities.xml',\n", " 'rel': 'wmts',\n", " 'title': 'WMTS',\n", " 'type': 'application/json'}]}\n" ] } ], "source": [ "mosaic_links = requests.post(\n", " url=f\"{titiler_endpoint}/mosaics\",\n", " headers={\n", " \"Content-Type\": \"application/vnd.titiler.mosaicjson+json\",\n", " },\n", " json=mosaicdata.model_dump(exclude_none=True),\n", ").json()\n", "\n", "pprint(mosaic_links)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The response from the post request gives endpoints for different services (eg. mosaicjson, tilejson, tiles, wmts, etc). We're fetching the `tilejson` endpoint. " ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "tags": [] }, "outputs": [ { "data": { "text/plain": [ "[{'href': 'https://titiler.maap-project.org/mosaics/a6739144-b9ab-4ef6-82c3-f092100858cb/tilejson.json',\n", " 'rel': 'tilejson',\n", " 'type': 'application/json',\n", " 'title': 'TileJSON'}]" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "tilejson_endpoint = list(\n", " filter(lambda x: x.get(\"rel\") == \"tilejson\", dict(mosaic_links)[\"links\"])\n", ")\n", "tilejson_endpoint" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Display Tiles\n", "\n", "NOTE: The important bit is the \"tiles\" endpoint returned from `f\"{titiler_endpoint}/mosaics`. This endpoint (e.g. `https://titiler.maap-project.org/mosaics/{mosaic_id}/tiles/{z}/{x}/{y}`) could be used in any map client which can tile `x,y,z` layers.\n", "\n", "Read more about this endpoint at [MAAP's TiTiler docs](https://titiler.maap-project.org/docs#/MosaicJSON)." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "tags": [] }, "outputs": [ { "data": { "text/html": [ "
Make this Notebook Trusted to load map: File -> Trust Notebook
" ], "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "params = {\"colormap_name\": \"viridis\", \"rescale\": \"0,4000\"}\n", "\n", "r_te = requests.get(tilejson_endpoint[0][\"href\"], params=params).json()\n", "\n", "tiles = TileLayer(tiles=f\"{r_te['tiles'][0]}\", opacity=1, attr=\"USGS\")\n", "\n", "tiles.add_to(m)\n", "m" ] }, { "cell_type": "code", "execution_count": null, "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" } }, "nbformat": 4, "nbformat_minor": 4 }