Irrigation circles detection prototype

This script is support material for the blog post: Searching for aliens

It is based on Landsat 8 Surface Reflectance High Level Data Products.

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/

For the blog post example, I downloaded the scene: LC82290822016035LGN00

In [36]:
import cv2
import math
import numpy as np

from osgeo import gdal

from affine import Affine
from collections import OrderedDict, namedtuple
In [37]:
DATA_FILES = [
    'LC82290822016035LGN00/LC82290822016035LGN00_sr_band1.tif',
    'LC82290822016035LGN00/LC82290822016035LGN00_sr_band2.tif',
    'LC82290822016035LGN00/LC82290822016035LGN00_sr_band3.tif',
    'LC82290822016035LGN00/LC82290822016035LGN00_sr_band4.tif',
    'LC82290822016035LGN00/LC82290822016035LGN00_sr_band5.tif',
    'LC82290822016035LGN00/LC82290822016035LGN00_sr_band6.tif',
    'LC82290822016035LGN00/LC82290822016035LGN00_sr_band7.tif',
]


# Minimum possible radius to look for (in pixels)
MIN_RADIUS = 14

# Maximum possible radius to look for (in pixels)
MAX_RADIUS = 23

# Two circles closer that this (in pixels) will be treated as one
MIN_DIST = 14

# Configuration for the circles detection algorithm (cv2.HoughCircles)
DETECTOR_CONF = OrderedDict()
DETECTOR_CONF['method'] = cv2.HOUGH_GRADIENT
DETECTOR_CONF['dp'] = 0.5
DETECTOR_CONF['minDist'] = 23
DETECTOR_CONF['param1'] = 45
DETECTOR_CONF['param2'] = 20
DETECTOR_CONF['minRadius'] = MIN_RADIUS
DETECTOR_CONF['maxRadius'] = MAX_RADIUS

Circle = namedtuple('Circle', ('x', 'y', 'radius', 'lon', 'lat', 'source'))
In [38]:
def read_geotiff(data_file_path):
    """
    Utility to read the first band of a Geotiff image.

    :param data_file_path: Path to a Geotiff file.
    :return: (Numpy Array, geo projection, geo transform)

    """
    raster_dataset = gdal.Open(data_file_path, gdal.GA_ReadOnly)
    band = raster_dataset.GetRasterBand(1)
    band_data = band.ReadAsArray()
    proj = raster_dataset.GetProjection()
    geo = raster_dataset.GetGeoTransform()
    raster_dataset = None
    return band_data, proj, geo


def cast_l8sr_data_to_byte(image_data):
    """
    Transform Landsat 8 Surface Reflectance data to 8-bit.
    
    WARNING: Modifies the input data!

    :param image_data: Numpy Array
    :param fname_suffix: String with the data source filename suffix.
    :return: Numpy Array of np.ubyte dtype

    """
    # Reflectance data. INT16. Valid range = (-2000, 16000). Fill=-9999. Saturate=20000
    min_value, max_value = (-2000, 16000)
    fill_value = -9999

    # Loose the no-data values.
    image_data[image_data == fill_value] = min_value
    offset = abs(min_value)
    image_data = ((image_data + offset) / (max_value + offset)) * 255.0
    return image_data.round().astype(np.ubyte)



def detect_circles(data_file_path):
    """
    Detect circles in the given image, using cv2.HoughCircles.
    :param data_file_path: str. Path to Geotiff file.
    :return: dict with 'fname' and 'circles' keys.

    """
    image, _, geo_transform = read_geotiff(data_file_path)
    ubyte_image = cast_l8sr_data_to_byte(image)

    detection_results = cv2.HoughCircles(ubyte_image, **DETECTOR_CONF)

    circles = []
    if detection_results.shape[0]:
        coordinates_transformation = Affine.from_gdal(*geo_transform)
        for (col, row, radius) in np.round(detection_results[0, :, :3]).astype("int"):
            lon, lat = coordinates_transformation * (col, row)
            circles.append(
                Circle(col, row, radius, lon, lat, source=data_file_path)
            )
    return circles


def merge_circles(current_circles, candidate_circles):
    """
    Update (and return) the current_circles list with the elements in detection_results that are not
    too close to already existing circles.

    :param current_circles: list of current circles. Might be updated.
    :param detection_results: dict, with 'circles' and 'fname' keys (output from <_detect_circles>)
    :return: list of filtered circles.

    """
    filtered_circles_data = []
    for new_point in candidate_circles:
        is_new = True
        for point in current_circles:
            if distance(new_point, point) < MIN_DIST:
                is_new = False
                break
        if is_new:
            filtered_circles_data.append(new_point)
            
    current_circles = current_circles + filtered_circles_data
    return current_circles


def distance(point_a, point_b):
    """
    Euclidean distance between point_a and point_b.

    :param point_a: List or tuple with two elements, (x, y).
    :param point_b: List or tuple with two elements, (x, y).
    :return: float

    """
    return math.sqrt((point_a[0] - point_b[0])**2 + (point_a[1] - point_b[1])**2)

Detect

In [39]:
detection_results = []
for band_fname in DATA_FILES:
    circles_in_band = detect_circles(band_fname)
    print(len(circles_in_band), 'circles detected in', band_fname)
    detection_results.append(circles_in_band)
46 circles detected in LC82290822016035LGN00/LC82290822016035LGN00_sr_band1.tif
68 circles detected in LC82290822016035LGN00/LC82290822016035LGN00_sr_band2.tif
88 circles detected in LC82290822016035LGN00/LC82290822016035LGN00_sr_band3.tif
140 circles detected in LC82290822016035LGN00/LC82290822016035LGN00_sr_band4.tif
469 circles detected in LC82290822016035LGN00/LC82290822016035LGN00_sr_band5.tif
327 circles detected in LC82290822016035LGN00/LC82290822016035LGN00_sr_band6.tif
293 circles detected in LC82290822016035LGN00/LC82290822016035LGN00_sr_band7.tif
In [40]:
print("Found %i circles in total" % sum(map(len, detection_results)))
Found 1431 circles in total

Merge

In [41]:
# Merge circles whose centers are too close and transform image (x, y) pixel location to map
# (lon, lat) geographic coordinates.
circles = []
for candidate_circles in detection_results:
    circles = merge_circles(circles, candidate_circles)
In [42]:
print("After merging, it ended up with %i circles" % len(circles))
After merging, it ended up with 1103 circles
In [ ]: