LANDSAT_download_indices_from_polygon_aoi.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# LANDSAT8 - Downloading band data for multiple dates\n",
    "This code has been specifically written for Landsat ARD data. The idea is that you provide the specific dates that you wish to process over a given area (defined by the polygon), that you have already sussed for cloud QA. Be warned, datasets are big.  Always try to stay about 1gb under your threshold to ensure your work gets processed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import datacube\n",
    "import geopandas as gpd\n",
    "import numpy as np\n",
    "from datacube.utils import geometry\n",
    "from datacube.utils.cog import write_cog\n",
    "\n",
    "import sys\n",
    "sys.path.insert(1, '../Tools/')\n",
    "from dea_tools.datahandling import load_ard\n",
    "from dea_tools.plotting import map_shapefile\n",
    "from dea_tools.bandindices import calculate_indices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Connect to the datacube"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Temporary solution to account for Collection 3 data being in a different\n",
    "# database on the NCI\n",
    "try:\n",
    "    dc = datacube.Datacube(app='Analyse_multiple_polygons', env='c3-samples')\n",
    "except:\n",
    "    dc = datacube.Datacube(app='Analyse_multiple_polygons')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis parameters\n",
    "\n",
    "* `time_list` : Enter dates, in units YYYY-MM-DD, for each target date  e.g. `'2019-01-01'`\n",
    "* `vector_file` : A path to a vector file (ESRI Shapefile or GeoJSON). Code below presumes your vector to live in a directory called 'vectors' located in the same folder as this notebook. Please ensure it is projected to Albers Equal Area projection\n",
    "* `attribute_col` : A column in the vector file used to label the output `xarray` datasets containing satellite images. Each row of this column should have a unique identifier\n",
    "* `products` : Here we are using Landsat8 only\n",
    "* `measurements` : A list of band names to load from the satellite product e.g. `['nbart_red', 'nbart_green']`\n",
    "* `resolution` : The spatial resolution of the loaded satellite data e.g. for Landsat, this is `(-30, 30)`\n",
    "* `output_crs` : The coordinate reference system/map projection to load data into, e.g. `'EPSG:3577'` to load data in the Albers Equal Area projection\n",
    "* `align` : How to align the x, y coordinates respect to each pixel. Landsat Collection 3 should be centre aligned `align = (15, 15)` if data is loaded in its native UTM zone projection, e.g. `'EPSG:32756'` "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Enter specific dates\n",
    "time_list = (\n",
    "\"2022-05-13\",\n",
    "\"2022-05-13\"\n",
    "            )\n",
    "\n",
    "# Converts the list to a list of tuples (which the code needs)\n",
    "time_range = [(time_list[i],time_list[i]) for i in range(0,len(time_list))]\n",
    "\n",
    "\n",
    "# Extents of the images to download (note: Albers 3577). These vector files need to be first uploaded to the sandbox directory\n",
    "vector_file = './vectors/ningaloo_processing_aoi.shp'\n",
    "\n",
    "# give attribute name to id outputs\n",
    "attribute_col = 'aoi' \n",
    "\n",
    "# These are the Landsat 8 product descriptions, do not change\n",
    "products = ['ga_ls8c_ard_3'] \n",
    "\n",
    "measurements = ['nbart_blue', 'nbart_green', 'nbart_red', 'nbart_nir'] # choose only the bands you need to save space\n",
    "resolution = (-30, 30)\n",
    "output_crs = 'EPSG:3577'\n",
    "align = (0, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/jovyan/Bart/../Tools/dea_tools/plotting.py:418: FutureWarning: The `map_shapefile` function is deprecated, and will be removed from future versions of `dea-tools`. Please use Geopanda's built-in `.explore` functionality instead.\n",
      "  warnings.warn(\"The `map_shapefile` function is deprecated, and will \"\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ef62ef77f5164f5baaf574d3fe9d2fa9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Label(value='')"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "476789795d564de18e6470ca1f9b635c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(center=[-21.960797363881863, 113.93438124230774], controls=(ZoomControl(options=['position', 'zoom_in_text…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# read in original vectors\n",
    "gdf = gpd.read_file(vector_file)\n",
    "geom = geometry.Geometry(gdf.iloc[0].geometry.__geo_interface__, \n",
    "                          crs=geometry.CRS(output_crs))\n",
    "\n",
    "# un-comment if you wish to visualise the polygon\n",
    "map_shapefile(gdf, attribute=attribute_col)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('2022-05-13', '2022-05-13')\n",
      "Finding datasets\n",
      "    ga_ls8c_ard_3\n",
      "Loading 1 time steps\n",
      "\n",
      "Writing image: ...\n",
      "ningaloo_2022-05-13_L8_multiband.tif\n",
      "\n",
      "Next image...\n",
      "\n",
      "('2022-05-13', '2022-05-13')\n",
      "Finding datasets\n",
      "    ga_ls8c_ard_3\n",
      "Loading 1 time steps\n",
      "\n",
      "Writing image: ...\n",
      "ningaloo_2022-05-13_L8_multiband.tif\n",
      "\n",
      "Next image...\n",
      "\n",
      "Fin!\n"
     ]
    }
   ],
   "source": [
    "# iterating over tuple elements\n",
    "for time_step in time_range:\n",
    "    \n",
    "    print(time_step)\n",
    "      \n",
    "    # Create a datacube dictionary query object    \n",
    "    query = {'time': time_step,\n",
    "        'geopolygon': geom,\n",
    "        'measurements': measurements,\n",
    "        'resolution': resolution,\n",
    "        'output_crs': output_crs,\n",
    "        'align': align,\n",
    "        }\n",
    "    \n",
    "    try: # The try except clauses stop the script failing if the image does not excist in the sandbox but it does in glovis for example, if goes to the next image\n",
    "        # Load the image(s)\n",
    "        ds = load_ard(dc=dc, \n",
    "                       products=products,\n",
    "                       mask_pixel_quality=False,            \n",
    "                       group_by='solar_day',\n",
    "                       **query)\n",
    "        \n",
    " \n",
    "        ### Create output file variables and write to GeoTiff\n",
    "\n",
    "        # Get time and site variable for naming\n",
    "        time = ds.groupby('time')\n",
    "        site = gdf.loc[0 ,attribute_col]\n",
    "        \n",
    "        # Select a single time-slice of each new index band\n",
    "        ds = ds.isel(time=0).to_array()\n",
    "\n",
    "        # Write GeoTIFF\n",
    "        name = f'{site}_{str(time)[57:67]}_L8_multiband.tif'\n",
    "        \n",
    "        print(\"\\nWriting image: ...\")\n",
    "        write_cog(geo_im=ds,\n",
    "                  fname=name,\n",
    "                  overwrite=True) # outputs in ipynp dir\n",
    "        \n",
    "        print(name)\n",
    "        \n",
    "\n",
    "    except (ValueError, KeyError): # I do not know how to use error exception, but this prints the details of the image that the sandbox does not have.\n",
    "        \n",
    "        print(\"\\nThe image combination of below variables not accessable:\",\n",
    "              \"----------------------------------------------------------\",\n",
    "              \"Time_step: {}\".format(time_step),\n",
    "              \"Product: {}\".format(products),\n",
    "              \"ValueError: {}\".format(ValueError),\n",
    "              \"KeyError: {}\".format(KeyError),\n",
    "              sep='\\n')\n",
    "        \n",
    "        pass\n",
    "    finally:\n",
    "        print(\"\\nNext image...\\n\")\n",
    "        \n",
    "print(\"Fin!\")       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T51HUE_test_2020-05-29_L8_NBR_10m.tif\n",
      "1T51HUE_test_2020-05-29_L8_NBR_10m\n",
      "tar: Removing leading `/' from member names\n",
      "/home/jovyan/testing/\n",
      "/home/jovyan/testing/T51HUE_test_2020-05-13_L8_NBR_10m.tif\n",
      "/home/jovyan/testing/T51HUE_test_2020-05-13_L8_NDVI_10m.tif\n",
      "/home/jovyan/testing/T51HUE_test_2020-05-29_L8_NBR_10m.tif\n",
      "/home/jovyan/testing/SENTINEL_select_dates_index_download.ipynb\n",
      "/home/jovyan/testing/GL_Landsat_scene_download_20210701_Loop_working.ipynb\n",
      "/home/jovyan/testing/vectors/\n",
      "/home/jovyan/testing/vectors/AOI_alb.prj\n",
      "/home/jovyan/testing/vectors/AOI_alb.shp\n",
      "/home/jovyan/testing/vectors/AOI_alb.shx\n",
      "/home/jovyan/testing/vectors/.ipynb_checkpoints/\n",
      "/home/jovyan/testing/vectors/AOI_alb.dbf\n",
      "/home/jovyan/testing/Landsat_select_dates_index_download-Copy1.ipynb\n",
      "/home/jovyan/testing/.ipynb_checkpoints/\n",
      "/home/jovyan/testing/.ipynb_checkpoints/Landsat_select_dates_index_download-Copy1-checkpoint.ipynb\n",
      "/home/jovyan/testing/.ipynb_checkpoints/SENTINEL_select_dates_index_download-checkpoint.ipynb\n",
      "/home/jovyan/testing/.ipynb_checkpoints/GL_Landsat_scene_download_20210701_Loop_working-checkpoint.ipynb\n",
      "/home/jovyan/testing/T51HUE_test_2020-05-29_L8_NDVI_10m.tif\n",
      "tar: /home/jovyan/testing: file changed as we read it\n"
     ]
    }
   ],
   "source": [
    "# Tarball - CHANGE directory name to yours\n",
    "# remove block quotes (''') from beginning and end and 'play' like all cells above if they are there\n",
    "# Tar ball the whole directory with a file name variable of the name of the last file created\n",
    "\n",
    "# print(name2)\n",
    "# tarName = \"1\" + name2[:-4] # in the shell the first char is drop off!!!\n",
    "# print(tarName)\n",
    "# ! tar -czvf ${tarName}.tar.gz '/home/jovyan/testing/'"
   ]
  },
  {
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}