python4oceanographers

Turning ripples into waves

Finding oceanic fronts with edge detection techniques

After a friend pointed me to this nice blog post about images derivatives with openCV I had to try it with an Sea Surface Temperature (SST) to see if we can obtain something useful.

First let's download some data. In the first cell I use iris to load the SST data from the SABGOM model results for 2014-8-22. Next I plot the data using matplotlib imshow.

In [3]:
import warnings
from datetime import datetime

import iris
iris.FUTURE.cell_datetime_objects = True

url = 'http://omgsrv1.meas.ncsu.edu:8080/thredds/dodsC/fmrc/sabgom/'
url += 'SABGOM_Forecast_Model_Run_Collection_best.ncd'

cube = iris.load_cube(url, 'sea_water_potential_temperature')

date = datetime(2014, 8, 22)
constraint = iris.Constraint(time=lambda cell: cell == datetime(2014, 8, 15, 12, 0, 0))
cube = cube.extract(constraint)
temp = cube.data[-1, ...]  # Surface.
In [4]:
%matplotlib inline
import matplotlib.pyplot as plt

figsize = (7, 7)
cbarkw = dict(shrink=0.6, extend='both')

fig, ax = plt.subplots(figsize=figsize)
i = ax.imshow(temp, origin='lower')
cbar = fig.colorbar(i, **cbarkw)
_ = ax.axis('off')

matplotlib's imshow has a neat trick to return the image data: the method _rgbacache. Note also that openCV expects the image in Blue Green Red Alpha instead of the usual Red Green Blue Alpha.

In the next cell I just convert the image to openCV BGRA, use a Gaussian blur to smooth things out and convert the image to gray.

In [5]:
import cv2
import numpy as np

def matplotlib_to_opencv(i):
    image = i._rgbacache
    r, g, b, a = cv2.split(image)
    return np.flipud(cv2.merge([b, g, r, a]))

image = matplotlib_to_opencv(i)
img = cv2.GaussianBlur(image, (3, 3), 0)
gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)

fig, ax = plt.subplots(figsize=figsize)
ax.imshow(gray, cmap=plt.cm.gray)
_ = ax.axis('off')

Like in the original post I will try Sobel and Scharr methods by using the x and y derivatives and displaying them by blending with equal weights and by just adding them up.

In [6]:
# Sobel.

ddepth = cv2.CV_16S
kw = dict(ksize=3, scale=1, delta=0, borderType=cv2.BORDER_DEFAULT)

# Gradient-X.
grad_x = cv2.Sobel(gray, ddepth, 1, 0, **kw)

# Gradient-Y.
grad_y = cv2.Sobel(gray, ddepth, 0, 1, **kw)

# Converting back to uint8.
abs_grad_x = cv2.convertScaleAbs(grad_x)
abs_grad_y = cv2.convertScaleAbs(grad_y)

sobel = cv2.addWeighted(abs_grad_x, 0.5, abs_grad_y, 0.5, 0)
sobel_no_blend = cv2.add(abs_grad_x, abs_grad_y)
In [7]:
# Scharr.

# Gradient-X
grad_x = cv2.Scharr(gray, ddepth, 1, 0)

# Gradient-Y.
grad_y = cv2.Scharr(gray, ddepth, 0, 1)

# Converting back to uint8.
abs_grad_x = cv2.convertScaleAbs(grad_x)
abs_grad_y = cv2.convertScaleAbs(grad_y)

scharr = cv2.addWeighted(abs_grad_x, 0.5, abs_grad_y, 0.5, 0)
scharr_no_blend = cv2.add(abs_grad_x, abs_grad_y)
In [8]:
def turn_off_labels(ax):
    ax.axes.get_xaxis().set_visible(False)
    ax.axes.get_yaxis().set_visible(False)


kw = dict(cmap=plt.cm.Greys_r)

fig, ((ax0, ax1), (ax2, ax3)) = plt.subplots(ncols=2, nrows=2, figsize=(9, 8))
ax0.imshow(sobel, **kw)
ax0.set_title('Sobel')

ax1.imshow(scharr, **kw)
ax1.set_title('Scharr')

ax2.imshow(sobel_no_blend, **kw)
ax2.set_title('Sobel no weights')

ax3.imshow(scharr_no_blend, **kw)
ax3.set_title('Scharr no weights')

map(turn_off_labels, (ax0, ax1, ax2, ax3))
fig.tight_layout(h_pad=0.1, w_pad=0.5)

I liked the results from Scharr with the weights. The results seem sharper then Sobel and the Gulf stream shows up nicely.

The original blog post also tries a Laplacian operator and a Canny edge detector.

In [9]:
# Laplacian.

ddepth = cv2.CV_16S
kw = dict(ksize=3, scale=1, delta=0, borderType=cv2.BORDER_DEFAULT)

gray_lap = cv2.Laplacian(gray, ddepth, **kw)
laplacian = cv2.convertScaleAbs(gray_lap)

fig, ax = plt.subplots(figsize=figsize)
ax.imshow(laplacian, cmap=plt.cm.Greys_r)
_ = ax.axis('off')