Detect_Trucks_Sentinel2_Upscaling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Detect Trucks Sentinel-2 - Europe\n",
    "________________                                            "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This script is a processing chain for Sentinel-2 truck detection on large scale. It usses EuroDataCube and Sentinel Hub resources. \n",
    "\n",
    "The script covers roughly the following steps:\n",
    "- Generate a processing grid\n",
    "- Get OSM road data through Overpass API\n",
    "- Mask clouds in acquisition\n",
    "- Mask bands to OSM road raster\n",
    "- Detect trucks\n",
    "- Calculate detected trucks per observation in a period\n",
    "\n",
    "All processing parameters can be adjusted from the top in the 'General Parameters' section.\n",
    "\n",
    "Author: Henrik Fisser, 2020"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load creds\n",
    "%load_ext dotenv\n",
    "%dotenv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_______________\n",
    "## 1 | Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# general\n",
    "import os\n",
    "import datetime\n",
    "import pandas as pd\n",
    "import subprocess\n",
    "import sys\n",
    "import shutil\n",
    "import bottleneck\n",
    "import time\n",
    "from datetime import date, datetime, timedelta\n",
    "from time import sleep\n",
    "from glob import glob\n",
    "\n",
    "# installations\n",
    "def install_package(pkg):\n",
    "    subprocess.check_call([sys.executable, \"-m\", \"pip\", \"install\", pkg])\n",
    "install_package(\"OSMPythonTools\")\n",
    "install_package(\"geocube\")\n",
    "\n",
    "# OSM API\n",
    "from OSMPythonTools.overpass import overpassQueryBuilder, Overpass\n",
    "\n",
    "# AWS Bucket\n",
    "import s3fs\n",
    "import requests\n",
    "import boto3\n",
    "import botocore\n",
    "bucket = \"edc-b62330c3-ed93-4b41-a90f-28f3a73107e9\"\n",
    "\n",
    "# xcube\n",
    "from xcube_sh.cube import open_cube\n",
    "from xcube_sh.config import CubeConfig\n",
    "from xcube.core.maskset import MaskSet\n",
    "\n",
    "# numba\n",
    "from numba import jit\n",
    "\n",
    "# spatial\n",
    "import xarray as xr\n",
    "import numpy as np\n",
    "import geopandas as gpd\n",
    "import geocube\n",
    "from rasterio import features\n",
    "from affine import Affine\n",
    "from shapely import geometry, coords\n",
    "from shapely.geometry import Polygon, Point\n",
    "#from osgeo import gdal #, gdal_array, ogr\n",
    "\n",
    "# plotting\n",
    "import matplotlib as plt\n",
    "import IPython.display\n",
    "%matplotlib inline\n",
    "figsize = [10,10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "____\n",
    "## 2 | General Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "overwrite_results = True\n",
    "country = \"Sweden\"\n",
    "obs_full_threshold = 98 # consider mosaic as full when reached\n",
    "n_obs = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Cube"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = \"S2L2A\"\n",
    "spatial_res = 0.00009 # approx. 10m\n",
    "bands = [\"B02\", \"B03\", \"B04\", \"B08\", \"B11\", \"SCL\"]\n",
    "day_bins = \"1D\" # for cube\n",
    "tile_size = [512, 512]\n",
    "minimum_valid_observations = 15 # percent\n",
    "grid_spacing = 1. # processing grid box size [degree]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Directories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_main = os.getcwd()\n",
    "dir_not_commit = os.path.join(dir_main, \"not_commit\")\n",
    "dir_ancil = os.path.join(dir_not_commit, \"ancillary_data\")\n",
    "dirs = {\"dir\":dir_main, \"dir_not_commit\":dir_not_commit, \"ancil\":dir_ancil, \"processing\":os.path.join(dir_main, \"processing\"),\n",
    "       \"processed\":os.path.join(dir_not_commit, \"processed\"), \"ancil_roads\":os.path.join(dir_ancil, \"roads\"), \"ancil_gadm\":os.path.join(dir_ancil, \"gadm\")}\n",
    "for directory in list(dirs.values()):\n",
    "    if not os.path.exists(directory): os.mkdir(directory)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "files = {\"gadm_efta\":os.path.join(dirs[\"ancil_gadm\"], \"gadm36_0_EU_EFTA.gpkg\"), \n",
    "         \"gadm_europe\":os.path.join(dirs[\"ancil_gadm\"], \"gadm36_0_europe.gpkg\"), \n",
    "         \"gadm_europe_union\":os.path.join(dirs[\"ancil_gadm\"], \"gadm0_europe_union.gpkg\"),\n",
    "         \"proc_grid\":os.path.join(dirs[\"processing\"], \"processing_grid_gadm_eu27.geojson\")}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Temporal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "weekdays = [\"tuesday\", \"wednesday\", \"thursday\"] \n",
    "target = {\"first\":datetime(2020, 4, 1), \"last\":datetime(2020, 6, 21)}\n",
    "timestamps_sub_period = 1 # how many aggregated timestamps per sub-period\n",
    "baseline_years = [2017, 2018, 2019] # if processing for a yearly period as in target"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### OSM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "osm_key = \"highway\"\n",
    "osm_values = [\"motorway\", \"trunk\", \"primary\"]\n",
    "roads_buffer = 0.00022 # degree, for motorway, the others lower (see OSM methods below)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Cloud masking"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "cloud_masking_thresholds = {\"rgb\":0.25,\n",
    "                            \"blue_green\":0.2,\n",
    "                            \"blue_red\":0.2}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "______\n",
    "## 3 | Detection Parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thresholds were determined based on a manually selected dataset of 150 truth points of a test acquisition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "thresholds = {\"min_blue\":0.06, \n",
    "              \"min_green\":0.04, \n",
    "              \"min_red\":0.04,\n",
    "              \"max_red\":0.15,\n",
    "              \"max_green\":0.15,\n",
    "              \"max_blue\":0.2,\n",
    "              \"max_ndvi\":0.5,\n",
    "              \"max_ndwi\":0.0001,\n",
    "              \"max_ndsi\":0.0001,\n",
    "              \"min_blue_green_ratio\":0.03,\n",
    "              \"min_blue_red_ratio\":0.05, \n",
    "              \"max_blue_green_ratio\":0.17, \n",
    "              \"max_blue_red_ratio\":0.2}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "________________\n",
    "## 4 | Utils"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def EPSG_4326(): return \"EPSG:4326\"\n",
    "def EPSG_3857(): return \"EPSG:3857\"\n",
    "def GEOJSON(): return \"GeoJSON\"\n",
    "def GPKG(): return \"GPKG\"\n",
    "def GEOJSON_EXT(): return \".geojson\"\n",
    "def GPKG_EXT(): return \".gpkg\"\n",
    "def NC_EXT(): return \".nc\"\n",
    "def BBOX_ID(): return \"bbox_id\"\n",
    "def ABS_N_TRUCKS(): return \"acquisitions_trucks\"\n",
    "def TRUCKS_VEC(): return \"trucks_points\"\n",
    "def TRUCKS_VEC_PHR(): return \"trucks_points_placeholder\"\n",
    "def TRUCKS_VEC_MERGED(): return \"trucks_points_merged\"\n",
    "def N_OBS(): return \"n_observations\"\n",
    "def VALID_PIXELS(): return \"valid_pixels\"\n",
    "def TRUCK_COUNT_NORM(): return \"truck_count_normalized\"\n",
    "def SUM_TRUCKS_SUB(): return \"sum_trucks_sub_period\"\n",
    "def MEAN_TRUCKS(): return \"mean_trucks\"\n",
    "\n",
    "def name_bucket(): return \"edc-b62330c3-ed93-4b41-a90f-28f3a73107e9\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### File names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fname_osm(directory, bbox_id, osm_key, ext = GPKG_EXT()): return os.path.join(directory, str(bbox_id) + \"_\" + osm_key + ext)\n",
    "def construct_fname(dirs_ts, dir_ts_key, bbox_id, ext):\n",
    "    return os.path.join(dirs_ts[dir_ts_key], dir_ts_key + \"_\" + BBOX_ID() + str(bbox_id) + ext)\n",
    "def fname_acquisition_trucks(dirs_ts, bbox_id, date, ext = NC_EXT()): \n",
    "    return os.path.join(dirs_ts[ABS_N_TRUCKS()], str(date) + \"_\" + ABS_N_TRUCKS() + \"_\" + BBOX_ID() + str(bbox_id) + ext)\n",
    "def fname_sum_obs(dirs_ts, bbox_id, date, ext = NC_EXT()): \n",
    "    return construct_fname(dirs_ts, N_OBS(), str(bbox_id) + \"_\" + str(date), ext)\n",
    "def fname_trucks_vec(dirs_ts, bbox_id, date, ext):\n",
    "    return construct_fname(dirs_ts, TRUCKS_VEC(), str(bbox_id) + \"_\" + str(date), ext)\n",
    "def fname_trucks_vec_placeholder(dirs_ts, bbox_id, date, ext = \".txt\"):\n",
    "    return construct_fname(dirs_ts, TRUCKS_VEC(), str(bbox_id) + \"_\" + str(date), ext)\n",
    "\n",
    "def sub_period_fnames(dirs_ts, bbox_id, date, ext_arr, ext_vec):\n",
    "    files = {N_OBS():fname_sum_obs(dirs_ts, bbox_id, date, ext_arr),\n",
    "             TRUCKS_VEC():fname_trucks_vec(dirs_ts, bbox_id, date, ext_vec),\n",
    "             TRUCKS_VEC_PHR():fname_trucks_vec_placeholder(dirs_ts, bbox_id, date)}\n",
    "    return files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Directory structure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_dirs_ts(dir_ts):\n",
    "    dir_ts_overall = os.path.join(dir_ts, \"overall\")\n",
    "    if not os.path.exists(dir_ts_overall): os.mkdir(dir_ts_overall)\n",
    "    dirs_ts = {ABS_N_TRUCKS():os.path.join(dir_ts, ABS_N_TRUCKS()),\n",
    "               N_OBS():os.path.join(dir_ts_overall, N_OBS()),\n",
    "               TRUCKS_VEC():os.path.join(dir_ts_overall, TRUCKS_VEC()),\n",
    "               TRUCKS_VEC_MERGED():os.path.join(dir_ts_overall, TRUCKS_VEC_MERGED())}\n",
    "    for direc in dirs_ts.values():\n",
    "        if not os.path.exists(direc): os.mkdir(direc)\n",
    "    return dirs_ts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def delete_recurs(directory, validation):\n",
    "    if validation == \"YES_DELETE_THAT\":\n",
    "        shutil.rmtree(directory)\n",
    "    else:\n",
    "        raise Exception(\"Validation wrong\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Processing messages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def start_msg(weekdays, n_boxes, periods, timestamps_sub_period, osm_values, baseline_years, minimum_valid_observations, overwrite_results):\n",
    "    print(\"%s\\n\\nStarting truck detection processing \\n%s\\n\\nWeekdays: \\\n",
    "%s\\nNumber of grid cells to process: %s\\nNumber of periods: %s\\nTimestamps sub-period: \\\n",
    "%s\\nOSM roads: %s\\nBaseline years: %s\\nMinimum valid observations to consider acquisition: \\\n",
    "%s %%\\nOverwrite results: %s\\n%s\\n\\n%s\" %(\"=\"*100, \".\"*50, weekdays, str(n_boxes), str(len(periods[\"first\"])), \n",
    "                                          str(timestamps_sub_period), str(osm_values), str(baseline_years), \n",
    "                                          str(minimum_valid_observations), str(overwrite_results), \".\"*50, \"=\"*100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Processing grid\n",
    "Create a grid covering Europe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_grid(grid_spacing, files):\n",
    "    crs = EPSG_4326()\n",
    "    gadm0_eu_efta = gpd.read_file(files[\"gadm_efta\"])\n",
    "    xmin,ymin,xmax,ymax = gadm0_eu_efta.total_bounds\n",
    "    width = xmax - xmin\n",
    "    height = ymax - ymin\n",
    "    cols = int((width) / grid_spacing)\n",
    "    rows = int((height) / grid_spacing)\n",
    "    box_width = width / cols\n",
    "    box_height = height / rows\n",
    "    boxes = []\n",
    "    for row in range(rows):\n",
    "        for col in range(cols):\n",
    "            col_right = col + 1\n",
    "            row_lower = row + 1\n",
    "            y_up = ymax-row*box_height\n",
    "            y_low = ymax-row_lower*box_height\n",
    "            x_left = xmin+col*box_width\n",
    "            x_right = xmin+col_right*box_width\n",
    "            boxes.append(Polygon([(x_left, y_up), \n",
    "                                  (x_right, y_up), \n",
    "                                  (x_right, y_low), \n",
    "                                  (x_left, y_low)]))\n",
    "    grid = gpd.GeoDataFrame({\"geometry\":boxes})\n",
    "    grid.crs = crs\n",
    "    gadm0_europe = gpd.read_file(files[\"gadm_europe\"])\n",
    "    gadm0_europe.crs = crs\n",
    "    file_union = files[\"gadm_europe_union\"]\n",
    "    if os.path.exists(file_union):\n",
    "        gadm0_europe_clip_union = gpd.read_file(file_union)\n",
    "    else:\n",
    "        gadm0_europe[\"continent\"] = [\"EUROPE\"] * len(gadm0_europe)\n",
    "        gadm0_europe_clip = gpd.overlay(gadm0_europe, gpd.GeoDataFrame({\"a\":[1],\"geometry\":test}), how = \"intersection\")\n",
    "        gadm0_europe_clip_union = gadm0_europe_clip.dissolve(by = \"continent\")\n",
    "        gadm0_europe_clip_union.to_file(file_union, driver = GPKG())\n",
    "    # intersect with gadm for cutting boxes at coast line\n",
    "    grid_intersect = gpd.overlay(grid, gadm0_europe_clip_union, how = \"intersection\")\n",
    "    # get country attributes\n",
    "    grid_gadm = gpd.sjoin(grid_intersect, gadm0_europe, how = \"left\", op = \"intersects\")\n",
    "    grid_gadm[BBOX_ID()] = range(len(grid_gadm))\n",
    "    boxes = grid_gadm.geometry.apply(lambda geom : geom.bounds)\n",
    "    grid_gadm.geometry = [Polygon(geometry.box(b[0], b[1], b[2], b[3])) for b in boxes] # use the light bboxes\n",
    "    delete_cols = [\"GID_0_left\", \"NAME_0_left\", \"a\", \"continent\", \"index_right\"]\n",
    "    for column in delete_cols:\n",
    "        if column in grid_gadm.columns: \n",
    "            grid_gadm = grid_gadm.drop(column, 1)\n",
    "    grid_gadm = grid_gadm.rename(columns = {\"GID_0_right\":\"GID_0\", \"NAME_0_right\":\"NAME_0\"})\n",
    "    grid_gadm = grid_gadm[grid_gadm[\"GID_0\"] != \"RUS\"] # do not include\n",
    "    grid_gadm.to_file(files[\"proc_grid\"], driver = GEOJSON())\n",
    "    return grid_gadm\n",
    "\n",
    "def unique_grid(file_grid_gadm):\n",
    "    grid_gadm = gpd.read_file(file_grid_gadm)\n",
    "    unique = []\n",
    "    unique_geoms = []\n",
    "    for i, geom in enumerate(grid_gadm.geometry):\n",
    "        if grid_gadm.GID_0[i] == \"GER\":\n",
    "            unique.append(i)\n",
    "            unique_geoms.append(geom)\n",
    "    for i, geom in enumerate(grid_gadm.geometry):\n",
    "        if geom not in unique_geoms:\n",
    "            unique.append(i)\n",
    "            unique_geoms.append(geom)\n",
    "    grid_gadm = grid_gadm.iloc[unique,:]\n",
    "    grid_gadm.to_file(file_grid_gadm, driver=GEOJSON())\n",
    "    return grid_gadm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Vectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def transform_lat_lon(lat, lon):\n",
    "    lat = np.asarray(lat)\n",
    "    lon = np.asarray(lon)\n",
    "    trans = Affine.translation(lon[0], lat[0])\n",
    "    scale = Affine.scale(lon[1] - lon[0], lat[1] - lat[0])\n",
    "    return trans * scale\n",
    "\n",
    "def rasterize(polygons, lat, lon, fill=np.nan):\n",
    "    transform = transform_lat_lon(lat, lon)\n",
    "    out_shape = (len(lat), len(lon))\n",
    "    raster = features.rasterize(polygons.geometry, out_shape=out_shape,\n",
    "                                fill=fill, transform=transform,\n",
    "                                dtype=float)\n",
    "    return xr.DataArray(raster, coords={\"lat\":lat, \"lon\":lon}, dims=(\"lat\", \"lon\"))\n",
    "\n",
    "def vec_driver_from_ext(ext):\n",
    "    return {GPKG_EXT():GPKG(), GEOJSON_EXT():GEOJSON()}[ext]\n",
    "\n",
    "def points_to_raster(points, reference_raster):\n",
    "    arr = np.zeros((reference_raster.shape[0], reference_raster.shape[1]))\n",
    "    for i in range(len(points)):\n",
    "        geom = points[i:i+1].geometry\n",
    "        y,x = float(geom.y), float(geom.x)\n",
    "        index_y = np.where(reference_raster.lat == y)\n",
    "        index_x = np.where(reference_raster.lon == x)\n",
    "        arr[index_y, index_x] = 1\n",
    "    return arr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data np array\n",
    "# lat_lon dict of \"lat\" and \"lon\" holding np arrays of coordinates\n",
    "def create_xr_dataset(data, lat_lon, name):\n",
    "    lat, lon = \"lat\", \"lon\"\n",
    "    data_array = xr.DataArray(data, coords={lat:lat_lon[lat], lon:lat_lon[lon]}, dims=(lat, lon))\n",
    "    return xr.Dataset({name:data_array})\n",
    "\n",
    "# extracts coordinates at value in np array and returns points as GeoDataFrame\n",
    "# data 2d np array\n",
    "# match_value Float value in data where point coordinates are extracted\n",
    "# lon_lat dict of:\n",
    "### \"lon\": np array longitude values\"\n",
    "### \"lat\": np array latitude values\"\n",
    "# crs String EPSG:XXXX\n",
    "def points_from_np(data, match_value, lon_lat, crs):\n",
    "    indices = np.argwhere(data == match_value)\n",
    "    if len(indices) > 0:\n",
    "        lat_indices = indices[:,[0]]\n",
    "        lon_indices = indices[:,[1]]\n",
    "        lat_coords = lon_lat[\"lat\"][lat_indices]\n",
    "        lon_coords = lon_lat[\"lon\"][lon_indices]\n",
    "        points = gpd.GeoDataFrame(geometry = gpd.points_from_xy(lon_coords, lat_coords))\n",
    "        points.crs = crs\n",
    "        return points\n",
    "    \n",
    "def raster_to_points(raster, lon_lat, field_name, crs):\n",
    "    points_list = []\n",
    "    match_values = np.unique(raster[(raster != 0) * ~np.isnan(raster)]) # by pixel value\n",
    "    for x in match_values:\n",
    "        points = points_from_np(raster, x, lon_lat, crs=crs)\n",
    "        points[field_name] = [x] * len(points)\n",
    "        points_list.append(points)\n",
    "    return gpd.GeoDataFrame(pd.concat(points_list, ignore_index=True))\n",
    "\n",
    "# take xarray and ensure each value with 1 in data has no neighbor with 1 in an extended 3x3 block. Extended means: horizontally and vertically\n",
    "# it is also checked for the second-next pixel\n",
    "# Method checks only surrounding of values equal 1\n",
    "# arr xarray DataArray with values and lat lon\n",
    "def filter_spatial_3x3_extended(arr):\n",
    "    values = arr.values\n",
    "    lon = arr.lon\n",
    "    lat = arr.lat\n",
    "    valid = np.where(arr == 1)\n",
    "    for y,x in zip(valid[0], valid[1]):\n",
    "        y_above = y - 1\n",
    "        y_above_next = y - 2\n",
    "        x_left = x - 1\n",
    "        x_right = x + 1\n",
    "        x_left_next = x - 2\n",
    "        space_left = x_left >= 0\n",
    "        space_right = x_right >= 0 and x_right < len(lon)\n",
    "        space_above = y_above >= 0\n",
    "        val_left_above = values[y_above, x_left] if space_left and space_above else 0\n",
    "        val_right_above = values[y_above, x_right] if space_right and space_above else 0\n",
    "        val_left = values[y, x_left] if space_left else 0\n",
    "        val_above = values[y_above, x] if space_above else 0\n",
    "        val_left_next = values[y, x_left_next] if x_left_next >= 0 else 0\n",
    "        val_above_next = values[y_above_next, x] if y_above_next >= 0 else 0\n",
    "        # if any of the values left, above and left above has 1 set current value 0\n",
    "        if (val_left_above + val_right_above + val_left + val_above + val_left_next + val_above_next) >= 1:\n",
    "            values[y,x] = 0\n",
    "    arr.values = values\n",
    "    return arr\n",
    "\n",
    "def calc_rgb_cloud_mask(band_stack, cloud_masking_thresholds):\n",
    "    B02, B03, B04 = band_stack.B02, band_stack.B03, band_stack.B04\n",
    "    c = cloud_masking_thresholds[\"rgb\"]\n",
    "    clouds_rgb = ((B02 > c) + (B03 > c) + (B04 > c)) >= 1\n",
    "    # attempt to mask haze without masking out truck pixels (similar! higher blue than red and green)\n",
    "    blue_green_ratio = (B02-B03) / (B02+B03)\n",
    "    blue_red_ratio = (B02-B04) / (B02+B04)\n",
    "    clouds_blue_green = blue_green_ratio > cloud_masking_thresholds[\"blue_green\"]\n",
    "    clouds_blue_red = blue_red_ratio > cloud_masking_thresholds[\"blue_red\"]\n",
    "    clouds = (clouds_rgb + clouds_blue_green + clouds_blue_red) >= 1\n",
    "    return clouds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Temporal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# date Datetime object to be checked\n",
    "# weekday String weekday to be checked against\n",
    "# returns Boolean if date is weekday\n",
    "def is_weekday(date_x, weekday):\n",
    "    if not isinstance(date_x, type(datetime.date)): date_x = pd.to_datetime(date_x).date()\n",
    "    weekday = weekday.lower()[0:3]\n",
    "    y, m, d = 2000, 1, 3\n",
    "    wd = [\"mon\", \"tue\", \"wed\", \"thu\", \"fri\", \"sat\", \"sun\"]\n",
    "    ref = {}\n",
    "    for i in range(len(wd)): ref[wd[i]] = datetime(y, m, d + i).date()\n",
    "    return (date_x - ref[weekday]).days % 7 == 0 # check if date is in a sequence of 7\n",
    "\n",
    "# Calculates the dates of the periods representing a timestamp each\n",
    "# n_days_sub Integer how many dates shall one timestamp (sub-period) cover\n",
    "# baseline Dict start and end dates (datetime) overall baseline period\n",
    "# target Dict start and end dates (datetime) target period\n",
    "def calc_periods(n_days_sub, baseline, target):\n",
    "    fst, lst = \"first\", \"last\"\n",
    "    base_first = baseline[fst]\n",
    "    n_days = target[fst] - base_first\n",
    "    n_sub = int(n_days.days / n_days_sub) # number of timestamps\n",
    "    start = [base_first] * n_sub\n",
    "    start = [start[i] + timedelta(n_days_sub * i) for i in range(len(start))]\n",
    "    end = [start[i+1] - timedelta(1) for i in range(len(start)-1)]\n",
    "    end.append(baseline[lst])\n",
    "    periods = {\"ts\":range(len(start)), fst:start, lst:end}\n",
    "    periods[\"first\"].append(target[\"first\"]) # append start date of target period\n",
    "    periods[\"last\"].append(target[\"last\"]) # append end date of target period\n",
    "    return periods\n",
    "\n",
    "# Returns periods equivalent to target period for other years\n",
    "# target Dict start and end dates (datetime) target period\n",
    "# years List of int years\n",
    "# timestamps_sub_period int\n",
    "def yearly_period_from_target(target, years = [2017, 2018, 2019], timestamps_sub_period = 2):\n",
    "    n_days_target = (target[\"last\"]-target[\"first\"]).days\n",
    "    n_days_timestamp = int(n_days_target / timestamps_sub_period)\n",
    "    trgt_fst = []\n",
    "    trgt_lst = []\n",
    "    for i in range(timestamps_sub_period):\n",
    "        start = target[\"first\"] + timedelta(days = i * (n_days_timestamp + 1))\n",
    "        trgt_fst.append(start)\n",
    "        trgt_lst.append(start + timedelta(days = n_days_timestamp))\n",
    "    start = []\n",
    "    end = []\n",
    "    for year in years:\n",
    "        for fst, lst in zip(trgt_fst, trgt_lst):\n",
    "            start.append(datetime(year, fst.month, fst.day))\n",
    "            end.append(datetime(year, lst.month, lst.day))\n",
    "    appended = [start.append(x) for x in trgt_fst]\n",
    "    appended = [end.append(x) for x in trgt_lst]\n",
    "    periods = {\"ts\":range(len(start)), \"first\":start, \"last\":end}\n",
    "    return periods\n",
    "\n",
    "# date_str String YYYY-MM-DD\n",
    "def string_to_datetime(date_str):\n",
    "    y, m, d = date_str[0:4], date_str[5:7], date_str[-2:]\n",
    "    if m[0] == \"0\": m = m[1]\n",
    "    if d[0] == \"0\": d = d[1]\n",
    "    return datetime(int(y), int(m), int(d))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### OSM data\n",
    "Utils for retrieving road data from OSM through its API"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# re-order bbox from W,S,E,N to S,W,N,E\n",
    "def convert_bbox_osm(bbox):\n",
    "    offset = 0.05 # add a buffer to bbox in order to be sure cube is entirely covered\n",
    "    bbox_osm = [bbox[1], bbox[0], bbox[3], bbox[2]]\n",
    "    bbox_osm[0] -= offset # min lat\n",
    "    bbox_osm[1] -= offset # min lon\n",
    "    bbox_osm[2] += offset # max lat\n",
    "    bbox_osm[3] += offset # max lon\n",
    "    return bbox_osm\n",
    "\n",
    "# bbox List of four coords\n",
    "# bbox_id Integer processing id of bbox\n",
    "# osm_value String OSM value\n",
    "# osm_key String OSM key\n",
    "# element_type List of String\n",
    "# returns GeoPandasDataFrame\n",
    "def get_osm(bbox, \n",
    "            bbox_id,\n",
    "            osm_value = \"motorway\",\n",
    "            osm_key = \"highway\", # in OSM 'highway' contains several road types: https://wiki.openstreetmap.org/wiki/Key:highway\n",
    "            element_type = [\"way\", \"relation\"]):\n",
    "    \n",
    "    bbox_osm = convert_bbox_osm(bbox)\n",
    "    quot = '\"'\n",
    "    select = quot+osm_key+quot + '=' + quot+osm_value+quot\n",
    "    select_link = select.replace(osm_value, osm_value + \"_link\") # also get road links\n",
    "    select_junction = select.replace(osm_value, osm_value + \"_junction\")\n",
    "    geoms = []\n",
    "    for selector in [select, select_link, select_junction]:  \n",
    "        try:\n",
    "            query = overpassQueryBuilder(bbox=bbox_osm, \n",
    "                                         elementType=element_type, \n",
    "                                         selector=selector, \n",
    "                                         out='body',\n",
    "                                         includeGeometry=True)\n",
    "            elements = Overpass().query(query, timeout=60).elements()\n",
    "            # create multiline of all elements\n",
    "            if len(elements) > 0:\n",
    "                for i in range(len(elements)):\n",
    "                    elem = elements[i]\n",
    "                    geoms.append(elem.geometry())\n",
    "        except:\n",
    "            Warning(\"Could not retrieve \" + select)\n",
    "    try:\n",
    "        lines = gpd.GeoDataFrame(crs = EPSG_4326(), geometry = geoms)\n",
    "        n = len(geoms)\n",
    "        lines[BBOX_ID()] = [bbox_id]*n\n",
    "        lines[\"osm_value\"] = [osm_value]*n # add road type\n",
    "        return lines\n",
    "    except:\n",
    "        Warning(\"Could not merge \" + osm_value)\n",
    "        \n",
    "# buffer Float road buffer distance [m]\n",
    "# bbox List of four coords\n",
    "# bbox_id Integer processing id of bbox\n",
    "# osm_values List of String OSM values\n",
    "# osm_key String OSM key\n",
    "# roads_buffer Float buffer width\n",
    "# dir_write\n",
    "def get_roads(bbox, bbox_id, osm_values, osm_key, roads_buffer, dir_write):\n",
    "    fwrite = fname_osm(dir_write, bbox_id, osm_key)\n",
    "    if not os.path.exists(fwrite):\n",
    "        roads = []\n",
    "        has_error = []\n",
    "        offset = 0.00002\n",
    "        buffer_dist = \"buffer_distance\"\n",
    "        # buffer according to road type\n",
    "        m,t,p,s,ter = \"motorway\", \"trunk\", \"primary\", \"secondary\", \"tertiary\"\n",
    "        buffers = {m:roads_buffer, t:roads_buffer-offset, p:roads_buffer-(2*offset), s:roads_buffer-(3*offset), ter:roads_buffer-(4*offset)}\n",
    "        osm_values_int = {m:1, t:2, p:3, s:4, ter:5}\n",
    "        for osm_value in osm_values:\n",
    "            try:\n",
    "                roads_osm = get_osm(bbox = bbox, bbox_id = bbox_id, osm_value = osm_value)\n",
    "                roads_osm[buffer_dist] = [buffers[osm_value]] * len(roads_osm)\n",
    "                roads_osm[\"osm_value_int\"] = osm_values_int[osm_value]\n",
    "                roads.append(roads_osm)\n",
    "            except:\n",
    "                has_error.append(1)\n",
    "                Warning(\"'get_osm'\" + \"failed for bbox_id \"+ str(bbox_id) + \"osm_value \" + osm_value + \"osm_key\" + osm_key)\n",
    "        if len(roads) > len(has_error):\n",
    "            roads_merge = gpd.GeoDataFrame(pd.concat(roads, ignore_index=True), crs=roads[0].crs)\n",
    "            buffered = roads_merge.buffer(distance=roads_merge[buffer_dist])\n",
    "            roads_merge.geometry = buffered\n",
    "            roads_merge.to_file(fwrite, driver = GPKG())\n",
    "    return fwrite\n",
    "\n",
    "# osm geodataframe of polygons\n",
    "# reference_raster xarray with lat and lon\n",
    "def rasterize_osm(osm, reference_raster):\n",
    "    osm_values = list(set(osm[\"osm_value\"]))\n",
    "    nan_placeholder = 100\n",
    "    road_rasters = []\n",
    "    for osm_value in osm_values:\n",
    "        osm_subset = osm[osm[\"osm_value\"] == osm_value]\n",
    "        raster = rasterize(osm_subset, reference_raster.lat, reference_raster.lon)\n",
    "        cond = np.isfinite(raster)\n",
    "        raster_osm = np.where(cond, list(osm_subset.osm_value_int)[0], nan_placeholder) # use placeholder instead of nan first\n",
    "        raster_osm = raster_osm.astype(np.float)\n",
    "        road_rasters.append(raster_osm)        \n",
    "    # merge road types in one layer\n",
    "    road_raster_np = np.array(road_rasters).min(axis=0) # now use the lowest value (highest road level) because some intersect\n",
    "    road_raster_np[road_raster_np == nan_placeholder] = 0\n",
    "    return road_raster_np # 0=no_road 1=motorway, 2=trunk, ...\n",
    "\n",
    "def osm_values_to_name(values):\n",
    "    mapping = {1:\"Motorway\", 2:\"Trunk\", 3:\"Primary\", 4:\"Secondary\", 5:\"Tertiary\"}\n",
    "    return [mapping[value] for value in values]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### AWS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# copy file to AWS bucket\n",
    "def copy_to_bucket(file):\n",
    "    s3 = s3fs.S3FileSystem()\n",
    "    fname = os.path.basename(file)\n",
    "    s3.put(file, name_bucket()+os.sep+fname)\n",
    "    \n",
    "def get_from_bucket(file):\n",
    "    if file_exists_bucket(file):\n",
    "        s3 = s3fs.S3FileSystem()\n",
    "        s3.get(name_bucket()+os.sep+os.path.basename(file), file)\n",
    "    \n",
    "def file_exists_bucket(file):\n",
    "    s3 = boto3.resource('s3')\n",
    "    try:\n",
    "        s3.Object(name_bucket(), os.path.basename(file)).load()\n",
    "        return True\n",
    "    except botocore.exceptions.ClientError as e:\n",
    "        if e.response['Error']['Code'] == \"404\":\n",
    "            return False\n",
    "\n",
    "def copy_all_to_cache(folder):\n",
    "    for file in os.listdir(folder):\n",
    "        if file[-3:] == \"txt\" or file[-2:] == \"nc\":\n",
    "            src = os.path.join(folder, file)\n",
    "            dst = os.path.join(dirs[\"cache\"], os.path.basename(src))\n",
    "            if not os.path.exists(dst):\n",
    "                shutil.copyfile(src, dst)\n",
    "\n",
    "def copy_all_to_aws(folder):\n",
    "    files = os.listdir(folder)\n",
    "    for file in files:\n",
    "        is_targeted = file[-3:] == \"txt\" or file[-2:] == \"nc\"\n",
    "        if is_targeted and not file_exists_bucket(file):\n",
    "            copy_to_bucket(os.path.join(folder, file))\n",
    "        \n",
    "def clean_up_procedure(dirs_ts):\n",
    "    copy_all_to_aws(dirs_ts[N_OBS()])\n",
    "    \n",
    "def delete_procedure(dirs_ts):\n",
    "    files_acq = os.listdir(dirs_ts[ABS_N_TRUCKS()])\n",
    "    files_obs = os.listdir(dirs_ts[N_OBS()])\n",
    "    msg = \"Does not exist in AWS bucket: \"\n",
    "    for file in files_obs:\n",
    "        path = os.path.join(dirs_ts[N_OBS()], file)\n",
    "        if file_exists_bucket(file):\n",
    "            os.remove(path)\n",
    "        else:\n",
    "            print(\"%s%s\" %(msg, path))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_________\n",
    "## 5 | Truck detection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### TruckDetector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TruckDetector detects trucks at acquisition-level\n",
    "class TruckDetector():   \n",
    "    def __init__(self, band_stack, thresholds):\n",
    "        is_none = band_stack is None\n",
    "        self.band_stack = band_stack.chunk(band_stack.dims[\"lat\"], band_stack.dims[\"lon\"]) if not is_none else None\n",
    "        self.thresholds = thresholds\n",
    "        self.no_truck_mask = None\n",
    "        self.trucks = None\n",
    "    \n",
    "    def detect(self):\n",
    "        B02 = self.band_stack.B02.persist().values\n",
    "        B03 = self.band_stack.B03.persist().values\n",
    "        B04 = self.band_stack.B04.persist().values\n",
    "        B08 = self.band_stack.B08.persist().values\n",
    "        B11 = self.band_stack.B11.persist().values\n",
    "        no_truck_mask = calc_no_trucks(B02, B03, B04, B08, B11)\n",
    "        trucks = detect_trucks(B02, B03, B04, no_truck_mask)\n",
    "        lon_lat = {\"lat\":self.band_stack.lat, \"lon\":self.band_stack.lon}\n",
    "        self.no_truck_mask = xr.DataArray(no_truck_mask, coords=lon_lat, dims=(\"lat\", \"lon\"))\n",
    "        self.trucks = xr.DataArray(trucks, coords=lon_lat, dims=(\"lat\", \"lon\"))\n",
    "        self.filter_trucks()\n",
    "                            \n",
    "    def filter_trucks(self):\n",
    "        self.trucks = filter_spatial_3x3_extended(self.trucks)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate a binary mask where pixels that are definitely no trucks are represented as 0.\n",
    "# thresholds Dict with at least:\n",
    "### max_ndvi Float above this val: no trucks. For Vegetation\n",
    "### max_ndwi Float above this val: no trucks. For Water\n",
    "### max_ndsi Float above this val: no_trucks. For Snow\n",
    "### min_rgb Float above this val: no_trucks. For dark surfaces, e.g. shadows\n",
    "### max_blue Float above this val: no_trucks\n",
    "### max_green Float above this val: no trucks\n",
    "### max_red Float above this val: no trucks\n",
    "### min_b11 Float below this val: no trucks. For dark surfaces, e.g. shadows\n",
    "### max_b11 Float below this val: no trucks. For bright (sealed) surfaces, e.g. buildings\n",
    "@jit(nopython=True, parallel=True)\n",
    "def calc_no_trucks(B02, B03, B04, B08, B11):\n",
    "    th = {\"min_blue\":0.06, \n",
    "          \"min_green\":0.04, \n",
    "          \"min_red\":0.04,\n",
    "          \"max_red\":0.15,\n",
    "          \"max_green\":0.15,\n",
    "          \"max_blue\":0.2,\n",
    "          \"max_ndvi\":0.5,\n",
    "          \"max_ndwi\":0.0001,\n",
    "          \"max_ndsi\":0.0001,\n",
    "          \"min_blue_green_ratio\":0.03,\n",
    "          \"min_blue_red_ratio\":0.05, \n",
    "          \"max_blue_green_ratio\":0.17, \n",
    "          \"max_blue_red_ratio\":0.2}\n",
    "    ndvi_mask = ratio(B08, B04) < th[\"max_ndvi\"]\n",
    "    ndwi_mask = ratio(B02, B11) < th[\"max_ndwi\"]\n",
    "    ndsi_mask = ratio(B03, B11) < th[\"max_ndsi\"]\n",
    "    low_rgb_mask = (B02 > th[\"min_blue\"]) * (B03 > th[\"min_green\"]) * (B04 > th[\"min_red\"])\n",
    "    high_rgb_mask = (B02 < th[\"max_blue\"]) * (B03 < th[\"max_green\"]) * (B04 < th[\"max_red\"])\n",
    "    no_truck_mask = ndvi_mask * ndwi_mask * ndsi_mask * low_rgb_mask * high_rgb_mask\n",
    "    return no_truck_mask\n",
    "\n",
    "# Calculate a binary mask where trucks are represented as 1 and no trucks as 0.\n",
    "# thresholds Dict with at least:\n",
    "### min_green_ratio Float, minimum value of blue-green ratio\n",
    "### min_red_ratio Float, minimum value of blue-red ratio\n",
    "@jit(nopython=True, parallel=True)\n",
    "def detect_trucks(B02, B03, B04, no_truck_mask):\n",
    "    th = {\"min_blue\":0.06, \n",
    "          \"min_green\":0.04, \n",
    "          \"min_red\":0.04,\n",
    "          \"max_red\":0.15,\n",
    "          \"max_green\":0.15,\n",
    "          \"max_blue\":0.2,\n",
    "          \"max_ndvi\":0.5,\n",
    "          \"max_ndwi\":0.0001,\n",
    "          \"max_ndsi\":0.0001,\n",
    "          \"min_blue_green_ratio\":0.03,\n",
    "          \"min_blue_red_ratio\":0.05, \n",
    "          \"max_blue_green_ratio\":0.17, \n",
    "          \"max_blue_red_ratio\":0.2}\n",
    "    bg_ratio = ratio(B02, B03)\n",
    "    br_ratio = ratio(B02, B04)\n",
    "    bg = bg_ratio > th[\"min_blue_green_ratio\"]\n",
    "    br = br_ratio > th[\"min_blue_red_ratio\"]\n",
    "    bg_max = bg_ratio < th[\"max_blue_green_ratio\"]\n",
    "    br_max = br_ratio < th[\"max_blue_red_ratio\"]\n",
    "    trucks = (bg * br * bg_max * br_max) * no_truck_mask\n",
    "    return trucks\n",
    "\n",
    "@jit(nopython=True, parallel=True)\n",
    "def ratio(a, b):\n",
    "    return (a-b)/(a+b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_________\n",
    "## 6 | Processing classes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### PeriodProcessor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PeriodProcessor processes a period of dates, represented in one cube\n",
    "class PeriodProcessor():\n",
    "    def __init__(self, start, end, bbox, bbox_id):\n",
    "        self.start = start\n",
    "        self.end = end\n",
    "        self.bbox = bbox\n",
    "        self.bbox_id = bbox_id\n",
    "        self.cube = None\n",
    "        self.dates = None\n",
    "        self.lon_lat = None\n",
    "        self.mean_trucks = None\n",
    "        self.sum_trucks = None\n",
    "        self.sum_obs = None\n",
    "        self.trucks_points_result = None\n",
    "        self.osm_mask = None\n",
    "        self.used_dates = []\n",
    "        self.percentage_valid = 0\n",
    "    \n",
    "    def get_cube(self, dataset, bands, tile_size, spatial_res, day_bins):\n",
    "        config = CubeConfig(dataset_name = dataset,\n",
    "                            band_names = bands,\n",
    "                            tile_size = tile_size,\n",
    "                            geometry = self.bbox,\n",
    "                            spatial_res = spatial_res,\n",
    "                            time_range = [self.start, self.end])\n",
    "        worked, i = False, 0\n",
    "        while i < 50 and not worked:            \n",
    "            try:\n",
    "                cube = open_cube(config)\n",
    "                worked = True\n",
    "            except:\n",
    "                print(\"Holding\")\n",
    "                sleep(120)\n",
    "                i += 1\n",
    "        self.cube = cube\n",
    "        self.dates = cube.time.values\n",
    "        self.lon_lat = {\"lon\":cube.lon.values, \"lat\":cube.lat.values}\n",
    "    \n",
    "    def calc_osm_mask(self, osm):\n",
    "        osm_raster = rasterize_osm(osm, self.cube.B02.sel(time = self.dates[0]))\n",
    "        self.osm_mask = create_xr_dataset(osm_raster, self.lon_lat, \"roadmask\")\n",
    "    \n",
    "    # Add from acquisition methods\n",
    "    def add_n_observations(self, no_clouds):\n",
    "        obs = np.where(no_clouds==1, 1, 0)\n",
    "        obs[np.isnan(obs)] = 0\n",
    "        if self.sum_obs is None:\n",
    "            self.sum_obs = obs\n",
    "        else:\n",
    "            self.sum_obs += obs\n",
    "        \n",
    "    def percentage_covered(self):\n",
    "        self.sum_obs[self.osm_mask.roadmask.values == 0] = -1\n",
    "        n_valid = len(self.sum_obs[self.sum_obs >= 1].ravel())\n",
    "        n_all = len(self.sum_obs[self.sum_obs >= 0].ravel())\n",
    "        self.percentage_valid = (n_valid / n_all) * 100\n",
    "        \n",
    "    def obs_time(self, n_obs):\n",
    "        n_pixels_roads = len(period.sum_obs[period.sum_obs >= 0])\n",
    "        n_pixels_obs = len(period.sum_obs[period.sum_obs >= n_obs])\n",
    "        return (n_pixels_obs / n_pixels_roads) * 100\n",
    "    \n",
    "    def add_detections(self, trucks):\n",
    "        if self.sum_trucks is None:\n",
    "            self.sum_trucks = trucks\n",
    "        else:\n",
    "            self.sum_trucks += trucks\n",
    "    \n",
    "    def add_acquisition(self, acquisition, date_str):\n",
    "        self.add_n_observations(acquisition.no_clouds)\n",
    "        self.add_detections(acquisition.detector.trucks.values)\n",
    "        self.used_dates.append(date_str)\n",
    "    \n",
    "    # Temporal summary methods\n",
    "    def sum_truck_n(self):\n",
    "        self.sum_trucks = parallel_sum(np.array(self.detections))\n",
    "    \n",
    "    def sum_observations(self):\n",
    "        self.sum_obs = parallel_sum(np.array(self.n_observations))\n",
    "    \n",
    "    def mask_sum_obs_to_trucks(self):\n",
    "        self.sum_obs[self.sum_trucks == 0] = 0\n",
    "    \n",
    "    def mean_truck_n(self):\n",
    "        self.mean_trucks = np.divide(self.sum_trucks, self.sum_obs)\n",
    "        self.mean_trucks[np.isnan(self.mean_trucks)] = 0\n",
    "    \n",
    "    def mask_osm_to_trucks(self, osm_mask):\n",
    "        self.osm_mask = self.osm_mask.roadmask.values\n",
    "        self.osm_mask[self.sum_trucks == 0] = 0\n",
    "            \n",
    "    # Write methods\n",
    "    def write_n_observations(self, fname):\n",
    "        sum_obs_xr = create_xr_dataset(self.sum_obs, self.lon_lat, os.path.basename(fname))\n",
    "        sum_obs_xr.to_netcdf(fname)\n",
    "        \n",
    "    def write_sum_trucks(self, fname):\n",
    "        sum_xr = create_xr_dataset(period.sum_trucks, period.lon_lat, os.path.basename(fname))\n",
    "        sum_xr.to_netcdf(fname)\n",
    "    \n",
    "    def write_mean_trucks(self, fname):\n",
    "        mean_xr = create_xr_dataset(self.mean_trucks, self.lon_lat, os.path.basename(fname))\n",
    "        mean_xr.to_netcdf(fname)\n",
    "        move_to_bucket(fname)\n",
    "        \n",
    "    def write_trucks_vec(self, points, fname, fname_placeholder):\n",
    "        got_points = points is not None and len(points) > 0\n",
    "        if got_points:\n",
    "            points.to_file(fname, driver = vec_driver_from_ext(\".\" + fname.split(\".\")[1]))\n",
    "        else:\n",
    "            # write txt as placeholder\n",
    "            with open(fname_placeholder, \"w\") as file:\n",
    "                file.write(\"%s has length: %s. Nothing to write\" %(os.path.basename(fname_placeholder), str(len(points))))\n",
    "        \n",
    "    def wrap_period(self, fnames):\n",
    "        crs = EPSG_4326()\n",
    "        ind_right = \"index_right\"\n",
    "        fname_obs = fnames[N_OBS()]\n",
    "        self.write_n_observations(fname_obs)\n",
    "        copy_to_bucket(fname_obs)\n",
    "        self.mean_truck_n()\n",
    "        self.mask_osm_to_trucks(self.osm_mask) # mask to trucks\n",
    "        self.mask_sum_obs_to_trucks() # mask to trucks\n",
    "        # merge n observations, n trucks and mean trucks into single points layer\n",
    "        try:\n",
    "            n_obs_points = raster_to_points(self.sum_obs, self.lon_lat, \"sum_observations\", crs)\n",
    "            sum_trucks_points = raster_to_points(self.sum_trucks, self.lon_lat, \"sum_trucks_sub_period\", crs)\n",
    "            sum_trucks_obs = gpd.sjoin(sum_trucks_points, n_obs_points, how=\"inner\", op=\"intersects\")\n",
    "            sum_trucks_obs = sum_trucks_obs.drop([ind_right], axis=1)\n",
    "            osm_points = raster_to_points(self.osm_mask, self.lon_lat, \"osm_value\", crs)\n",
    "            osm_points[\"osm_name\"] = osm_values_to_name(osm_points[\"osm_value\"])\n",
    "            sum_trucks_obs_osm = gpd.sjoin(sum_trucks_obs, osm_points, how=\"inner\", op=\"intersects\")\n",
    "            sum_trucks_obs_osm = sum_trucks_obs_osm.drop([ind_right], axis=1)\n",
    "            mean_trucks_points = raster_to_points(self.mean_trucks, self.lon_lat, \"mean_trucks\", crs)\n",
    "            self.trucks_points_result = gpd.sjoin(sum_trucks_obs_osm, mean_trucks_points, how=\"inner\", op=\"intersects\")\n",
    "            self.trucks_points_result = self.trucks_points_result.drop([ind_right], axis=1)\n",
    "            self.trucks_points_result[VALID_PIXELS()] = [self.percentage_valid] * len(self.trucks_points_result)\n",
    "            self.write_trucks_vec(self.trucks_points_result, fnames[TRUCKS_VEC()], fnames[TRUCKS_VEC_PHR()])\n",
    "        except:\n",
    "            Warning(\"No results could be written, detections could be empty\")\n",
    "        if file_exists_bucket(fname_obs) and os.path.exists(fname_obs):\n",
    "            os.remove(fname_obs)\n",
    "                    \n",
    "    def cleanup(self, period_files):\n",
    "        for file in period_files:\n",
    "            if os.path.exists(file):\n",
    "                os.remove(file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "@jit(nopython=True, parallel=True)\n",
    "def parallel_sum(arr):\n",
    "    return arr.sum(axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### AcquisitionProcessor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# AcquisitionProcessor processes all valid pixels of a single acquisition in cube\n",
    "class AcquisitionProcessor():\n",
    "    def __init__(self, date_np64, cube):\n",
    "        self.date_np64 = date_np64\n",
    "        self.cube = cube\n",
    "        self.band_stack = cube.sel(time = date_np64)\n",
    "        self.detector = None\n",
    "        self.no_clouds = None\n",
    "        self.osm_mask = None\n",
    "        self.has_obs = None\n",
    "        self.improves_cov = False\n",
    "           \n",
    "    def mask_clouds(self, cloud_masking_thresholds):       \n",
    "        scl = MaskSet(self.band_stack.SCL)\n",
    "        high_prob = scl.clouds_high_probability\n",
    "        med_prob = scl.clouds_medium_probability\n",
    "        cirrus = scl.cirrus\n",
    "        no_data = scl.no_data\n",
    "        rgb_cloud_mask = calc_rgb_cloud_mask(self.band_stack, cloud_masking_thresholds)\n",
    "        self.no_clouds = (high_prob + med_prob + cirrus + no_data + rgb_cloud_mask) == 0\n",
    "        self.band_stack = self.band_stack.where(self.no_clouds)\n",
    "        \n",
    "    # percentage Float secifies the percentage of valid observations\n",
    "    # that is needed to be considered as valid acquisition\n",
    "    def has_observations(self, minimum_valid_percentage):\n",
    "        values = self.no_clouds.values\n",
    "        # mask valid pixel mask to OSM roads\n",
    "        values[self.osm_mask.roadmask == 0] = 0\n",
    "        n_vals = np.count_nonzero(self.osm_mask.roadmask)\n",
    "        n_valid = np.count_nonzero(values)\n",
    "        percent_valid = (n_valid / n_vals) * 100\n",
    "        self.has_obs = percent_valid >= minimum_valid_percentage\n",
    "        \n",
    "    def improves_coverage(self, period, min_improvement):\n",
    "        first_acq = period.sum_obs is None\n",
    "        if first_acq:\n",
    "            period.sum_obs = np.zeros((self.band_stack.B02.shape[0], self.band_stack.B02.shape[1]))\n",
    "        sum_obs_copy = period.sum_obs.copy() # get copy because original will be altered temporarilly\n",
    "        cov_before = period.percentage_valid\n",
    "        period.add_n_observations(self.no_clouds)\n",
    "        period.percentage_covered()\n",
    "        cov_after = period.percentage_valid\n",
    "        period.sum_obs = sum_obs_copy # re-assign non-altered sum_obs array\n",
    "        self.improves_cov = (cov_after - cov_before) > min_improvement\n",
    "    \n",
    "    def mask_with_osm(self, osm_mask):\n",
    "        self.osm_mask = osm_mask\n",
    "        self.band_stack = self.band_stack.where(self.osm_mask.roadmask != 0)\n",
    "                \n",
    "    def do_detection(self, thresholds):\n",
    "        self.detector = TruckDetector(self.band_stack, thresholds)\n",
    "        self.detector.detect()\n",
    "        \n",
    "    def write_detections(self, fname_vec):\n",
    "        lon_lat = {\"lat\":self.cube.lat.values, \"lon\":self.cube.lon.values}\n",
    "        points = raster_to_points(self.detector.trucks.values, lon_lat, \"trucks\", EPSG_4326())\n",
    "        points.to_file(fname_vec, driver=GPKG())\n",
    "        \n",
    "    def read_detections(self, fname, fname_vec):\n",
    "        self.detector = TruckDetector(None, None)\n",
    "        if os.path.exists(fname):\n",
    "            dataset = xr.open_dataset(fname)\n",
    "            self.no_clouds = dataset[\"no_clouds\"]\n",
    "            trucks = dataset[os.path.basename(fname)]\n",
    "            self.detector.trucks = trucks\n",
    "        elif os.path.exists(fname_vec):\n",
    "            points = gpd.read_file(fname_vec)\n",
    "            raster = points_to_raster(points, self.band_stack.B02).astype(np.bool)\n",
    "            self.detector.trucks = xr.DataArray(raster, coords={\"lat\":self.cube.lat.values, \"lon\":self.cube.lon.values}, dims=[\"lat\", \"lon\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__________\n",
    "____________\n",
    "## 7 | Execute processing\n",
    "Process by grid __box__ (bbox_id), __sub-period__, __acquisition__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Main process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--------------------------------------------------------------------------------\n",
      "====================================================================================================\n",
      "\n",
      "Starting truck detection processing \n",
      "..................................................\n",
      "\n",
      "Weekdays: ['tuesday', 'wednesday', 'thursday']\n",
      "Number of grid cells to process: 97\n",
      "Number of periods: 4\n",
      "Timestamps sub-period: 1\n",
      "OSM roads: ['motorway', 'trunk', 'primary']\n",
      "Baseline years: [2017, 2018, 2019]\n",
      "Minimum valid observations to consider acquisition: 15 %\n",
      "Overwrite results: True\n",
      "..................................................\n",
      "\n",
      "====================================================================================================\n",
      "--------------------------------------------------------------------------------\n",
      "\n",
      "Processing: 68\n",
      "bbox_id: 398\n",
      "\n",
      "--------------------------------------------------------------------------------\n",
      "TS: 0 Period Start: 2017-04-01 00:00:00 End: 2017-06-21 00:00:00\n",
      "Getting cube\n",
      "Calculating OSM mask\n",
      "--------------------------------------------------------------------------------\n",
      "Date: 2017-04-04T10_11_19\n",
      "Calculating cloud masks\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/lib/python3.7/site-packages/dask/core.py:121: RuntimeWarning: invalid value encountered in greater\n",
      "  return func(*(_execute_task(a, cache) for a in args))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Skipping, valid obs. below 15 %\n",
      "--------------------------------------------------------------------------------\n",
      "Date: 2017-05-30T10_30_24\n",
      "Calculating cloud masks\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/lib/python3.7/site-packages/dask/core.py:121: RuntimeWarning: invalid value encountered in greater\n",
      "  return func(*(_execute_task(a, cache) for a in args))\n"
     ]
    }
   ],
   "source": [
    "# make or read processing grid\n",
    "if os.path.exists(files[\"proc_grid\"]):\n",
    "    try:\n",
    "        grid_gadm = gpd.read_file(files[\"proc_grid\"])\n",
    "        grid_gadm = unique_grid(files[\"proc_grid\"])\n",
    "    except:\n",
    "        raise Exception(\"Failed reading proc grid from: \" + files[\"proc_grid\"])\n",
    "else:\n",
    "    grid_gadm = make_grid(grid_spacing, files)\n",
    "    grid_gadm = unique_grid(files[\"proc_grid\"])\n",
    "if country is not None:\n",
    "    grid_gadm = grid_gadm[grid_gadm[\"NAME_0\"]==country]\n",
    "# calc temporal bounds of baseline sub-periods\n",
    "if len(baseline_years) > 0:\n",
    "    periods = yearly_period_from_target(target, baseline_years, timestamps_sub_period)\n",
    "else:\n",
    "    periods = calc_periods(n_days_sub, baseline, target)\n",
    "sep = \"-\" * 80\n",
    "print(sep)\n",
    "trace = []\n",
    "ext_arr = NC_EXT()\n",
    "ext_vec = GPKG_EXT()\n",
    "n_boxes = len(grid_gadm)\n",
    "start_msg(weekdays, n_boxes, periods, timestamps_sub_period, osm_values, baseline_years, minimum_valid_observations, overwrite_results)\n",
    "for i in range(n_boxes):\n",
    "    bbox = list(grid_gadm.geometry)[i].bounds\n",
    "    bbox_id = list(grid_gadm[BBOX_ID()])[i]\n",
    "    bbox_id_str = str(bbox_id)\n",
    "    first, last = periods[\"first\"], periods[\"last\"]\n",
    "    print(\"%s\\n\\nProcessing: %s\\nbbox_id: %s\\n\\n%s\" %(sep, str(i), bbox_id_str, sep))\n",
    "    try:\n",
    "        file_osm = get_roads(bbox, bbox_id, osm_values, osm_key, roads_buffer, dirs[\"ancil_roads\"])\n",
    "        # retry\n",
    "        if not os.path.exists(file_osm): file_osm = get_roads(bbox, bbox_id, osm_values, osm_key, roads_buffer, dirs[\"ancil_roads\"])\n",
    "        osm = gpd.read_file(file_osm)\n",
    "        # subset according to desired osm road levels (file might have been written with other values)\n",
    "        value_not_desired = [v not in osm_values for v in list(osm[\"osm_value\"])]\n",
    "        osm = osm.drop(osm[value_not_desired].index)\n",
    "    except:\n",
    "        msg = \"Could not get OSM roads: \" + bbox_id_str\n",
    "        Warning(msg)\n",
    "        trace.append(msg)\n",
    "        #continue\n",
    "    osm_mask = None\n",
    "    for start, end in zip(first, last):\n",
    "        period_files = []\n",
    "        tracker_start = datetime.now()\n",
    "        ts = first.index(start)\n",
    "        ts_str = str(ts)\n",
    "        dir_ts = dirs[\"processed\"] # write into output pool\n",
    "        if not os.path.exists(dir_ts): os.mkdir(dir_ts)\n",
    "        dirs_ts = make_dirs_ts(dir_ts)\n",
    "        print(\"TS: %s Period Start: %s End: %s\" %(ts_str, str(start), str(end)))\n",
    "        # check if yet processed\n",
    "        fnames = sub_period_fnames(dirs_ts, bbox_id, str(start.date()) + \"_\" + str(end.date()), ext_arr, ext_vec)\n",
    "        exists = {}\n",
    "        for key, f in fnames.items():\n",
    "            exists[key] = os.path.exists(f) if not f[-3:]==\"txt\" else True\n",
    "        exists[N_OBS()] = not overwrite_results            \n",
    "        if not exists[TRUCKS_VEC()]: exists[TRUCKS_VEC()] = os.path.exists(fnames[TRUCKS_VEC_PHR()]) # placeholder might exist instead\n",
    "        already_proc = \"because already processed\"\n",
    "        osm_msg = \"Calculating OSM mask\\n\" + sep\n",
    "        skip = all(exists.values()) and not overwrite_results\n",
    "        if skip:\n",
    "            print(\"Skipping \" + already_proc)\n",
    "        else:\n",
    "            print(\"Getting cube\")\n",
    "            period = PeriodProcessor(start, end, bbox, bbox_id)\n",
    "            period.get_cube(dataset, bands, tile_size, spatial_res, day_bins)\n",
    "            t1 = datetime.now() # track time cube is open to prevent timeout\n",
    "            print(osm_msg)\n",
    "            period.calc_osm_mask(osm)\n",
    "            p = 0\n",
    "            obs_full = False\n",
    "            # order by weekday of target weekdays\n",
    "            period_dates_sorted = []\n",
    "            for w in weekdays:\n",
    "                dates_w = period.dates[[is_weekday(period_date, w) for period_date in period.dates]].flatten()\n",
    "                for d in range(len(dates_w)): period_dates_sorted.append(dates_w[d])\n",
    "            period.dates = period_dates_sorted\n",
    "            while p < len(period.dates) and not obs_full:\n",
    "                period_date = period.dates[p]\n",
    "                if (datetime.now()-t1).seconds > 1800: # if 30 minutes exceeded get cube again to prevent timeout\n",
    "                    print(\"Getting cube again to prevent timeout\")\n",
    "                    period.get_cube(dataset, bands, tile_size, spatial_res, day_bins)\n",
    "                    t1 = datetime.now()\n",
    "                    period.calc_osm_mask(osm)\n",
    "                acquisition = AcquisitionProcessor(period_date, period.cube)\n",
    "                date_str = str(period_date).replace(\":\",\"_\")[0:-10]\n",
    "                print(\"Date: %s\" %(date_str))\n",
    "                fname = fname_acquisition_trucks(dirs_ts, bbox_id, date_str, ext_arr)\n",
    "                fname_vec = fname_acquisition_trucks(dirs_ts, bbox_id, date_str, ext_vec)\n",
    "                acq_processed = [os.path.exists(fname_vec), os.path.exists(fname), file_exists_bucket(fname)]\n",
    "                if any(acq_processed):\n",
    "                    print(\"File already processed\\n%s\" %(sep))\n",
    "                    if not acq_processed[0] and not acq_processed[1] and acq_processed[2]:\n",
    "                        get_from_bucket(fname)\n",
    "                    acquisition.read_detections(fname, fname_vec)\n",
    "                    if not os.path.exists(fname_vec):\n",
    "                        acquisition.write_detections(fname_vec)\n",
    "                        os.remove(fname)\n",
    "                        acquisition.read_detections(fname, fname_vec)\n",
    "                    n_try, max_try = 0, 2\n",
    "                    while n_try <= max_try:\n",
    "                        try:\n",
    "                            acquisition.mask_clouds(cloud_masking_thresholds)\n",
    "                            period.add_acquisition(acquisition, date_str)\n",
    "                            n_try = max_try+1\n",
    "                        except:\n",
    "                            n_try += 1\n",
    "                    period.used_dates.append(period_date)\n",
    "                    period_files.append(fname)\n",
    "                else:\n",
    "                    try:\n",
    "                        min_valid_str = str(minimum_valid_observations)\n",
    "                        fname_placeholder = os.path.join(dirs_ts[ABS_N_TRUCKS()], \n",
    "                                                         date_str + min_valid_str + \"_\" + bbox_id_str + \n",
    "                                                         \"_placeholder.txt\").replace(\"-\",\"_\")\n",
    "                        can_have_obs = not os.path.exists(fname_placeholder)\n",
    "                        if can_have_obs:\n",
    "                            print(\"Calculating cloud masks\")\n",
    "                            acquisition.mask_clouds(cloud_masking_thresholds)\n",
    "                            acquisition.mask_with_osm(period.osm_mask) # mask to OSM roads\n",
    "                            # check if enough observations to be considered\n",
    "                            acquisition.has_observations(minimum_valid_observations)\n",
    "                            if not acquisition.has_obs:\n",
    "                                with open(fname_placeholder, \"w\") as f:\n",
    "                                    f.write(\"%s not fulfilled for date %s\" %(min_valid_str, date_str))\n",
    "                        else:\n",
    "                            has_obs = can_have_obs\n",
    "                            print(sep)\n",
    "                    except:\n",
    "                        print(\"Cloud masking failed. Assuming not enough observations\")\n",
    "                        acquisition.has_obs = False\n",
    "                    if acquisition.has_obs:\n",
    "                        print(\"Processing acquisition: %s\" %(date_str))\n",
    "                        try:\n",
    "                            acquisition.do_detection(thresholds)\n",
    "                            print(\"Writing\")\n",
    "                            acquisition.write_detections(fname_vec)\n",
    "                            print(\"Number of trucks: %s\" %(str(np.count_nonzero(acquisition.detector.trucks))))\n",
    "                            period.add_acquisition(acquisition, date_str)\n",
    "                            period.used_dates.append(period_date)\n",
    "                            period_files.append(fname)\n",
    "                            print(\"Done\")\n",
    "                        except:\n",
    "                            print(\"!!Failed: %s!!\" %(date_str))\n",
    "                        print(\"Done with acquisition: %s\\n%s\" %(date_str, sep))\n",
    "                    else:\n",
    "                        print(\"Skipping, valid obs. below %s %%\\n%s\" %(minimum_valid_observations, sep))\n",
    "                if period.sum_obs is not None:\n",
    "                    period.percentage_covered()\n",
    "                    print(\"Percentage valid in mosaic: %s\" %(str(period.percentage_valid)))\n",
    "                    n_obs_time = period.obs_time(n_obs)\n",
    "                    print(\"Percentage covered by at least %s observations: %s\" %(str(n_obs), str(n_obs_time)))\n",
    "                    obs_full = period.percentage_valid >= obs_full_threshold and n_obs_time >= obs_full_threshold\n",
    "                p += 1\n",
    "            if len(period_files) > 0:\n",
    "                print(\"%s\\nNumber of files: %s\" %(sep, str(len(period_files))))\n",
    "                # read all processed files from disk in order not to hold in memory during processing\n",
    "                print(\"Aggregating sub-period\")\n",
    "                period.wrap_period(fnames)        \n",
    "                fname_used_acq = os.path.basename(fnames[TRUCKS_VEC()])[0:-3]+\"txt\"\n",
    "                with open(os.path.join(dirs_ts[ABS_N_TRUCKS()], fname_used_acq), \"w\") as f:\n",
    "                    for date in period.used_dates:\n",
    "                        f.write(str(date))\n",
    "                period.cleanup(period_files) # delete rasters, vector points are still there & can be reread\n",
    "            else:\n",
    "                msg = \"No acquisitions in period %s to %s. In bbox_id: %s\" %(str(start), str(end), bbox_id_str)\n",
    "                Warning(msg)\n",
    "                trace.append(msg)\n",
    "            print(\"%s\\nProcessing sub-period took: %s minutes\" %(sep, str((datetime.now()-tracker_start).seconds/60)))\n",
    "        print(\"%s\\nDone with period %s of bbox_id %s\\n%s\" %(sep, ts_str, bbox_id_str, sep))    \n",
    "    print(\"%s\\nDone with bbox_id: %s\\n%s\\n%s\" %(sep, bbox_id_str, sep, sep))\n",
    "    print(\"%s\\nMoving files to AWS bucket..\\n%s\" %(sep, sep))\n",
    "    clean_up_procedure(dirs_ts)\n",
    "print(\"%s\\nDone with all requested bboxes\\n%s\\n%s\" %(sep, sep, sep*2))\n",
    "print(\"%s\\nDeleting files on local disk that have been moved to AWS bucket\\n%s\" %(sep, sep))\n",
    "delete_procedure(dirs_ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Merge detections of different processing boxes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_ts = dirs[\"processed\"] # write into output pool\n",
    "if not os.path.exists(dir_ts): os.mkdir(dir_ts)\n",
    "dirs_ts = make_dirs_ts(dir_ts)\n",
    "all_files = glob(dirs_ts[TRUCKS_VEC()]+os.sep+\"*.gpkg\")\n",
    "first = periods[\"first\"]\n",
    "last = periods[\"last\"]\n",
    "for start, end in zip(periods[\"first\"], periods[\"last\"]):\n",
    "    print(\"Merging from %s to %s\" %(start, end))\n",
    "    points_period = []\n",
    "    for f in all_files:\n",
    "        if string_to_datetime(f.split(\"_\")[9]) == start:\n",
    "            try:\n",
    "                points = gpd.read_file(f)\n",
    "                points_period.append(points)\n",
    "            except:\n",
    "                Warning(\"Could not read: %s\" %(f))\n",
    "    if len(points_period) > 0:\n",
    "        merged = pd.concat(points_period)\n",
    "        fsep = \"_\"\n",
    "        fname = TRUCKS_VEC_MERGED() + fsep + fsep + str(start.date()) + fsep + str(end.date()) + GPKG_EXT()\n",
    "        merged[TRUCK_COUNT_NORM()] = merged[MEAN_TRUCKS()] * (100/merged[VALID_PIXELS()])\n",
    "        merged.to_file(os.path.join(dirs_ts[TRUCKS_VEC_MERGED()], fname), driver = GPKG())\n",
    "    else:\n",
    "        print(\"No files. \\nStart: %s\\nEnd: %s\" %(start, end))\n",
    "\n",
    "results_years = baseline_years.copy()\n",
    "results_years.append(2020)\n",
    "merged_files = glob(dirs_ts[TRUCKS_VEC_MERGED()]+os.sep+\"*.gpkg\")\n",
    "median_values = []\n",
    "for year in results_years:\n",
    "    points_year = []\n",
    "    for file in merged_files:\n",
    "        if \"__\" + str(year) in os.path.basename(file):\n",
    "            points = gpd.read_file(file)\n",
    "            points_year.append(points)\n",
    "    n_elements = np.array([np.array(x[x[\"osm_value\"]==1][TRUCK_COUNT_NORM()]).sum() for x in points_year])\n",
    "    merged = pd.concat(points_year)\n",
    "    file_str = TRUCKS_VEC_MERGED() + \"_\" + str(year)\n",
    "    merged.to_file(os.path.join(dirs_ts[TRUCKS_VEC_MERGED()], file_str + GPKG_EXT()), driver=GPKG())\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}