<commit_before><commit_msg>Kmeans light, NaN are simply converted to -1. Hence, no need to indexes and joins to reconstruct the GeoTiff with clusters Id.<commit_after>{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kmeans over a set of GeoTiffs\n",
    "\n",
    "This notebook loads a set of GeoTiffs into a **RDD** of Tiles, with each Tile being a band in the GeoTiff. Each GeoTiff file contains **SpringIndex-** or **LastFreeze-** value for one year over the entire USA.\n",
    "\n",
    "Kmeans takes years as dimensions. Hence, the matrix has cells as rows and the years as columns. To cluster on all years, the matrix needs to be transposed. The notebook has two flavors of matrix transpose, locally by the Spark-driver or distributed using the Spark-workers. Once transposed the matrix is converted to a **RDD** of dense vectors to be used by **Kmeans** algorithm from **Spark-MLlib**. The end result is a grid where each cell has a cluster ID which is then saved into a SingleBand GeoTiff. By saving the result into a GeoTiff, the reader can plot it using a Python notebook as the one defined in the [python examples](../examples/python).\n",
    "\n",
    "<span style=\"color:red\">In this notebook the reader only needs to modify the variables in **Mode of Operation Setup**</span>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dependencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "...\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import sys.process._\n",
    "\n",
    "import java.io.{ByteArrayInputStream, ByteArrayOutputStream, ObjectInputStream, ObjectOutputStream}\n",
    "\n",
    "\n",
    "import geotrellis.proj4.CRS\n",
    "import geotrellis.raster.{CellType, ArrayTile, DoubleArrayTile, Tile, UByteCellType}\n",
    "import geotrellis.raster.io.geotiff._\n",
    "import geotrellis.raster.io.geotiff.writer.GeoTiffWriter\n",
    "import geotrellis.raster.io.geotiff.{GeoTiff, SinglebandGeoTiff}\n",
    "import geotrellis.spark.io.hadoop._\n",
    "import org.apache.hadoop.io._\n",
    "import geotrellis.vector.{Extent, ProjectedExtent}\n",
    "import org.apache.spark.broadcast.Broadcast\n",
    "import org.apache.spark.mllib.clustering.{KMeans, KMeansModel}\n",
    "import org.apache.spark.mllib.linalg.distributed.{CoordinateMatrix, MatrixEntry, RowMatrix}\n",
    "import org.apache.spark.mllib.linalg.{Vector, Vectors}\n",
    "import org.apache.spark.rdd.RDD\n",
    "import org.apache.spark.{SparkConf, SparkContext}\n",
    "\n",
    "import org.apache.hadoop.io.{IOUtils, SequenceFile}\n",
    "import org.apache.hadoop.io.SequenceFile.Writer\n",
    "\n",
    "//Spire is a numeric library for Scala which is intended to be generic, fast, and precise.\n",
    "import spire.syntax.cfor._"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mode of operation\n",
    "\n",
    "Here the user can define the mode of operation.\n",
    "* **rdd_offline_mode**: If false it means the notebook will create all data from scratch and store protected_extent and num_cols_rows into HDFS. Otherwise, these data structures are read from HDFS.\n",
    "* **matrix_offline_mode**: If false it means the notebook will create a mtrix,  transposed it and save it to HDFS. Otherwise, these data structures are read from HDFS.\n",
    "* **kmeans_offline_mode**: If false it means the notebook will train kmeans and run kemans and store kmeans model into HDFS. Otherwise, these data structures are read from HDFS.\n",
    "\n",
    "It is also possible to define which directory of GeoTiffs is to be used and on which **band** to run Kmeans. The options are\n",
    "* **BloomFinal** or **LeafFinal** which are multi-band (**4 bands**)\n",
    "* **DamageIndex** and **LastFreeze** which are single-band and if set band_num higher, it will reset to 0\n",
    "\n",
    "For kmeans the user can define the **number of iterations** and **number of clusters** as an inclusive range. Such range is defined using **minClusters**, **maxClusters**, and **stepClusters**. These variables will set a loop starting at **minClusters** and stopping at **maxClusters** (inclusive), iterating **stepClusters** at the time. <span style=\"color:red\">Note that when using a range **kemans offline mode** is not possible and it will be reset to **online mode**</span>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mode of Operation setup\n",
    "<a id='mode_of_operation_setup'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "rdd_offline_mode = true\n",
       "matrix_offline_mode = true\n",
       "kmeans_offline_mode = true\n",
       "dir_path = hdfs:///user/hadoop/spring-index/\n",
       "offline_dir_path = hdfs:///user/emma/spring-index/\n",
       "geoTiff_dir = BloomFinal\n",
       "band_num = 4\n",
       "toBeMasked = true\n",
       "mask_path = hdfs:///user/hadoop/usa_mask.tif\n",
       "numIterations = 35\n",
       "minClusters = 3\n",
       "maxClusters = 16\n",
       "stepClusters = 1\n",
       "save_kmeans_model = false\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "false"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "//Operation mode\n",
    "var rdd_offline_mode = true\n",
    "var matrix_offline_mode = true\n",
    "var kmeans_offline_mode = true\n",
    "\n",
    "//GeoTiffs to be read from \"hdfs:///user/hadoop/spring-index/\"\n",
    "var dir_path = \"hdfs:///user/hadoop/spring-index/\"\n",
    "var offline_dir_path = \"hdfs:///user/emma/spring-index/\"\n",
    "var geoTiff_dir = \"BloomFinal\"\n",
    "var band_num = 4\n",
    "\n",
    "//Mask\n",
    "val toBeMasked = true\n",
    "val mask_path = \"hdfs:///user/hadoop/usa_mask.tif\"\n",
    "\n",
    "//Kmeans number of iterations and clusters\n",
    "var numIterations = 35\n",
    "var minClusters = 13\n",
    "var maxClusters = 13\n",
    "var stepClusters = 1\n",
    "var save_kmeans_model = false"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<span style=\"color:red\">DON'T MODIFY ANY PIECE OF CODE FROM HERE ON!!!</span>.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mode of operation validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Waiting for a Spark session to start..."
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\"Load GeoTiffs\" offline mode is not set properly, i.e., either it was set to false and the required file does not exist or vice-versa. We will reset it to false\n",
      "\"Matrix\" offline mode is not set properly, i.e., either it was set to false and the required file does not exist or vice-versa. We will reset it to false\n",
      "14\n",
      "There is already a GeoTiff with the path: hdfs:///user/emma/spring-index/BloomFinal/clusters_3_35.tif. Please make either a copy or move it to another location, otherwise, it will be over-written.\n",
      "There is already a GeoTiff with the path: hdfs:///user/emma/spring-index/BloomFinal/clusters_4_35.tif. Please make either a copy or move it to another location, otherwise, it will be over-written.\n",
      "There is already a GeoTiff with the path: hdfs:///user/emma/spring-index/BloomFinal/clusters_5_35.tif. Please make either a copy or move it to another location, otherwise, it will be over-written.\n",
      "There is already a GeoTiff with the path: hdfs:///user/emma/spring-index/BloomFinal/clusters_6_35.tif. Please make either a copy or move it to another location, otherwise, it will be over-written.\n",
      "There is already a GeoTiff with the path: hdfs:///user/emma/spring-index/BloomFinal/clusters_7_35.tif. Please make either a copy or move it to another location, otherwise, it will be over-written.\n",
      "There is already a GeoTiff with the path: hdfs:///user/emma/spring-index/BloomFinal/clusters_8_35.tif. Please make either a copy or move it to another location, otherwise, it will be over-written.\n",
      "There is already a GeoTiff with the path: hdfs:///user/emma/spring-index/BloomFinal/clusters_9_35.tif. Please make either a copy or move it to another location, otherwise, it will be over-written.\n",
      "There is already a GeoTiff with the path: hdfs:///user/emma/spring-index/BloomFinal/clusters_10_35.tif. Please make either a copy or move it to another location, otherwise, it will be over-written.\n",
      "There is already a GeoTiff with the path: hdfs:///user/emma/spring-index/BloomFinal/clusters_11_35.tif. Please make either a copy or move it to another location, otherwise, it will be over-written.\n",
      "There is already a GeoTiff with the path: hdfs:///user/emma/spring-index/BloomFinal/clusters_12_35.tif. Please make either a copy or move it to another location, otherwise, it will be over-written.\n",
      "There is already a GeoTiff with the path: hdfs:///user/emma/spring-index/BloomFinal/clusters_13_35.tif. Please make either a copy or move it to another location, otherwise, it will be over-written.\n",
      "There is already a GeoTiff with the path: hdfs:///user/emma/spring-index/BloomFinal/clusters_14_35.tif. Please make either a copy or move it to another location, otherwise, it will be over-written.\n",
      "There is already a GeoTiff with the path: hdfs:///user/emma/spring-index/BloomFinal/clusters_15_35.tif. Please make either a copy or move it to another location, otherwise, it will be over-written.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "single_band = false\n",
       "mask_str = _mask\n",
       "grids_noNaN_path = hdfs:///user/emma/spring-index/BloomFinal/grids_noNaN_mask\n",
       "metadata_path = hdfs:///user/emma/spring-index/BloomFinal/metadata_mask\n",
       "grids_matrix_path = hdfs:///user/emma/spring-index/BloomFinal/grids_matrix_mask\n",
       "conf = Configuration: core-default.xml, core-site.xml, mapred-default.xml, mapred-site.xml, yarn-default.xml, yarn-site.xml, hdfs-default.xml, hdfs-site.xml, file:/usr/lib/spark-2.1.1-bin-without-hadoop/conf/hive-site.xml\n",
       "fs = DFS[DFSClient[clientName=DFSClient_NONMAPREDUCE_934661335_39, ugi=emma (auth:SIMPLE)]]\n",
       "rdd_offline_exists = false\n",
       "matrix_offline_exists = false\n",
       "num_kmeans = 14\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "kmeans_m...\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "14"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "//Validation, do not modify these lines.\n",
    "var single_band = false\n",
    "if (geoTiff_dir == \"BloomFinal\" || geoTiff_dir == \"LeafFinal\") {\n",
    "    single_band = false\n",
    "} else if (geoTiff_dir == \"LastFreeze\" || geoTiff_dir == \"DamageIndex\") {\n",
    "    single_band = true\n",
    "    if (band_num > 0) {\n",
    "        println(\"Since LastFreezze and DamageIndex are single band, we will use band 0!!!\")\n",
    "        band_num  = 0\n",
    "    }\n",
    "} else {\n",
    "    println(\"Directory unknown, please set either BloomFinal, LeafFinal, LastFreeze or DamageIndex!!!\")\n",
    "}\n",
    "\n",
    "if (minClusters > maxClusters) {\n",
    "    maxClusters = minClusters\n",
    "    stepClusters = 1\n",
    "}\n",
    "if (stepClusters < 1) {\n",
    "    stepClusters = 1\n",
    "}\n",
    "\n",
    "//Paths to store data structures for Offline runs\n",
    "var mask_str = \"\"\n",
    "if (toBeMasked)\n",
    "    mask_str = \"_mask\"\n",
    "var grids_noNaN_path = offline_dir_path + geoTiff_dir + \"/grids_noNaN\" + mask_str\n",
    "var metadata_path = offline_dir_path + geoTiff_dir + \"/metadata\" + mask_str\n",
    "var grids_matrix_path = offline_dir_path + geoTiff_dir + \"/grids_matrix\" + mask_str\n",
    "\n",
    "//Check offline modes\n",
    "var conf = sc.hadoopConfiguration\n",
    "var fs = org.apache.hadoop.fs.FileSystem.get(conf)\n",
    "\n",
    "val rdd_offline_exists = fs.exists(new org.apache.hadoop.fs.Path(grids_noNaN_path))\n",
    "val matrix_offline_exists = fs.exists(new org.apache.hadoop.fs.Path(grids_matrix_path))\n",
    "                                      \n",
    "if (rdd_offline_mode != rdd_offline_exists) {\n",
    "    println(\"\\\"Load GeoTiffs\\\" offline mode is not set properly, i.e., either it was set to false and the required file does not exist or vice-versa. We will reset it to \" + rdd_offline_exists.toString())\n",
    "    rdd_offline_mode = rdd_offline_exists\n",
    "} \n",
    "if (matrix_offline_mode != matrix_offline_exists) {\n",
    "    println(\"\\\"Matrix\\\" offline mode is not set properly, i.e., either it was set to false and the required file does not exist or vice-versa. We will reset it to \" + matrix_offline_exists.toString())\n",
    "    matrix_offline_mode = matrix_offline_exists\n",
    "}\n",
    "\n",
    "if (!fs.exists(new org.apache.hadoop.fs.Path(mask_path))) {\n",
    "    println(\"The mask path: \" + mask_path + \" is invalid!!!\")\n",
    "}\n",
    "\n",
    "var num_kmeans :Int  = 1\n",
    "if (minClusters != maxClusters) {\n",
    "    num_kmeans = ((maxClusters - minClusters) / stepClusters) + 1\n",
    "}\n",
    "println(num_kmeans)\n",
    "var kmeans_model_paths :Array[String] = Array.fill[String](num_kmeans)(\"\")\n",
    "var wssse_path :String = offline_dir_path + geoTiff_dir + \"/wssse\"\n",
    "var geotiff_hdfs_paths :Array[String] = Array.fill[String](num_kmeans)(\"\")\n",
    "var geotiff_tmp_paths :Array[String] = Array.fill[String](num_kmeans)(\"\")\n",
    "\n",
    "if (num_kmeans > 1) {\n",
    "    var numClusters_id = 0\n",
    "    cfor(minClusters)(_ <= maxClusters, _ + stepClusters) { numClusters =>\n",
    "        kmeans_model_paths(numClusters_id) = offline_dir_path + geoTiff_dir + \"/kmeans_model_\" + numClusters + \"_\" + numIterations\n",
    "        \n",
    "        //Check if the file exists\n",
    "        val kmeans_exist = fs.exists(new org.apache.hadoop.fs.Path(kmeans_model_paths(numClusters_id)))\n",
    "        if (kmeans_exist && !kmeans_offline_mode) {\n",
    "            println(\"The kmeans model path \" + kmeans_model_paths(numClusters_id) + \" exists, please remove it.\")\n",
    "        } else if (!kmeans_exist && kmeans_offline_mode) {\n",
    "            kmeans_offline_mode = false\n",
    "        }\n",
    "        \n",
    "        geotiff_hdfs_paths(numClusters_id) = offline_dir_path + geoTiff_dir + \"/clusters_\" + numClusters + \"_\" + numIterations + \".tif\"\n",
    "        geotiff_tmp_paths(numClusters_id) = \"/tmp/clusters_\" + geoTiff_dir + \"_\" + numClusters + \"_\" + numIterations + \".tif\"\n",
    "        if (fs.exists(new org.apache.hadoop.fs.Path(geotiff_hdfs_paths(numClusters_id)))) {\n",
    "            println(\"There is already a GeoTiff with the path: \" + geotiff_hdfs_paths(numClusters_id) + \". Please make either a copy or move it to another location, otherwise, it will be over-written.\")\n",
    "        }\n",
    "        numClusters_id += 1\n",
    "    }\n",
    "    kmeans_offline_mode = false\n",
    "} else { \n",
    "    kmeans_model_paths(0) = offline_dir_path + geoTiff_dir + \"/kmeans_model_\" + minClusters + \"_\" + numIterations\n",
    "    val kmeans_offline_exists = fs.exists(new org.apache.hadoop.fs.Path(kmeans_model_paths(0)))\n",
    "    if (kmeans_offline_mode != kmeans_offline_exists) {\n",
    "        println(\"\\\"Kmeans\\\" offline mode is not set properly, i.e., either it was set to false and the required file does not exist or vice-versa. We will reset it to \" + kmeans_offline_exists.toString())\n",
    "        kmeans_offline_mode = kmeans_offline_exists\n",
    "    }\n",
    "    geotiff_hdfs_paths(0) = offline_dir_path + geoTiff_dir + \"/clusters_\" + minClusters + \"_\" + numIterations + \".tif\"\n",
    "    geotiff_tmp_paths(0) = \"/tmp/clusters_\" + geoTiff_dir + \"_\" + minClusters + \"_\" + numIterations + \".tif\"\n",
    "    if (fs.exists(new org.apache.hadoop.fs.Path(geotiff_hdfs_paths(0)))) {\n",
    "        println(\"There is already a GeoTiff with the path: \" + geotiff_hdfs_paths(0) + \". Please make either a copy or move it to another location, otherwise, it will be over-written.\")\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Functions to (de)serialize any structure into Array[Byte]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "serialize: (value: Any)Array[Byte]\n",
       "deserialize: (bytes: Array[Byte])Any\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def serialize(value: Any): Array[Byte] = {\n",
    "    val out_stream: ByteArrayOutputStream = new ByteArrayOutputStream()\n",
    "    val obj_out_stream = new ObjectOutputStream(out_stream)\n",
    "    obj_out_stream.writeObject(value)\n",
    "    obj_out_stream.close\n",
    "    out_stream.toByteArray\n",
    "}\n",
    "\n",
    "def deserialize(bytes: Array[Byte]): Any = {\n",
    "    val obj_in_stream = new ObjectInputStream(new ByteArrayInputStream(bytes))\n",
    "    val value = obj_in_stream.readObject\n",
    "    obj_in_stream.close\n",
    "    value\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load GeoTiffs\n",
    "\n",
    "Using GeoTrellis all GeoTiffs of a directory will be loaded into a RDD. Using the RDD, we extract a grid from the first file to lated store the Kmeans cluster_IDS, we build an Index for populate such grid and we convert all NaN values to **-1**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Stage 5:=======================================================> (35 + 1) / 36]Elapsed time: 720700223515ns\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "t0 = 116568168752299\n",
       "projected_extent = ProjectedExtent(Extent(-126.30312894720473, 14.29219617034159, -56.162671563152486, 49.25462702827337),geotrellis.proj4.CRS$$anon$3@41d0d1b7)\n",
       "grids_noNaN_RDD = MapPartitionsRDD[21] at map at <console>:138\n",
       "num_cols_rows = (7808,3892)\n",
       "cellT = float64raw\n",
       "mask_tile0 = FloatRawArrayTile([F@c44c716,7808,3892)\n",
       "pattern = tif\n",
       "filepath = hdfs:///user/hadoop/spring-index/BloomFinal\n",
       "t1 = 117288868975814\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "117288868975814"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val t0 = System.nanoTime()\n",
    "//Global variables\n",
    "var projected_extent = new ProjectedExtent(new Extent(0,0,0,0), CRS.fromName(\"EPSG:3857\"))\n",
    "var grids_noNaN_RDD: RDD[Array[Double]] = sc.emptyRDD\n",
    "var num_cols_rows :(Int, Int) = (0, 0)\n",
    "var cellT :CellType = UByteCellType\n",
    "var mask_tile0 :Tile = new SinglebandGeoTiff(geotrellis.raster.ArrayTile.empty(cellT, num_cols_rows._1, num_cols_rows._2), projected_extent.extent, projected_extent.crs, Tags.empty, GeoTiffOptions.DEFAULT).tile\n",
    "\n",
    "//Load Mask\n",
    "if (toBeMasked) {\n",
    "    val mask_tiles_RDD = sc.hadoopGeoTiffRDD(mask_path).values\n",
    "    val mask_tiles_withIndex = mask_tiles_RDD.zipWithIndex().map{case (e,v) => (v,e)}\n",
    "    mask_tile0 = (mask_tiles_withIndex.filter(m => m._1==0).filter(m => !m._1.isNaN).values.collect())(0)\n",
    "}\n",
    "\n",
    "//Local variables\n",
    "val pattern: String = \"tif\"\n",
    "val filepath: String = dir_path + geoTiff_dir\n",
    "\n",
    "if (rdd_offline_mode) {\n",
    "    grids_noNaN_RDD = sc.objectFile(grids_noNaN_path)\n",
    "\n",
    "    val metadata = sc.sequenceFile(metadata_path, classOf[IntWritable], classOf[BytesWritable]).map(_._2.copyBytes()).collect()\n",
    "    projected_extent = deserialize(metadata(0)).asInstanceOf[ProjectedExtent]\n",
    "    num_cols_rows = (deserialize(metadata(1)).asInstanceOf[Int], deserialize(metadata(2)).asInstanceOf[Int])\n",
    "} else {\n",
    "    if (single_band) {\n",
    "        //Lets load a Singleband GeoTiffs and return RDD just with the tiles.\n",
    "        val geos_RDD = sc.hadoopGeoTiffRDD(filepath, pattern)\n",
    "        val tiles_RDD = geos_RDD.values\n",
    "            \n",
    "        //Retrive the numbre of cols and rows of the Tile's grid\n",
    "        val tiles_withIndex = tiles_RDD.zipWithIndex().map{case (e,v) => (v,e)}\n",
    "        val tile0 = (tiles_withIndex.filter(m => m._1==0).values.collect())(0)\n",
    "        num_cols_rows = (tile0.cols,tile0.rows)\n",
    "        cellT = tile0.cellType\n",
    "    \n",
    "        //Retrieve the ProjectExtent which contains metadata such as CRS and bounding box\n",
    "        val projected_extents_withIndex = geos_RDD.keys.zipWithIndex().map{case (e,v) => (v,e)}\n",
    "        projected_extent = (projected_extents_withIndex.filter(m => m._1 == 0).values.collect())(0)\n",
    "        \n",
    "        if (toBeMasked) {\n",
    "            val mask_tile_broad :Broadcast[Tile] = sc.broadcast(mask_tile0)\n",
    "            grids_noNaN_RDD = tiles_RDD.map(m => m.localInverseMask(mask_tile_broad.value, 1, 0).toArrayDouble().map(m => if ((m == 0.0) || (m.isNaN)) -1.0 else m))\n",
    "        } else {\n",
    "            grids_noNaN_RDD = tiles_RDD.map(m => m.toArrayDouble().map(m => if (m.isNaN) -1.0 else m))\n",
    "        }\n",
    "    } else {\n",
    "        //Lets load Multiband GeoTiffs and return RDD just with the tiles.\n",
    "        val geos_RDD = sc.hadoopMultibandGeoTiffRDD(filepath, pattern)\n",
    "        val tiles_RDD = geos_RDD.values\n",
    "            \n",
    "        //Retrive the numbre of cols and rows of the Tile's grid\n",
    "        val tiles_withIndex = tiles_RDD.zipWithIndex().map{case (e,v) => (v,e)}\n",
    "        val tile0 = (tiles_withIndex.filter(m => m._1==0).values.collect())(0)\n",
    "\n",
    "        num_cols_rows = (tile0.cols,tile0.rows)\n",
    "        cellT = tile0.cellType\n",
    "    \n",
    "        //Retrieve the ProjectExtent which contains metadata such as CRS and bounding box\n",
    "        val projected_extents_withIndex = geos_RDD.keys.zipWithIndex().map{case (e,v) => (v,e)}\n",
    "        projected_extent = (projected_extents_withIndex.filter(m => m._1 == 0).values.collect())(0)\n",
    "    \n",
    "        //Lets read the average of the Spring-Index which is stored in the 4th band\n",
    "        if (toBeMasked) {\n",
    "            val mask_tile_broad :Broadcast[Tile] = sc.broadcast(mask_tile0)\n",
    "            grids_noNaN_RDD = tiles_RDD.map(m => m.band(3).localInverseMask(mask_tile_broad.value, 1, 0).toArrayDouble().map(m => if ((m == 0.0) || (m.isNaN)) -1.0 else m))\n",
    "        } else {\n",
    "            grids_noNaN_RDD = tiles_RDD.map(m => m.band(3).toArrayDouble().map(m => if (m.isNaN) -1.0 else m))\n",
    "        }\n",
    "    }  \n",
    "    //Store data in HDFS\n",
    "    grids_noNaN_RDD.saveAsObjectFile(grids_noNaN_path)\n",
    "    \n",
    "    val writer: SequenceFile.Writer = SequenceFile.createWriter(conf,\n",
    "        Writer.file(metadata_path),\n",
    "        Writer.keyClass(classOf[IntWritable]),\n",
    "        Writer.valueClass(classOf[BytesWritable])\n",
    "    )\n",
    "\n",
    "    writer.append(new IntWritable(1), new BytesWritable(serialize(projected_extent)))\n",
    "    writer.append(new IntWritable(2), new BytesWritable(serialize(num_cols_rows._1)))\n",
    "    writer.append(new IntWritable(3), new BytesWritable(serialize(num_cols_rows._2)))\n",
    "    writer.hflush()\n",
    "    writer.close()\n",
    "}\n",
    "val t1 = System.nanoTime()\n",
    "println(\"Elapsed time: \" + (t1 - t0) + \"ns\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Matrix\n",
    "\n",
    "We need to do a Matrix transpose to have clusters per cell and not per year. With a GeoTiff representing a single year, the loaded data looks liks this:\n",
    "```\n",
    "bands_RDD.map(s => Vectors.dense(s)).cache()\n",
    "\n",
    "//The vectors are rows and therefore the matrix will look like this:\n",
    "[\n",
    "Vectors.dense(0.0, 1.0, 2.0),\n",
    "Vectors.dense(3.0, 4.0, 5.0),\n",
    "Vectors.dense(6.0, 7.0, 8.0),\n",
    "Vectors.dense(9.0, 0.0, 1.0)\n",
    "]\n",
    "```\n",
    "\n",
    "To achieve that we convert the **RDD[Vector]** into a distributed Matrix, a [**CoordinateMatrix**](https://spark.apache.org/docs/latest/mllib-data-types.html#coordinatematrix), which as a **transpose** method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Stage 12:======================================================> (31 + 1) / 32]Elapsed time: 1628018879622ns\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "t0 = 117294627315186\n",
       "grids_matrix = MapPartitionsRDD[39] at map at <console>:86\n",
       "t1 = 118922646194808\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "118922646194808"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val t0 = System.nanoTime()\n",
    "//Global variables\n",
    "var grids_matrix: RDD[Vector] = sc.emptyRDD\n",
    "\n",
    "if (matrix_offline_mode) {\n",
    "    grids_matrix = sc.objectFile(grids_matrix_path)\n",
    "} else {\n",
    "    val mat :RowMatrix = new RowMatrix(grids_noNaN_RDD.map(m => Vectors.dense(m)))\n",
    "\n",
    "    // Split the matrix into one number per line.\n",
    "    val byColumnAndRow = mat.rows.zipWithIndex.map {\n",
    "        case (row, rowIndex) => row.toArray.zipWithIndex.map {\n",
    "            case (number, columnIndex) => new MatrixEntry(rowIndex, columnIndex, number)\n",
    "        }   \n",
    "    }.flatMap(x => x)\n",
    "    \n",
    "    val matt: CoordinateMatrix = new CoordinateMatrix(byColumnAndRow)\n",
    "    val matt_T = matt.transpose()\n",
    "    //grids_matrix = matt_T.toRowMatrix().rows\n",
    "    grids_matrix = matt_T.toIndexedRowMatrix().rows.sortBy(_.index).map(_.vector)\n",
    "    grids_matrix.saveAsObjectFile(grids_matrix_path)\n",
    "}\n",
    "val t1 = System.nanoTime()\n",
    "println(\"Elapsed time: \" + (t1 - t0) + \"ns\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Kmeans\n",
    "\n",
    "We use Kmeans from Sparl-MLlib. The user should only modify the variables on Kmeans setup."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Kmeans Training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3\n",
      "Within Set Sum of Squared Errors = 1.2960778421007756E11                        \n",
      "4\n",
      "Within Set Sum of Squared Errors = 7.512159562280254E10                         \n",
      "5\n",
      "Within Set Sum of Squared Errors = 5.123675278927168E10                         \n",
      "6\n",
      "Within Set Sum of Squared Errors = 3.846708033994191E10                         \n",
      "7\n",
      "Within Set Sum of Squared Errors = 3.105469424051701E10                         \n",
      "8\n",
      "Within Set Sum of Squared Errors = 2.657219149361745E10                         \n",
      "9\n",
      "Within Set Sum of Squared Errors = 2.3102424425223885E10                        \n",
      "10\n",
      "Within Set Sum of Squared Errors = 2.0896273140021423E10                        \n",
      "11\n",
      "Within Set Sum of Squared Errors = 1.935491109936614E10                         \n",
      "12\n",
      "Within Set Sum of Squared Errors = 1.8039011596280605E10                        \n",
      "13\n",
      "Within Set Sum of Squared Errors = 1.7999790559466667E10                        \n",
      "14\n",
      "Within Set Sum of Squared Errors = 1.63046263164963E10                          \n",
      "15\n",
      "Within Set Sum of Squared Errors = 1.574263468598316E10                         \n",
      "16\n",
      "Within Set Sum of Squared Errors = 1.4976494885726418E10                        \n",
      "We will delete the wssse file\n",
      "Lets create it with the new data\n",
      "Elapsed time: 3218492730934ns\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "t0 = 118941683543930\n",
       "kmeans_models = Array(org.apache.spark.mllib.clustering.KMeansModel@5f6f6941, org.apache.spark.mllib.clustering.KMeansModel@73e7d1ee, org.apache.spark.mllib.clustering.KMeansModel@7ef9f3e2, org.apache.spark.mllib.clustering.KMeansModel@5f941c34, org.apache.spark.mllib.clustering.KMeansModel@7b1553e6, org.apache.spark.mllib.clustering.KMeansModel@36b085cb, org.apache.spark.mllib.clustering.KMeansModel@7e7061c, org.apache.spark.mllib.clustering.KMeansModel@786dcf68, org.apache.spark.mllib.clustering.KMeansModel@29b639e0, org.apache.spark.mllib.clustering.KMeansModel@1b9a6fc0, org.apache.spark.mllib.clustering.KMeansModel@24725d49, org.apache.spark.mllib.clustering.KMeansModel@32f2c020, org.apache.spark.mllib....\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[org.apache.spark.mllib.clustering.KMeansModel@5f6f6941, org.apache.spark.mllib.clustering.KMeansModel@73e7d1ee, org.apache.spark.mllib.clustering.KMeansModel@7ef9f3e2, org.apache.spark.mllib.clustering.KMeansModel@5f941c34, org.apache.spark.mllib.clustering.KMeansModel@7b1553e6, org.apache.spark.mllib.clustering.KMeansModel@36b085cb, org.apache.spark.mllib.clustering.KMeansModel@7e7061c, org.apache.spark.mllib.clustering.KMeansModel@786dcf68, org.apache.spark.mllib.clustering.KMeansModel@29b639e0, org.apache.spark.mllib.clustering.KMeansModel@1b9a6fc0, org.apache.spark.mllib.clustering.KMeansModel@24725d49, org.apache.spark.mllib.clustering.KMeansModel@32f2c020, org.apache.spark.mllib.clustering.KMeansModel@2a9396a7, org.apache.spark.mllib.clustering.KMeansModel@62bf085a]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val t0 = System.nanoTime()\n",
    "//Global variables\n",
    "var kmeans_models :Array[KMeansModel] = new Array[KMeansModel](num_kmeans)\n",
    "var wssse_data :List[(Int, Int, Double)] = List.empty\n",
    "\n",
    "if (kmeans_offline_mode) {\n",
    "    var numClusters_id = 0\n",
    "    cfor(minClusters)(_ <= maxClusters, _ + stepClusters) { numClusters =>\n",
    "        if (!fs.exists(new org.apache.hadoop.fs.Path(kmeans_model_paths(numClusters_id)))) {\n",
    "            println(\"One of the files does not exist, we will abort!!!\")\n",
    "            System.exit(0)\n",
    "        } else {\n",
    "            kmeans_models(numClusters_id) = KMeansModel.load(sc, kmeans_model_paths(numClusters_id))\n",
    "        }\n",
    "        numClusters_id += 1\n",
    "    }\n",
    "    val wssse_data_RDD :RDD[(Int, Int, Double)]  = sc.objectFile(wssse_path)\n",
    "    wssse_data  = wssse_data_RDD.collect().toList\n",
    "} else {\n",
    "    var numClusters_id = 0\n",
    "    if (fs.exists(new org.apache.hadoop.fs.Path(wssse_path))) {\n",
    "        val wssse_data_RDD :RDD[(Int, Int, Double)]  = sc.objectFile(wssse_path)\n",
    "        wssse_data  = wssse_data_RDD.collect().toList\n",
    "    }\n",
    "    grids_matrix.cache()\n",
    "    cfor(minClusters)(_ <= maxClusters, _ + stepClusters) { numClusters =>\n",
    "        println(numClusters)\n",
    "        kmeans_models(numClusters_id) = {\n",
    "            KMeans.train(grids_matrix, numClusters, numIterations)\n",
    "        }\n",
    "\n",
    "        // Evaluate clustering by computing Within Set Sum of Squared Errors\n",
    "        val WSSSE = kmeans_models(numClusters_id).computeCost(grids_matrix)\n",
    "        println(\"Within Set Sum of Squared Errors = \" + WSSSE)\n",
    "                \n",
    "        wssse_data = wssse_data :+ (numClusters, numIterations, WSSSE)\n",
    "        \n",
    "        //Save kmeans model\n",
    "        if (save_kmeans_model) {\n",
    "            if (!fs.exists(new org.apache.hadoop.fs.Path(kmeans_model_paths(numClusters_id)))) {\n",
    "                kmeans_models(numClusters_id).save(sc, kmeans_model_paths(numClusters_id))\n",
    "            }\n",
    "        }\n",
    "        numClusters_id += 1\n",
    "    }\n",
    "\n",
    "    //Un-persist it to save memory\n",
    "    grids_matrix.unpersist()\n",
    "    \n",
    "    if (fs.exists(new org.apache.hadoop.fs.Path(wssse_path))) {\n",
    "        println(\"We will delete the wssse file\")\n",
    "        try { fs.delete(new org.apache.hadoop.fs.Path(wssse_path), true) } catch { case _ : Throwable => { } }\n",
    "    }\n",
    "    \n",
    "    println(\"Lets create it with the new data\")\n",
    "    sc.parallelize(wssse_data, 1).saveAsObjectFile(wssse_path)\n",
    "}\n",
    "val t1 = System.nanoTime()\n",
    "println(\"Elapsed time: \" + (t1 - t0) + \"ns\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inspect WSSSE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Name: Compile Error\n",
       "Message: <console>:10: error: eof expected but '}' found.\n",
       "}\n",
       "^\n",
       "\n",
       "StackTrace: "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val t0 = System.nanoTime()\n",
    "//current\n",
    "println(wssse_data)\n",
    "\n",
    "//from disk\n",
    "if (fs.exists(new org.apache.hadoop.fs.Path(wssse_path))) {\n",
    "    var wssse_data_tmp :RDD[(Int, Int, Double)] = sc.objectFile(wssse_path)//.collect()//.toList\n",
    "    println(wssse_data_tmp.collect().toList)    \n",
    "}\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run Kmeans clustering\n",
    "\n",
    "Run Kmeans and obtain the clusters per each cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elapsed time: 142030419ns\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "t0 = 122162075645709\n",
       "kmeans_res = Array(MapPartitionsRDD[1257] at map at KMeansModel.scala:69, MapPartitionsRDD[1258] at map at KMeansModel.scala:69, MapPartitionsRDD[1259] at map at KMeansModel.scala:69, MapPartitionsRDD[1260] at map at KMeansModel.scala:69, MapPartitionsRDD[1261] at map at KMeansModel.scala:69, MapPartitionsRDD[1262] at map at KMeansModel.scala:69, MapPartitionsRDD[1263] at map at KMeansModel.scala:69, MapPartitionsRDD[1264] at map at KMeansModel.scala:69, MapPartitionsRDD[1265] at map at KMeansModel.scala:69, MapPartitionsRDD[1266] at map at KMeansModel.scala:69, MapPartitionsRDD[1267] at map at KMeansModel.scala:69, MapPartitionsRDD[1268] at map at KMeansModel.scala:69, MapPartitionsRDD[1269] at map at KMeansModel.scala:69...\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[MapPartitionsRDD[1257] at map at KMeansModel.scala:69, MapPartitionsRDD[1258] at map at KMeansModel.scala:69, MapPartitionsRDD[1259] at map at KMeansModel.scala:69, MapPartitionsRDD[1260] at map at KMeansModel.scala:69, MapPartitionsRDD[1261] at map at KMeansModel.scala:69, MapPartitionsRDD[1262] at map at KMeansModel.scala:69, MapPartitionsRDD[1263] at map at KMeansModel.scala:69, MapPartitionsRDD[1264] at map at KMeansModel.scala:69, MapPartitionsRDD[1265] at map at KMeansModel.scala:69, MapPartitionsRDD[1266] at map at KMeansModel.scala:69, MapPartitionsRDD[1267] at map at KMeansModel.scala:69, MapPartitionsRDD[1268] at map at KMeansModel.scala:69, MapPartitionsRDD[1269] at map at KMeansModel.scala:69, MapPartitionsRDD[1270] at map at KMeansModel.scala:69]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val t0 = System.nanoTime()\n",
    "//Cache it so kmeans is more efficient\n",
    "grids_matrix.cache()\n",
    "\n",
    "var kmeans_res: Array[RDD[Int]] = Array.fill(num_kmeans)(sc.emptyRDD)\n",
    "var numClusters_id = 0\n",
    "cfor(minClusters)(_ <= maxClusters, _ + stepClusters) { numClusters =>\n",
    "    kmeans_res(numClusters_id) = kmeans_models(numClusters_id).predict(grids_matrix)\n",
    "    numClusters_id += 1\n",
    "}\n",
    "\n",
    "//Un-persist it to save memory\n",
    "grids_matrix.unpersist()\n",
    "val t1 = System.nanoTime()\n",
    "println(\"Elapsed time: \" + (t1 - t0) + \"ns\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Sanity test\n",
    "\n",
    "It can be skipped, it only shows the cluster ID for the first 50 cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111150\n",
      "Elapsed time: 3689349210ns\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "t0 = 122163379270287\n",
       "kmeans_res_out = Array(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)\n",
       "t1 = 122167068619497\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "122167068619497"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val t0 = System.nanoTime()\n",
    "val kmeans_res_out = kmeans_res(0).take(150)\n",
    "kmeans_res_out.foreach(print)\n",
    "\n",
    "println(kmeans_res_out.size)\n",
    "val t1 = System.nanoTime()\n",
    "println(\"Elapsed time: \" + (t1 - t0) + \"ns\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build GeoTiff with Kmeans cluster_IDs\n",
    "\n",
    "The Grid with the cluster IDs is stored in a SingleBand GeoTiff and uploaded to HDFS."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assign cluster ID to each grid cell and save the grid as SingleBand GeoTiff\n",
    "\n",
    "To assign the clusterID to each grid cell it is necessary to get the indices of gird cells they belong to. To join the two RDDS the knowledge was obtaing from a stackoverflow post on [how to perform basic joins of two rdd tables in spark using python](https://stackoverflow.com/questions/31257077/how-do-you-perform-basic-joins-of-two-rdd-tables-in-spark-using-python)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Stage 2258:====================================================> (31 + 1) / 32]DEPRECATED: Use of this script to execute hdfs command is deprecated.\n",
      "Instead use the hdfs command for it.\n",
      "\n",
      "[Stage 2261:====================================================> (31 + 1) / 32]DEPRECATED: Use of this script to execute hdfs command is deprecated.\n",
      "Instead use the hdfs command for it.\n",
      "\n",
      "[Stage 2264:====================================================> (31 + 1) / 32]DEPRECATED: Use of this script to execute hdfs command is deprecated.\n",
      "Instead use the hdfs command for it.\n",
      "\n",
      "[Stage 2267:==================================================>   (30 + 2) / 32]DEPRECATED: Use of this script to execute hdfs command is deprecated.\n",
      "Instead use the hdfs command for it.\n",
      "\n",
      "[Stage 2270:====================================================> (31 + 1) / 32]DEPRECATED: Use of this script to execute hdfs command is deprecated.\n",
      "Instead use the hdfs command for it.\n",
      "\n",
      "[Stage 2273:====================================================> (31 + 1) / 32]DEPRECATED: Use of this script to execute hdfs command is deprecated.\n",
      "Instead use the hdfs command for it.\n",
      "\n",
      "[Stage 2276:====================================================> (31 + 1) / 32]DEPRECATED: Use of this script to execute hdfs command is deprecated.\n",
      "Instead use the hdfs command for it.\n",
      "\n",
      "[Stage 2279:====================================================> (31 + 1) / 32]DEPRECATED: Use of this script to execute hdfs command is deprecated.\n",
      "Instead use the hdfs command for it.\n",
      "\n",
      "[Stage 2282:====================================================> (31 + 1) / 32]DEPRECATED: Use of this script to execute hdfs command is deprecated.\n",
      "Instead use the hdfs command for it.\n",
      "\n",
      "[Stage 2285:====================================================> (31 + 1) / 32]DEPRECATED: Use of this script to execute hdfs command is deprecated.\n",
      "Instead use the hdfs command for it.\n",
      "\n",
      "[Stage 2288:====================================================> (31 + 1) / 32][Stage 2288:=====================================>                (22 + 5) / 32]DEPRECATED: Use of this script to execute hdfs command is deprecated.\n",
      "Instead use the hdfs command for it.\n",
      "\n",
      "[Stage 2291:====================================================> (31 + 1) / 32]DEPRECATED: Use of this script to execute hdfs command is deprecated.\n",
      "Instead use the hdfs command for it.\n",
      "\n",
      "[Stage 2294:====================================================> (31 + 1) / 32]DEPRECATED: Use of this script to execute hdfs command is deprecated.\n",
      "Instead use the hdfs command for it.\n",
      "\n",
      "[Stage 2297:====================================================> (31 + 1) / 32]DEPRECATED: Use of this script to execute hdfs command is deprecated.\n",
      "Instead use the hdfs command for it.\n",
      "\n",
      "Elapsed time: 645597426996ns\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "t0 = 122168532275082\n",
       "numClusters_id = 14\n",
       "t1 = 122814129702078\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "warning: there were two feature warnings; re-run with -feature for details\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "122814129702078"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val t0 = System.nanoTime()\n",
    "var numClusters_id = 0\n",
    "\n",
    "cfor(minClusters)(_ <= maxClusters, _ + stepClusters) { numClusters =>\n",
    "    val grid_clusters_res = kmeans_res(numClusters_id)\n",
    "        \n",
    "    //Define a Tile\n",
    "    val cluster_cells :Array[Double] = grid_clusters_res.map(m => m.toDouble).collect()\n",
    "    val cluster_cellsD = DoubleArrayTile(cluster_cells, num_cols_rows._1, num_cols_rows._2)\n",
    "    val cluster_tile = geotrellis.raster.DoubleArrayTile.empty(num_cols_rows._1, num_cols_rows._2)\n",
    "    cfor(0)(_ < num_cols_rows._1, _ + 1) { col =>\n",
    "        cfor(0)(_ < num_cols_rows._2, _ + 1) { row =>\n",
    "            val v = cluster_cellsD.getDouble(col, row)\n",
    "            if (v == -1.0)\n",
    "                cluster_tile.setDouble(col, row, 0.0)\n",
    "            else\n",
    "                cluster_tile.setDouble(col, row, v)\n",
    "        }\n",
    "    }\n",
    "\n",
    "    val geoTif = new SinglebandGeoTiff(cluster_tile, projected_extent.extent, projected_extent.crs, Tags.empty, GeoTiffOptions(compression.DeflateCompression))\n",
    "    \n",
    "    //Save to /tmp/\n",
    "    GeoTiffWriter.write(geoTif, geotiff_tmp_paths(numClusters_id))\n",
    "\n",
    "    //Upload to HDFS\n",
    "    var cmd = \"hadoop dfs -copyFromLocal -f \" + geotiff_tmp_paths(numClusters_id) + \" \" + geotiff_hdfs_paths(numClusters_id)\n",
    "    Process(cmd)!\n",
    "    \n",
    "    //Remove from /tmp/\n",
    "    cmd = \"rm -fr \" + geotiff_tmp_paths(numClusters_id)\n",
    "    Process(cmd)!\n",
    "    \n",
    "    numClusters_id += 1\n",
    "}\n",
    "val t1 = System.nanoTime()\n",
    "println(\"Elapsed time: \" + (t1 - t0) + \"ns\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [Visualize results](../stable/plot_kmeans_clusters.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Apache Toree - Scala",
   "language": "scala",
   "name": "apache_toree_scala"
  },
  "language_info": {
   "codemirror_mode": "text/x-scala",
   "file_extension": ".scala",
   "mimetype": "text/x-scala",
   "name": "scala",
   "pygments_lexer": "scala",
   "version": "2.11.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
{'fraction_non_alphanumeric': 0.18032545312927606, 'fraction_numerical': 0.036124034271519695, 'mean_word_length': 2.6725164282953227, 'pattern_counts': {'":': 252, '<': 22, '<?xml version=': 0, '>': 67, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 9}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3427018', 'n_tokens_mistral': 18649, 'n_tokens_neox': 17065, 'n_words': 3795}