step3_get_10m_dems.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import geojson\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as colors\n",
    "import seaborn as sns\n",
    "import geopandas as gpd\n",
    "from shapely import geometry, ops\n",
    "import rasterio\n",
    "from rasterio import features\n",
    "import geemap\n",
    "import ee\n",
    "\n",
    "from pathlib import Path\n",
    "\n",
    "import csv\n",
    "\n",
    "import requests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ee.Authenticate()\n",
    "ee.Initialize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shed_list = pd.read_csv(\"shed_list_for_dd_new.csv\", dtype = {'HYBAS_ID':'int', 'DD':'float64'},\n",
    "                        #  index_col='HYBAS_ID'\n",
    "                         )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "big_file_list = []\n",
    "\n",
    "for row in range(len(shed_list)):\n",
    "\n",
    "  HYBAS_ID = int(shed_list.iloc[row]['HYBAS_ID']) # Needs to be int here\n",
    "  # print(f'Now looking at {HYBAS_ID}')\n",
    "  \n",
    "  hybas_path = Path(f'./lsdtt/{HYBAS_ID}/')\n",
    "\n",
    "  if (hybas_path / f'{HYBAS_ID}.tif').exists() == True: # In case the script breaks while running\n",
    "    # print(f'{HYBAS_ID}.tif already exists')\n",
    "    continue\n",
    "  else:\n",
    "    \n",
    "    shed = ee.FeatureCollection(\"WWF/HydroSHEDS/v1/Basins/hybas_10\").filter(ee.Filter.eq('HYBAS_ID', HYBAS_ID))\n",
    "\n",
    "    long = shed.geometry().centroid().coordinates().get(0).getInfo()\n",
    "    lat = shed.geometry().centroid().coordinates().get(1).getInfo()\n",
    "    epsg = int(32700-(np.round((45+lat)/90,0)*100)+np.round((183+long)/6,0))\n",
    "\n",
    "    # if epsg < 32610:\n",
    "    #   print(\"Whoops this is a single digit UTM zone so it'll break lsdtt, skipping\")\n",
    "    #   continue\n",
    "\n",
    "    # else:\n",
    "    if lat < 59.0:\n",
    "      elevation_0 = ee.Image('USGS/3DEP/10m').select('elevation').clip(shed)\n",
    "      elevation = elevation_0.resample('bilinear').reproject(crs=f'EPSG:{epsg}', scale=10) # It is not clear why I have to do this again, but the sheds weren't clipping otherwise\n",
    "    else:\n",
    "      elevation_0 = ee.Image('UMN/PGC/ArcticDEM/V3/2m_mosaic').select('elevation').clip(shed)\n",
    "      elevation = elevation_0.resample('bilinear').reproject(crs=f'EPSG:{epsg}', scale=10)\n",
    "      # elevation = elevation_0.resample('bilinear').reproject(crs='EPSG:4326', scale=10)\n",
    "\n",
    "    try: \n",
    "      url = elevation.getDownloadUrl({\n",
    "          'region': shed.geometry(),\n",
    "          'scale': 10,\n",
    "          'crs': f'EPSG:{epsg}',\n",
    "          # 'crs': 'EPSG:4326',\n",
    "          'maxPixels': 10063164208,\n",
    "          'format': 'GEO_TIFF'\n",
    "      }) \n",
    "      print(f'Downloading {HYBAS_ID}.tif')\n",
    "      response = requests.get(url)\n",
    "      Path(f'./lsdtt/{HYBAS_ID}/').mkdir(parents=True, exist_ok=True)\n",
    "      with open(hybas_path / f'{HYBAS_ID}.tif', 'wb') as fd:\n",
    "        fd.write(response.content)\n",
    "    except KeyboardInterrupt:\n",
    "        print('Interrupted')\n",
    "        break\n",
    "    except Exception: # I can't get the EEException to work?? so now it's generic error which is not great\n",
    "      # print(f'{HYBAS_ID} is too big! :(')\n",
    "      big_file_list.append([f'{HYBAS_ID}'])\n",
    "      # print(big_file_list)\n",
    "\n",
    "# And now record the too big sheds\n",
    "\n",
    "fields = ['HYBAS_ID']\n",
    "\n",
    "with open('big_sheds.csv', 'w') as f:\n",
    "  # using csv.writer method from CSV package\n",
    "  write = csv.writer(f)\n",
    "    \n",
    "  write.writerow(fields)\n",
    "  write.writerows(big_file_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = Path('./lsdtt/')\n",
    "downloaded = pd.DataFrame({'HYBAS_ID':[f.name for f in p.iterdir() if f.is_dir()]}).astype(int)\n",
    "shed_list = pd.read_csv(\"shed_list_for_dd_new.csv\", dtype = {'HYBAS_ID':'int', 'DD':'float64'},\n",
    "                        #  index_col='HYBAS_ID'\n",
    "                         )\n",
    "merged = shed_list.merge(downloaded, on=\"HYBAS_ID\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "merged['EXTENT'].value_counts()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "gee",
   "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.9.13"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}