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
import cv2
import math
import numpy as np
from osgeo import gdal
from affine import Affine
from collections import OrderedDict, namedtuple
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'))
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)
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
print("Found %i circles in total" % sum(map(len, detection_results)))
Found 1431 circles in total
# 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)
print("After merging, it ended up with %i circles" % len(circles))
After merging, it ended up with 1103 circles