Counting features in satellite images using scikit-image

The example below uses scikit-image library to detect circular features in farms using center pivot irrigation in Saudi Arabia. It then counts and reports the number of farms. This is one of the ways in which libraries from the scientific Python ecosystem can be integrated with the ArcGIS platform.

It uses the Multispectral Landsat imagery available at ArcGIS Online.

Note: to run this sample, you need a few extra libraries in your conda environment. If you don't have the libraries, install them by running the following commands from cmd.exe or your shell

conda install scipy
conda install matplotlib
conda install scikit-image
In [1]:
from arcgis.gis import GIS
agol = GIS()
In [2]:
l8 = agol.content.search('"Multispectral Landsat"', 'Imagery Layer')[0]
l8
Out[2]:
Multispectral Landsat
Landsat 8 OLI, 30m multispectral and multitemporal 8-band imagery, with on-the-fly renderings and indices. This imagery layer is sourced from the Landsat on AWS collections and is updated daily with new imagery.Imagery Layer by esri
Last Modified: May 04, 2018
0 comments, 1,19,600 views
In [3]:
l8lyr = l8.layers[0]

The code below sets the extent of the layer, to a known extent of farms in Saudi Arabia and then visualizes the landsat layer:

In [4]:
l8lyr.extent = {'spatialReference': {'latestWkid': 3857, 'wkid': 102100},
                                         'type': 'extent',
                                         'xmax': 4296559.143733407,
                                         'xmin': 4219969.241391764,
                                         'ymax': 3522726.823081019,
                                         'ymin': 3492152.0117669892}
l8lyr
Out[4]:

We can preprocess the imagery using raster functions. The code below uses the ndvi raster function to identify areas that have healthy vegetation. This preprocessing step makes the scikit-image blob detection algorithm work better.

In [5]:
from arcgis.raster.functions import ndvi, stretch
In [6]:
stretch(ndvi(l8lyr), stretch_type='PercentClip', min_percent=30, max_percent=70, dra=True)
Out[6]:

The code below exports the imagery to a file from which we read it using matplotlib's image API and plot it:

In [7]:
img = stretch(ndvi(l8lyr), stretch_type='PercentClip', min_percent=30, max_percent=70, dra=True).export_image(bbox=l8lyr.extent, bbox_sr=102100, size=[1200, 450],
                       export_format='jpeg', save_folder='.', save_file='centerpivotfarms.jpg', f='image')
In [9]:
import matplotlib.pyplot as plt
import matplotlib.image as mpimg

img = mpimg.imread('centerpivotfarms.jpg')

# what does it look like?
plt.imshow(img)
plt.show()

Blob detection using scikit-image

The code below uses scikit-image library to find blobs in the given grayscale image, and reports the number of farms thus detected. It also plots them for visualization using matplotlib.

Blobs are found using the Difference of Gaussian (DoG) method.

In [10]:
from skimage import feature, color
import matplotlib.pyplot as plt
bw = img.mean(axis=2)

fig = plt.figure(figsize = (15,15))
ax = fig.add_subplot(1,1,1)

blobs_dog = [(x[0],x[1],x[2]) for x in feature.blob_dog(-bw, 
                                                        min_sigma=4, 
                                                        max_sigma=8,
                                                        threshold=0.1,
                                                        overlap=0.6)]

#remove duplicates
blobs_dog = set(blobs_dog)

img_blobs = color.gray2rgb(img)

for blob in blobs_dog:
    y, x, r = blob
    c = plt.Circle((x, y), r+1, color='red', linewidth=2, fill=False)
    ax.add_patch(c)

plt.imshow(img_blobs)
plt.title('Center Pivot Farms')

plt.show()
print('Number of center pivot farms detected: ' + str(len(blobs_dog)))
Number of center pivot farms detected: 964

Feedback on this topic?