0crop_S2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import gdal\n",
    "from os import listdir, makedirs\n",
    "from os.path import isdir, join, isfile\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "FOLDER_PATH = (\n",
    "    \"C:\\\\Users\\\\mracic\\\\Univerza v Ljubljani\\\\Fetai, Bujar - satimagets-master\\\\\"\n",
    ")\n",
    "dfs = pd.read_excel(FOLDER_PATH + \"Harmonisation_PS_S2.xlsx\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "locations = [\"Radenci\", \"Izola\", \"Jesenice\", \"Kranj\"]\n",
    "INPUT_PATH = \"Y:\\\\S2_10m\\\\\"\n",
    "OUTPUT_PATH = \"Y:\\\\Harmonisation_PS_S2\\\\M3sat_S2output\\\\\""
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "sub_dfs.iloc[0][\"PS_Path\"]"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "for loc_name in locations:\n",
    "    print(loc_name)\n",
    "    location_path = OUTPUT_PATH + loc_name + \"5\\\\\"\n",
    "    sub_dfs = dfs[dfs[\"Location\"] == loc_name]\n",
    "    # Open the image of the selected area\n",
    "    PS_PATH = sub_dfs.iloc[0][\"PS_Path\"]\n",
    "    if isfile(PS_PATH):\n",
    "        PS = gdal.Open(PS_PATH)\n",
    "    else:\n",
    "        print(\"Can't find file in path.\")\n",
    "\n",
    "    # Get bbox of AOI\n",
    "    geoTransform = PS.GetGeoTransform()\n",
    "    minx = geoTransform[0]\n",
    "    maxy = geoTransform[3]\n",
    "    maxx = minx + geoTransform[1] * PS.RasterXSize\n",
    "    miny = maxy + geoTransform[5] * PS.RasterYSize\n",
    "    data = None\n",
    "    extent = [minx, miny, maxx, maxy]\n",
    "\n",
    "    for year in range(2017, 2021):\n",
    "        list_S2 = [\n",
    "            o\n",
    "            for o in listdir(INPUT_PATH + str(year))\n",
    "            if isdir(join(INPUT_PATH + str(year), o))\n",
    "        ]\n",
    "\n",
    "        for observation in list_S2:\n",
    "            # check if image is in dataset\n",
    "            if sub_dfs[\"S2_filename\"].str.contains(observation).any():\n",
    "                # print(observation)\n",
    "                tiffs = [\n",
    "                    o\n",
    "                    for o in listdir(INPUT_PATH + str(year) + \"\\\\\" + observation)\n",
    "                    if \".tif\" in o\n",
    "                ]\n",
    "                data = [o for o in tiffs if \"p2atm_d96tm\" in o]\n",
    "                mask_data = [o for o in tiffs if \"p2atm_mask_d96tm\" in o]\n",
    "                # crop the raster data\n",
    "                gdal.Warp(\n",
    "                    location_path + \"\\\\\" + str(year) + \"\\\\\" + observation + \"_S2.tif\",\n",
    "                    INPUT_PATH + \"\\\\\" + str(year) + \"\\\\\" + observation + \"\\\\\" + data[0],\n",
    "                    outputBounds=extent,\n",
    "                    dstSRS=\"EPSG:32633\",\n",
    "                )\n",
    "                gdal.Warp(\n",
    "                    location_path\n",
    "                    + \"\\\\\"\n",
    "                    + str(year)\n",
    "                    + \"\\\\\"\n",
    "                    + observation\n",
    "                    + \"_mask_S2.tif\",\n",
    "                    INPUT_PATH\n",
    "                    + \"\\\\\"\n",
    "                    + str(year)\n",
    "                    + \"\\\\\"\n",
    "                    + observation\n",
    "                    + \"\\\\\"\n",
    "                    + mask_data[0],\n",
    "                    outputBounds=extent,\n",
    "                    dstSRS=\"EPSG:32633\",\n",
    "                )  #  srcSRS=\"EPSG:3794\","
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# check for number of unique files\n",
    "dfs[(dfs[\"Location\"] == \"Radenci\") & (dfs[\"Year\"] == 2020)][\n",
    "    \"S2_filename\"\n",
    "].unique().shape"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  }
 ],
 "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}