00_data_downloads.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8bce0490-1004-4cee-8d65-ca565e5fc8ea",
   "metadata": {},
   "source": [
    "# Data downloads\n",
    "\n",
    "Download required data not included in the repo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "d7de0c34-22c3-4e7d-b2b8-192c748ee6cb",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Note: you may need to restart the kernel to use updated packages.\n"
     ]
    }
   ],
   "source": [
    "#pip -q install earthaccess"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d6d801a2-ab9f-4736-9520-819aa014801a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import requests\n",
    "import zipfile\n",
    "import geopandas as gpd\n",
    "import earthaccess\n",
    "from osgeo import gdal\n",
    "import rioxarray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "1a54b1df-266d-4647-b5fb-ffbc639f0993",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "data_dir = \"data/external/\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf156e2d-55e4-456a-a604-c3f1b311f61c",
   "metadata": {},
   "source": [
    "# GRWL\n",
    "Global river widths from Landsat. Vector file with rivers visible from Landsat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "3c941e5a-0ef2-4104-a514-1d79af625e2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# GRWL summary stats\n",
    "# https://zenodo.org/records/1297434\n",
    "url = 'https://zenodo.org/records/1297434/files/GRWL_summaryStats_V01.01.zip?download=1'\n",
    "out_fn = data_dir + 'GRWL_summaryStats_V01.01.zip'\n",
    "response = requests.get(url)\n",
    "if response.status_code == 200:\n",
    "    open(out_fn, 'wb').write(response.content)\n",
    "    \n",
    "# unzip\n",
    "with zipfile.ZipFile(out_fn, \"r\") as zip_ref:\n",
    "    zip_ref.extractall(data_dir + 'GRWL_summaryStats_V01')\n",
    "\n",
    "# remove zipfile\n",
    "os.remove(out_fn)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c422cd72-7f40-412d-ac2b-7674ecbe2349",
   "metadata": {},
   "source": [
    "# HydroRIVERS\n",
    "River reach vector file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "80afa37f-946d-4a4f-86dd-6788553c4925",
   "metadata": {},
   "outputs": [],
   "source": [
    "# HydroRIVERS\n",
    "# https://www.hydrosheds.org/products/hydrorivers\n",
    "url = 'https://data.hydrosheds.org/file/HydroRIVERS/HydroRIVERS_v10_as_shp.zip'\n",
    "out_fn = data_dir + 'HydroRIVERS_v10_as_shp.zip'\n",
    "response = requests.get(url)\n",
    "if response.status_code == 200:\n",
    "    open(out_fn, 'wb').write(response.content)\n",
    "    \n",
    "# unzip\n",
    "with zipfile.ZipFile(out_fn, \"r\") as zip_ref:\n",
    "    zip_ref.extractall(data_dir) # it will unzip into a directory\n",
    "\n",
    "# remove zipfile\n",
    "os.remove(out_fn)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2904bad5-b674-4275-a88d-adc1ca3da424",
   "metadata": {},
   "source": [
    "# Country boundaries\n",
    "For all of Asia"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "026c8fd8-7689-4f5c-8b3e-48230251d5b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# asia polygons (for visualization)\n",
    "# in this case, just load directly into memory and only \n",
    "# save a simplified version\n",
    "url = 'https://stacks.stanford.edu/file/druid:yg089df0008/data.zip'\n",
    "out_fn = data_dir + 'asia_polygons.gpkg'\n",
    "\n",
    "asia_polygons = gpd.read_file(url)\n",
    "asia = asia_polygons.to_crs('EPSG:32648')\n",
    "asia['geometry'] = asia.simplify(1000)\n",
    "asia = asia.to_crs(asia_polygons.crs)\n",
    "asia = asia[asia['REGION'] == 'ASIA']\n",
    "asia.to_file(out_fn)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "345b50f5-f0a5-41bd-bd74-35154804da52",
   "metadata": {},
   "source": [
    "# Water classification\n",
    "Water classification raster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "bd0f36c0-439d-4398-b074-026cd679b29f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# JRC yearly water classification\n",
    "url = 'https://storage.googleapis.com/global-surface-water/downloads2021/occurrence/occurrence_100E_20Nv1_4_2021.tif'\n",
    "out_fn = data_dir + 'occurrence_100E_20Nv1_4_2021.tif'\n",
    "response = requests.get(url)\n",
    "if response.status_code == 200:\n",
    "            #Write to disk\n",
    "            open(out_fn, 'wb').write(response.content)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8abeb89-2f53-4e18-a33d-1f1bbf6148c6",
   "metadata": {},
   "source": [
    "# Tonle Sap Lake\n",
    "Shapefile of TSL for plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "aec92826-d917-42e4-a845-8f9774135f39",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Tonle Sap Lake\n",
    "# from https://data.opendevelopmentmekong.net/dataset/water-bodies-in-cambodia\n",
    "url = 'https://data.opendevelopmentcambodia.net/en/dataset/3c74812b-0a78-4407-a072-6f91c58de45b/resource/a76793df-76aa-469a-97bc-ae3a04b6faad/download/water-bodies.geojson'\n",
    "out_fn = 'data/' + 'tonle_sap_lake.geojson'\n",
    "lakes = gpd.read_file(url)\n",
    "tsl = lakes.loc[lakes['name'] == 'Boeung Tonle Sap']\n",
    "tsl.to_file(out_fn)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "392a3d94-68f9-4d8f-bbaf-7c0fe51060f6",
   "metadata": {},
   "source": [
    "# Country boundaries\n",
    "Shapefile with country boundaries for plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "id": "cea3b292-5b0d-4416-b393-8365cdc53d64",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# LMB (just remove China from Mekong shapefile)\n",
    "countries = gpd.read_file(data_dir + 'asia_polygons.gpkg')\n",
    "lmb_countries = [\n",
    "    'VIETNAM',\n",
    "    'THAILAND',\n",
    "    'CAMBODIA', \n",
    "    'LAOS',\n",
    "    'BURMA'\n",
    "]\n",
    "mekong = gpd.read_file('data/mekong.geojson')\n",
    "lmb = mekong.clip(countries[countries['NAME'].isin(lmb_countries)])\n",
    "lmb.to_file('data/lmb.geojson')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e14e973-fafe-488d-9e35-a57f9c1ecad7",
   "metadata": {},
   "source": [
    "## Load HLS dataset over the extent needed\n",
    "Use this as a template for setting up masks etc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "7f2078c9-7d15-45f3-9307-2874cb16a1bf",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# open HLS dataset using rioxarray\n",
    "# EarthAccess setup for HLS\n",
    "# https://github.com/nasa/HLS-Data-Resources/blob/main/python/tutorials/HLS_Tutorial.ipynb\n",
    "\n",
    "earthaccess.login(persist=True)\n",
    "\n",
    "# GDAL configurations used to successfully access LP DAAC Cloud Assets via vsicurl \n",
    "gdal.SetConfigOption('GDAL_HTTP_COOKIEFILE','~/cookies.txt')\n",
    "gdal.SetConfigOption('GDAL_HTTP_COOKIEJAR', '~/cookies.txt')\n",
    "gdal.SetConfigOption('GDAL_DISABLE_READDIR_ON_OPEN','EMPTY_DIR')\n",
    "gdal.SetConfigOption('CPL_VSIL_CURL_ALLOWED_EXTENSIONS','TIF')\n",
    "gdal.SetConfigOption('GDAL_HTTP_UNSAFESSL', 'YES')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "37921097-eec8-42ea-8e10-a99e7c4f0956",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# hls url\n",
    "# see 00-1_pull_hls_data to see where this url comes from\n",
    "url_pxa = 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSS30.020/HLS.S30.T48PXA.2023325T032029.v2.0/HLS.S30.T48PXA.2023325T032029.v2.0.B04.tif'\n",
    "url_pwv = 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSS30.020/HLS.S30.T48PWV.2023325T032029.v2.0/HLS.S30.T48PWV.2023325T032029.v2.0.B04.tif'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "3f3ebb10-59e2-411e-b50f-bb9eb33c70d9",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The COGs have been loaded into memory!\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# Use vsicurl to load the data directly into memory (be patient, may take a few seconds)\n",
    "chunk_size = dict(band=1, x=512, y=512) # Tiles have 1 band and are divided into 512x512 pixel chunks\n",
    "# Sometimes a vsi curl error occurs so we need to retry if it does\n",
    "max_retries = 10\n",
    "\n",
    "# Load PXA\n",
    "for _i in range(max_retries):\n",
    "    try:\n",
    "        # pull hls (red band)\n",
    "        hls_red_pxa = rioxarray.open_rasterio(url_pxa, chunks=chunk_size, masked=True).squeeze('band')\n",
    "        hls_red_pxa.attrs['scale_factor'] = 0.0001 # hard coded the scale_factor attribute\n",
    "        break # Break out of the retry loop\n",
    "    except Exception as ex:\n",
    "        print(f\"vsi curl error: {ex}. Retrying...\")\n",
    "else:\n",
    "    print(f\"Failed to process {url_pxa} after {max_retries} retries. Please check to see you're authenticated with earthaccess.\")\n",
    "    \n",
    "# Load PWV \n",
    "for _i in range(max_retries):\n",
    "    try:\n",
    "        # pull hls (red band)\n",
    "        hls_red_pwv = rioxarray.open_rasterio(url_pwv, chunks=chunk_size, masked=True).squeeze('band')\n",
    "        hls_red_pwv.attrs['scale_factor'] = 0.0001 # hard coded the scale_factor attribute\n",
    "        break # Break out of the retry loop\n",
    "    except Exception as ex:\n",
    "        print(f\"vsi curl error: {ex}. Retrying...\")\n",
    "else:\n",
    "    print(f\"Failed to process {url_pwv} after {max_retries} retries. Please check to see you're authenticated with earthaccess.\")\n",
    "print(\"The COGs have been loaded into memory!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "8eac8a53-9507-4428-a175-6114847a3ef2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Create a window for HLS data\n",
    "# Use utm rivers because that's what HLS is in\n",
    "window_bounds = (560172.6343010075, 1450022.6611399346, 654025.1724448078, 1542844.933869364)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "8956b420-8366-440a-a211-3daf341237ab",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# merge the two images\n",
    "# This takes kind of a long time\n",
    "# but it creates a dataset that spans all the sample sites. \n",
    "hls_fn = data_dir + 'hls_region.nc'\n",
    "if not os.path.exists(hls_fn):\n",
    "    import rioxarray.merge\n",
    "    # need to pass in list of xarray.DataSet so cast DataArrays to DataSets\n",
    "    hls_merged = rioxarray.merge.merge_datasets([hls_red_pwv.to_dataset(name='red'), \n",
    "                                                 hls_red_pxa.to_dataset(name='red')], \n",
    "                                                bounds = window_bounds) \n",
    "    hls_merged.to_netcdf(hls_fn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4dc72a1e-8a53-4cc1-9098-d353815ea6f6",
   "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}