exercise.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Algal Bloom Detection\n",
    "## or\n",
    "# Earth Observation and Data Science: data preparation\n",
    "\n",
    "### Example how to:\n",
    "\n",
    "* Reproject spatial data;\n",
    "* Clean scenes from the clouds;\n",
    "* Detect harmful algal blooms in the Landsat 8 scenes by band indices taken from the literature review."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import these modules and functions: os, numpy, pandas, rasterio, fiona\n",
    "# from_epsg from fiona.crs, Proj from pyproj, transform from pyproj,\n",
    "# pyplot from matplotlib, clip_area from scripts.clip_area,\n",
    "# read_landsat_images from scripts.read_landsat_images, show_band from scripts.show_band\n",
    "# set %matplotlib inline or %matplotlib notebook\n",
    "# optionally set np.seterr(divide='ignore', invalid='ignore')\n",
    "\n",
    "# Runetime Warnings Disabling\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Class for preparation of Landsat bands is given for you. If you would like to implement it\n",
    "# by yourself then remember to preapre data structure which contains paths to the parent\n",
    "# folders with Landsat 8 bands and has field dedicated for the timestamp\n",
    "\n",
    "\n",
    "class RasterTimeSeries:\n",
    "    \"\"\"Class prepares dataframe object with addresses to the folders with Landsat bands. Dataframe index is a time when each scene was retrieved.\n",
    "    It is initialized with the folder which contains all landsat bands.\n",
    "    \n",
    "    Class methods:\n",
    "    prepare_directories(main_folder=None) - method prepares a list of folders with satellite images. Method performs search in the parent directory\n",
    "    given by the user or if it is None then it is taken from the class itself.\"\"\"\n",
    "    \n",
    "    def __init__(self, main_folder):\n",
    "        self.main_folder = main_folder\n",
    "        self.landsat_folders = None\n",
    "        self.landsat_dataframe = None\n",
    "    \n",
    "    def prepare_directories(self, main_folder=None):\n",
    "        \"\"\"Function prepares list with Landsat images directories in the given folder.\n",
    "        input 1: main_folder (string)\n",
    "        If None then self.main_folder is used as the parent directory of the folders with satellite images.\n",
    "        output: prepared_list - list with all folders with satellite images.\"\"\"\n",
    "        \n",
    "        if main_folder:\n",
    "            folder_name = os.listdir(main_folder)\n",
    "        else:\n",
    "            folder_name = self.main_folder\n",
    "        raw_list = os.listdir(folder_name)\n",
    "        prepared_list = []\n",
    "        for folder in raw_list:\n",
    "            if ((folder.startswith('LC')) and not (folder.endswith('.tar.gz'))):\n",
    "                prepared_list.append(folder_name + '/' + folder + '/')\n",
    "        self.landsat_folders = prepared_list\n",
    "        return prepared_list\n",
    "    \n",
    "    @staticmethod\n",
    "    def read_timestamp(folder):\n",
    "        \"\"\"Function reads date of acquisition of each Landsat 8 scene based on the folder with MTL file.\n",
    "        input 1: folder - string with Landsat folder (MTL file must be placed in it).\n",
    "        output: string with the scene acquisition time\"\"\"\n",
    "    \n",
    "        files = os.listdir(folder)\n",
    "        for file in files:\n",
    "            if file.startswith('LC') and file.endswith('MTL.txt'):\n",
    "                address = folder + '/' + file\n",
    "                with open(address, 'r') as f:\n",
    "                    for line in f:\n",
    "                        line = ''.join(line.split())\n",
    "                        if line.startswith('DATE_ACQUIRED='):\n",
    "                            line = line.replace('DATE_ACQUIRED=', '')\n",
    "                            return line\n",
    "    \n",
    "    def prepare_dataframe(self):\n",
    "        \"\"\"Function prepares dataframe with index column as a date of acquisition and the main column DIR which points\n",
    "        to the directories with Landsat 8 images. Function takes one argument: the list of folders with Landsat 8 datasets.\n",
    "        output: pandas dataframe with indexes set to the dates of acquisition and column indices folders with Landsat images.\"\"\"\n",
    "        \n",
    "        # DataFrame preparation\n",
    "        time_series = pd.DataFrame(self.landsat_folders, columns=['DIR'])\n",
    "        time_series.index = time_series['DIR'].apply(lambda x: str(pd.to_datetime(self.read_timestamp(x)).date()))\n",
    "        time_series.index.name = 'DATE'\n",
    "        time_series = time_series.sort_index()\n",
    "        self.landsat_dataframe = time_series\n",
    "        return time_series\n",
    "    \n",
    "    def __str__(self):\n",
    "        if self.landsat_folders is None:\n",
    "            return('Parent folder of all scenes: {}.\\\n",
    "            \\nList of subfolders is not created.\\\n",
    "            \\nDataframe with dates is not created.'.format(\n",
    "            self.main_folder))\n",
    "        elif self.landsat_dataframe is None:\n",
    "            return('Parent folder of all scenes: {}\\\n",
    "            \\nList of subfolders has {} directories.\\\n",
    "            \\nDataframe with dates is not created.'.format(\n",
    "            self.main_folder, len(self.landsat_folders)))\n",
    "        else:\n",
    "            return('Parent folder of all scenes: {}\\\n",
    "            \\nList of subfolders has {} directories.\\\n",
    "            \\nDataframe is created. List of all dates: {}.'.format(\n",
    "            self.main_folder,\n",
    "            len(self.landsat_folders),\n",
    "            self.landsat_dataframe.index.values))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create RasterTimeSeries instance with directory which contains Landsat 8 folders\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare list of Landsat 8 directories and timestamp of each scene (it could be read from the\n",
    "# directory name). Use prepare_directories() and then prepare_dataframe() methods of \n",
    "# the RasterTimeSeries class.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Again, class for Vector Files reading and reprojecting is given here. It is hard to prepare\n",
    "# it if you've never worked with the polygon files.\n",
    "\n",
    "# The main class for vector file handling\n",
    "\n",
    "class VectorData:\n",
    "    \"\"\"Class reads, stores, shows and reprojects vector shapefiles.\n",
    "    It is initialized with filename.\n",
    "    \n",
    "    Class methods:\n",
    "    get_vector_data(parameters='all', write_params=True) - method reads vector data and parse its parameters such as geometry,\n",
    "    properties, crs and schema into dictionaries. Method may store parameters inside the object if write_params is set\n",
    "    to True. Method returns desired parameters.\n",
    "    reproject_geometry(destination_epsg, update=False) - method changes projection of each point inside the Polygon into\n",
    "    projection given as an epsg number. It may update object geometry if the update parameter is set to True. Method returns\n",
    "    reprojected geometry dictionary.\n",
    "    show_vector_data(geometry=None) - method shows polygon. If geometry is not given then method uses class geometry.\"\"\"\n",
    "    \n",
    "    def __init__(self, filename):\n",
    "        self.v_file = filename\n",
    "        self.geometry = {}\n",
    "        self.properties = None\n",
    "        self.vec_crs = None\n",
    "        self.vec_schema = None\n",
    "        \n",
    "    def get_vector_data(self, parameters='all', write_params=True):\n",
    "        \"\"\"\n",
    "        Function reads vector and its parameters and return them and / or store them inside the objects instance.\n",
    "        input 1: parameters (string)\n",
    "        parameters: 'all', 'none', 'geometry', 'properties', 'crs', 'schema'\n",
    "        'all': returns tuple with geometry, properties, crs and schema;\n",
    "        'none': does not return anything;\n",
    "        'geometry', 'properties', 'crs' or 'schema': returns chosen parameter.\n",
    "        input 2: write_params (bool)\n",
    "        write_params: True, False\n",
    "        True: store all parameters in the object\n",
    "        False: do not store anything in the object\n",
    "        output: chosen parameters as a dict\"\"\"\n",
    "        \n",
    "        with fio.open(self.v_file, 'r') as masking_region:\n",
    "            geometry = [feature[\"geometry\"] for feature in masking_region]\n",
    "            properties = [feature['properties'] for feature in masking_region]\n",
    "            vec_crs = masking_region.crs\n",
    "            vec_schema = masking_region.schema\n",
    "        \n",
    "        if write_params:\n",
    "            print('--- Object geometry, properties, crs and schema updated ---')\n",
    "            self.geometry = geometry\n",
    "            self.properties = properties\n",
    "            self.vec_crs = vec_crs\n",
    "            self.vec_schema = vec_schema\n",
    "        else:\n",
    "            print('--- Object parameters not updated ---')\n",
    "            \n",
    "        output_dict = {'all': ({'geometry': geometry, 'properties': properties, 'crs': vec_crs, 'scheme': vec_schema}),\n",
    "                      'none': 0,\n",
    "                      'geometry': geometry,\n",
    "                      'properties': properties,\n",
    "                      'crs': vec_crs,\n",
    "                      'schema': vec_schema}\n",
    "        try:\n",
    "            return output_dict[parameters]\n",
    "        except KeyError:\n",
    "            raise KeyError(\"Parameter not available. Available parameters: 'all', 'none', 'geometry', 'properties', 'crs', 'schema'\")\n",
    "      \n",
    "    def reproject_geometry(self, destination_epsg, update=False):\n",
    "        \"\"\"\n",
    "        Function reprojects vector geometry and updates it.\n",
    "        input 1: destination_epsg (string or int)\n",
    "        input 2: update (bool)\n",
    "        True: update object's geometry\n",
    "        False: do not update object's geometry\n",
    "        output: reprojected geometry as a dict\"\"\"\n",
    "        \n",
    "        try:\n",
    "            destination_crs = from_epsg(destination_epsg)\n",
    "            proj_crs_in = Proj(self.vec_crs)\n",
    "            coordinates_list = self.geometry[0]['coordinates'][0]\n",
    "        except RuntimeError:\n",
    "            raise ValueError('Given EPSG is wrong or it is not stored in fiona to run from_epsg method')\n",
    "        except KeyError:\n",
    "            raise ValueError('You must update your object parameters by get_vector_data method with write_params set to True.')\n",
    "        proj_crs_out = Proj(init = destination_crs['init'])\n",
    "        coordinates_array = np.asarray(coordinates_list)\n",
    "        projected_g = []\n",
    "        \n",
    "        for g in coordinates_array:\n",
    "            transformed = transform(proj_crs_in, proj_crs_out, g[0], g[1])\n",
    "            points = (transformed[0], transformed[1],)\n",
    "            projected_g.append(points)\n",
    "        geometry_dict = {'coordinates': [projected_g], 'type': 'Polygon'}\n",
    "        \n",
    "        if update:\n",
    "            self.geometry = [geometry_dict]\n",
    "            self.vec_crs = destination_crs\n",
    "            \n",
    "        return [geometry_dict]\n",
    "    \n",
    "    def show_vector_data(self, geometry=None):\n",
    "        \"\"\"\n",
    "        Function shows vector geometry.\n",
    "        input 1: geometry (dict, geometry retrieved as a parameter from the vector data type).\n",
    "        If geometry is not given then method takes it from the objects instance.\"\"\"\n",
    "        \n",
    "        if geometry:\n",
    "            g = geometry\n",
    "        else:\n",
    "            g = self.geometry\n",
    "            \n",
    "        try:\n",
    "            coordinates_list = g[0]['coordinates'][0]\n",
    "        except KeyError:\n",
    "            raise ValueError(\n",
    "                'Geometry is not defined. Please, update geometry with get_vector_data method or provide valid geometry object to the method.')\n",
    "            \n",
    "        coordinates_array = np.asarray(coordinates_list)\n",
    "        plt.figure()\n",
    "        plt.plot(coordinates_array[:, 0], coordinates_array[:, 1])\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You may test this class and its methods with vector file of Poland. It is stored in\n",
    "# data/vector/Polska.shp file\n",
    "\n",
    "# First read file by the class constructor and get all parameters of this file \n",
    "# by get_vector_data() method and then show it by show_vector_data() method.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read CRS of this data. (VectorData class has vec_crs object, use it)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Change crs and projection to the epsg:32634 and save it to the new variable\n",
    "# Use reproject_geometry method. Do not overwrite actual CRS and set update\n",
    "# parameter to False.\n",
    "# Show the new geometry. Compare axes of this plot and plot in original CRS.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Additionally: plot both projections on the same figure to see that they cover different\n",
    "# geolocations. Use show_multiple_vectors function and pass into it list of your vectors.\n",
    "\n",
    "def show_multiple_vectors(geometries):\n",
    "    \"\"\"\n",
    "    Function shows vector geometry.\n",
    "    input 1: list of geometries \n",
    "    (dictionaries, geometry retrieved as a parameter from the vector data type).\"\"\"\n",
    "\n",
    "    plt.figure()\n",
    "    for geometry in geometries:\n",
    "        try:\n",
    "            coordinates_list = np.asarray(geometry[0]['coordinates'][0])\n",
    "            plt.scatter(coordinates_list[:, 0], coordinates_list[:, 1], s=0.1)\n",
    "        except KeyError:\n",
    "            raise ValueError(\n",
    "                'Geometry is not defined. Please, update geometry with get_vector_data method or provide valid geometry object to the method.')\n",
    "    plt.show()\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare new function to get the epsg number from Landsat raster. Use rasterio.open method.\n",
    "# You shouldn't read whole raster, if you use rasterio.open to get the Landsat file\n",
    "# then you may read only its crs. You will get dict CRS({'init': 'epsg:32634'})\n",
    "# EPSG is a number after epsg: string (32634) and this should be returned by the function\n",
    "\n",
    "def get_epsg_from_raster(raster_address):\n",
    "    \"\"\"Function reads raster data and gets its coordinate reference system.\n",
    "    Then it is converted to the European Petroleum Survey Group (EPSG) reference number.\n",
    "    Function was tested only with the Landsat rasters.\"\"\"\n",
    "    pass\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now you have the epsg number and you know how to reproject vector geometry.\n",
    "# It is the time to prepare area of interest for the further study. Read file zatoka_gdanska.shp \n",
    "# and prepare its geometry. Do not reproject this file yet. Try to use function clip_area:\n",
    "# clip_area(vector_data.geometry, address_to_landsat_image, output_filename.tif)\n",
    "# Read the error returned by rasterio. It is important to understand these errors.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reproject vector geometry and try to clip raster area to the new file with updated geometry\n",
    "# You may overwrite data (vector file)\n",
    "# Save clipped file as a \"test.tif\"\n",
    "# Show \"test.tif\". Read file with rasterio.open and show band with show_band() function.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# From this point you know how to prepare single band for analysis.\n",
    "# The next step is to make process automatic so you cannot write functions for each raster.\n",
    "# Maybe you've noticed that Landsat has many bands - we will use only few of them so your\n",
    "# method should include this fact.\n",
    "# Below is an example of implementation. Explanations for this implementation are presented\n",
    "# in the article.\n",
    "\n",
    "def clip_all(dataset, clipping_vector, bands=(2, 3, 4), folder_name='clipped/'):\n",
    "    \"\"\"Function performs clipping of all Landsat 8 bands given in the DIR column of the df parameter.\n",
    "    input 1: dataset - dataframe with DIR column where each record points to the directory with a Landsat bands\n",
    "    input 2: clipping_vector - vector geometry which is used as a area of interest\n",
    "    input 3: bands=(2, 3, 4) - tuple of bands which should be clipped. Bands 2, 3, 4 are used for the BGR image creation.\n",
    "    input 4: folder_name='clipped/' - directory where all clipped bands should be stored.\n",
    "    \n",
    "    output: [correct, damaged] - list of dicts with keys (dates) and values - correct images\n",
    "    and damaged images. Damaged images were not process and they not exist in the destination folder.\"\"\"\n",
    "    \n",
    "    data_range = len(dataset)\n",
    "    damaged = {}\n",
    "    correct = {}\n",
    "    for i in range(0, data_range):\n",
    "        date_str = dataset.index[i]\n",
    "        images = read_landsat_images(dataset['DIR'][i])\n",
    "        for band_number in bands:\n",
    "            band_to_clip = images[band_number]\n",
    "            filename = folder_name + 'LC_band' + str(band_number) + '_' + date_str + '.tif'\n",
    "            try:\n",
    "                clip_area(clipping_vector, band_to_clip, filename)\n",
    "            except (rio.errors.WindowError, ValueError):\n",
    "                print('\\nPolygon does not cover raster area. Error occurs when projections are different or raster extent is translated and it does not overlap polygon extent.\\n')\n",
    "                try: \n",
    "                    damaged[date_str].append(filename)\n",
    "                except KeyError:\n",
    "                    damaged[date_str] = [filename]\n",
    "                pass\n",
    "            \n",
    "            # Check file\n",
    "            try: \n",
    "                damaged[date_str]\n",
    "                print('File: {} not saved. Raster has not been clipped.'.format(\n",
    "                filename))\n",
    "            except KeyError:\n",
    "                try:\n",
    "                    with rio.open(filename) as tf:\n",
    "                         pass\n",
    "                except rio.RasterioIOError:\n",
    "                    raise IOError('File was: {} not saved successfully'.format(\n",
    "                    filename))\n",
    "                else:\n",
    "                    print('File: {} saved successfully'.format(\n",
    "                    filename))\n",
    "                    try: \n",
    "                        correct[date_str].append(filename)\n",
    "                    except KeyError:\n",
    "                        correct[date_str] = [filename]\n",
    "    return [correct, damaged]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare tuple of bands (1, 2, 3, 4, 5) and clip multiple files with this tuple.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare pandas DataFrame from generated dict - take into account only correct values.\n",
    "# Use method pandas.DataFrame.from_dict(). Set orient parameter to 'index'. Show head() of this\n",
    "# DataFrame.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try to show band from the date '2018-08-03'. Use DataFrame generated in the previous step.\n",
    "# Use df.loc[row_indexer, column_indexer] to get the proper value\n",
    "# There's a little problem in there...\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Clouds are common problem with remote sensing data. We must get rid of them!\n",
    "# It is somewhat complicated - you must dig in the Landsat 8 documentation to prepare\n",
    "# functions for this task. I've prepared it for you as a class PixelQA. We will work with it\n",
    "# for a while. You may see that methods of this class have many parameters. It is good to know\n",
    "# them because sometimes your algorithms must be very accurate and sometimes not.\n",
    "\n",
    "class PixelQA:\n",
    "    \"\"\"Class removes clouds and cloud shadows from images based on the quality band given with Landsat's scene.\n",
    "    Class is initialized with the folder with all Landsat scenes and the QA band.\n",
    "    \n",
    "    Class methods:\n",
    "    prepare values(clear_key='all') - method prepares values to remove with accordance to the Landsat documentation:\n",
    "    https://landsat.usgs.gov/landsat-surface-reflectance-quality-assessment.\n",
    "    clear_pixels(band, qa_input = None, clear_by='all', nofill_value=0.0) - method clears pixels of a given band.\n",
    "    Method takes band for cleaning (band) multiply it with (qa_input) where all values to remove are zeros and\n",
    "    other values are ones. Removed pixels are filled by nofill_value.\n",
    "    IMPORTANT:\n",
    "    Script works only with the Landsat 8 Level-2 data products. To process Level 1 products you should use different\n",
    "    values for cloud and cloud confidence estimation.\"\"\"\n",
    "    \n",
    "    \n",
    "    def __init__(self, qa_file=None, data_folder=None):\n",
    "        \n",
    "        # prepare qa file address\n",
    "        if qa_file is None:\n",
    "            file_list = os.listdir(data_folder)\n",
    "            for f in file_list:\n",
    "                if (f.startswith('LC') and f.endswith('pixel_qa.tif')):\n",
    "                    self.bqa_file = data_folder + f\n",
    "        else:\n",
    "            self.bqa_file = qa_file\n",
    "        \n",
    "        # prepare qa band\n",
    "        with rio.open(self.bqa_file) as f:\n",
    "            self.qa_band = f.read(1)\n",
    "            \n",
    "        # Dictionary with cloud pixel values for pixel_qa band \n",
    "        # (https://landsat.usgs.gov/landsat-surface-reflectance-quality-assessment)\n",
    "        self.clouds = {'Cloud Shadow': [328, 392, 840, 904, 1350],\n",
    "                      'Cloud': [352, 368, 416, 432, 480, 864, 880, 928, 944, 992],\n",
    "                      'Cloud confidence': {'Low': [322, 324, 328, 336, 352, 368, 834, 836, 840, 848, 864, 880],\n",
    "                                           'Medium': [386, 388, 392, 400, 416, 432, 900, 904, 928, 944],\n",
    "                                           'High': [480, 992]},\n",
    "                      'Cirrus confidence': {'Low': [322, 324, 328, 336, 352, 368, 386, 388, 392, 400, 416, 432, 480],\n",
    "                                           'High': [834, 836, 840, 848, 864, 880, 898, 900, 904, 912, 928, 944, 992]}}\n",
    "        \n",
    "    def prepare_values(self, clear_key):\n",
    "        \"\"\"\n",
    "        Function prepares values to remove (cloud pixels) based on the input.\n",
    "        input 1: clear_key (string)\n",
    "        parameters:\n",
    "        'all' - returns values for Clouds and Cloud Shadows,\n",
    "        'Cloud Shadow' - returns values for Cloud Shadows,\n",
    "        'Cloud' - returns values for Clouds,\n",
    "        'Cloud low' - returns pixel values which could be a cloud with low, medium and high probability,\n",
    "        'Cloud medium' - returns pixel values which could be a cloud with medium and high probability,\n",
    "        'Cloud high' - returns pixel values which could be a cloud with high probability,\n",
    "        'Cirrus low' - returns pixel values which could be a cirrus with low and high probability,\n",
    "        'Cirrus high' - returns pixel values which could be a cirrus with high probability.\n",
    "        output: python list with pixel values to remove from the scene\"\"\"\n",
    "        \n",
    "        values = {'all': [self.clouds['Cloud Shadow'], self.clouds['Cloud']],\n",
    "                 'Cloud Shadow': [self.clouds['Cloud Shadow']],\n",
    "                 'Cloud': [self.clouds['Cloud']],\n",
    "                 'Cloud low': [self.clouds['Cloud confidence']['Low'],\n",
    "                               self.clouds['Cloud confidence']['Medium'],\n",
    "                               self.clouds['Cloud confidence']['High']],\n",
    "                 'Cloud medium': [self.clouds['Cloud confidence']['Medium'],\n",
    "                                  self.clouds['Cloud confidence']['High']],\n",
    "                 'Cloud high': [self.clouds['Cloud confidence']['High']],\n",
    "                 'Cirrus low': [self.clouds['Cirrus confidence']['Low'],\n",
    "                                self.clouds['Cirrus confidence']['High']],\n",
    "                 'Cirrus high': [self.clouds['Cirrus confidence']['High']]}\n",
    "        values_to_remove = values[clear_key]\n",
    "        \n",
    "        if len(values_to_remove) > 1:\n",
    "            vals = []\n",
    "            for values_list in values_to_remove:\n",
    "                vals = vals + values_list\n",
    "        else:\n",
    "            vals = values_to_remove[0]\n",
    "        return vals\n",
    "    \n",
    "    def clear_pixels(self, band, qa_input = None, clear_by='all', nofill_value=0.0):\n",
    "        \"\"\"\n",
    "        Function removes cloud pixels from the chosen scene based on the scene quality band.\n",
    "        input 1: band (numpy array) - band for the cloud removal.\n",
    "        input 2: qa_input - if None then self.qa_band is used as the mask for cloud pixels. Else mask should be\n",
    "        given as a numpy array.\n",
    "        input 3: clear_all - parameter to obtain values of cloud pixels:\n",
    "        'all' - returns values for Clouds and Cloud Shadows,\n",
    "        'Cloud Shadow' - returns values for Cloud Shadows,\n",
    "        'Cloud' - returns values for Clouds,\n",
    "        'Cloud low' - returns pixel values which could be a cloud with low, medium and high probability,\n",
    "        'Cloud medium' - returns pixel values which could be a cloud with medium and high probability,\n",
    "        'Cloud high' - returns pixel values which could be a cloud with high probability,\n",
    "        'Cirrus low' - returns pixel values which could be a cirrus with low and high probability,\n",
    "        'Cirrus high' - returns pixel values which could be a cirrus with high probability.\n",
    "        input 4: nofill_value - value which replaces cloud pixels.\n",
    "        output: numpy array (scene) without cloud pixels.\"\"\"\n",
    "        \n",
    "        values_to_remove = self.prepare_values(clear_by)\n",
    "        if qa_input is None:\n",
    "            qa_channel = self.qa_band\n",
    "        else:\n",
    "            qa_channel = qa_input\n",
    "        for val in values_to_remove:\n",
    "            qa_channel[qa_channel==val] = nofill_value\n",
    "        qa_channel[qa_channel > 0] = 1\n",
    "        output_band = band * qa_channel\n",
    "        return output_band"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Class gives you opportunity to load QA file initialized only by providing the directory\n",
    "# with it. Try it - for the same date as for the file with clouds. Then show qa_band instance.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now you know where clouds are and what is the name of the quality band. It will be good to\n",
    "# update function clip_all() to handle qa band everytime when you run it. I've added this\n",
    "# function below with new parameter: qa_band = True. (Parameter is set to True. We imply that\n",
    "# qa_band is always needed for the analysis).\n",
    "\n",
    "def clip_all(dataset, clipping_vector, bands=(2, 3, 4), qa_band=True, folder_name='data/clipped/'):\n",
    "    \"\"\"Function performs clipping of all Landsat 8 bands given in the DIR column of the df parameter.\n",
    "    input 1: dataset - dataframe with DIR column where each record points to the directory with a Landsat bands\n",
    "    input 2: clipping_vector - vector geometry which is used as a area of interest\n",
    "    input 3: bands=(2, 3, 4) - tuple of bands which should be clipped. Bands 2, 3, 4 are used for the BGR image creation.\n",
    "    input 4: folder_name='clipped/' - directory where all clipped bands should be stored.\n",
    "    \n",
    "    output: [correct, damaged] - list of dicts with keys (dates) and values - correct images\n",
    "    and damaged images. Damaged images were not process and they not exist in the destination folder.\"\"\"\n",
    "    \n",
    "    data_range = len(dataset)\n",
    "    damaged = {}\n",
    "    correct = {}\n",
    "    for i in range(0, data_range):\n",
    "        date_str = dataset.index[i]\n",
    "        images = read_landsat_images(dataset['DIR'][i])\n",
    "        \n",
    "        # Cut bands\n",
    "        for band_number in bands:\n",
    "            band_to_clip = images[band_number]\n",
    "            filename = folder_name + 'LC_band' + str(band_number) + '_' + date_str + '.tif'\n",
    "            try:\n",
    "                clip_area(clipping_vector, band_to_clip, filename)\n",
    "            except (rio.errors.WindowError, ValueError):\n",
    "                print('\\nPolygon does not cover raster area. Error occurs when projections are different or raster extent is translated and it does not overlap polygon extent.\\n')\n",
    "                try: \n",
    "                    damaged[date_str].append(filename)\n",
    "                except KeyError:\n",
    "                    damaged[date_str] = [filename]\n",
    "                pass\n",
    "            \n",
    "            # Check file\n",
    "            try: \n",
    "                damaged[date_str]\n",
    "                print('File: {} not saved. Raster has not been clipped.'.format(\n",
    "                filename))\n",
    "            except KeyError:\n",
    "                try:\n",
    "                    with rio.open(filename) as tf:\n",
    "                         pass\n",
    "                except rio.RasterioIOError:\n",
    "                    raise IOError('File was: {} not saved successfully'.format(\n",
    "                    filename))\n",
    "                else:\n",
    "                    print('File: {} saved successfully'.format(\n",
    "                    filename))\n",
    "                    try: \n",
    "                        correct[date_str].append(filename)\n",
    "                    except KeyError:\n",
    "                        correct[date_str] = [filename]\n",
    "        \n",
    "        # qa_band\n",
    "        if qa_band:\n",
    "            file_list = os.listdir(dataset['DIR'][i])\n",
    "            for f in file_list:\n",
    "                if (f.startswith('LC') and f.endswith('pixel_qa.tif')):\n",
    "                    filename = folder_name + 'LC_band_' + date_str + '_pixel_qa.tif'\n",
    "                    f = dataset['DIR'][i] + f\n",
    "                    try:\n",
    "                        clip_area(clipping_vector, f, filename)\n",
    "                    except (rio.errors.WindowError, ValueError):\n",
    "                        print('\\nPolygon does not cover raster area. Error occurs when projections are different or raster extent is translated and it does not overlap polygon extent.\\n')\n",
    "                        try: \n",
    "                            damaged[date_str].append(filename)\n",
    "                        except KeyError:\n",
    "                            damaged[date_str] = [filename]\n",
    "                        pass\n",
    "            \n",
    "                    # Check file\n",
    "                    try: \n",
    "                        damaged[date_str]\n",
    "                        print('File: {} not saved. Raster has not been clipped.'.format(filename))\n",
    "                    except KeyError:\n",
    "                        try:\n",
    "                            with rio.open(filename) as tf:\n",
    "                                 pass\n",
    "                        except rio.RasterioIOError:\n",
    "                            raise IOError('File was: {} not saved successfully'.format(filename))\n",
    "                        else:\n",
    "                            print('QA file: {} saved successfully'.format(filename))\n",
    "                            try: \n",
    "                                correct[date_str].append(filename)\n",
    "                            except KeyError:\n",
    "                                correct[date_str] = [filename]\n",
    "            \n",
    "    return [correct, damaged]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# What do we have now? It is a good time for little programming summary, so try below:\n",
    "# a) Prepare dataframe of Landsat 8 directories with RasterTimeSeries class and its methods\n",
    "\n",
    "\n",
    "# b) Show first five records of a DataFrame created in the previous step. Use pandas' df.head()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# c) Prepare new vector file (this file is correct for the task, explanation in the article).\n",
    "# Use correct_vec_address variable and VectorData class\n",
    "\n",
    "\n",
    "# d) Prepare 'all' vector parameters with get_vector_data() method and show it with\n",
    "# show_vector_data() function\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# e) Read vector's crs (vec_crs attribute)\n",
    "\n",
    "\n",
    "# f) And reproject vector geometry to the EPSG: 32634, do not update vector's instance\n",
    "# Instead create variable with vector geometry\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# g) Create geo_dict variable and clip all rasters with qa_band set to True. Use a geometry \n",
    "# from the previous point as a mask for clipping. Save files wherever you want, \n",
    "# but rather not to the folder with the jupyter notebooks.\n",
    "# Remember that your:\n",
    "# - dataset is a DataFrame with Landsat 8 directories,\n",
    "# - clipping_vector is a geometry from the previous steps\n",
    "# - bands are (1, 2, 3, 4, 5)\n",
    "# - qa_band must be written (we want to remove clouds and cloud shadows)\n",
    "# - destination folder must be there (you have two simple choices: update clip_all() \n",
    "# function with os.makedirs method - with checking if folder exists by os.path.exists \n",
    "# method or create directory manually)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# h) Prepare DataFrame with all of these values. Use method:\n",
    "# Use method pandas.DataFrame.from_dict(dictonary, orient='index', columns = [1, 2, 3, 4, 5, 'qa'])\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# i) Show the first ten records of this dataframe\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# At this point we have everything to perform analysis. We will test three different indices\n",
    "# to find the dates when algal bloom exists: Enhanced Vegetation Index, kab1, kab2. I've implemented\n",
    "# EVI and kab1. Your role is to implement function to calculate kab2.\n",
    "\n",
    "# Prepare function to calculate EVI index, kab1 index and kab2 index\n",
    "\n",
    "# EVI is calculated as (bands numbers in square brackets):\n",
    "# evi = 2.5 * ([B5] - [B4]) / ([B5] + 6 * [B4] - 7.5 * [B2] + 1)\n",
    "# equation above derived from: \n",
    "# https://landsat.usgs.gov/sites/default/files/documents/si_product_guide.pdf\n",
    "# EVI should be scaled by the factor of 0.00001\n",
    "# EVI must be in a range between -1 and 1. To avoid problems with further analysis EVI values\n",
    "# below are set to 0\n",
    "\n",
    "def calculate_evi(band2, band4, band5):\n",
    "    with rio.open(band5) as a:\n",
    "        nir_band = (a.read()[0]/10000).astype(np.float)\n",
    "    with rio.open(band4) as b:\n",
    "        r_band = (b.read()[0]/10000).astype(np.float)\n",
    "    with rio.open(band2) as c:\n",
    "        b_band = (c.read()[0]/10000).astype(np.float)\n",
    "        \n",
    "    numerator = nir_band - r_band\n",
    "    g = 2.5\n",
    "    c1 = 6\n",
    "    c2 = 7.5\n",
    "    l = 1\n",
    "    denominator = nir_band + c1 * r_band - c2 * b_band + l\n",
    "    evi = g * (numerator / denominator)\n",
    "    evi[evi > 1] = 1\n",
    "    return evi\n",
    "\n",
    "# Kab1 is calculated as (bands numbers in square brackets):\n",
    "# kab1 = 1.67 - 3.94 * ln[B2] + 3.78 * ln[B3]\n",
    "# You should be prepared for zeros and negative values (you are using logarithm here) so use\n",
    "# numpy.log1p() insted numpy.log()\n",
    "\n",
    "def calculate_kab1(band2, band3):\n",
    "    pass\n",
    "\n",
    "# Kab2 is calculated as (bands numbers in square brackets):\n",
    "# kab2 = 6.92274 - 5.7581 * (ln[B1] / ln[B3])\n",
    "# You should be prepared for zeros and negative values (you are using logarithm here)\n",
    "\n",
    "def calculate_kab2(band1, band3):\n",
    "    pass\n",
    "\n",
    "# Kab1 and Kab2 are derived from: Kabbara, N., J. Benkhelil, M. Awad, and V. Barale. 2008. \n",
    "# Monitoring water quality in the coastal area of Tripoli (Lebanon) using high-resolution \n",
    "# satellite data. ISPRS Journal of Photogrammetry and Remote Sensing 63:488–495"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now you have all indices. Try to calculate one of it for the date 2015-08-11.\n",
    "# Remember that your output is a numpy array. You could and should show it. Use 'winter'\n",
    "# as a color_map.\n",
    "# To select proper bands use DataFrame.loc[index_name, column_name]\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# It will be good to observe how changes mean of each value with the time. Prepare new columns\n",
    "# in data frame for the mean of each index. Name columns 'EVI', 'KAB1' and 'KAB2'.\n",
    "# Below is an example how to perform this task for the 'EVI' column. Your role is to do it\n",
    "# for 'KAB1' and 'KAB2'.\n",
    "\n",
    "# EVI\n",
    "df_with_bands['EVI'] = df_with_bands.apply(\n",
    "    lambda x: np.nanmean(calculate_evi(x[2], x[4], x[5])), axis=1)\n",
    "\n",
    "# Kab1\n",
    "\n",
    "\n",
    "# Kab2\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print kab1 and kab2 values using matplotlib.pyplot. Plot values as a points not as a lines.\n",
    "# Use example below (for evi). Plot also line for all values mean (you must calculate this mean).\n",
    "# The easiest way is to use numpy.mean() function.\n",
    "\n",
    "# evi\n",
    "evi_mean = [np.mean(df_with_bands['EVI']) for x in range(0, len(df_with_bands['EVI']))]\n",
    "\n",
    "plt.figure(figsize=(15, 3))\n",
    "plt.plot(df_with_bands['EVI'], 'ro')\n",
    "plt.plot(evi_mean, 'b--')\n",
    "plt.title('EVI over time')\n",
    "plt.legend(['evi', 'mean of evi'])\n",
    "plt.show()\n",
    "\n",
    "# kab1\n",
    "\n",
    "# kab2\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fast look into a plots shows that Kab1 and Kab2 are very similar indicators.\n",
    "# In that case we will use only one of them. Which one would you like to use?\n",
    "# My answer is in the article with explanation why.\n",
    "\n",
    "# We forget about very important processing step! We should remove clouds from which satellite\n",
    "# image. Clouds are affecting our calculations because cloud pixels are very bright and\n",
    "# in the other hand cloud shadows make pixels darker than they are in reality.\n",
    "\n",
    "# To fix this we will use PixelQA class and before we calculate mean we will delete cloud and\n",
    "# cloud shadow pixels from a scene. You may use remove_cloud() function with two parameters:\n",
    "# output band (evi, kab1 or kab2) and quality map for a given index. The example how to use\n",
    "# it is given for a kab1. Try to do it for evi! (Set remove_negative parameter in show_band\n",
    "# function to False)\n",
    "\n",
    "def remove_clouds(path_to_cloud_band, band_to_clean):\n",
    "    qa_band = PixelQA(path_to_cloud_band)\n",
    "    cleaned_band = qa_band.clear_pixels(band_to_clean)\n",
    "    return cleaned_band\n",
    "\n",
    "# Kab1\n",
    "\n",
    "test_kab1 = calculate_kab1(\n",
    "    df_with_bands.loc['2013-05-17', 2], df_with_bands.loc['2013-05-17', 3])\n",
    "test_kab1 = remove_clouds(\n",
    "    df_with_bands.loc['2013-05-17', 'qa'], test_kab1)\n",
    "show_band(test_kab1, color_map='winter')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# evi\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now you know how to remove clouds and finally you may calculate mean for each index.\n",
    "# Calculate mean for EVI and kab1 and name new columns as 'EVI-cloudless' and 'KAB1-cloudless'.\n",
    "# Use examples above as the reference.\n",
    "# The important point is that we must count mean of values greater than 0 because after\n",
    "# cloud removing operation we set all unwanted pixels to 0. Use for it\n",
    "# calculate_nanmean_of_positive() function given below.\n",
    "\n",
    "def calculate_nanmean_of_positive(idx_func, idx_bands, qa_band):\n",
    "    retrieved_band = idx_func(*idx_bands)\n",
    "    cloudless = remove_clouds(qa_band, retrieved_band)\n",
    "    value = np.nanmean(cloudless[cloudless > 0])\n",
    "    return value\n",
    "\n",
    "\n",
    "# EVI-cloudless\n",
    "\n",
    "# KAB1-cloudless\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show EVI and EVI-cloudless on the same plot as points with red and green colors. On the same\n",
    "# plot draw yellow line represnting mean of EVI and blue line representing mean of EVI-cloudless.\n",
    "# Show KAB1 and KAB1-cloudless on the same plot as points with red and green colors. On the same\n",
    "# plot draw yellow line represnting mean of KAB1 and blue line representing mean of KAB1-cloudless.\n",
    "\n",
    "# evi\n",
    "\n",
    "# kab1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Analysis:\n",
    "# what is the simplest method of classification of this dataset? What do you think?\n",
    "# We use simple threshold to find a time when algal blooms have occured.\n",
    "# But this is a classic example of problem which can be solved which machine learning.\n",
    "# At the end of this course we will tackle it with a simplest approach and in the future you may\n",
    "# expect more complicated methods.\n",
    "\n",
    "# We will set a threshold to the mean of the bands means. The first example is for EVI.\n",
    "# Try to do the same for the kab1.\n",
    "\n",
    "# EVI\n",
    "\n",
    "threshold = np.mean(df_with_bands['EVI-cloudless'])\n",
    "algal_blooms_by_evi = df_with_bands[df_with_bands['EVI-cloudless'] > threshold]\n",
    "algal_blooms_by_evi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check... the best option is to look into a EVI image itself\n",
    "\n",
    "def perform_visual_test_evi(df_of_interest):\n",
    "    test_band = df_of_interest.apply(lambda x: calculate_evi(x[2],\n",
    "                                                            x[4],\n",
    "                                                            x[5]), axis=1)\n",
    "    indexes = test_band.index.values\n",
    "    for i in indexes:\n",
    "        print('Date:', i)\n",
    "        show_band(test_band[i], color_map='viridis', remove_negative=False)\n",
    "    \n",
    "\n",
    "perform_visual_test_evi(algal_blooms_by_evi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# And kab1 results\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Checking... the best option is to look into a kab1 image itself\n",
    "\n",
    "def perform_visual_test_kab1(df_of_interest):\n",
    "    test_band = df_of_interest.apply(lambda x: calculate_kab1(x[2],\n",
    "                                                            x[3]), axis=1)\n",
    "    indexes = test_band.index.values\n",
    "    for i in indexes:\n",
    "        print('Date:', i)\n",
    "        show_band(test_band[i], color_map='viridis', remove_negative=False)\n",
    "    \n",
    "\n",
    "perform_visual_test_kab1()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# And what do you think about it?"
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}