Searching for aliens.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Irrigation circles detection prototype\n",
    "\n",
    "This script is support material for the blog post: [Searching for aliens](http://www.machinalis.com/blog/searching-for-aliens/)\n",
    "\n",
    "It is based on [Landsat 8 Surface Reflectance High Level Data Products](http://landsat.usgs.gov/CDR_LSR.php). \n",
    "\n",
    "In particular, I'm using Landsat 8 Surface Reflectance data (bands 1 to 7). It can be freely obtained from: http://espa.cr.usgs.gov/index/\n",
    "\n",
    "For the blog post example, I downloaded the scene: [LC82290822016035LGN00](http://earthexplorer.usgs.gov/metadata/8704/LC82290822016035LGN00/) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import cv2\n",
    "import math\n",
    "import numpy as np\n",
    "\n",
    "from osgeo import gdal\n",
    "\n",
    "from affine import Affine\n",
    "from collections import OrderedDict, namedtuple"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "DATA_FILES = [\n",
    "    'LC82290822016035LGN00/LC82290822016035LGN00_sr_band1.tif',\n",
    "    'LC82290822016035LGN00/LC82290822016035LGN00_sr_band2.tif',\n",
    "    'LC82290822016035LGN00/LC82290822016035LGN00_sr_band3.tif',\n",
    "    'LC82290822016035LGN00/LC82290822016035LGN00_sr_band4.tif',\n",
    "    'LC82290822016035LGN00/LC82290822016035LGN00_sr_band5.tif',\n",
    "    'LC82290822016035LGN00/LC82290822016035LGN00_sr_band6.tif',\n",
    "    'LC82290822016035LGN00/LC82290822016035LGN00_sr_band7.tif',\n",
    "]\n",
    "\n",
    "\n",
    "# Minimum possible radius to look for (in pixels)\n",
    "MIN_RADIUS = 14\n",
    "\n",
    "# Maximum possible radius to look for (in pixels)\n",
    "MAX_RADIUS = 23\n",
    "\n",
    "# Two circles closer that this (in pixels) will be treated as one\n",
    "MIN_DIST = 14\n",
    "\n",
    "# Configuration for the circles detection algorithm (cv2.HoughCircles)\n",
    "DETECTOR_CONF = OrderedDict()\n",
    "DETECTOR_CONF['method'] = cv2.HOUGH_GRADIENT\n",
    "DETECTOR_CONF['dp'] = 0.5\n",
    "DETECTOR_CONF['minDist'] = 23\n",
    "DETECTOR_CONF['param1'] = 45\n",
    "DETECTOR_CONF['param2'] = 20\n",
    "DETECTOR_CONF['minRadius'] = MIN_RADIUS\n",
    "DETECTOR_CONF['maxRadius'] = MAX_RADIUS\n",
    "\n",
    "Circle = namedtuple('Circle', ('x', 'y', 'radius', 'lon', 'lat', 'source'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def read_geotiff(data_file_path):\n",
    "    \"\"\"\n",
    "    Utility to read the first band of a Geotiff image.\n",
    "\n",
    "    :param data_file_path: Path to a Geotiff file.\n",
    "    :return: (Numpy Array, geo projection, geo transform)\n",
    "\n",
    "    \"\"\"\n",
    "    raster_dataset = gdal.Open(data_file_path, gdal.GA_ReadOnly)\n",
    "    band = raster_dataset.GetRasterBand(1)\n",
    "    band_data = band.ReadAsArray()\n",
    "    proj = raster_dataset.GetProjection()\n",
    "    geo = raster_dataset.GetGeoTransform()\n",
    "    raster_dataset = None\n",
    "    return band_data, proj, geo\n",
    "\n",
    "\n",
    "def cast_l8sr_data_to_byte(image_data):\n",
    "    \"\"\"\n",
    "    Transform Landsat 8 Surface Reflectance data to 8-bit.\n",
    "    \n",
    "    WARNING: Modifies the input data!\n",
    "\n",
    "    :param image_data: Numpy Array\n",
    "    :param fname_suffix: String with the data source filename suffix.\n",
    "    :return: Numpy Array of np.ubyte dtype\n",
    "\n",
    "    \"\"\"\n",
    "    # Reflectance data. INT16. Valid range = (-2000, 16000). Fill=-9999. Saturate=20000\n",
    "    min_value, max_value = (-2000, 16000)\n",
    "    fill_value = -9999\n",
    "\n",
    "    # Loose the no-data values.\n",
    "    image_data[image_data == fill_value] = min_value\n",
    "    offset = abs(min_value)\n",
    "    image_data = ((image_data + offset) / (max_value + offset)) * 255.0\n",
    "    return image_data.round().astype(np.ubyte)\n",
    "\n",
    "\n",
    "\n",
    "def detect_circles(data_file_path):\n",
    "    \"\"\"\n",
    "    Detect circles in the given image, using cv2.HoughCircles.\n",
    "    :param data_file_path: str. Path to Geotiff file.\n",
    "    :return: dict with 'fname' and 'circles' keys.\n",
    "\n",
    "    \"\"\"\n",
    "    image, _, geo_transform = read_geotiff(data_file_path)\n",
    "    ubyte_image = cast_l8sr_data_to_byte(image)\n",
    "\n",
    "    detection_results = cv2.HoughCircles(ubyte_image, **DETECTOR_CONF)\n",
    "\n",
    "    circles = []\n",
    "    if detection_results.shape[0]:\n",
    "        coordinates_transformation = Affine.from_gdal(*geo_transform)\n",
    "        for (col, row, radius) in np.round(detection_results[0, :, :3]).astype(\"int\"):\n",
    "            lon, lat = coordinates_transformation * (col, row)\n",
    "            circles.append(\n",
    "                Circle(col, row, radius, lon, lat, source=data_file_path)\n",
    "            )\n",
    "    return circles\n",
    "\n",
    "\n",
    "def merge_circles(current_circles, candidate_circles):\n",
    "    \"\"\"\n",
    "    Update (and return) the current_circles list with the elements in detection_results that are not\n",
    "    too close to already existing circles.\n",
    "\n",
    "    :param current_circles: list of current circles. Might be updated.\n",
    "    :param detection_results: dict, with 'circles' and 'fname' keys (output from <_detect_circles>)\n",
    "    :return: list of filtered circles.\n",
    "\n",
    "    \"\"\"\n",
    "    filtered_circles_data = []\n",
    "    for new_point in candidate_circles:\n",
    "        is_new = True\n",
    "        for point in current_circles:\n",
    "            if distance(new_point, point) < MIN_DIST:\n",
    "                is_new = False\n",
    "                break\n",
    "        if is_new:\n",
    "            filtered_circles_data.append(new_point)\n",
    "            \n",
    "    current_circles = current_circles + filtered_circles_data\n",
    "    return current_circles\n",
    "\n",
    "\n",
    "def distance(point_a, point_b):\n",
    "    \"\"\"\n",
    "    Euclidean distance between point_a and point_b.\n",
    "\n",
    "    :param point_a: List or tuple with two elements, (x, y).\n",
    "    :param point_b: List or tuple with two elements, (x, y).\n",
    "    :return: float\n",
    "\n",
    "    \"\"\"\n",
    "    return math.sqrt((point_a[0] - point_b[0])**2 + (point_a[1] - point_b[1])**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Detect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "46 circles detected in LC82290822016035LGN00/LC82290822016035LGN00_sr_band1.tif\n",
      "68 circles detected in LC82290822016035LGN00/LC82290822016035LGN00_sr_band2.tif\n",
      "88 circles detected in LC82290822016035LGN00/LC82290822016035LGN00_sr_band3.tif\n",
      "140 circles detected in LC82290822016035LGN00/LC82290822016035LGN00_sr_band4.tif\n",
      "469 circles detected in LC82290822016035LGN00/LC82290822016035LGN00_sr_band5.tif\n",
      "327 circles detected in LC82290822016035LGN00/LC82290822016035LGN00_sr_band6.tif\n",
      "293 circles detected in LC82290822016035LGN00/LC82290822016035LGN00_sr_band7.tif\n"
     ]
    }
   ],
   "source": [
    "detection_results = []\n",
    "for band_fname in DATA_FILES:\n",
    "    circles_in_band = detect_circles(band_fname)\n",
    "    print(len(circles_in_band), 'circles detected in', band_fname)\n",
    "    detection_results.append(circles_in_band)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found 1431 circles in total\n"
     ]
    }
   ],
   "source": [
    "print(\"Found %i circles in total\" % sum(map(len, detection_results)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Merge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Merge circles whose centers are too close and transform image (x, y) pixel location to map\n",
    "# (lon, lat) geographic coordinates.\n",
    "circles = []\n",
    "for candidate_circles in detection_results:\n",
    "    circles = merge_circles(circles, candidate_circles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "After merging, it ended up with 1103 circles\n"
     ]
    }
   ],
   "source": [
    "print(\"After merging, it ended up with %i circles\" % len(circles))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1+"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}