OT_BulkAccessCOGs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "unknown-witch",
   "metadata": {},
   "source": [
    "# Programmatic Access to OpenTopography’s Cloud Optimized GeoTIFF (COG) Global Datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "indian-silence",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Authors\" data-toc-modified-id=\"Authors-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Authors</a></span></li><li><span><a href=\"#Purpose\" data-toc-modified-id=\"Purpose-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>Purpose</a></span></li><li><span><a href=\"#Funding\" data-toc-modified-id=\"Funding-3\"><span class=\"toc-item-num\">3&nbsp;&nbsp;</span>Funding</a></span></li><li><span><a href=\"#Keywords\" data-toc-modified-id=\"Keywords-4\"><span class=\"toc-item-num\">4&nbsp;&nbsp;</span>Keywords</a></span></li><li><span><a href=\"#Citation\" data-toc-modified-id=\"Citation-5\"><span class=\"toc-item-num\">5&nbsp;&nbsp;</span>Citation</a></span></li><li><span><a href=\"#Setup\" data-toc-modified-id=\"Setup-6\"><span class=\"toc-item-num\">6&nbsp;&nbsp;</span>Setup</a></span><ul class=\"toc-item\"><li><span><a href=\"#Library-Imports\" data-toc-modified-id=\"Library-Imports-6.1\"><span class=\"toc-item-num\">6.1&nbsp;&nbsp;</span>Library Imports</a></span></li></ul></li><li><span><a href=\"#Data-Access-and-Visualization\" data-toc-modified-id=\"Data-Access-and-Visualization-7\"><span class=\"toc-item-num\">7&nbsp;&nbsp;</span>Data Access and Visualization</a></span><ul class=\"toc-item\"><li><span><a href=\"#Access-SRTM-Tiles-with-AWS-CLI\" data-toc-modified-id=\"Access-SRTM-Tiles-with-AWS-CLI-7.1\"><span class=\"toc-item-num\">7.1&nbsp;&nbsp;</span>Access SRTM Tiles with AWS CLI</a></span></li><li><span><a href=\"#Plot-SRTM-Tiles\" data-toc-modified-id=\"Plot-SRTM-Tiles-7.2\"><span class=\"toc-item-num\">7.2&nbsp;&nbsp;</span>Plot SRTM Tiles</a></span></li><li><span><a href=\"#SRTM-Bulk-Access-via-GDAL-VRTs\" data-toc-modified-id=\"SRTM-Bulk-Access-via-GDAL-VRTs-7.3\"><span class=\"toc-item-num\">7.3&nbsp;&nbsp;</span>SRTM Bulk Access via GDAL VRTs</a></span></li><li><span><a href=\"#Plot-SRTM-Mosaic\" data-toc-modified-id=\"Plot-SRTM-Mosaic-7.4\"><span class=\"toc-item-num\">7.4&nbsp;&nbsp;</span>Plot SRTM Mosaic</a></span></li><li><span><a href=\"#SRTM-Access-via-OpenTopography's-API\" data-toc-modified-id=\"SRTM-Access-via-OpenTopography's-API-7.5\"><span class=\"toc-item-num\">7.5&nbsp;&nbsp;</span>SRTM Access via OpenTopography's API</a></span></li></ul></li><li><span><a href=\"#Data-Processing\" data-toc-modified-id=\"Data-Processing-8\"><span class=\"toc-item-num\">8&nbsp;&nbsp;</span>Data Processing</a></span><ul class=\"toc-item\"><li><span><a href=\"#Download-Colorado-River-Vector-Data\" data-toc-modified-id=\"Download-Colorado-River-Vector-Data-8.1\"><span class=\"toc-item-num\">8.1&nbsp;&nbsp;</span>Download Colorado River Vector Data</a></span></li><li><span><a href=\"#Clip-River-Data\" data-toc-modified-id=\"Clip-River-Data-8.2\"><span class=\"toc-item-num\">8.2&nbsp;&nbsp;</span>Clip River Data</a></span></li><li><span><a href=\"#Buffer-River-Centerline\" data-toc-modified-id=\"Buffer-River-Centerline-8.3\"><span class=\"toc-item-num\">8.3&nbsp;&nbsp;</span>Buffer River Centerline</a></span></li><li><span><a href=\"#Dissolve-River-Buffer-Polygon\" data-toc-modified-id=\"Dissolve-River-Buffer-Polygon-8.4\"><span class=\"toc-item-num\">8.4&nbsp;&nbsp;</span>Dissolve River Buffer Polygon</a></span></li><li><span><a href=\"#Overplot-River-Shapefiles-on-DEM\" data-toc-modified-id=\"Overplot-River-Shapefiles-on-DEM-8.5\"><span class=\"toc-item-num\">8.5&nbsp;&nbsp;</span>Overplot River Shapefiles on DEM</a></span></li><li><span><a href=\"#Calculate-Projected-Bounds-of-River-Polygon\" data-toc-modified-id=\"Calculate-Projected-Bounds-of-River-Polygon-8.6\"><span class=\"toc-item-num\">8.6&nbsp;&nbsp;</span>Calculate Projected Bounds of River Polygon</a></span></li><li><span><a href=\"#Clip-Projected-River-Polygon\" data-toc-modified-id=\"Clip-Projected-River-Polygon-8.7\"><span class=\"toc-item-num\">8.7&nbsp;&nbsp;</span>Clip Projected River Polygon</a></span></li><li><span><a href=\"#Plot-Clipped-River-Buffer-on-DEM\" data-toc-modified-id=\"Plot-Clipped-River-Buffer-on-DEM-8.8\"><span class=\"toc-item-num\">8.8&nbsp;&nbsp;</span>Plot Clipped River Buffer on DEM</a></span></li><li><span><a href=\"#Extracting-and-Masking-SRTM-Data-in-One-Step\" data-toc-modified-id=\"Extracting-and-Masking-SRTM-Data-in-One-Step-8.9\"><span class=\"toc-item-num\">8.9&nbsp;&nbsp;</span>Extracting and Masking SRTM Data in One Step</a></span></li></ul></li><li><span><a href=\"#Conclusion\" data-toc-modified-id=\"Conclusion-9\"><span class=\"toc-item-num\">9&nbsp;&nbsp;</span>Conclusion</a></span></li><li><span><a href=\"#Resources\" data-toc-modified-id=\"Resources-10\"><span class=\"toc-item-num\">10&nbsp;&nbsp;</span>Resources</a></span></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "preceding-threshold",
   "metadata": {},
   "source": [
    "## Authors\n",
    "-  Author1 = {\"name\": \"Matthew Beckley\", \"affiliation\": \"UNAVCO, Inc.\", \"email\": \"beckley@unavco.org\", \"orcid\":\"\"}\n",
    "-  Author2 = {\"name\": \"Christopher Crosby\", \"affiliation\": \"UNAVCO, Inc.\", \"email\": \"crosby@unavco.org\", \"orcid\": \"https://orcid.org/0000-0003-2522-4193\"}\n",
    "-  Author3 = {\"name\": \"Viswanath Nandigam\", \"affiliation\": \"San Diego Supercomputer Center, UC San Diego\", \"email\": \"viswanat@sdsc.edu\", \"orcid\": \"https://orcid.org/0000-0003-0928-9851\"}\n",
    "-  Author4 = {\"name\": \"Minh Phan\", \"affiliation\": \"San Diego Supercomputer Center, UC San Diego\", \"email\": \"mnphan@sdsc.edu\",\"orcid\":\"https://orcid.org/0000-0002-2806-5450\"}\n",
    "\n",
    "## Purpose\n",
    "[OpenTopography](https://opentopography.org/) strives to provide a variety of access methods to our datasets to address the varying needs of our users. While most users access data through the web map interface, others prefer to programmatically access data from our cloud-based storage. In this tutorial we will explore a couple of different methods of accessing data over the Grand Canyon from the cloud-based [NASA Shuttle Radar Topography Mission (SRTM) catalog in OpenTopography](https://doi.org/10.5069/G9445JDF).  We will also demonstrate how we can use open-source, command line tools to programmatically access data without the need for a web-map interface. \n",
    "\n",
    "In addition, this tutorial should highlight the power of [Cloud Optimized GeoTIFFs (COGs),](https://www.cogeo.org/) and how they can be used to reduce download file sizes as well as increase data access speeds.  COGs are quickly becoming a standard format for hosting raster data in the cloud due to its optimized file structure to enable faster access and subsetting. In order to take advantage of this benefit, [OpenTopography](https://opentopography.org/) has recently converted its entire global dataset collection to COGs. With the COG format, programs and web requests can read the top-level header and quickly extract data for a region of interest as opposed to having to read the entire file.  This drastically reduces file read-time as well as reducing the number of http GET requests to access the data.  For a more detailed discussion on COGs and their benefits see the resources section at the bottom of this page.\n",
    "\n",
    "## Funding\n",
    "OpenTopography is supported by the National Science Foundation under Award Numbers 1948997, 1948994 & 1948857\n",
    "\n",
    "## Keywords\n",
    "keywords=[\"COG\", \"OpenTopography\", \"SRTM\", \"GDAL\", \"VRT\"]\n",
    "\n",
    "## Citation\n",
    "To cite this notebook:  Beckley, M., Crosby, C., Nandigam, V. & Phan, M. (2021). Programmatic Access to OpenTopography’s Cloud Optimized GeoTIFF (COG) Global Datasets (Version v1.0). https://doi.org/10.5069/G9BC3WRK Accessed: MM/DD/YYYY"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "boxed-island",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "If running this tutorial locally: set up a new conda environment with GDAL 3.2 and AWS CLI tools installed:\n",
    "\n",
    "`# Create a conda environment called, 'COGS', and install GDAL v 3.2:\n",
    " conda create -n COGS -c conda-forge gdal==3.2.0`\n",
    "\n",
    " `#Activate your new environment\n",
    " conda activate COGS` \n",
    " \n",
    " `#install AWS command line tools in your environment\n",
    " conda install --yes -c conda-forge awscli`\n",
    "\n",
    " We will also be using the following python libraries:\n",
    " -  rasterio\n",
    " -  fiona\n",
    " -  wget\n",
    " -  pyproj\n",
    " -  geopandas\n",
    " \n",
    "Within this conda environment, install these libraries with pip:\n",
    "`pip install pip install rasterio fiona wget geopandas pyproj matplotlib`\n",
    " \n",
    "If you are having difficulties getting the environment set up, you can always run the notebook from [Binder](https://mybinder.org/v2/gh/OpenTopography/Jupyter/main?filepath=OT_01_BulkAccessCOGs.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "material-powder",
   "metadata": {},
   "source": [
    "### Library Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fantastic-eugene",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Visualization\n",
    "import rasterio as rio                                                                                                                  \n",
    "from rasterio.plot import show  \n",
    "import rasterio.mask\n",
    "import matplotlib.pyplot as plt                                                                                                         \n",
    "\n",
    "#Working with Shapefiles\n",
    "import geopandas as gpd \n",
    "import fiona\n",
    "\n",
    "#Projections and Transformations\n",
    "from pyproj import Transformer\n",
    "\n",
    "#Downloads\n",
    "import wget\n",
    "\n",
    "print(\"Success! Imported all necessary libraries.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "lightweight-notice",
   "metadata": {},
   "source": [
    "## Data Access and Visualization\n",
    "\n",
    "Below is a Google Earth image of a section of the Grand Canyon, AZ where we will extract data and do some data processing:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "continental-benefit",
   "metadata": {},
   "source": [
    "![GrandCanyon](img/GC_Basemap.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "raising-cleaning",
   "metadata": {},
   "source": [
    "---\n",
    "### Access SRTM Tiles with AWS CLI\n",
    "The [SRTM catalog in OpenTopography](https://doi.org/10.5069/G9445JDF) is organized into tiled COGs where each tile is named after its lower left corner (e.g. N35W112.tif is the tile that has its South-West corner at 35 Degrees North, 112 degrees West). \n",
    "\n",
    "For this exercise we will gather data around the Grand Canyon with a lat/lon bounds of [-114, 35, -111, 37].  \n",
    "\n",
    "Using this Area of Interest (AOI) information we can grab the tiles that cover this region.  OpenTopography users can access all hosted data via AWS Command Line Tools.  So, knowing the lat/lon bounds of our area of interest, we can download just the tiles for this region using [aws command line tools](https://anaconda.org/conda-forge/awscli).  \n",
    "\n",
    "It is useful to incorporate wildcards to limit the data you need to download.  The example belows downloads all files that match: N35W112.tif, N35W113.tif, N35W114.tif, N36W114.tif, N36W112.tif, & N36W113.tif \n",
    "\n",
    "Within Jupyter Notebooks, shell commands can be executed by starting the command with an exclamation mark (\"!\"). In this notebook, we will explore a variety of tools - using both shell and python commands to demonstrate different methods to complete common geoprocessing tasks.\n",
    "\n",
    "**NOTE:  The AWS copy command may take up to a minute to complete depending on your internet connection.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dominant-snapshot",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "print(\"Downloading files from AWS S3...\")\n",
    "!aws s3 cp s3://raster/SRTM_GL1/SRTM_GL1_srtm/ . --recursive --exclude \"*\" --include \"N3[56]W11[234].tif\" --endpoint-url https://opentopography.s3.sdsc.edu --no-sign-request\n",
    "print(\"Success!  Download complete.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "induced-germany",
   "metadata": {},
   "source": [
    "### Plot SRTM Tiles\n",
    "Now that we have downloaded the SRTM tiles that cover the Grand Canyon, we can plot them using the Python Library, [Rasterio:](https://rasterio.readthedocs.io/en/latest/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "finnish-witch",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "images =['N35W112.tif','N35W113.tif','N35W114.tif','N36W112.tif','N36W113.tif','N36W114.tif'] \n",
    "fig, ax = plt.subplots(figsize=(10,10))\n",
    "plt.xlim([-114, -111])\n",
    "plt.ylim([35, 37])\n",
    "\n",
    "for im in images:                                                                                                                       \n",
    "    print('loading: ',im)\n",
    "    raster = rio.open(im)\n",
    "    show(raster,ax=ax,cmap='gray')\n",
    "print(\"Success! Plotted SRTM Tiles.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "secondary-image",
   "metadata": {},
   "source": [
    "Downloading the data in this manner, we end up with 6 tifs that broadly cover our area. But what if we wanted to confine our search area around the Colorado River, and reduce time spent on post processing of the data?  GDAL's [Virtual Raster Datasets (VRT)](https://gdal.org/programs/gdalbuildvrt.html) can help in this regard."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "foreign-plenty",
   "metadata": {},
   "source": [
    "### SRTM Bulk Access via GDAL VRTs\n",
    "As part of our ingest process, OpenTopography builds GDAL [Virtual Raster Datasets (VRTs)](https://gdal.org/programs/gdalbuildvrt.html) for all its rasters to enable faster and more efficient data access. With GDAL and VRTs, users can extract a specific area of interest and obtain a merged, clipped raster all in one operation. For example, by using the vsicurl command with a targeted bounding box, we can extract a raster that covers the area of the Grand Canyon around the Colorado River."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "british-lawyer",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!time gdal_translate /vsicurl/https://opentopography.s3.sdsc.edu/raster/SRTM_GL1/SRTM_GL1_srtm.vrt -projwin -113.98 36.563 -111.777 35.7 -co\\\n",
    " COMPRESS=deflate srtm_GC.tif\n",
    "print(\"Success!  Downloaded SRTM from VRT.\")\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "documented-brisbane",
   "metadata": {},
   "source": [
    "So the above command runs pretty quickly, and if using GDAL version 3.2, which has native COG support, we are able to output a COG-formatted mosaic of SRTM data over our area of interest.  Note the '-co COMPRESS-deflate' command additionally compresses the output file using the 'deflate' algorithm to reduce the output filesize.  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "weighted-flash",
   "metadata": {},
   "source": [
    "### Plot SRTM Mosaic\n",
    "Using the same plotting commands as earlier, we can plot the single, mosaiced image to see that we now have a single image that is clipped to the portion of the Grand Canyon that is of interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "neural-fence",
   "metadata": {},
   "outputs": [],
   "source": [
    "image ='srtm_GC.tif' \n",
    "fig, ax = plt.subplots(figsize=(10,10))\n",
    "plt.xlim([-114, -111])\n",
    "plt.ylim([35, 37])\n",
    "\n",
    "print('loading: ',image)\n",
    "raster = rio.open(image)\n",
    "show(raster,ax=ax,cmap='gray')\n",
    "print(\"Success! Plotted SRTM Mosaic.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "developmental-forward",
   "metadata": {},
   "source": [
    "---\n",
    "### SRTM Access via OpenTopography's API\n",
    "\n",
    "Users can get the same data product by using [OpenTopography's API](https://opentopography.org/developers).  Users can [build the REST API call](https://portal.opentopography.org/apidocs/#/Public/getGlobalDem) and it will output the appropriate request to use in a curl command, or a http GET request.  Here is the URL request that can be used via https for SRTM GL1 data over the GrandCanyon:\n",
    "\n",
    "`https://portal.opentopography.org/API/globaldem?demtype=SRTMGL1&south=35.7&north=36.54&west=-114&east=-111.7&outputFormat=GTiff`\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "enclosed-stability",
   "metadata": {},
   "source": [
    "---\n",
    "## Data Processing\n",
    "### Download Colorado River Vector Data\n",
    "\n",
    "By using the SRTM VRT, we have now have a single, mosaiced GeoTIFF that is much smaller in size than the original tiled images. But we can take this further by using the power of COGs and other Python libraries. So, let's say we really only want data around the Colorado River. To do this, we first need to get a vector representation of the Colorado River. Luckily, there is a fantastic site called, [Natural Earth Data](https://www.naturalearthdata.com/) that has a wealth of vector and raster map data at a variety of scales.  For this example, we'll use the rivers and lakes centerlines file. You can download the data from the [Natural Earth Data Download site](https://www.naturalearthdata.com/downloads/10m-physical-vectors/) or you can use wget."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "overhead-patrick",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "wget.download('https://naturalearth.s3.amazonaws.com/10m_physical/ne_10m_rivers_lake_centerlines.zip','rivers.zip')\n",
    "print(\"Success! Downloaded river vector data.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "velvet-packaging",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "!unzip -o rivers.zip\n",
    "print(\"Success! Extracted river shapefiles from archive.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "alpha-fleet",
   "metadata": {},
   "source": [
    "---\n",
    "### Clip River Data\n",
    "Now we have river data, but it is for the entire world. Next we'll clip the river dataset to our area of interest using [GDAL's vector toolkit](https://gdal.org/programs/ogr2ogr.html#ogr2ogr). We also want to create a buffer around the river, but in order to do that we need to project the data. As downloaded, the river data is in a Geographic Coordinate System, so we will need to output the data in a projected coordinate system in order to do the buffering operation. Arizona is in UTM zone 12N, so we'll use [EPSG:32612](https://epsg.io/32612) for the output projection. Using the ogr2ogr tool we can output a projected, clipped version of the river shapefile with the following command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "intensive-jackson",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!ogr2ogr ColoradoRiver_AOI.shp ne_10m_rivers_lake_centerlines.shp -t_srs EPSG:32612 -spat -113.98 35.7 -111.777 36.563 -clipsrc spat_extent\n",
    "print(\"Success! Clip global river data to Grand Canyon.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cordless-wesley",
   "metadata": {},
   "source": [
    "### Buffer River Centerline\n",
    "Next we need to create a buffer around the Colorado River.  We'll make a buffer that is 10km around either side so that we include the interesting features of the Grand Canyon.  With OGR you can pass in [SQLite commands](https://gdal.org/user/sql_sqlite_dialect.html). In this example we will use the ST_Buffer command to create a 10km buffer around the Colorado River."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "lyric-measure",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!ogr2ogr ColoradoRiver_AOI_10kmBuffer.shp ColoradoRiver_AOI.shp -dialect sqlite -sql \"select ST_buffer(geometry,10000) as geometry from ColoradoRiver_AOI\"\n",
    "print(\"Success! Created buffer from river centerline.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sitting-carry",
   "metadata": {},
   "source": [
    "### Dissolve River Buffer Polygon\n",
    "Since the river shapefile had segments, we should dissolve the buffer so that we end up with a single dissolved polygon that covers the buffer region of the river. To do this, we can once again use the \"ST_Union\" command with sqlite:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "iraqi-original",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!ogr2ogr ColoradoRiver_AOI_Buffer_Dissolved.shp ColoradoRiver_AOI_10kmBuffer.shp -dialect sqlite -sql \"select ST_Union(geometry) as geometry from ColoradoRiver_AOI_10kmBuffer\"\n",
    "print(\"Success! Dissolved buffer polygon.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "after-disease",
   "metadata": {},
   "source": [
    "### Overplot River Shapefiles on DEM\n",
    "Using the GeoPandas library, we can easily read in the buffered river shapefile.  Note we are plotting in lat/lon space, so we will need to reproject the shapefile coordinates to latitude and longitude. This code shows yet another way of how you could reproject your coordinates - this time using the GeoPandas library:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "flying-robin",
   "metadata": {},
   "outputs": [],
   "source": [
    "image ='srtm_GC.tif'                                                                                                                    \n",
    "river_buffered = 'ColoradoRiver_AOI_Buffer_Dissolved.shp'                                                                               \n",
    "river = 'ColoradoRiver_AOI.shp'                                                                                                         \n",
    "fig, ax = plt.subplots(figsize=(10,10))                                                                                                   \n",
    "plt.xlim([-114, -111])                                                                                                                  \n",
    "plt.ylim([35, 37])                                                                                                                      \n",
    "                                                                                                                                        \n",
    "print('loading: ',image)                                                                                                                \n",
    "raster = rio.open(image)                                                                                                                \n",
    "show(raster,ax=ax,cmap='gray')                                                                                                                                                                                  \n",
    "\n",
    "#get the coordinate system of the raster\n",
    "raster_CRS = raster.crs                 \n",
    "\n",
    "#Read in the river centerline shapefile with geopandas\n",
    "river = gpd.read_file(river)                          \n",
    "\n",
    "#Convert the river centerline to the same coordinate system as the DEM (Lat/Lon)\n",
    "river_LL = river.to_crs(raster_CRS)                                            \n",
    "\n",
    "\n",
    "#perform the same operations for the buffered river shapefile                                                                                                                                        \n",
    "river_buffered = gpd.read_file(river_buffered)                                                                                          \n",
    "river_buffered_LL = river_buffered.to_crs(raster_CRS)                                                                                   \n",
    "                                                                                                                                        \n",
    "#Overplot the river centerline and the river buffer on the SRTM DEM\n",
    "river_LL.plot(ax=ax)                                                                                                                    \n",
    "river_buffered_LL.plot(ax=ax,alpha=0.4)  \n",
    "\n",
    "print(\"Success! Plotted river data on SRTM DEM.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "gorgeous-consumption",
   "metadata": {},
   "source": [
    "---\n",
    "### Calculate Projected Bounds of River Polygon\n",
    "Since we buffered the data, it now extends outside our AOI. We can clip it to the bounds, but the buffered shapefile is now in UTM 12N projected coordinates, so we need to find what the bounding box coordinates are in UTM 12N coordinates. There are many different libraries that can perform this operation, but for this step, we'll highlight using pyproj Python library and OGR to do this tranformation and clipping. First we will get the UTM bounds of our area of interest by converting the lat/lon bounds to UTM coordinates with a simple pyproj script: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "neither-detective",
   "metadata": {},
   "outputs": [],
   "source": [
    "#set up the transformation in pyproj\n",
    "transformer = Transformer.from_crs(\"epsg:4326\", \"epsg:32612\")\n",
    "\n",
    "#Convert lower-left corner of AOI\n",
    "LL_lon,LL_lat = -113.98,35.7\n",
    "LL_x,LL_y = transformer.transform(LL_lat,LL_lon)\n",
    "\n",
    "#Convert upper-right corner of AOI\n",
    "UR_lon,UR_lat = -111.777,36.563\n",
    "UR_x,UR_y = transformer.transform(UR_lat,UR_lon)\n",
    "\n",
    "print('Lower Left Corner in UTM X,Y is: {},{}\\nUpper Right Corner in UTM X,Y is: {},{}'.format(LL_x,LL_y,UR_x,UR_y))\n",
    "print(\"Success! Extracted UTM coordinates for river bounds.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "mechanical-portuguese",
   "metadata": {},
   "source": [
    "---\n",
    "### Clip Projected River Polygon\n",
    "Now we can feed these UTM bounding coordinates back into a ogr to clip the river polygon to the area of interest:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "interesting-jamaica",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!ogr2ogr ColoradoRiver_AOI_Final.shp ColoradoRiver_AOI_Buffer_Dissolved.shp -clipdst 230367.047752107 3954768.801329 430471.347916117 4046677.17926595\n",
    "print(\"Success! Clipped river buffer to area of interest.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "spoken-drill",
   "metadata": {},
   "source": [
    "### Plot Clipped River Buffer on DEM\n",
    "Now that we have a clipped river buffer, we can plot it, and the river centerline on top of the SRTM DEM to verify that the polygon has been clipped:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ethical-falls",
   "metadata": {},
   "outputs": [],
   "source": [
    "image ='srtm_GC.tif' \n",
    "river = 'ColoradoRiver_AOI.shp'\n",
    "river_buffered = 'ColoradoRiver_AOI_Final.shp'\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10,10))\n",
    "plt.xlim([-114, -111])\n",
    "plt.ylim([35, 37])\n",
    "\n",
    "print('loading: ',image)\n",
    "raster = rio.open(image)\n",
    "show(raster,ax=ax,cmap='gray')\n",
    "\n",
    "raster_CRS = raster.crs                                                                                                                 \n",
    "river = gpd.read_file(river)                                                                                                            \n",
    "river_LL = river.to_crs(raster_CRS)                                                                                                     \n",
    "                                                                                                                                        \n",
    "river_buffered = gpd.read_file(river_buffered)                                                                                          \n",
    "river_buffered_LL = river_buffered.to_crs(raster_CRS)                                                                                   \n",
    "                                                                                                                                        \n",
    "                                                                                                                                        \n",
    "river_LL.plot(ax=ax)                                                                                                                    \n",
    "river_buffered_LL.plot(ax=ax,alpha=0.4) \n",
    "\n",
    "print(\"Success! Plotted clipped river data on SRTM DEM..\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "rental-costa",
   "metadata": {},
   "source": [
    "### Extracting and Masking SRTM Data in One Step \n",
    "Many applications can now take advantage of COGs. [Rasterio](https://rasterio.readthedocs.io/en/latest/) is a Python library that takes advantage of COGs ability to use the vsis3 and vsicurl file system handlers to perform HTTP read operations in a more efficient way. Ultimately, this benefit is passed on to the user by providing faster downloads that take up less disk space. As an example, we will use [Fiona](https://pypi.org/project/Fiona/) and [Rasterio](https://rasterio.readthedocs.io/en/latest/) to download only the SRTM data that fits within our 10km Colorado River buffer. \n",
    "\n",
    "Because we are ultimately going to intersect this boundary with the SRTM VRT, they both need to be in the same coordinate system.  So, we need to re-project this final boundary back to Geographic Coordinates so that it is compatible with the SRTM dataset. We can do this with a ogr2ogr command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "growing-stomach",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!ogr2ogr ColoradoRiver_AOI_Final_WGSLatLon.shp ColoradoRiver_AOI_Final.shp -s_srs \"EPSG:32612\" -t_srs \"EPSG:4326\"\n",
    "print(\"Success! Reprojected river buffer to WGS84 Lat/Lon\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "minimal-tiffany",
   "metadata": {},
   "source": [
    "Now that we have a mask in the same coordinate system as the SRTM VRT, we can read, mask, and plot the data with a few lines of code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "hidden-manhattan",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "inshp = 'ColoradoRiver_AOI_Final_WGSLatLon.shp'                                                         \n",
    "fp = 'https://opentopography.s3.sdsc.edu/raster/SRTM_GL1/SRTM_GL1_srtm.vrt'                                                                                     \n",
    "\n",
    "print(\"Loading Mask...\")\n",
    "#read in the buffered shapefile using fiona\n",
    "with fiona.open(inshp, \"r\") as shapefile:                                                                                                                       \n",
    "    shapes = [feature[\"geometry\"] for feature in shapefile]                                                                                                     \n",
    "                                                                                                                                                                \n",
    "#Using rasterio, read directly from the SRTM VRT of COGs, and crop\n",
    "# data to the buffer polygon\n",
    "print(\"Extracting SRTM to mask...\")\n",
    "with rasterio.open(fp) as src:                                                                                                                                  \n",
    "    out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)                                                                                       \n",
    "    out_meta = src.meta                                                                                                                                         \n",
    "\n",
    "#set up the output parameters\n",
    "out_meta.update({\"driver\": \"GTiff\",                                                                                                                             \n",
    "                 \"height\": out_image.shape[1],                                                                                                                  \n",
    "                 \"width\": out_image.shape[2],\n",
    "                 \"compress\": \"LZW\",\n",
    "                 \"tiled\":\"True\", \n",
    "                 \"blockxsize\":\"256\",\n",
    "                 \"blockysize\":\"256\",\n",
    "                 \"transform\": out_transform})                                                                                                                   \n",
    "                                 \n",
    "\n",
    "#Write out just the SRTM data that is within the buffer polygon\n",
    "print(\"Writing out Final SRTM DEM\")\n",
    "with rasterio.open(\"GrandCanyon_masked.tif\", \"w\", **out_meta) as dest:                                                                                                   \n",
    "    dest.write(out_image)     \n",
    "    \n",
    "fig, ax = plt.subplots(figsize=(10,10))\n",
    "plt.xlim([-114, -111])\n",
    "plt.ylim([35, 37])\n",
    "\n",
    "print('Plotting: ',\"GrandCanyon_masked.tif\")\n",
    "raster = rio.open(\"GrandCanyon_masked.tif\")\n",
    "show(raster,ax=ax,cmap='gray')\n",
    "\n",
    "print(\"Success! Plotted final result\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "massive-folder",
   "metadata": {},
   "source": [
    "---\n",
    "Using Rasterio with OpenTopography's VRT of COGs, we are able to download only the portion of the SRTM data we need. In this example, we downloaded only data that is within our 10 km polygon buffer of the portion of Colorado River around the Grand Canyon. **Our final output raster is only 12MB, versus the 86MB of tiffs we started with at the beginning!!** By putting all these tools together, users can create powerful scripts that automate data access, and ultimately reduce time on data downloads and post-processing.  \n",
    "\n",
    "\n",
    "## Conclusion\n",
    "Hopefully this tutorial has been helpful, and has demonstrated alternative methods to accessing cloud-based data. There are often many different ways to accomplish the same task, and we highlighted some of the more common tools and methods to perform some basic geospatial operations. We encourage users to explore datasets hosted at OpenTopography with these powerful tools. OpenTopography has various global datasets that are available through our API and/or bulk access as described in this article:\n",
    "* [SRTM](https://doi.org/10.5069/G9445JDF)\n",
    "* [ALOS World 3D - 30m](https://doi.org/10.5069/G94M92HB)\n",
    "* [SRTM15+](https://doi.org/10.5069/G92R3PT9)\n",
    "* [NASADEM](https://doi.org/10.5069/G93T9FD9)\n",
    "* [Copernicus](https://doi.org/10.5069/G9028PQB)\n",
    "\n",
    "\n",
    "\n",
    "## Resources\n",
    "Check out some of these excellent resources on GDAL, Rasterio, and COGs\n",
    "* [A good introduction to COGs](https://www.cogeo.org/ )\n",
    "* [Another good intro to COGs](https://medium.com/planet-stories/a-handy-introduction-to-cloud-optimized-geotiffs-1f2c9e716ec3)                    \n",
    "* [Detailed info about the COG format in GDAL](https://trac.osgeo.org/gdal/wiki/CloudOptimizedGeoTIFF)                 \n",
    "* [Python notebook on how to use Rasterio](https://gist.github.com/sgillies/7e5cd548110a5b4d45ac1a1d93cb17a3)\n",
    "* [An excellent explanation of COG file structure](https://medium.com/element84/cloud-optimized-geotiff-vs-the-meta-raster-format-d24c1a77dc2e)\n",
    "* [Reading COGs with Rasterio (and much more)](https://automating-gis-processes.github.io/CSC/notebooks/L5/read-cogs.html)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.2"
  },
  "toc": {
   "base_numbering": "1",
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "223px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}