calculate-biomass/CalcBiomass.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import os\n",
    "from osgeo import gdal, osr\n",
    "import matplotlib.pyplot as plt\n",
    "import sys\n",
    "import tempfile\n",
    "import os\n",
    "import rasterio\n",
    "import math\n",
    "import pickle\n",
    "import lithops\n",
    "from lithops import Storage\n",
    "from scipy import ndimage as ndi\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "from skimage.segmentation import watershed\n",
    "from skimage.feature import peak_local_max\n",
    "from skimage.measure import regionprops\n",
    "from sklearn.ensemble import RandomForestRegressor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "BUCKET_NAME = 'YOUR_BUCKET_NAME'"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "storage = Storage()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Define functions"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "* `plot_band_array`: function to plot NEON spatial data."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "def plot_band_array(band_array, image_extent, title, cmap_title, colormap, colormap_limits):\n",
    "    plt.imshow(band_array, extent=image_extent)\n",
    "    cbar = plt.colorbar()\n",
    "    plt.set_cmap(colormap)\n",
    "    plt.clim(colormap_limits)\n",
    "    cbar.set_label(cmap_title, rotation=270, labelpad=20)\n",
    "    plt.title(title)\n",
    "    ax = plt.gca()\n",
    "    ax.ticklabel_format(useOffset=False, style='plain')\n",
    "    rotatexlabels = plt.setp(ax.get_xticklabels(),rotation=90)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "* `array2raster`: function to output geotiff files."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "def array2raster(newRasterfn, rasterOrigin_file, pixelWidth, pixelHeight, array_file, epsg):\n",
    "    storage = Storage()\n",
    "    array_byte = storage.get_cloudobject(array_file)\n",
    "    array_smooth = pickle.loads(array_byte)\n",
    "    array = np.array(array_smooth,dtype=float)\n",
    "\n",
    "    metadata_byte = storage.get_cloudobject(rasterOrigin_file)\n",
    "    metadata = pickle.loads(metadata_byte)\n",
    "    rasterOrigin = (metadata['ext_dict']['xMin'],metadata['ext_dict']['yMax'])\n",
    "\n",
    "    cols = array.shape[1]\n",
    "    rows = array.shape[0]\n",
    "    originX = rasterOrigin[0]\n",
    "    originY = rasterOrigin[1]\n",
    "\n",
    "    driver = gdal.GetDriverByName('GTiff')\n",
    "    outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Float32)\n",
    "    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))\n",
    "    outband = outRaster.GetRasterBand(1)\n",
    "    outband.WriteArray(array)\n",
    "    outRasterSRS = osr.SpatialReference()\n",
    "    outRasterSRS.ImportFromEPSG(epsg)\n",
    "    outRaster.SetProjection(outRasterSRS.ExportToWkt())\n",
    "    outband.FlushCache()\n",
    "\n",
    "    with open(newRasterfn, 'rb') as tif_temp:\n",
    "        storage.put_cloudobject(tif_temp.read(), key=(\"geotiff/\" + newRasterfn))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "* `raster2array`: function to conver rasters to an array."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "def raster2array(geotif_file, bucket_name, ibm_cos):\n",
    "    storage = Storage()\n",
    "\n",
    "    ibm_cos.download_file(bucket_name, geotif_file, \"/tmp/\" + geotif_file.replace(\"chunks/\", \"\"))\n",
    "\n",
    "    geotif_file = geotif_file.replace(\"chunks/\", \"\")\n",
    "\n",
    "    metadata = {}\n",
    "    dataset = gdal.Open(\"/tmp/\" + geotif_file)\n",
    "    metadata['array_rows'] = dataset.RasterYSize\n",
    "    metadata['array_cols'] = dataset.RasterXSize\n",
    "    metadata['bands'] = dataset.RasterCount\n",
    "    metadata['driver'] = dataset.GetDriver().LongName\n",
    "    metadata['projection'] = dataset.GetProjection()\n",
    "    metadata['geotransform'] = dataset.GetGeoTransform()\n",
    "\n",
    "    mapinfo = dataset.GetGeoTransform()\n",
    "    metadata['pixelWidth'] = mapinfo[1]\n",
    "    metadata['pixelHeight'] = mapinfo[5]\n",
    "\n",
    "    metadata['ext_dict'] = {}\n",
    "    metadata['ext_dict']['xMin'] = mapinfo[0]\n",
    "    metadata['ext_dict']['xMax'] = mapinfo[0] + dataset.RasterXSize/mapinfo[1]\n",
    "    metadata['ext_dict']['yMin'] = mapinfo[3] + dataset.RasterYSize/mapinfo[5]\n",
    "    metadata['ext_dict']['yMax'] = mapinfo[3]\n",
    "\n",
    "    metadata['extent'] = (metadata['ext_dict']['xMin'],metadata['ext_dict']['xMax'],\n",
    "                          metadata['ext_dict']['yMin'],metadata['ext_dict']['yMax'])\n",
    "\n",
    "    if metadata['bands'] == 1:\n",
    "        raster = dataset.GetRasterBand(1)\n",
    "        metadata['noDataValue'] = raster.GetNoDataValue()\n",
    "        metadata['scaleFactor'] = raster.GetScale()\n",
    "\n",
    "        # band statistics\n",
    "        metadata['bandstats'] = {} # make a nested dictionary to store band stats in same\n",
    "        stats = raster.GetStatistics(True,True)\n",
    "        metadata['bandstats']['min'] = round(stats[0],2)\n",
    "        metadata['bandstats']['max'] = round(stats[1],2)\n",
    "        metadata['bandstats']['mean'] = round(stats[2],2)\n",
    "        metadata['bandstats']['stdev'] = round(stats[3],2)\n",
    "\n",
    "        array = dataset.GetRasterBand(1).ReadAsArray(0,0,\n",
    "                                                     metadata['array_cols'],\n",
    "                                                     metadata['array_rows']).astype(np.float)\n",
    "        array[array == int(metadata['noDataValue'])] = np.nan\n",
    "        if metadata['scaleFactor'] is not None:\n",
    "            array = array / metadata['scaleFactor']\n",
    "\n",
    "        ref_array = storage.put_cloudobject(pickle.dumps(array), key=('arrayAux/array_object_' + (geotif_file.split(\"COB1_\")[1]).replace(\".tif\", \"\")))\n",
    "        ref_metadata = storage.put_cloudobject(pickle.dumps(metadata), key=('arrayAux/metadata_object_' + (geotif_file.split(\"COB1_\")[1]).replace(\".tif\", \"\")))\n",
    "\n",
    "        return ref_array, ref_metadata\n",
    "\n",
    "    elif metadata['bands'] > 1:\n",
    "        print('More than one band ... need to modify function for case of multiple bands')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "* `crown_geometric_volume_pth`: function to get tree crown volumn."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "def crown_geometric_volume_pth(tree_data,min_tree_height,pth):\n",
    "    p = np.percentile(tree_data, pth)\n",
    "    tree_data_pth = [v if v < p else p for v in tree_data]\n",
    "    crown_geometric_volume_pth = np.sum(tree_data_pth - min_tree_height)\n",
    "    return crown_geometric_volume_pth, p"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "* `get_predictors`: function to get the trees from the biomass data."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "def get_predictors(tree,chm_array, labels):\n",
    "    indexes_of_tree = np.asarray(np.where(labels==tree.label)).T\n",
    "    tree_crown_heights = chm_array[indexes_of_tree[:,0],indexes_of_tree[:,1]]\n",
    "\n",
    "    full_crown = np.sum(tree_crown_heights - np.min(tree_crown_heights))\n",
    "\n",
    "    crown50, p50 = crown_geometric_volume_pth(tree_crown_heights,tree.min_intensity,50)\n",
    "    crown60, p60 = crown_geometric_volume_pth(tree_crown_heights,tree.min_intensity,60)\n",
    "    crown70, p70 = crown_geometric_volume_pth(tree_crown_heights,tree.min_intensity,70)\n",
    "\n",
    "    return [tree.label,\n",
    "            np.float(tree.area),\n",
    "            tree.major_axis_length,\n",
    "            tree.max_intensity,\n",
    "            tree.min_intensity,\n",
    "            p50, p60, p70,\n",
    "            full_crown, crown50, crown60, crown70]"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "* `asc_to_geotiff`: function to convert asc file to geotiff file."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "def asc_to_geotiff(ref_file, bucket_name, ibm_cos):\n",
    "    storage = Storage()\n",
    "    ibm_cos.download_file(bucket_name, ref_file, \"/tmp/\" + ref_file)\n",
    "    file = \"/tmp/\" + ref_file\n",
    "\n",
    "    asc_file_name = os.path.basename(file)\n",
    "    tile_id, _ = os.path.splitext(asc_file_name)\n",
    "    out_path = os.path.join(tile_id + '.tiff')\n",
    "    out_key = os.path.join(tile_id + '.tiff')\n",
    "\n",
    "    print(f'Converting {tile_id} to GeoTIFF...')\n",
    "    with rasterio.open(file, 'r') as src:\n",
    "        profile = src.profile\n",
    "        # Cloud optimized GeoTiff parameters\n",
    "        profile.update(driver='GTiff')\n",
    "        profile.update(blockxsize=256)\n",
    "        profile.update(blockysize=256)\n",
    "        profile.update(tiled=True)\n",
    "        profile.update(compress='deflate')\n",
    "        profile.update(interleave='band')\n",
    "\n",
    "        with rasterio.open(out_path, 'w', **profile) as dest:\n",
    "            dest.write(src.read())\n",
    "\n",
    "        with open(out_path, 'rb') as tif_temp:\n",
    "            storage.put_cloudobject(tif_temp.read(), key=out_path)\n",
    "\n",
    "    return out_key"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "* `data_chunker`: function to split the geotiff file to data chunks."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "def data_chunker(file, dst, n_splits, block_x, block_y, bucket_name, ibm_cos):\n",
    "    storage = Storage()\n",
    "\n",
    "    ibm_cos.download_file(bucket_name, file, \"/tmp/\" + file)\n",
    "\n",
    "    tile_key = os.path.basename(file)\n",
    "    tile_id, _ = os.path.splitext(tile_key)\n",
    "\n",
    "    with rasterio.open(\"/tmp/\" + file) as src:\n",
    "        transform = src.transform\n",
    "\n",
    "        # Compute working window\n",
    "        step_w = src.width / n_splits\n",
    "        step_h = src.height / n_splits\n",
    "\n",
    "        offset_h = round(step_h * block_x)\n",
    "        offset_w = round(step_w * block_y)\n",
    "\n",
    "        profile = src.profile\n",
    "\n",
    "        width = math.ceil(step_w * (block_y + 1) - offset_w)\n",
    "        height = math.ceil(step_h * (block_x + 1) - offset_h)\n",
    "\n",
    "        profile.update(width=width)\n",
    "        profile.update(height=height)\n",
    "\n",
    "        window = rasterio.windows.Window(offset_w, offset_h, width, height)\n",
    "\n",
    "        chunk_file = os.path.join(dst, \"/tmp/\" + tile_id + '_' + str(block_x) + '_' + str(block_y) + '.tif')\n",
    "\n",
    "        with rasterio.open(chunk_file, 'w', **profile) as dest:\n",
    "            dest.write(src.read(window=window))\n",
    "\n",
    "        with open(chunk_file, 'rb') as tif_temp:\n",
    "            storage.put_cloudobject(tif_temp.read(), key=chunk_file.replace('/tmp', 'chunks'))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "* `make_plots`: function to make plots of the obtained data."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "def make_plots(array_file, metadata_file, j):\n",
    "    storage = Storage()\n",
    "\n",
    "    # Dowloading data\n",
    "    array_byte = storage.get_cloudobject(array_file)\n",
    "    array = pickle.loads(array_byte)\n",
    "\n",
    "    metadata_byte = storage.get_cloudobject(metadata_file)\n",
    "    metadata = pickle.loads(metadata_byte)\n",
    "\n",
    "    # Making plots\n",
    "    plt.clf()\n",
    "\n",
    "    plt.figure(1)\n",
    "\n",
    "    plot_band_array(array,metadata['extent'],\n",
    "                'Canopy height Model',\n",
    "                'Canopy height (m)',\n",
    "                'Greens',[0, 10])\n",
    "\n",
    "    # Path to save the plots\n",
    "    path = just_chm_file[0:-5]+'_CHM_' + str(j) + '.png'\n",
    "\n",
    "    plt.savefig(path,dpi=300,orientation='landscape',\n",
    "                bbox_inches='tight',\n",
    "                pad_inches=0.1)\n",
    "\n",
    "    # Updating plots to cloud storage\n",
    "    with open(path, 'rb') as png_temp:\n",
    "        storage.put_cloudobject(png_temp.read(), key=(\"plots/CHM/\" + path))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "* `smooth_CHM`: function smooth the CHM."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "#Smooth the CHM using a gaussian filter to remove spurious points\n",
    "def smooth_CHM(array_file, i):\n",
    "    # Dowloading data\n",
    "    storage = Storage()\n",
    "    array_byte = storage.get_cloudobject(array_file)\n",
    "\n",
    "    array = pickle.loads(array_byte)\n",
    "\n",
    "    array_smooth = ndi.gaussian_filter(array,2,mode='constant',cval=0,truncate=2.0)\n",
    "    array_smooth[array==0] = 0\n",
    "\n",
    "    # Updating plots to cloud storage\n",
    "    ref_array = storage.put_cloudobject(pickle.dumps(array_smooth), key=(\"arrayAux/smooth_array_object_\" + str(i)))\n",
    "\n",
    "    return ref_array"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "* `nan_to_num`: function to convert all nans position to zero in an array."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Converting nan position to zero\n",
    "def nan_to_num(array_file, i):\n",
    "    # Dowloading data\n",
    "    storage = Storage()\n",
    "\n",
    "    array_byte = storage.get_cloudobject(array_file)\n",
    "    array = pickle.loads(array_byte)\n",
    "\n",
    "    nan_array = np.nan_to_num(array)\n",
    "\n",
    "    # Updating plots to cloud storage\n",
    "    ref_array = storage.put_cloudobject(pickle.dumps(nan_array), key=(\"arrayAux/nan_array_object_\" + str(i)))\n",
    "\n",
    "    return ref_array"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "* `peak_max`: function to peak the position on an array with the max value."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "def peak_max(array_file, i):\n",
    "    # Dowloading data\n",
    "    storage = Storage()\n",
    "    array_byte = storage.get_cloudobject(array_file)\n",
    "    array = pickle.loads(array_byte)\n",
    "\n",
    "    array_peak = peak_local_max(array, indices=False, footprint=np.ones((5, 5)))\n",
    "\n",
    "    # Updating plots to cloud storage\n",
    "    ref_array = storage.put_cloudobject(pickle.dumps(array_peak), key=(\"arrayAux/local_maxi_array_object_\" + str(i)))\n",
    "\n",
    "    return ref_array"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "* `make_plots_local_maximus`: function to make plots of the local maximus data."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "#Plot the local maximums\n",
    "def make_plots_local_maximus(maxi_file, metadata_file, j):\n",
    "    # Dowloading data\n",
    "    storage = Storage()\n",
    "\n",
    "    maxi_byte = storage.get_cloudobject(maxi_file)\n",
    "    maxi = pickle.loads(maxi_byte)\n",
    "\n",
    "    metadata_byte = storage.get_cloudobject(metadata_file)\n",
    "    metadata = pickle.loads(metadata_byte)\n",
    "\n",
    "    # Making the plots\n",
    "    plt.clf()\n",
    "\n",
    "    plt.figure(2)\n",
    "    plot_band_array(maxi.astype(int),metadata['extent'],\n",
    "                    'Maximum',\n",
    "                    'Maxi',\n",
    "                    'Greys',\n",
    "                    [0, 1])\n",
    "\n",
    "    # Path to save the plots\n",
    "    path = just_chm_file[0:-4]+'_Maximums_' + str(j) + '.png'\n",
    "\n",
    "    plt.savefig(path,\n",
    "                dpi=300,orientation='landscape',\n",
    "                bbox_inches='tight',pad_inches=0.1)\n",
    "\n",
    "    # Updating plots to cloud storage\n",
    "    with open(path, 'rb') as png_temp:\n",
    "        storage.put_cloudobject(png_temp.read(), key=(\"plots/Maximums/\" + path))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "* `markers`: function to get the markers."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "def markers(maxi_file, i):\n",
    "    # Dowloading data\n",
    "    storage = Storage()\n",
    "\n",
    "    maxi_byte = storage.get_cloudobject(maxi_file)\n",
    "    maxi = pickle.loads(maxi_byte)\n",
    "\n",
    "    markers = ndi.label(maxi)[0]\n",
    "\n",
    "    # Updating plots to cloud storage\n",
    "    ref_marker = storage.put_cloudobject(pickle.dumps(markers), key=(\"arrayAux/marker_object_\" + str(i)))\n",
    "\n",
    "    return ref_marker"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "* `perform_watershed`: function to perform the watershed of our data."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "def perform_watershed(smooth_file, marker_file, mask_file, i):\n",
    "    # Dowloading data\n",
    "    storage = Storage()\n",
    "\n",
    "    smooth_byte = storage.get_cloudobject(smooth_file)\n",
    "    smooth = pickle.loads(smooth_byte)\n",
    "\n",
    "    marker_byte = storage.get_cloudobject(marker_file)\n",
    "    marker = pickle.loads(marker_byte)\n",
    "\n",
    "    mask_byte = storage.get_cloudobject(mask_file)\n",
    "    array_mask = pickle.loads(mask_byte)\n",
    "\n",
    "\n",
    "    labels = watershed(smooth, marker, mask=array_mask)\n",
    "\n",
    "    labels_for_plot = labels.copy()\n",
    "    labels_for_plot = np.array(labels_for_plot,dtype = np.float32)\n",
    "    labels_for_plot[labels_for_plot==0] = np.nan\n",
    "\n",
    "    max_labels = np.max(labels)\n",
    "\n",
    "    # Updating plots to cloud storage\n",
    "    ref_labels = storage.put_cloudobject(pickle.dumps(labels), key=(\"arrayAux/labels_object_\" + str(i)))\n",
    "\n",
    "    ref_labels_for_plot = storage.put_cloudobject(pickle.dumps(labels_for_plot), key=(\"arrayAux/labels_for_plot_object_\" + str(i)))\n",
    "\n",
    "    ref_max_labels = storage.put_cloudobject(pickle.dumps(max_labels), key=(\"arrayAux/max_labels_object_\" + str(i)))\n",
    "\n",
    "    return ref_labels, ref_labels_for_plot, ref_max_labels"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "* `create_mask`: function to create a mask with the smooth array."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "def create_mask(smooth_file, i):\n",
    "    # Dowloading data\n",
    "    storage = Storage()\n",
    "\n",
    "    smooth_byte = storage.get_cloudobject(smooth_file)\n",
    "    smooth = pickle.loads(smooth_byte)\n",
    "\n",
    "    chm_mask = smooth\n",
    "    chm_mask[smooth != 0] = 1\n",
    "\n",
    "    # Updating plots to cloud storage\n",
    "    ref_mask = storage.put_cloudobject(pickle.dumps(chm_mask), key=(\"arrayAux/mask_object_\" + str(i)))\n",
    "\n",
    "    return ref_mask"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "* `make_plots_segments`: function to make plots of the segments."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "def make_plots_segments(labels_for_plots_file, metadata_file, max_labels_file, j):\n",
    "    # Dowloading data\n",
    "    storage = Storage()\n",
    "\n",
    "    labels_for_plots_byte = storage.get_cloudobject(labels_for_plots_file)\n",
    "    labels_for_plots = pickle.loads(labels_for_plots_byte)\n",
    "\n",
    "    metadata_byte = storage.get_cloudobject(metadata_file)\n",
    "    metadata = pickle.loads(metadata_byte)\n",
    "\n",
    "    max_labels_byte = storage.get_cloudobject(max_labels_file)\n",
    "    max_labels = pickle.loads(max_labels_byte)\n",
    "\n",
    "\n",
    "    plt.clf()\n",
    "\n",
    "    plt.figure(2)\n",
    "    plot_band_array(labels_for_plots,metadata['extent'],\n",
    "                'Crown Segmentation','Tree Crown Number',\n",
    "                'Spectral',[0, max_labels])\n",
    "\n",
    "    # Path to save the plots\n",
    "    path = just_chm_file[0:-4]+'_Segmentation_' + str(j) + '.png'\n",
    "\n",
    "    plt.savefig(path,\n",
    "                dpi=300,orientation='landscape',\n",
    "                bbox_inches='tight',pad_inches=0.1)\n",
    "\n",
    "    # Updating plots to cloud storage\n",
    "    with open(path, 'rb') as png_temp:\n",
    "        storage.put_cloudobject(png_temp.read(), key=(\"plots/Segmentation/\" + path))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "* `region_props`: function to get the properties of each segment. "
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "def region_props(labels_file, array_file, i):\n",
    "    # Dowloading data\n",
    "    storage = Storage()\n",
    "\n",
    "    labels_byte = storage.get_cloudobject(labels_file)\n",
    "    labels = pickle.loads(labels_byte)\n",
    "\n",
    "    array_byte = storage.get_cloudobject(array_file)\n",
    "    array = pickle.loads(array_byte)\n",
    "\n",
    "\n",
    "    tree_properties = regionprops(labels,array)\n",
    "\n",
    "    predictors_chm = np.array([get_predictors(tree, array, labels) for tree in tree_properties])\n",
    "\n",
    "\n",
    "    if (not predictors_chm.shape[0] == 0):\n",
    "        X = np.nan_to_num(predictors_chm[:,1:])\n",
    "        tree_ids = np.nan_to_num(predictors_chm[:,0])\n",
    "    else:\n",
    "        X = np.resize(predictors_chm,(1,11))\n",
    "        tree_ids = np.resize(predictors_chm,(1,11))\n",
    "\n",
    "\n",
    "    # Updating plots to cloud storage\n",
    "    ref_tree = storage.put_cloudobject(pickle.dumps(tree_properties), key=(\"arrayAux/tree_properties_object_\" + str(i)))\n",
    "\n",
    "    ref_X = storage.put_cloudobject(pickle.dumps(X), key=(\"arrayAux/X_object_\" + str(i)))\n",
    "\n",
    "    ref_tree_ids = storage.put_cloudobject(pickle.dumps(tree_ids), key=(\"arrayAux/tree_ids_object_\" + str(i)))\n",
    "\n",
    "    return ref_tree, ref_X, ref_tree_ids"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "* `calculate_estimated_biomass`: function to apply the model to the predictors."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "def calculate_estimated_biomass(X_file, biomass_predictors_file, biomass_file, i):\n",
    "    # Dowloading data\n",
    "    storage = Storage()\n",
    "\n",
    "    X_byte = storage.get_cloudobject(X_file)\n",
    "    X = pickle.loads(X_byte)\n",
    "\n",
    "    biomass_predictors_byte = storage.get_cloudobject(biomass_predictors_file)\n",
    "    biomass_predictors = pickle.loads(biomass_predictors_byte)\n",
    "\n",
    "    biomass_byte = storage.get_cloudobject(biomass_file)\n",
    "    biomass = pickle.loads(biomass_byte)\n",
    "\n",
    "    regr_rf.fit(biomass_predictors,biomass)\n",
    "\n",
    "    estimated_biomass = regr_rf.predict(X)\n",
    "\n",
    "    # Updating plots to cloud storage\n",
    "    ref_estimated_biomass = storage.put_cloudobject(pickle.dumps(estimated_biomass), key=(\"arrayAux/estimated_biomass_object_\" + str(i)))\n",
    "\n",
    "    return ref_estimated_biomass"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "* `out_raster`: function to set an out raster with the same size as the labels."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "def out_raster(labels_file, tree_ids_file, estimated_biomass_file, i):\n",
    "    # Dowloading data\n",
    "    storage = Storage()\n",
    "\n",
    "    labels_byte = storage.get_cloudobject(labels_file)\n",
    "    labels = pickle.loads(labels_byte)\n",
    "\n",
    "    tree_ids_byte = storage.get_cloudobject(tree_ids_file)\n",
    "    tree_ids = pickle.loads(tree_ids_byte)\n",
    "\n",
    "    estimated_biomass_byte = storage.get_cloudobject(estimated_biomass_file)\n",
    "    estimated_biomass = pickle.loads(estimated_biomass_byte)\n",
    "\n",
    "\n",
    "    biomass_map =  np.array((labels),dtype=float)\n",
    "\n",
    "    #Assign the appropriate biomass to the labels\n",
    "    biomass_map[biomass_map==0] = np.nan\n",
    "    for tree_id, biomass_of_tree_id in zip(tree_ids, estimated_biomass):\n",
    "        biomass_map[biomass_map == tree_id] = biomass_of_tree_id\n",
    "\n",
    "    # Updating plots to cloud storage\n",
    "    ref_biomass_map = storage.put_cloudobject(pickle.dumps(biomass_map), key=(\"arrayAux/biomass_map_object_\" + str(i)))\n",
    "\n",
    "    return ref_biomass_map"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "* `biomass_stats`: function to get biomass stats for plotting."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "def biomass_stats(estimated_biomass_file, i):\n",
    "    # Dowloading data\n",
    "    storage = Storage()\n",
    "\n",
    "    estimated_biomass_byte = storage.get_cloudobject(estimated_biomass_file)\n",
    "    estimated_biomass = pickle.loads(estimated_biomass_byte)\n",
    "\n",
    "\n",
    "    mean_biomass = np.mean(estimated_biomass)\n",
    "    std_biomass = np.std(estimated_biomass)\n",
    "    min_biomass = np.min(estimated_biomass)\n",
    "    sum_biomass = np.sum(estimated_biomass)\n",
    "\n",
    "\n",
    "    # Updating plots to cloud storage\n",
    "    ref_mean_biomass = storage.put_cloudobject(pickle.dumps(mean_biomass), key=(\"arrayAux/mean_biomass_object_\" + str(i)))\n",
    "    ref_std_biomass = storage.put_cloudobject(pickle.dumps(std_biomass), key=(\"arrayAux/std_biomass_object_\" + str(i)))\n",
    "    ref_min_biomass = storage.put_cloudobject(pickle.dumps(min_biomass), key=(\"arrayAux/min_biomass_object_\" + str(i)))\n",
    "    ref_sum_biomass = storage.put_cloudobject(pickle.dumps(sum_biomass), key=(\"arrayAux/sum_biomass_object_\" + str(i)))\n",
    "\n",
    "    return ref_mean_biomass, ref_std_biomass, ref_min_biomass, ref_sum_biomass"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "* `make_plots_biomass`: function to make plots of the biomass data."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "def make_plots_biomass(biomass_map_file, metadata_file, mean_biomass_file, std_biomass_file, min_biomass_file, sum_biomass_file, i):\n",
    "    # Dowloading data\n",
    "    storage = Storage()\n",
    "\n",
    "    biomass_map_byte = storage.get_cloudobject(biomass_map_file)\n",
    "    biomass_map = pickle.loads(biomass_map_byte)\n",
    "\n",
    "    metadata_byte = storage.get_cloudobject(metadata_file)\n",
    "    metadata = pickle.loads(metadata_byte)\n",
    "\n",
    "    mean_biomass_byte = storage.get_cloudobject(mean_biomass_file)\n",
    "    mean_biomass = pickle.loads(mean_biomass_byte)\n",
    "\n",
    "    std_biomass_byte = storage.get_cloudobject(std_biomass_file)\n",
    "    std_biomass = pickle.loads(std_biomass_byte)\n",
    "\n",
    "    min_biomass_byte = storage.get_cloudobject(min_biomass_file)\n",
    "    min_biomass = pickle.loads(min_biomass_byte)\n",
    "\n",
    "    sum_biomass_byte = storage.get_cloudobject(sum_biomass_file)\n",
    "    sum_biomass = pickle.loads(sum_biomass_byte)\n",
    "\n",
    "\n",
    "    plt.clf()\n",
    "\n",
    "    plt.figure(5)\n",
    "    plot_band_array(biomass_map,metadata['extent'],\n",
    "                'Biomass (kg)','Biomass (kg)',\n",
    "                'winter',\n",
    "                [min_biomass+std_biomass, mean_biomass+std_biomass*3])\n",
    "\n",
    "    # Path to save the plots\n",
    "    path = 'ResultBiomass_' + str(i) + '.png'\n",
    "\n",
    "    plt.savefig(path,\n",
    "            dpi=300,orientation='landscape',\n",
    "            bbox_inches='tight',\n",
    "            pad_inches=0.1)\n",
    "\n",
    "\n",
    "    # Updating plots to cloud storage\n",
    "    with open(path, 'rb') as png_temp:\n",
    "        storage.put_cloudobject(png_temp.read(), key=(\"plots/resultBiomass/\" + path))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Data Preparation: Canopy Height Models"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "SPLITS = 50"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "chm_file = 'NDSM-Vegetacion-ETRS89-H31-0473-COB1.asc'"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "with open(chm_file, 'rb') as asc_temp:\n",
    "    storage.put_cloudobject(asc_temp.read(), key=chm_file)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Convert asc file to geotiff\n",
    "fexec = lithops.FunctionExecutor(runtime_memory=2048)\n",
    "\n",
    "fexec.map(asc_to_geotiff, (chm_file, BUCKET_NAME))\n",
    "\n",
    "fexec.get_result()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Preparing data\n",
    "chm_chunk_files = []\n",
    "iterdata = []\n",
    "for i in range(SPLITS):\n",
    "     for j in range(SPLITS):\n",
    "            file_name = 'chunks/' + chm_file.replace('.asc', '') + '_' + str(i) + '_' + str(j) + '.tif'\n",
    "            iterdata.append((chm_file,'chunks',SPLITS,i,j, BUCKET_NAME))\n",
    "            chm_chunk_files.append((file_name , BUCKET_NAME))\n",
    "\n",
    "\n",
    "# Transform the geotiff file to little data chunks\n",
    "fexec = lithops.FunctionExecutor(runtime_memory=1024)\n",
    "\n",
    "fexec.map(data_chunker, iterdata)\n",
    "\n",
    "fexec.get_result()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "When we output the results, we will want to include the same file information as the input, so we will gather the file name information."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "#Get info from chm file for outputting results\n",
    "just_chm_file = os.path.basename(chm_file)\n",
    "print(just_chm_file)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now we will get the CHM data..."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# chm_chunk_files"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Converting chunk file to numpy array\n",
    "fexec = lithops.FunctionExecutor()\n",
    "\n",
    "fexec.map(raster2array, chm_chunk_files)\n",
    "\n",
    "results = fexec.get_result()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "..., plot it, and save the figure."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Preparing data\n",
    "i = 0\n",
    "iterdata = []\n",
    "array_only = []\n",
    "metadata_only = []\n",
    "\n",
    "for array, metadata in results:\n",
    "    array_only.append(array)\n",
    "    metadata_only.append(metadata)\n",
    "    iterdata.append((array, metadata, i))\n",
    "    i = i + 1\n",
    "\n",
    "#Plot the original CHM\n",
    "fexec = lithops.FunctionExecutor()\n",
    "\n",
    "fexec.map(make_plots, iterdata)\n",
    "\n",
    "results = fexec.get_result()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "It looks like SJER primarily has low vegetation with scattered taller trees.\n",
    "\n",
    "## Create Filtered CHM\n",
    "\n",
    "Now we will use a Gaussian smoothing kernal (convolution) across the data set to remove spurious high vegetation points. This will help ensure we are finding the treetops properly before running the watershed segmentation algorithm.\n",
    "\n",
    "For different forest types it may be necessary to change the input parameters. Information on the function can be found in the <a href=\"https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.ndimage.filters.gaussian_filter.html\" target=\"_blank\">SciPy documentation</a>.\n",
    "\n",
    "Of most importance are the second and fifth inputs. The second input defines the standard deviation of the Gaussian smoothing kernal. Too large a value will apply too much smoothing, too small and some spurious high points may be left behind. The fifth, the truncate value, controls after how many standard deviations the Gaussian kernal will get cut off (since it theoretically goes to infinity)."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Preparing data\n",
    "iterdata = []\n",
    "i = 0\n",
    "for array in array_only:\n",
    "    iterdata.append((array, i))\n",
    "    i = i + 1\n",
    "\n",
    "# Remove spurious points from array\n",
    "fexec = lithops.FunctionExecutor()\n",
    "\n",
    "fexec.map(smooth_CHM, iterdata)\n",
    "\n",
    "results = fexec.get_result()\n",
    "\n",
    "array_smooth_only = results"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now save a copy of filtered CHM. We will later use this in our code, so we'll output it into our data directory."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Preparing data\n",
    "\n",
    "iterdata = []\n",
    "i = 0\n",
    "\n",
    "for metadata, smooth in zip(metadata_only, array_smooth_only):\n",
    "    iterdata.append(('chm_filter_' + str(i) + '.tif',metadata,1,-1,smooth,32611))\n",
    "    i = i + 1\n",
    "\n",
    "#Save the smoothed CHM\n",
    "fexec = lithops.FunctionExecutor()\n",
    "\n",
    "fexec.map(array2raster, iterdata)\n",
    "\n",
    "results = fexec.get_result()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "#Plot the filtered CHM\n",
    "i = len(metadata_only) + 1\n",
    "iterdata = []\n",
    "\n",
    "for metadata, smooth in zip(metadata_only, array_smooth_only):\n",
    "    iterdata.append((smooth, metadata, i))\n",
    "    i = i + 1\n",
    "\n",
    "fexec = lithops.FunctionExecutor()\n",
    "\n",
    "fexec.map(make_plots, iterdata)\n",
    "\n",
    "results = fexec.get_result()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Preparing data\n",
    "iterdata = []\n",
    "i = 0\n",
    "\n",
    "for array in array_only:\n",
    "    iterdata.append((array, i))\n",
    "    i = i + 1\n",
    "\n",
    "fexec = lithops.FunctionExecutor()\n",
    "\n",
    "fexec.map(nan_to_num, iterdata)\n",
    "\n",
    "array_smooth_only = fexec.get_result()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# array_smooth_only"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Determine local maximums\n",
    "\n",
    "Now we will run an algorithm to determine local maximums within the image. Setting indices to 'False' returns a raster of the maximum points, as opposed to a list of coordinates. The footprint parameter is an area where only a single peak can be found. This should be approximately the size of the smallest tree. Information on more sophisticated methods to define the window can be found in Chen (2006)."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Calculate local maximum points in the smoothed CHM\n",
    "# Preparing data\n",
    "iterdata = []\n",
    "i = 0\n",
    "\n",
    "for array in array_smooth_only:\n",
    "    iterdata.append((array, i))\n",
    "    i = i + 1\n",
    "\n",
    "fexec = lithops.FunctionExecutor(runtime_memory=2048)\n",
    "\n",
    "fexec.map(peak_max, iterdata)\n",
    "\n",
    "results = fexec.get_result()\n",
    "\n",
    "local_maxi = results"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Our new object `local_maxi` is an array of boolean values where each pixel is identified as either being the local maximum (`True`) or not being the local maximum (`False`)."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# local_maxi"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "This is very helpful, but it can be difficult to visualizee boolean values using our typical numeric plotting procedures as defined in the `plot_band_array` function above. Therefore, we will need to convert this boolean array to an numeric format to use this function. Booleans convert easily to integers with values of `False=0` and `True=1` using the `.astype(int)` method."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# local_maxi.astype(int)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Next ,we can plot the raster of local maximums bo coercing the boolean array into an array ofintegers inline. The following figure shows the difference in finding local maximums for a filtered vs. non-filtered CHM.\n",
    "\n",
    "We will save the graphics (.png) in an outputs folder sister to our working directory and data outputs (.tif) to our data directory."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Preparing data\n",
    "iterdata = []\n",
    "j = 0\n",
    "\n",
    "for maxi, metadata in zip(local_maxi, metadata_only):\n",
    "    iterdata.append((maxi, metadata, j))\n",
    "    j = j + 1\n",
    "\n",
    "fexec = lithops.FunctionExecutor()\n",
    "\n",
    "# Making plots of local_maximus\n",
    "fexec.map(make_plots_local_maximus, iterdata)\n",
    "\n",
    "fexec.get_result()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Preparing data\n",
    "iterdata = []\n",
    "i = 0\n",
    "\n",
    "for metadata, maxi in zip(metadata_only, local_maxi):\n",
    "    iterdata.append(('maximum_' + str(i) + '.tif',metadata,1,-1,maxi,32611))\n",
    "    i = i + 1\n",
    "\n",
    "# Converting array to geotiff file\n",
    "fexec = lithops.FunctionExecutor()\n",
    "\n",
    "fexec.map(array2raster, iterdata)\n",
    "\n",
    "fexec.get_result()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "If we were to look at the overlap between the tree crowns and the local maxima from each method, it would appear a bit like this raster.\n",
    "\n",
    " <figure>\n",
    "\t<a href=\"https://raw.githubusercontent.com/NEONScience/NEON-Data-Skills/main/graphics/raster-general/raster-classification-filter-vs-nonfilter.jpg\">\n",
    "\t<img src=\"https://raw.githubusercontent.com/NEONScience/NEON-Data-Skills/main/graphics/raster-general/raster-classification-filter-vs-nonfilter.jpg\"></a>\n",
    "\t<figcaption> The difference in finding local maximums for a filtered vs.\n",
    "\tnon-filtered CHM.\n",
    "\tSource: National Ecological Observatory Network (NEON)\n",
    "\t</figcaption>\n",
    "</figure>\n",
    "\n",
    "\n",
    "Apply labels to all of the local maximum points"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "#Identify all the maximum points\n",
    "iterdata = []\n",
    "i = 0\n",
    "\n",
    "for maxi in local_maxi:\n",
    "    iterdata.append((maxi, i))\n",
    "    i = i + 1\n",
    "\n",
    "fexec = lithops.FunctionExecutor()\n",
    "\n",
    "fexec.map(markers, iterdata)\n",
    "\n",
    "markers_only = fexec.get_result()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Next we will create a mask layer of all of the vegetation points so that the watershed segmentation will only occur on the trees and not extend into the surrounding ground points. Since 0 represent ground points in the CHM, setting the mask to 1 where the CHM is not zero will define the mask"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "#Create a CHM mask so the segmentation will only occur on the trees\n",
    "fexec = lithops.FunctionExecutor()\n",
    "\n",
    "fexec.map(create_mask, iterdata)\n",
    "\n",
    "mask_only = fexec.get_result()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Watershed segmentation\n",
    "\n",
    "As in a river system, a watershed is divided by a ridge that divides areas. Here our watershed are the individual tree canopies and the ridge is the delineation between each one.\n",
    "\n",
    "<figure>\n",
    "\t<a href=\"https://raw.githubusercontent.com/NEONScience/NEON-Data-Skills/main/graphics/raster-general/raster-classification-watershed-segments.png\">\n",
    "\t<img src=\"https://raw.githubusercontent.com/NEONScience/NEON-Data-Skills/main/graphics/raster-general/raster-classification-watershed-segments.png\"></a>\n",
    "\t<figcaption> A raster classified based on watershed segmentation.\n",
    "\tSource: National Ecological Observatory Network (NEON)\n",
    "\t</figcaption>\n",
    "</figure>\n",
    "\n",
    "Next, we will perform the watershed segmentation which produces a raster of labels."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Preparing data\n",
    "iterdata = []\n",
    "i = 0\n",
    "\n",
    "for smooth, mark, mask in zip(array_smooth_only, markers_only, mask_only):\n",
    "    iterdata.append((smooth, mark, mask, i))\n",
    "    i = i + 1\n",
    "\n",
    "fexec = lithops.FunctionExecutor()\n",
    "\n",
    "fexec.map(perform_watershed, iterdata)\n",
    "\n",
    "labels_only = []\n",
    "labels_for_plot_only = []\n",
    "max_labels_only = []\n",
    "\n",
    "results = fexec.get_result()\n",
    "for labels, labels_for_plot, max_labels in results:\n",
    "    labels_only.append(labels)\n",
    "    labels_for_plot_only.append(labels_for_plot)\n",
    "    max_labels_only.append(max_labels)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Preparing data\n",
    "iterdata = []\n",
    "i = 0\n",
    "\n",
    "for labels_for_plots, metadata, max_labels in zip(labels_for_plot_only, metadata_only, max_labels_only):\n",
    "    iterdata.append((labels_for_plots, metadata, max_labels, i))\n",
    "    i = i + 1\n",
    "\n",
    "# Making plots of segments\n",
    "fexec = lithops.FunctionExecutor()\n",
    "\n",
    "fexec.map(make_plots_segments, iterdata)\n",
    "\n",
    "fexec.get_result()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Preparing data\n",
    "iterdata = []\n",
    "i = 0\n",
    "\n",
    "for labels, metadata in zip(labels_only, metadata_only):\n",
    "    iterdata.append(('labels' + str(i) + '.tif',metadata,1,-1,labels,32611))\n",
    "    i = i + 1\n",
    "\n",
    "# Converting labels to geotiff file\n",
    "fexec = lithops.FunctionExecutor()\n",
    "\n",
    "fexec.map(array2raster, iterdata)\n",
    "\n",
    "fexec.get_result()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now we will get several properties of the individual trees will be used as predictor variables."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Preparing data\n",
    "iterdata = []\n",
    "i = 0\n",
    "\n",
    "for labels, array in zip(labels_only, array_only):\n",
    "    iterdata.append((labels, array, i))\n",
    "    i = i + 1\n",
    "\n",
    "fexec = lithops.FunctionExecutor()\n",
    "\n",
    "fexec.map(region_props, iterdata)\n",
    "\n",
    "results = fexec.get_result()\n",
    "\n",
    "tree_properties_only = []\n",
    "X_only = []\n",
    "tree_ids_only = []\n",
    "\n",
    "for tree_properties, X, tree_ids in results:\n",
    "    tree_properties_only.append(tree_properties)\n",
    "    X_only.append(X)\n",
    "    tree_ids_only.append(tree_ids)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "cloudbutton-geospatial",
   "language": "python",
   "name": "cloudbutton-geospatial"
  },
  "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}