quickstarts/GDAL-and-STAC.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b1968d48-fe12-45bf-bc87-0075eefec408",
   "metadata": {},
   "source": [
    "## Example Notebook on getting data from STAC directly with GDAL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b86a756f-499f-4944-b63b-537a1d728291",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81c0f7bf-b23e-4a87-ac66-e54572e2d00c",
   "metadata": {},
   "source": [
    "Due to a Bug in GDAL, which did not accept local file:// URIS that terrabyte uses please make sure you have at minimum GDAL Version 3.8.0 installed. (If this is not available yet you can use the master branch installabe from the gdal-master conda channel "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "aa633c47-6a0d-4fe4-bed9-d0a260acef6e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GDAL 3.8.0dev-c95490018a-dirty, released 2023/07/21\n"
     ]
    }
   ],
   "source": [
    "!gdalinfo --version"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55e89101-110f-4b76-8b2d-d4ef84514be0",
   "metadata": {},
   "source": [
    "GDAL directly uses the search API of STAC. The resulting data of a search query can be used as input for a GDAL command. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "815cee1d-19d0-4f88-abb9-00ed741f6161",
   "metadata": {},
   "outputs": [],
   "source": [
    "# GDAL directly uses the search API of STAC\n",
    "stac_search_url=\"https://stac.terrabyte.lrz.de/public/api/search\"\n",
    "# STAC Collection to search\n",
    "collection=\"cop-dem-glo-90\"\n",
    "# Name of the Asset in the collection above we want do download\n",
    "asset=\"dem\"\n",
    "\n",
    "# Area of interest \n",
    "bboxMunich = [11.089325,47.893787,11.923599,48.329322 ]\n",
    "#so we can switch bounding box easily\n",
    "bbox=bboxMunich\n",
    "\n",
    "#extra gdal option to pass. Here we activate Bigtiff support and LZW Compression as well as give a name for the Band in the Format ArcGIS expects\n",
    "gdalOptions=f\"-co COMPRESS=LZW -co BIGTIFF=YES -mo band_1={asset}\"\n",
    "outputFilename=\"gdal_COPDEM_example.tif\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "93a7e57e-f1f3-418f-90b1-fa8a5b1e92d0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'STACIT:\"https://stac.terrabyte.lrz.de/public/api/search?collections=cop-dem-glo-90&bbox=11.089325,47.893787,11.923599,48.329322\":asset=data'"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# construct the gdal URL to access STAC search with our parameters\n",
    "gdal_inputURl=f\"STACIT:\\\"{stac_search_url}?collections={collection}&bbox={','.join(str(b) for b in bbox)}\\\":asset={asset}\"\n",
    "gdal_inputURl"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d712e5a-a7c6-45f2-9098-276e936c8824",
   "metadata": {},
   "source": [
    "Specifying the Bounding Box as as search Parameter to STAC only returns those results that intersect the bounding box"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "5530bb76-cc41-49fe-adb6-8d1e534d176a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gdalinfo 'STACIT:\"https://stac.terrabyte.lrz.de/public/api/search?collections=cop-dem-glo-30&bbox=11.089325,47.893787,11.923599,48.329322\":asset=data'\n"
     ]
    }
   ],
   "source": [
    "infoCmd=f\"gdalinfo '{gdal_inputURl}'\"\n",
    "print(infoCmd)\n",
    "#os.system(f\"echo {cmdUrl}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "49a6d43f-eee6-46f9-8981-944be89d19b2",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "256"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "os.system(infoCmd)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd17c534-2d1a-42f9-9832-dd663efbf354",
   "metadata": {},
   "source": [
    "We specify the Bounding Box again to gdal itself to clip the output exactly to the bounding box. Be aware that gdal expects a different Order of Coordinates "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5d56d8bf-3fa4-424c-98a3-074d72c51ba8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'-projwin 11.089325 48.329322 11.923599 47.893787'"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gdal_projwin=f\"-projwin {' '.join(str(bbox[i]) for i in [0,3,2,1])}\"\n",
    "gdal_projwin                                  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "9691eaf0-fd04-45c1-85a3-556d26fe225e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gdal_translate 'STACIT:\"https://stac.terrabyte.lrz.de/public/api/search?collections=cop-dem-glo-90&bbox=11.089325,47.893787,11.923599,48.329322\":asset=data' gdal_COPDEM_example.tif -projwin 11.089325 48.329322 11.923599 47.893787 -co COMPRESS=LZW -co BIGTIFF=YES -mo band_1=data\n"
     ]
    }
   ],
   "source": [
    "translateCmd=f\"gdal_translate '{gdal_inputURl}' {outputFilename} {gdal_projwin} {gdalOptions}\"\n",
    "print(translateCmd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "e8bc64b9-cf94-490b-8ef1-6f6e189c7269",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Input file size is 1201, 2401\n",
      "0...10...20...30...40...50...60...70...80...90...100 - done.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "os.system(translateCmd)"
   ]
  }
 ],
 "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}