multiTempDPM.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e15f0603-4ac4-4c8f-915f-1c0d005ab443",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, shutil\n",
    "import re\n",
    "import time\n",
    "\n",
    "from datetime import datetime as dt\n",
    "from datetime import timedelta\n",
    "from pathlib import Path\n",
    "from zipfile import ZipFile\n",
    "\n",
    "import requests\n",
    "import numpy as np\n",
    "import pyproj\n",
    "import geopandas as gpd\n",
    "import fiona\n",
    "from osgeo import gdal\n",
    "import rasterio as rio\n",
    "from rasterio.merge import merge\n",
    "from rasterio.features import shapes\n",
    "from ost import Sentinel1Batch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4775d350-e46b-45b0-acea-ce00a26f277e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_product_on_asf(identifier, uname, pword):\n",
    "    \n",
    "    url = f\"https://datapool.asf.alaska.edu/SLC/SA/{identifier}.zip\"\n",
    "    with requests.Session() as session:\n",
    "        session.auth = (uname, pword)\n",
    "        request = session.request(\"get\", url)\n",
    "        response = session.get(request.url, auth=(uname, pword), stream=True)\n",
    "        return response.status_code\n",
    "\n",
    "\n",
    "def create_dmp(aoi, start_date, end_date, scihub_uname, scihub_pword, asf_uname, asf_pword):\n",
    "    \n",
    "    # get start and end date into dt format\n",
    "    event_start = dt.strptime(start_date, '%Y-%m-%d')\n",
    "    event_end = dt.strptime(end_date, '%Y-%m-%d')\n",
    "    \n",
    "    # set a period around start and end date for data search\n",
    "    search_start = dt.strftime(event_start + timedelta(days=-60), '%Y-%m-%d')\n",
    "    search_end = dt.strftime(event_end + timedelta(days=30), '%Y-%m-%d')\n",
    "    \n",
    "    # define project dir \n",
    "    aoi_name = Path(aoi).stem\n",
    "    if event_start == event_end:\n",
    "        project_dir = Path().home().joinpath(f'module_results/Damage_Proxy_Maps/{aoi_name}_{start_date}')\n",
    "    else:\n",
    "        project_dir = Path().home().joinpath(f'module_results/Damage_Proxy_Maps/{aoi_name}_{start_date}_{end_date}')\n",
    "    \n",
    "    print(' Setting up project')\n",
    "    s1_slc = Sentinel1Batch(\n",
    "        project_dir=project_dir,\n",
    "        aoi = aoi,\n",
    "        start = search_start,\n",
    "        end = search_end,\n",
    "        product_type='SLC',\n",
    "        ard_type='OST-RTC'\n",
    "    )\n",
    "\n",
    "    from ost.helpers.settings import HERBERT_USER\n",
    "    \n",
    "    s1_slc.scihub_uname = scihub_uname if scihub_uname else HERBERT_USER['uname'] \n",
    "    s1_slc.scihub_pword = scihub_pword if scihub_pword else HERBERT_USER['pword']\n",
    "    s1_slc.asf_uname = asf_uname if asf_uname else HERBERT_USER['uname'] \n",
    "    s1_slc.asf_pword = asf_pword if asf_pword else HERBERT_USER['asf_pword']\n",
    "\n",
    "\n",
    "    print(' Searching for data')\n",
    "    s1_slc.search(base_url='https://scihub.copernicus.eu/dhus/')\n",
    "    \n",
    "    print(f' Found {len(s1_slc.inventory.relativeorbit.unique())} tracks to process')\n",
    "    for i, track in enumerate(s1_slc.inventory.relativeorbit.unique()):\n",
    "        \n",
    "        print(f' Starting to process track {track}.')\n",
    "        # filter by track\n",
    "        df = s1_slc.inventory[s1_slc.inventory.relativeorbit == track].copy()\n",
    "        \n",
    "        # make sure all products are on ASF\n",
    "        for i, row in df.iterrows():\n",
    "            status = check_product_on_asf(row.identifier, s1_slc.asf_uname, s1_slc.asf_pword)\n",
    "            if status != 200:\n",
    "                df = df[df.identifier != row.identifier]\n",
    "        \n",
    "        print(df)\n",
    "        break\n",
    "        # get all acquisitions dates for that track\n",
    "        datelist = sorted([dt.strptime(date, '%Y%m%d') for date in df.acquisitiondate.unique()])\n",
    "        to_process_list = []\n",
    "        for date in datelist:\n",
    "            \n",
    "            # cehck if we have an image after end date\n",
    "            if datelist[-1] < event_end:\n",
    "                print(f' No image available after the end date for track {track}', 'warning')\n",
    "                break\n",
    "                \n",
    "            # ignore dates before start of event\n",
    "            if date < event_start:\n",
    "                continue\n",
    "            \n",
    "            # get index of date in datelist\n",
    "            idx = datelist.index(date)\n",
    "            if idx < 2:\n",
    "                print(f' Not enough pre-event images available for track {track}', 'warning')\n",
    "                break\n",
    "            \n",
    "            # add dates to process list, if not already there\n",
    "            # we take care of the two images needed before the start date\n",
    "            [to_process_list.append(date) for date in datelist[idx-2:idx+1] if date not in to_process_list]\n",
    "    \n",
    "            # once we added the last image after the the end of the event we can stop\n",
    "            if date > event_end:\n",
    "                break\n",
    "        \n",
    "            if len(to_process_list) < 3:\n",
    "                print(f' Not enough images available for track {track}', 'warning')\n",
    "                break\n",
    "\n",
    "        # turn the dates back to strings so we can use to filter the data inventory\n",
    "        final_dates = [dt.strftime(date, '%Y%m%d') for date in to_process_list]\n",
    "        print(' We\\'ll download/process scenes for the following dates:')\n",
    "        print(final_dates)\n",
    "        final_df = df[df.acquisitiondate.isin(final_dates)]\n",
    "        \n",
    "        print(f' {len(final_df)} identfied for download')\n",
    "        print(' Downloading relevant Sentinel-1 SLC scenes ... (this may take a while)')\n",
    "        s1_slc.download(final_df, mirror=2, concurrent=10, uname=s1_slc.asf_uname, pword=s1_slc.asf_pword)\n",
    "            \n",
    "        print(' Create burst inventory')\n",
    "        s1_slc.create_burst_inventory(final_df)\n",
    "\n",
    "        # setting ARD parameters\n",
    "        print(' Setting processing parameters')\n",
    "        s1_slc.ard_parameters['single_ARD']['resolution'] = 30 # in metres\n",
    "        s1_slc.ard_parameters['single_ARD']['create_ls_mask'] = False\n",
    "        s1_slc.ard_parameters['single_ARD']['backscatter'] = False\n",
    "        s1_slc.ard_parameters['single_ARD']['coherence'] = True\n",
    "        s1_slc.ard_parameters['single_ARD']['coherence_bands'] = 'VV'  # 'VV, VH'\n",
    "\n",
    "        # production of polarimetric layers\n",
    "        s1_slc.ard_parameters['single_ARD']['H-A-Alpha'] = False # does not give a lot of additional information\n",
    "\n",
    "        # resampling of image (not so important)\n",
    "        s1_slc.ard_parameters['single_ARD']['dem']['dem_name'] = \"SRTM 1Sec HGT\"\n",
    "        s1_slc.ard_parameters['single_ARD']['dem']['image_resampling'] = 'BILINEAR_INTERPOLATION'  # 'BILINEAR_INTERPOLATION'\n",
    "\n",
    "        # multi-temporal speckle filtering is quite effective\n",
    "        s1_slc.ard_parameters['time-series_ARD']['mt_speckle_filter']['filter'] = 'Boxcar'\n",
    "        s1_slc.ard_parameters['time-series_ARD']['remove_mt_speckle'] = True\n",
    "        s1_slc.ard_parameters['time-scan_ARD']['metrics'] = ['min', 'max']\n",
    "        s1_slc.ard_parameters['time-scan_ARD']['remove_outliers'] = False\n",
    "        s1_slc.ard_parameters['mosaic']['cut_to_aoi'] = True\n",
    "\n",
    "        # set tmp_dir\n",
    "        s1_slc.config_dict['temp_dir'] = '/ram'\n",
    "\n",
    "        #\n",
    "        workers = int(4) if os.cpu_count()/4 > 4 else int(os.cpu_count()/4)\n",
    "        print(f' We process {workers} bursts in parallel.')\n",
    "        s1_slc.config_dict['max_workers'] = workers\n",
    "        s1_slc.config_dict['executor_type'] = 'concurrent_processes'\n",
    "\n",
    "        # process\n",
    "        print(f' {len(s1_slc.burst_inventory)} bursts to process in total')\n",
    "        print(\" Processing... (this may take a while)\")\n",
    "        s1_slc.bursts_to_ards(\n",
    "            timeseries=True, \n",
    "            timescan=True, \n",
    "            mosaic=False,\n",
    "            overwrite=False # TO BE CHAGNED\n",
    "        )\n",
    "        \n",
    "        \n",
    "        print(\"calculate change and merge results\")\n",
    "        bursts = list(s1_slc.processing_dir.glob(f'[A,D]*{track}*'))\n",
    "        \n",
    "        for burst in bursts:\n",
    "            print(f' Coherent change calculation for burst {burst}')\n",
    "            # get track\n",
    "            track_name = burst.name[:4]\n",
    "\n",
    "            # in and out files\n",
    "            coh_min = burst.joinpath('Timescan/01.coh.VV.min.tif')\n",
    "            coh_max = burst.joinpath('Timescan/02.coh.VV.max.tif')\n",
    "            dstnt_file = burst.joinpath(f\"Timescan/ccd_{burst.name}.tif\")\n",
    "            print(coh_min, coh_max)\n",
    "            with rio.open(coh_max) as pre_coh:\n",
    "                pre_arr = pre_coh.read()\n",
    "                \n",
    "                # get metadata for destination file\n",
    "                meta = pre_coh.meta\n",
    "                meta.update(dtype='uint8', nodata=0)\n",
    "            \n",
    "                with rio.open(coh_min) as post_coh:\n",
    "                    post_arr = post_coh.read()\n",
    "\n",
    "                    # calulate difference\n",
    "                    coh_diff = np.subtract(pre_arr, post_arr)\n",
    "                    coh_diff[coh_diff < 0.27] = 0\n",
    "                    coh_diff = coh_diff * 100\n",
    " \n",
    "                    with rio.open(dstnt_file, 'w', **meta) as dstnt:\n",
    "                        dstnt.write(coh_diff.astype('uint8'))  \n",
    "            \n",
    "        # -----------------------------------------\n",
    "        # and merge the result\n",
    "        print(' Mosaic CCD files')\n",
    "        src_files_to_mosaic = []\n",
    "        for file in s1_slc.processing_dir.glob(f\"*[A,D]*{track}_*/Timescan/ccd*tif\"):\n",
    "            src = rio.open(file)\n",
    "            src_files_to_mosaic.append(src)\n",
    "\n",
    "\n",
    "        mosaic, out_trans = merge(src_files_to_mosaic)\n",
    "        out_meta = src.profile.copy()\n",
    "\n",
    "        # Update the metadata\n",
    "        out_meta.update(\n",
    "            driver='GTiff',\n",
    "            height=mosaic.shape[1],\n",
    "            width= mosaic.shape[2],\n",
    "            transform=out_trans,\n",
    "            crs=src.crs,\n",
    "            tiled=True,\n",
    "            blockxsize=128,\n",
    "            blockysize=128,\n",
    "            compress='lzw'\n",
    "        )\n",
    "\n",
    "        tmp_dir = Path(s1_slc.config_dict['temp_dir'])\n",
    "        tmp_mrg = tmp_dir.joinpath(f\"ccd_{track_name}.tif\")\n",
    "        with rio.open(tmp_mrg, \"w\", **out_meta) as dest:\n",
    "            dest.write(mosaic)\n",
    "\n",
    "        # close datasets\n",
    "        [src.close for src in src_files_to_mosaic]\n",
    "\n",
    "        # crop to aoi (some ost routine)\n",
    "        with fiona.open(aoi, \"r\") as shapefile:\n",
    "            shapes_ = [feature[\"geometry\"] for feature in shapefile]\n",
    "\n",
    "        with rio.open(tmp_mrg) as src:\n",
    "            out_image, out_transform = rio.mask.mask(src, shapes_, crop=True)\n",
    "            out_meta = src.profile\n",
    "\n",
    "        out_meta.update({\n",
    "            \"driver\": \"GTiff\",\n",
    "            \"height\": out_image.shape[1],\n",
    "            \"width\": out_image.shape[2],\n",
    "            \"transform\": out_transform\n",
    "        })\n",
    "\n",
    "        # create final output directory\n",
    "        outdir = s1_slc.project_dir.joinpath('Damage_Proxy_Map')\n",
    "        outdir.mkdir(parents=True, exist_ok=True)\n",
    "        out_ds_tif = outdir.joinpath(f\"ccd_{track_name}.tif\")\n",
    "        with rio.open(out_ds_tif, \"w\", **out_meta) as dest:\n",
    "            dest.write(out_image)\n",
    "\n",
    "        # delete tmpmerge \n",
    "        tmp_mrg.unlink()\n",
    "        # -----------------------------------------\n",
    "\n",
    "        # -----------------------------------------\n",
    "        # kmz and dmp output\n",
    "        # write a color file to tmp\n",
    "        ctfile = tmp_dir.joinpath('colourtable.txt')\n",
    "        f = open(ctfile, \"w\")\n",
    "        ct = [\"0 0 0 0 0\\n\"\n",
    "            \"27 253 246 50 255\\n\"\n",
    "            \"35 253 169 50 255\\n\"\n",
    "            \"43 253 100 50 255\\n\"\n",
    "            \"51 253 50 50 255\\n\"\n",
    "            \"59 255 10 10 255\\n\"\n",
    "            \"255 253 0 0 255\"\n",
    "            ]\n",
    "        f.writelines(ct)\n",
    "        f.close()\n",
    "        \n",
    "        print(' Create DPM tif')\n",
    "        out_dmp_tif = outdir.joinpath(f\"dmp_{track_name}.tif\")\n",
    "        demopts = gdal.DEMProcessingOptions(colorFilename=str(ctfile), addAlpha=True)\n",
    "        gdal.DEMProcessing(str(out_dmp_tif), str(out_ds_tif), 'color-relief', options=demopts)         \n",
    "\n",
    "        print(' Create DPM kmz')\n",
    "        opts = gdal.TranslateOptions(format='KMLSUPEROVERLAY', creationOptions=[\"format=png\"])\n",
    "        gdal.Translate(str(out_dmp_tif.with_suffix('.kmz')), str(out_dmp_tif), options=opts)\n",
    "        # -----------------------------------------\n",
    "\n",
    "        # -----------------------------------------\n",
    "        # polygonize (to points)\n",
    "        print(' Create DPM geojson')\n",
    "        with rio.open(out_ds_tif) as src:\n",
    "            image = src.read()\n",
    "            mask = image != 0\n",
    "\n",
    "            results = (\n",
    "                    {\n",
    "                        'properties': {'raster_val': v},\n",
    "                        'geometry': s\n",
    "                    }\n",
    "                    for i, (s, v)\n",
    "                    in enumerate(\n",
    "                        shapes(image, mask=mask, transform=src.transform)\n",
    "                    )\n",
    "            )\n",
    "\n",
    "        geoms = list(results)  \n",
    "        gpd_polygonized_raster  = gpd.GeoDataFrame.from_features(geoms)\n",
    "        gpd_polygonized_raster['geometry'] = gpd_polygonized_raster['geometry'].centroid\n",
    "        gpd_polygonized_raster.to_file(out_dmp_tif.with_suffix('.geojson'), driver='GeoJSON')\n",
    "        # -----------------------------------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74eb8d76-c039-4733-ad13-2bae2ed58496",
   "metadata": {},
   "outputs": [],
   "source": [
    "start = '2022-02-24'\n",
    "end = '2022-03-14'\n",
    "aoi = '/home/vollrath/DPMs/eriks/eriks_aoi.gpkg'\n",
    "scihub_uname = ''\n",
    "scihub_pword = ''\n",
    "asf_uname = ''\n",
    "asf_pword = ''\n",
    "\n",
    "create_dmp(aoi, start, end, scihub_uname, scihub_pword, asf_uname, asf_pword)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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": 5
}