# 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'

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)

grad_x = cv2.Sobel(gray, ddepth, 1, 0, **kw)

grad_y = cv2.Sobel(gray, ddepth, 0, 1, **kw)

# Converting back to uint8.


In [7]:
# Scharr.

grad_x = cv2.Scharr(gray, ddepth, 1, 0)

grad_y = cv2.Scharr(gray, ddepth, 0, 1)

# Converting back to uint8.


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))


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')