PC-RTC_download_dev

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "15201120-9441-448a-94b2-ee4e97468e65",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# based on exmaples from\n",
    "# https://planetarycomputer.microsoft.com/docs/tutorials/cloudless-mosaic-sentinel2/\n",
    "# https://planetarycomputer.microsoft.com/dataset/sentinel-2-l2a#Example-Notebook\n",
    "from pystac.extensions.eo import EOExtension as eo\n",
    "import pystac_client\n",
    "import planetary_computer\n",
    "import glob\n",
    "import rioxarray as rxr\n",
    "import re\n",
    "import datetime\n",
    "import pandas as pd\n",
    "from shapely.geometry import box\n",
    "import odc.stac\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import xarray as xr\n",
    "import rasterio as rio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9fb141b-afac-4839-a1bc-3e055e86e164",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# def rtc_for_aso_snowon(aso_raster_fn):\n",
    "#     time = pd.to_datetime(re.search(\"(\\d{4}\\d{2}\\d{2})\", aso_raster_fn).group())\n",
    "#     week_before = (time - datetime.timedelta(weeks=1)).strftime('%Y-%m-%d')\n",
    "#     week_after = (time + datetime.timedelta(weeks=1)).strftime('%Y-%m-%d')\n",
    "#     time_of_interest = f'{week_before}/{week_after}'\n",
    "    \n",
    "#     aso_raster = rxr.open_rasterio(aso_raster_fn).squeeze()\n",
    "#     aso_raster = aso_raster.where(aso_raster>=0, drop=True)\n",
    "#     aso_raster = aso_raster.interpolate_na(dim='x')\n",
    "#     bounds_latlon = box(*aso_raster.rio.transform_bounds(\"EPSG:4326\"))\n",
    "    \n",
    "#     catalog = pystac_client.Client.open(\n",
    "#     \"https://planetarycomputer.microsoft.com/api/stac/v1\",\n",
    "#     modifier=planetary_computer.sign_inplace)\n",
    "\n",
    "#     search = catalog.search(\n",
    "#         collections=[\"sentinel-1-rtc\"],\n",
    "#         intersects=bounds_latlon,\n",
    "#         datetime=time_of_interest)\n",
    "    \n",
    "#     # Check how many items were returned\n",
    "#     items = search.item_collection()\n",
    "    \n",
    "#     rtc_stac = odc.stac.load(items,chunks={\"x\": 2048, \"y\": 2048},resolution=50, groupby='sat:absolute_orbit')\n",
    "#     print(f\"Returned {len(rtc_stac.time)} acquisitions\")\n",
    "#     rtc_stac_clipped = rtc_stac.rio.clip_box(*bounds_latlon.bounds,crs=\"EPSG:4326\")\n",
    "    \n",
    "#     rel_orbits = [scene.properties['sat:relative_orbit'] for scene in items.items]\n",
    "#     ac_times = [scene.properties['datetime'] for scene in items.items]\n",
    "#     ac_times = [np.datetime64(item) for item in ac_times]\n",
    "    \n",
    "#     # clip to ASO extent\n",
    "#     rtc_stac_clipped = rtc_stac_clipped.rio.reproject_match(aso_raster, resampling=rio.enums.Resampling.bilinear)\n",
    "\n",
    "#     # limit to morning acquisitions\n",
    "#     rtc_ds = rtc_stac_clipped.where(rtc_stac_clipped.time.dt.hour > 11, drop=True)\n",
    "#     if 'vv' not in list(rtc_ds.keys()) or 'vh' not in list(rtc_ds.keys()):\n",
    "#         print('missing polarization')\n",
    "#         return None\n",
    "    \n",
    "#     if len(rtc_ds.time) == 0:\n",
    "#         print('no morning acquisitions')\n",
    "#         return None\n",
    "        \n",
    "#     # calculate percent vh coverage of each acquisition\n",
    "#     perc_cover = (rtc_ds.vh > 0).sum(dim=['x', 'y'])/(rtc_ds.vh >= -1000000000).sum(dim=['x', 'y'])\n",
    "    \n",
    "#     # if multiple with full coverage, grab nearest in time with full coverage\n",
    "#     if perc_cover.values.tolist().count(1) > 1:\n",
    "#         print('total snow-on coverage available')\n",
    "#         rtc_ds = rtc_ds.where(perc_cover == 1, drop=True).sortby('time')\n",
    "#         rtc_ds = rtc_ds.sel(time=time, method='nearest')\n",
    "\n",
    "#     # exit if no rasters have good vh coverage\n",
    "#     elif perc_cover.max() < 0.1:\n",
    "#         print('max vh coverage is < 10 percent--recommend skipping ASO raster')\n",
    "#         return None\n",
    "\n",
    "#     # otherwise, grab max coverage \n",
    "#     else:\n",
    "#         if perc_cover.max() == 1:\n",
    "#             print('total snow-on coverage available')\n",
    "#         else: \n",
    "#             print(f'{perc_cover.max().item()} snow-on coverage')\n",
    "#         rtc_ds = rtc_ds.sel(time=perc_cover.idxmax())\n",
    "        \n",
    "#     # mask negative areas\n",
    "#     rtc_ds = rtc_ds.where(rtc_ds.vh > 0, drop=True)\n",
    "\n",
    "#     # get relative orbit of scene\n",
    "#     rel_orbit = rel_orbits[ac_times.index(rtc_ds.time)]\n",
    "    \n",
    "#     #rtc_ds.to_netcdf(f'../data/S1_rtc/S1_snow-on_{rtc_ds.time.dt.strftime(\"%Y%m%d\").item()}_for_{aso_raster_fn.split(\"/\")[-1][:-4]}.nc')\n",
    "    \n",
    "#     return rel_orbit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f6332bf-c57e-483e-80c0-7854ec0804c9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# def rtc_for_aso_snowoff(aso_raster_fn, orbit):\n",
    "#     year = pd.to_datetime(re.search(\"(\\d{4}\\d{2}\\d{2})\", aso_raster_fn).group()).year\n",
    "#     time = pd.to_datetime(f'{year-1}0910')\n",
    "#     week_before = (time - datetime.timedelta(weeks=2)).strftime('%Y-%m-%d')\n",
    "#     week_after = (time + datetime.timedelta(weeks=2)).strftime('%Y-%m-%d')\n",
    "#     time_of_interest = f'{week_before}/{week_after}'\n",
    "\n",
    "#     aso_raster = rxr.open_rasterio(aso_raster_fn).squeeze()\n",
    "#     aso_raster = aso_raster.where(aso_raster>=0, drop=True)\n",
    "#     aso_raster = aso_raster.interpolate_na(dim='x')\n",
    "#     bounds_latlon = box(*aso_raster.rio.transform_bounds(\"EPSG:4326\"))\n",
    "\n",
    "#     catalog = pystac_client.Client.open(\n",
    "#     \"https://planetarycomputer.microsoft.com/api/stac/v1\",\n",
    "#     modifier=planetary_computer.sign_inplace)\n",
    "\n",
    "#     search = catalog.search(\n",
    "#         collections=[\"sentinel-1-rtc\"],\n",
    "#         intersects=bounds_latlon,\n",
    "#         datetime=time_of_interest)\n",
    "\n",
    "#     # Check how many items were returned\n",
    "#     items = search.item_collection()\n",
    "\n",
    "#     rel_orbits = [scene.properties['sat:relative_orbit'] for scene in items.items]\n",
    "#     ac_times = [scene.properties['datetime'] for scene in items.items]\n",
    "#     ac_times = [np.datetime64(item) for item in ac_times]\n",
    "\n",
    "#     rtc_stac = odc.stac.load(items,chunks={\"x\": 2048, \"y\": 2048},resolution=50, groupby='sat:absolute_orbit')\n",
    "#     print(f\"Returned {len(rtc_stac.time)} acquisitions\")\n",
    "#     rtc_stac_clipped = rtc_stac.rio.clip_box(*bounds_latlon.bounds,crs=\"EPSG:4326\")\n",
    "\n",
    "#     orbit_dict = {}\n",
    "#     for i, orbit in enumerate(rel_orbits):\n",
    "#         if orbit not in orbit_dict.keys():\n",
    "#             orbit_dict[orbit] = [ac_times[i]]\n",
    "#         else:\n",
    "#             orbit_dict[orbit].append(ac_times[i])\n",
    "\n",
    "#     if rel_orbit not in orbit_dict.keys():\n",
    "#         print('no acquisitons from same orbit, skipping')\n",
    "#         return\n",
    "\n",
    "#     rtc_stac_clipped = rtc_stac_clipped.where(rtc_stac_clipped.time.isin(orbit_dict[rel_orbit]), drop=True)\n",
    "\n",
    "#     # clip to ASO extent\n",
    "#     rtc_stac_clipped = rtc_stac_clipped.rio.reproject_match(aso_raster, resampling=rio.enums.Resampling.bilinear)\n",
    "\n",
    "#     # limit to morning acquisitions\n",
    "#     rtc_ds = rtc_stac_clipped.where(rtc_stac_clipped.time.dt.hour > 11, drop=True)\n",
    "#     if 'vv' not in list(rtc_ds.keys()) or 'vh' not in list(rtc_ds.keys()):\n",
    "#             print('missing polarization, skipping')\n",
    "#             return\n",
    "\n",
    "\n",
    "#     # calculate percent vh coverage of each acquisition\n",
    "#     perc_cover = (rtc_ds.vh > 0).sum(dim=['x', 'y'])/(rtc_ds.vh >= -1000000000).sum(dim=['x', 'y'])\n",
    "\n",
    "#     # if multiple with full coverage, grab nearest in time with full coverage\n",
    "#     if perc_cover.values.tolist().count(1) > 1:\n",
    "#         print('total snow-off coverage available')\n",
    "#         rtc_ds = rtc_ds.where(perc_cover == 1, drop=True).sortby('time')\n",
    "#         rtc_ds = rtc_ds.sel(time=time, method='nearest')\n",
    "\n",
    "#     # exit if no rasters have good vh coverage\n",
    "#     elif perc_cover.max() < 0.1:\n",
    "#         print('max vh coverage is < 10 percent--recommend skipping ASO raster')\n",
    "#         return\n",
    "\n",
    "#     # otherwise, grab max coverage \n",
    "#     else:\n",
    "#         if perc_cover.max() == 1:\n",
    "#             print('total snow-off coverage available')\n",
    "#         else: \n",
    "#             print(f'{perc_cover.max().item()} snow-off coverage')\n",
    "#         rtc_ds = rtc_ds.sel(time=perc_cover.idxmax())\n",
    "\n",
    "#     # mask negative areas\n",
    "#     rtc_ds = rtc_ds.where(rtc_ds.vh > 0, drop=True)\n",
    "    \n",
    "#     rtc_ds.to_netcdf(f'../data/S1_rtc/S1_snow-off_{rtc_ds.time.dt.strftime(\"%Y%m%d\").item()}_for_{aso_raster_fn.split(\"/\")[-1][:-4]}.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "8b30962e-b6d4-4157-985f-db8a2581fa19",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def rtc_for_aso_snowon_mean(aso_raster_fn):\n",
    "    time = pd.to_datetime(re.search(\"(\\d{4}\\d{2}\\d{2})\", aso_raster_fn).group())\n",
    "    week_before = (time - datetime.timedelta(weeks=2)).strftime('%Y-%m-%d')\n",
    "    week_after = (time + datetime.timedelta(weeks=2)).strftime('%Y-%m-%d')\n",
    "    time_of_interest = f'{week_before}/{week_after}'\n",
    "\n",
    "    aso_raster = rxr.open_rasterio(aso_raster_fn).squeeze()\n",
    "    aso_raster = aso_raster.where(aso_raster>=0, drop=True)\n",
    "    aso_raster = aso_raster.interpolate_na(dim='x')\n",
    "    bounds_latlon = box(*aso_raster.rio.transform_bounds(\"EPSG:4326\"))\n",
    "\n",
    "    catalog = pystac_client.Client.open(\n",
    "    \"https://planetarycomputer.microsoft.com/api/stac/v1\",\n",
    "    modifier=planetary_computer.sign_inplace)\n",
    "\n",
    "    search = catalog.search(\n",
    "        collections=[\"sentinel-1-rtc\"],\n",
    "        intersects=bounds_latlon,\n",
    "        datetime=time_of_interest)\n",
    "\n",
    "    # Check how many items were returned\n",
    "    items = search.item_collection()\n",
    "\n",
    "    rtc_stac = odc.stac.load(items,chunks={\"x\": 2048, \"y\": 2048},resolution=50, groupby='sat:absolute_orbit')\n",
    "    print(f\"Returned {len(rtc_stac.time)} acquisitions\")\n",
    "    rtc_stac_clipped = rtc_stac.rio.clip_box(*bounds_latlon.bounds,crs=\"EPSG:4326\")\n",
    "\n",
    "    rel_orbits = [scene.properties['sat:relative_orbit'] for scene in items.items]\n",
    "    ac_times = [scene.properties['datetime'] for scene in items.items]\n",
    "    ac_times = [np.datetime64(item) for item in ac_times]\n",
    "\n",
    "    # clip to ASO extent\n",
    "    rtc_stac_clipped = rtc_stac_clipped.rio.reproject_match(aso_raster, resampling=rio.enums.Resampling.bilinear)\n",
    "\n",
    "    # limit to morning acquisitions\n",
    "    rtc_ds = rtc_stac_clipped.where(rtc_stac_clipped.time.dt.hour > 11, drop=True)\n",
    "    if 'vv' not in list(rtc_ds.keys()) or 'vh' not in list(rtc_ds.keys()):\n",
    "        print('missing polarization')\n",
    "        return None\n",
    "\n",
    "    if len(rtc_ds.time) == 0:\n",
    "        print('no morning acquisitions')\n",
    "        return None\n",
    "\n",
    "\n",
    "    # calculate percent vh coverage of each acquisition\n",
    "    perc_cover = (rtc_ds.vh > 0).sum(dim=['x', 'y'])/(rtc_ds.vh >= -1000000000).sum(dim=['x', 'y'])\n",
    "\n",
    "    # if multiple with full coverage, grab nearest in time with full coverage\n",
    "    if perc_cover.values.tolist().count(1) > 1:\n",
    "        print('total snow-on coverage available')\n",
    "        rtc_ds = rtc_ds.where(perc_cover == 1, drop=True).sortby('time')\n",
    "        rtc_ds = rtc_ds.sel(time=time, method='nearest')\n",
    "\n",
    "    # exit if no rasters have good vh coverage\n",
    "    elif perc_cover.max() < 0.01:\n",
    "        print('max vh coverage is < 1%--recommend skipping ASO raster')\n",
    "\n",
    "    # otherwise, grab max coverage \n",
    "    else:\n",
    "        if perc_cover.max() == 1:\n",
    "            print('total snow-on coverage available')\n",
    "        else: \n",
    "            print(f'{perc_cover.max().item()} snow-on coverage')\n",
    "        rtc_ds = rtc_ds.sel(time=perc_cover.idxmax())\n",
    "\n",
    "    # get relative orbit of scene\n",
    "    rel_orbit = rel_orbits[ac_times.index(rtc_ds.time)]\n",
    "\n",
    "    orbit_dict = {}\n",
    "    for i, orbit in enumerate(rel_orbits):\n",
    "        if orbit not in orbit_dict.keys():\n",
    "            orbit_dict[orbit] = [ac_times[i]]\n",
    "        else:\n",
    "            orbit_dict[orbit].append(ac_times[i])\n",
    "\n",
    "    rtc_ds = rtc_stac_clipped.where(rtc_stac_clipped.time.isin(orbit_dict[rel_orbit]), drop=True)\n",
    "\n",
    "    # take mean of all acquisitions\n",
    "    print(f'taking mean of {rtc_ds.time.size} snow-on rasters')\n",
    "    rtc_ds = rtc_ds.mean(dim='time', skipna=True)\n",
    "\n",
    "    rtc_ds = rtc_ds.compute()\n",
    "\n",
    "    # mask negative areas\n",
    "    rtc_ds = rtc_ds.where(rtc_ds.vh > 0, drop=True)\n",
    "    rtc_ds = rtc_ds.where(rtc_ds.vv > 0, drop=True)\n",
    "\n",
    "    #rtc_ds.to_netcdf(f'../data/S1_rtc_mean/S1_snow-on_orbit{rel_orbit}_for_{aso_raster_fn.split(\"/\")[-1][:-4]}.nc')\n",
    "    \n",
    "    return rel_orbit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "3873b05e-1a5e-44c5-8d60-7fb3bbfcf2f4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "aso_raster_fn = '/home/jovyan/deep-snow/data/ASO/ASO_50m_SD_cleaned/utm10n/ASO_50M_SD_SacramentoMcCloud_20230209_clean.tif'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "4ef240af-38c0-4620-bb83-c73c34cb44f1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def rtc_for_aso_snowoff_mean(aso_raster_fn, rel_orbit):\n",
    "    year = pd.to_datetime(re.search(\"(\\d{4}\\d{2}\\d{2})\", aso_raster_fn).group()).year\n",
    "    time = pd.to_datetime(f'{year-1}0910')\n",
    "    week_before = (time - datetime.timedelta(weeks=2)).strftime('%Y-%m-%d')\n",
    "    week_after = (time + datetime.timedelta(weeks=2)).strftime('%Y-%m-%d')\n",
    "    time_of_interest = f'{week_before}/{week_after}'\n",
    "\n",
    "    aso_raster = rxr.open_rasterio(aso_raster_fn).squeeze()\n",
    "    aso_raster = aso_raster.where(aso_raster>=0, drop=True)\n",
    "    aso_raster = aso_raster.interpolate_na(dim='x')\n",
    "    bounds_latlon = box(*aso_raster.rio.transform_bounds(\"EPSG:4326\"))\n",
    "\n",
    "    catalog = pystac_client.Client.open(\n",
    "    \"https://planetarycomputer.microsoft.com/api/stac/v1\",\n",
    "    modifier=planetary_computer.sign_inplace)\n",
    "\n",
    "    search = catalog.search(\n",
    "        collections=[\"sentinel-1-rtc\"],\n",
    "        intersects=bounds_latlon,\n",
    "        datetime=time_of_interest)\n",
    "\n",
    "    # Check how many items were returned\n",
    "    items = search.item_collection()\n",
    "\n",
    "    rel_orbits = [scene.properties['sat:relative_orbit'] for scene in items.items]\n",
    "    ac_times = [scene.properties['datetime'] for scene in items.items]\n",
    "    ac_times = [np.datetime64(item) for item in ac_times]\n",
    "\n",
    "    rtc_stac = odc.stac.load(items,chunks={\"x\": 2048, \"y\": 2048},resolution=50, groupby='sat:absolute_orbit')\n",
    "    print(f\"Returned {len(rtc_stac.time)} acquisitions\")\n",
    "    rtc_stac_clipped = rtc_stac.rio.clip_box(*bounds_latlon.bounds,crs=\"EPSG:4326\")\n",
    "\n",
    "    orbit_dict = {}\n",
    "    for i, orbit in enumerate(rel_orbits):\n",
    "        if orbit not in orbit_dict.keys():\n",
    "            orbit_dict[orbit] = [ac_times[i]]\n",
    "        else:\n",
    "            orbit_dict[orbit].append(ac_times[i])\n",
    "\n",
    "    if rel_orbit not in orbit_dict.keys():\n",
    "        print('no acquisitons from same orbit, skipping')\n",
    "        return\n",
    "\n",
    "    rtc_stac_clipped = rtc_stac_clipped.where(rtc_stac_clipped.time.isin(orbit_dict[rel_orbit]), drop=True)\n",
    "\n",
    "    # clip to ASO extent\n",
    "    rtc_ds = rtc_stac_clipped.rio.reproject_match(aso_raster, resampling=rio.enums.Resampling.bilinear)\n",
    "\n",
    "    if 'vv' not in list(rtc_ds.keys()) or 'vh' not in list(rtc_ds.keys()):\n",
    "            print('missing polarization, skipping')\n",
    "            return\n",
    "    \n",
    "    if len(rtc_ds.time) == 0:\n",
    "        print('no morning acquisitions')\n",
    "        return None\n",
    "\n",
    "    # calculate percent vh coverage of each acquisition\n",
    "    perc_cover = (rtc_ds.vh > 0).sum(dim=['x', 'y'])/(rtc_ds.vh >= -1000000000).sum(dim=['x', 'y'])\n",
    "\n",
    "    # exit if no rasters have good vh coverage\n",
    "    if perc_cover.max() < 0.01:\n",
    "        print('max vh coverage is < 1%--recommend skipping ASO raster')\n",
    "        return\n",
    "    \n",
    "    # take mean of all acquisitions\n",
    "    print(f'taking mean of {rtc_ds.time.size} snow-off rasters')\n",
    "    rtc_ds = rtc_ds.mean(dim='time', skipna=True)\n",
    "\n",
    "    rtc_ds = rtc_ds.compute()\n",
    "\n",
    "    # mask negative areas\n",
    "    rtc_ds = rtc_ds.where(rtc_ds.vh > 0, drop=True)\n",
    "    rtc_ds = rtc_ds.where(rtc_ds.vv > 0, drop=True)\n",
    "    \n",
    "    #rtc_ds.to_netcdf(f'../data/S1_rtc/S1_snow-off_{rtc_ds.time.dt.strftime(\"%Y%m%d\").item()}_for_{aso_raster_fn.split(\"/\")[-1][:-4]}.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "41ee6e97-92a8-4a70-a90a-feca5a7904d8",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Returned 10 acquisitions\n",
      "taking mean of 2 snow-off rasters\n"
     ]
    }
   ],
   "source": [
    "rtc_for_aso_snowoff_mean(aso_raster_fn, rel_orbit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83556330-97b1-47f3-bf7e-1dd8006fd979",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21dd58a5-b992-42c7-bf69-51b655f015b4",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4d80ee0-922e-4c06-b780-309455a7edf7",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b724b26-166d-4bf9-8afe-163d0b398a5a",
   "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
}