{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Nowcasting con imágenes satelitales"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vamos a presentar una posible aplicación de nowcasting usando imágenes satelitales. En este caso, buscaremos detectar asentamientos precarios usando solamente esa fuente de datos. De esta manera, podemos actualizar periódicamente la cantindad de km2 de asentamientos y así tener un indicador que pueda anticiparse a las estadísticas de hábitat.\n",
    "\n",
    "El objetivo será mostrar como elaborar el mapa para que luego pueda ser replicable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tipos de imágenes satelitales"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Existen distintos tipos de imágenes satelitales (ópticas), una clasificación que se puede hacer es de acuerdo a su resolución.\n",
    "\n",
    "Podemos distinguir entre (no exhaustivo):\n",
    "\n",
    "* Resolución media: [MODIS](https://modis.gsfc.nasa.gov/), 250-1000m por pixel.\n",
    "    * Utilizadas para investigación en cambio climático, recursos naturales y análisis de grandes superficies.\n",
    "    * 1 imagen nueva por día.\n",
    "    * Gratuitas, disponibles en varias fuentes. Por ejemplo [Google Earth Engine](https://earthengine.google.com/)\n",
    "    * Ejemplos de nowcasting: predicción de cosecha para el agro para grandes supericies. Indices de calidad del aire.\n",
    "\n",
    "\n",
    "* Resolución alta: Landsat (30m), Sentinel-2 (10m), RapidEye (5m).\n",
    "    * Utilizadas para análisis urbano y agro.\n",
    "    * 1 imagen cada 15 días (Landsat) y 5 días (Sentinel-2).\n",
    "    * Landsat y Sentinel-2 son gratuitas, disponibles en varias fuentes.\n",
    "    * Ejemplos de nowcasting: predicción de cosecha e indicadores relacionados con agricultura de precisión.\n",
    "\n",
    "\n",
    "* Resolución muy alta: WorldView-4 (0.3m de resolución)\n",
    "    * utilizadas para detección de objetos y monitoreo de gran precisión.\n",
    "    * Revisita variable según la zona.\n",
    "    * No son gratuitas (~ 24 USD/km2).\n",
    "    * Ejemplos de nowcasting: conteo de autos en parking lots de supermercados y contenedores en puertos."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ejemplo Museo Louvre (Paris, Francia)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Sentinel-2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![louvre_sen](img/louvre_sen.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### WorldView-4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![louvre_wv](img/louvre_wv.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Enfoques para resolución de problemas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Teledetección tradicional"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combinaciones de bandas en las imágenes para resaltar aspectos relevantes: por ejemplo banda roja e infrarroja para resaltar vegetación."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![caba_nc](img/caba_nc.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![caba_ndvi](img/caba_ndvi.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Machine Learning a nivel pixel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Se clasifican las imágenes de acuerdo a los valores de los píxeles. Esto puede traer problemas por la variabilidad que puede haber (ruido en la captura de la imagen y sombra entre otros).\n",
    "\n",
    "Se resuelven problemas de clasificación de usos del suelo, detección de cambios, detección de cursos de agua y detección de vegetación, entre otros.\n",
    "\n",
    "Ejemplo de detección de cambios para el municipio de Pilar."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![pilar_after](img/pilar_after.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![pilar_before](img/pilar_before.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![pilar_cd](img/pilar_cd.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### OBIA (Object Based Image Analysis)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Se segmenta (agrupamiento de pixeles) previamente la imagen. De esta manera, cada fila es un segmento y a partir de estos calculamos atributos. Por ejemplo, para cada segmento se puede calcular distancia a cursos de agua, distancia a avenidas, intensidad de brillo promedio, etc.\n",
    "\n",
    "Se clasifican las imágenes de acuerdo a los valores de los segmentos."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![lma_raw](img/lma_raw.png) ![lma_seg](img/lma_seg.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![lma_seg_class](img/lma_seg_class.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Redes Neuronales y Deep Learning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Un enfoque posible es usando redes neuronales convolucionales para hacer la clasificación. Se combina con técnicas de computer vision para filtrar y segmentar las imágenes de acuerdo a la clase que se busca detectar.\n",
    "\n",
    "A continuación va un ejemplo de cómo procesar imáges de un territorio para clasificar asentamientos precarios usando redes neuronales convolucionales. Vamos a hacer transfer learning, aprovechando de redes ya entrenadas."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Se particiona la imagen entera en tiles de 256x256 pixeles y se cruza con las localizaciones de asentamientos precarios para establecer cuales son True y cuales False."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ejemplos de True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![guatemala_vya_t](img/guatemala_vya_t.jpg) ![guatemala_vya_t_2](img/guatemala_vya_t_2.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ejemplos de False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![guatemala_no_vya_f](img/guatemala_no_vya_f.jpg) ![guatemala_no_vya_f_2](img/guatemala_no_vya_f_2.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dada la estructura de los directorios,\n",
    "\n",
    "```\n",
    "dir/\n",
    "  train/\n",
    "    true/\n",
    "    false/\n",
    "  validation/\n",
    "    true/\n",
    "    false/\n",
    "```\n",
    "\n",
    "levantamos los datos con una estructura que permita trabajar con Keras."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_files = glob.glob(os.path.join('data/','train/**/', '*.jpg'))\n",
    "validation_files = glob.glob(os.path.join('data/', 'validation/**/', '*.jpg'))\n",
    "img_width, img_height = 256, 256\n",
    "train_data_dir = \"data/train/\"\n",
    "validation_data_dir = \"data/validation/\"\n",
    "nb_train_samples = len(train_files)\n",
    "nb_validation_samples = len(validation_files)\n",
    "batch_size = 64\n",
    "epochs = 200"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Inicializamos una red ResNet50"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = applications.resnet50.ResNet50(\n",
    "    weights = \"imagenet\",\n",
    "    include_top = False,\n",
    "    input_shape = (img_width, img_height, 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Freezamos las primeras 50 capas, vamos a reentrenar las restantes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for layer in model.layers[:50]:\n",
    "    layer.trainable = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Agregamos una última capa al final para obtener las predicciones"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = model.output\n",
    "x = Flatten()(x)\n",
    "x = Dense(1024, activation = \"relu\")(x)\n",
    "x = Dropout(0.5)(x)\n",
    "predictions = Dense(1, activation = \"sigmoid\")(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Cerramos el modelo y lo compilamos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_final = Model(inputs = model.input, outputs = predictions)\n",
    "\n",
    "# compile the model\n",
    "model_final.compile(\n",
    "    loss = 'binary_crossentropy',\n",
    "    optimizer = optimizers.SGD(lr = 0.0001, momentum = 0.9),\n",
    "    metrics = ['accuracy'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Generamos los conjuntos de datos de entrenamiento y validación\n",
    "\n",
    "Usamos data agumentation para contemplar rotaciones y otras variaciones. El concepto aquí es que no nos importa en qué parte de la imagen está la clase que queremos detectar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_datagen = ImageDataGenerator(\n",
    "    rescale = 1./255,\n",
    "    horizontal_flip = True,\n",
    "    fill_mode = \"nearest\",\n",
    "    zoom_range = 0.3,\n",
    "    width_shift_range = 0.3,\n",
    "    height_shift_range = 0.3,\n",
    "    rotation_range = 30)\n",
    "\n",
    "test_datagen = ImageDataGenerator(\n",
    "    rescale = 1./255,\n",
    "    horizontal_flip = True,\n",
    "    fill_mode = \"nearest\",\n",
    "    zoom_range = 0.3,\n",
    "    width_shift_range = 0.3,\n",
    "    height_shift_range = 0.3,\n",
    "    rotation_range = 30)\n",
    "\n",
    "train_generator = train_datagen.flow_from_directory(\n",
    "    train_data_dir,\n",
    "    target_size = (img_height, img_width),\n",
    "    batch_size = batch_size,\n",
    "    class_mode = \"binary\")\n",
    "\n",
    "validation_generator = test_datagen.flow_from_directory(\n",
    "    validation_data_dir,\n",
    "    target_size = (img_height, img_width),\n",
    "    class_mode = \"binary\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finalmente, usamos *early stopping* y usamos checkpoints para ir guardando cada corrida. Al final guardamos los pesos de la red."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "checkpoint = ModelCheckpoint(\"./models/demo_weights.hdf5\".format(trainable_layers),\n",
    "    monitor = 'val_acc',\n",
    "    verbose = 1,\n",
    "    save_best_only = True,\n",
    "    save_weights_only = False,\n",
    "    mode = 'auto',\n",
    "    period = 1)\n",
    "\n",
    "early = EarlyStopping(\n",
    "    monitor = 'val_acc',\n",
    "    min_delta = 0,\n",
    "    patience = 10,\n",
    "    verbose = 1,\n",
    "    mode = 'auto')\n",
    "\n",
    "model_final.fit_generator(\n",
    "    train_generator,\n",
    "    steps_per_epoch = nb_train_samples // batch_size,\n",
    "    epochs = epochs, \n",
    "    validation_data = validation_generator,\n",
    "    validation_steps = nb_validation_samples // batch_size,\n",
    "    callbacks = [checkpoint, early])\n",
    "\n",
    "\n",
    "model_final.save('./models/demo.h5'.format(trainable_layers))\n",
    "print('Done with {} layers'.format(trainable_layers))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Luego de esto hay que proceder a filtrar los resultados y segmentar las imágenes para terminar de detectar los objetos relevantes. Cabe destacar que hay arquitecturas de redes que permiten hacer esto último, como el caso de [U-Net](https://blog.deepsense.ai/deep-learning-for-satellite-imagery-via-image-segmentation/).\n",
    "\n",
    "Para segmentar las imágenes se puede usar la librería scikit-image u OpenCV."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Resultado ([link al mapa](http://dymaxionlabs.com/ap-latam/en/map/?id=guatemala))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Guatemala](img/guatemala_map.jpg)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}