acquistion_s2_l3a_dlr.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Info\n",
    "\n",
    "Script to download Level3A (WASP-generated) Sentinel2A data processed by the German Aerospace Centre.\n",
    "This data contains low-cloud/cloudfree, monthly averages of Sentinel2A scenes.\n",
    "\n",
    "Processing steps:\n",
    "* download tiles based on specified aoi, timeframe & bands\n",
    "* clipping tiles to aoi extent & creating a mosaic for each timestep and band\n",
    "\n",
    "What is not covered so far & may be improved in upcoming releases\n",
    "* consideration of quality control mask provided along with data for two reasons\n",
    "  + inaccuracy of the mask\n",
    "  + low significance of mask for clouds as processing chain removes almost all cloudy pixels\n",
    "\n",
    "But mask can be downloaded & clipped for AoI and then statistics be calculated to decide if it should be excluded (if above threshold value).\n",
    "Also visual inspection possible directly in script via plotting rbg image for region. Print statistics for each layer (on proprotion of masked pixels). Further idea: Script for downloading L2A data for specific AoI with cloud thresholding possibility (first look what is already provided on github in this regard). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import libraries\n",
    "import geopandas as gpd\n",
    "import glob\n",
    "import minihit\n",
    "import numpy as np\n",
    "import os\n",
    "import rasterio\n",
    "import requests\n",
    "import shutil\n",
    "import wget\n",
    "import zipfile\n",
    "\n",
    "from bs4 import BeautifulSoup\n",
    "from pathos.multiprocessing import ProcessingPool as Pool\n",
    "from tqdm.notebook import tqdm\n",
    "from rasterio.mask import mask\n",
    "from rasterio.merge import merge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "### variables to set ###\n",
    "# download dir specifies existing folder to write scences to disk\n",
    "download_dir = os.path.join(os.getcwd(), \"test_data\", \"magdeburger_boerde\", \"s2_l3a_data\")\n",
    "\n",
    "# aoi has to be specified as geojson\n",
    "# must be in parent dir of s2_l3_data folder\n",
    "# can be created online, e.g. via https://geojson.io\n",
    "aoi_geojson = os.path.join(os.path.dirname(download_dir), \"aoi.geojson\")\n",
    "\n",
    "# timeframe\n",
    "years = [2021]\n",
    "months = [x for x in range(3,11)]\n",
    "\n",
    "# bands\n",
    "bands = [\"B2\", \"B3\", \"B4\", \"B8\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\felix\\AppData\\Roaming\\Python\\Python310\\site-packages\\geopandas\\geodataframe.py:2196: UserWarning: `keep_geom_type=True` in overlay resulted in 4 dropped geometries of different geometry types than df1 has. Set `keep_geom_type=False` to retain all geometries\n",
      "  return geopandas.overlay(\n",
      "C:\\Users\\felix\\AppData\\Roaming\\Python\\Python310\\site-packages\\geopandas\\geodataframe.py:2196: UserWarning: `keep_geom_type=True` in overlay resulted in 48 dropped geometries of different geometry types than df1 has. Set `keep_geom_type=False` to retain all geometries\n",
      "  return geopandas.overlay(\n"
     ]
    }
   ],
   "source": [
    "# create dataset with utm tiles\n",
    "# download utm tile shp\n",
    "url = \"https://codeload.github.com/justinelliotmeyers/Sentinel-2-Shapefile-Index/zip/refs/heads/master\"\n",
    "wget.download(url, download_dir)\n",
    "with zipfile.ZipFile(os.path.join(download_dir, \"Sentinel-2-Shapefile-Index-master.zip\"), \"r\") as zip_ref:\n",
    "    zip_ref.extractall(download_dir)\n",
    "\n",
    "# read utm & delete copies on disk\n",
    "utm = gpd.read_file(os.path.join(download_dir, \"Sentinel-2-Shapefile-Index-master\", \"sentinel_2_index_shapefile.shp\"))\n",
    "shutil.rmtree(os.path.join(download_dir, \"Sentinel-2-Shapefile-Index-master\"))\n",
    "os.remove(os.path.join(download_dir, \"Sentinel-2-Shapefile-Index-master.zip\"))\n",
    "\n",
    "# restrict utm dataset to germany\n",
    "world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))\n",
    "germany = world[world[\"name\"] == \"Germany\"]\n",
    "\n",
    "# split utm tiles in overlapping areas\n",
    "_utm = utm.sjoin(germany).iloc[:,0:2].overlay(utm.sjoin(germany).iloc[:,0:2], how='union')\n",
    "_utm = _utm[_utm[\"Name_1\"] != _utm[\"Name_2\"]].dropna()\n",
    "\n",
    "_utm2 = utm.sjoin(germany).iloc[:,0:2].overlay(_utm, how='symmetric_difference')\n",
    "_utm2 = _utm2.rename(columns={\"Name\":\"Name_3\"})\n",
    "utm = _utm.append(_utm2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['32UPC']"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# intersect aoi with utm to determine relevant utm tiles \n",
    "aoi = gpd.read_file(aoi_geojson)\n",
    "intersect = gpd.sjoin(aoi, utm)\n",
    "set_utm_tiles = [{x,y,z} for x,y,z in list(zip(intersect[\"Name_1\"], intersect[\"Name_2\"], intersect[\"Name_3\"]))]\n",
    "set_utm_tiles = [{x for x in tile if x==x} for tile in set_utm_tiles]\n",
    "\n",
    "# solve hitting set problem (a.k.a. set covering problem)\n",
    "# aim: filter absolute minimum of necessary utm tiles to cover aoi\n",
    "rc = minihit.RcTree(set_utm_tiles)\n",
    "rc.solve(prune=True, sort=False)\n",
    "utm_tiles = sorted(list(rc.generate_minimal_hitting_sets())[0])\n",
    "utm_tiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "73d7993c7db04e878c4fdf7ee73c1686",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "constructing download urls & checking url validity:   0%|          | 0/1 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- Number/Availability of Scenes ---\n",
      "The current specification of aoi, timeframe and bands corresponds to a total of 32 individual tiles to be downloaded.\n",
      "Out of this 32 tiles (100%) are available on the server.\n",
      "\n",
      "--- Download volume ---\n",
      "Estimated download volume: ~7GB\n"
     ]
    }
   ],
   "source": [
    "# construction of download urls following the schema baseurl/utm_tile/year/month/band \n",
    "baseurl = \"https://download.geoservice.dlr.de/S2_L3A_WASP/files/\"\n",
    "download_urls = []\n",
    "\n",
    "for tile in tqdm(utm_tiles, desc='constructing download urls & checking url validity'):\n",
    "    # construct tiles\n",
    "    tile_url = os.path.join(baseurl, tile[0:2], tile[2:3], tile[3:]).replace('\\\\','/')\n",
    "    if requests.get(tile_url).status_code == 200:\n",
    "        # construct years\n",
    "        for year in years:\n",
    "            year_url = os.path.join(tile_url, str(year)).replace('\\\\','/')\n",
    "            if requests.get(year_url).status_code == 200:\n",
    "                # construct months\n",
    "                for month in months:\n",
    "                    month = \"0\" + str(month) if len(str(month))==1 else str(month)\n",
    "                    path_spec_month = \"SENTINEL2X_\" + str(year) + month + \"15-000000-000_L3A_T\" + tile + \"_C_V1-2\"\n",
    "                    month_url = os.path.join(year_url, path_spec_month).replace('\\\\','/')\n",
    "                    if requests.get(month_url).status_code == 200:\n",
    "                        # construct bands\n",
    "                        for band in bands:\n",
    "                            page_bands = requests.get(month_url)\n",
    "                            page_bands = BeautifulSoup(page_bands.content, \"html.parser\")\n",
    "                            path_spec_band = path_spec_month + \"_FRC_\" + band + \".tif\"\n",
    "                            if page_bands.find_all(\"a\", href=lambda x: x.startswith(path_spec_band)):\n",
    "                                band_url = os.path.join(month_url, path_spec_band).replace('\\\\','/')\n",
    "                                download_urls.append(band_url)\n",
    "                            else:\n",
    "                                print(f\"URL {band_url} does not exist or cannot be queried\")\n",
    "                    else:\n",
    "                        print(f\"URL {month_url} does not exist or cannot be queried\")\n",
    "            else:\n",
    "                print(f\"URL {year_url} does not exist or cannot be queried\")\n",
    "    else:\n",
    "        print(f\"URL {tile_url} does not exist or cannot be queried\")\n",
    "\n",
    "# evaluation of the availability and memory consumption for the selected tiles\n",
    "n_tiles = len(utm_tiles) * len(years) * len(months) * len(bands)\n",
    "\n",
    "memory_bands = []\n",
    "memory_bands.extend([240 for x in bands if x=='B8'])\n",
    "memory_bands.extend([200 for x in bands if x in ['B2','B3','B4']])\n",
    "memory_bands.extend([55 for x in bands if x in ['B5','B6','B7', 'B8A', 'B11', 'B12']])\n",
    "\n",
    "memory_all = n_tiles * sum(memory_bands) / len(bands)\n",
    "\n",
    "print(\"--- Number/Availability of Scenes ---\")\n",
    "print(f\"The current specification of aoi, timeframe and bands corresponds to a total of {n_tiles} individual tiles to be downloaded.\")\n",
    "print(f\"Out of this {len(download_urls)} tiles ({len(download_urls)/n_tiles:.0%}) are available on the server.\")\n",
    "print(\"\")\n",
    "\n",
    "print(\"--- Download volume ---\")\n",
    "print(f\"Estimated download volume: ~{round(memory_all/1000) if memory_all>2500 else round(memory_all/1000,1)}GB\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c8feaf47bb5d402ca7bb3916a858856f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "downloading data:   0%|          | 0/32 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A total of 32 files is now present in specified download directory.\n"
     ]
    }
   ],
   "source": [
    "# download of data\n",
    "# using multithreading to accelerate download\n",
    "# note: p.starmap stops download after each submitted chunk of tasks\n",
    "# note: using p.imap works & progress bar (tqdm) can still be displayed\n",
    "def download_scenes(url, dir):\n",
    "    import os\n",
    "    import sys\n",
    "    import wget\n",
    "    if not os.path.isfile(os.path.join(dir, url.split(\"/\")[-1])):\n",
    "        try:\n",
    "            wget.download(url, dir)\n",
    "        except:\n",
    "            pass\n",
    "\n",
    "with Pool(min(os.cpu_count(), 4)) as p:\n",
    "    _ = list(tqdm(p.imap(download_scenes, download_urls, [download_dir for x in range(len(download_urls))]), \n",
    "                  total=len(download_urls),\n",
    "                  desc='downloading data', \n",
    "                  smoothing=0.2))\n",
    "\n",
    "# in case of issues: use single core download\n",
    "# for url in tqdm(download_urls, desc='downloading data', smoothing=0.2):\n",
    "#     if not os.path.isfile(os.path.join(download_dir, url.split(\"/\")[-1])):\n",
    "#         wget.download(url, download_dir)\n",
    "\n",
    "print(f\"A total of {len(glob.glob(os.path.join(download_dir, 'SENTINEL2X*.tif')))} files is now present in specified download directory.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "93f680667fd544148295e36ddd545ccc",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "merging & masking of bands:   0%|          | 0/32 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# procedure to mosaic and mask tiles per date & band\n",
    "downloaded_tiles = glob.glob(os.path.join(download_dir, 'SENTINEL2X*.tif'))\n",
    "tiles_per_date_band = [\"_\".join(x.split(\"\\\\\")[-1].split(\"_\")[:3]) + \"*\" + x.split(\"\\\\\")[-1].split(\"_\")[-1] for x in downloaded_tiles]\n",
    "tiles_per_date_band = list(set(tiles_per_date_band))\n",
    "\n",
    "for tiles in tqdm(tiles_per_date_band, desc=\"merging & masking of bands\"):\n",
    "\n",
    "    # open tiles\n",
    "    src_files_to_mosaic = []\n",
    "    for tile in glob.glob(os.path.join(download_dir, tiles)):\n",
    "        src = rasterio.open(tile)\n",
    "        src_files_to_mosaic.append(src)\n",
    "\n",
    "    # merge tiles\n",
    "    mosaic, out_trans = merge(src_files_to_mosaic)\n",
    "\n",
    "    out_meta = src.meta.copy()\n",
    "    out_meta.update({\"driver\": \"GTiff\",\n",
    "                    \"height\": mosaic.shape[1],\n",
    "                    \"width\": mosaic.shape[2],\n",
    "                    \"transform\": out_trans}\n",
    "                    )\n",
    "\n",
    "    out_dir = os.path.join(download_dir, tiles.replace(\"-000000-000\", \"\").replace(\"*\", \"_\"))\n",
    "    with rasterio.open(out_dir, \"w\", **out_meta) as dest:\n",
    "        dest.write(mosaic)\n",
    "\n",
    "    # mask tiles with aoi geojson\n",
    "    with rasterio.open(out_dir) as src:\n",
    "        masked_mosaic, out_trans2 = mask(src, aoi[\"geometry\"].to_crs(out_meta[\"crs\"].to_string()), crop=True)\n",
    "\n",
    "        out_meta = src.meta.copy()\n",
    "        out_meta.update({\"driver\": \"GTiff\",\n",
    "                        \"height\": masked_mosaic.shape[1],\n",
    "                        \"width\": masked_mosaic.shape[2],\n",
    "                        \"transform\": out_trans2}\n",
    "                        )\n",
    "\n",
    "    with rasterio.open(out_dir, \"w\", **out_meta) as dest:\n",
    "        dest.write(masked_mosaic)\n",
    "\n",
    "    #remove original data & keep only mosaic\n",
    "    for rasters in src_files_to_mosaic:\n",
    "        rasters.close()\n",
    "        \n",
    "    for tile in glob.glob(os.path.join(download_dir, tiles)):\n",
    "        os.remove(tile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compatibility with ecognition workflow\n",
    "# rename scenes\n",
    "downloaded_scenes = [x for x in glob.glob(os.path.join(download_dir, 'SENTINEL2X*.tif'))]\n",
    "\n",
    "for scene_path in downloaded_scenes:\n",
    "    scene = scene_path.split(\"\\\\\")[-1]\n",
    "    scene_date = scene.split(\"_\")[1]\n",
    "    scene_band = scene.split(\"_\")[-1].split(\".tif\")[0]\n",
    "    bandidx_to_name = {\"B2\": \"blue\", \"B3\": \"green\", \"B4\": \"red\", \"B8\": \"nir\"}\n",
    "    new_scene_name = f\"s2_l3a_{scene_date}_{bandidx_to_name[scene_band]}.tif\"\n",
    "    os.rename(scene_path, os.path.join(download_dir, new_scene_name))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# output buffered aoi\n",
    "# aim: avoid border effects\n",
    "# means: 20m inwards buffer via Lambert Azimuthal Equal Area projection \n",
    "# aoi_ecognition = aoi.to_crs(crs=9821).buffer(-2000).to_crs(crs=4326)\n",
    "# aoi_ecognition_path = os.path.join(os.path.dirname(download_dir), \"aoi_ecognition.geojson\")\n",
    "# aoi_ecognition.to_file(aoi_ecognition_path, driver='GeoJSON')"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "26de051ba29f2982a8de78e945f0abaf191376122a1563185a90213a26c5da77"
  },
  "kernelspec": {
   "display_name": "Python 3.10.1 64-bit",
   "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.10.1"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}