{ "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 }