Based on a solution that I read at How to define the markers for Watershed in OpenCV?, I am trying apply watershed to grayscale data (not very visible but not all black), extracted from netcdf (precipitation data).
Here is a black and white version of the data (threshold at 0) so that you can see more easily, and the markers I want to use to define the different basins (basically just another threshold where precipitation is more intense).
The code I'm running is as follows:
import os,sys,string
from netCDF4 import Dataset as nc
import cv2
import numpy as np
import matplotlib.pyplot as mpl
import scipy.ndimage as ndimage
import scipy.spatial as spatial
from skimage import filter
from skimage.morphology import watershed
from scipy import ndimage
filename=["Cmorph-1999_01_03.nc"]
nc_data=nc(filename[0])
data=nc_data.variables["CMORPH"][23,0:250,250:750]
new_data=np.flipud(data)
ma_data=np.ma.masked_where(new_data<=0,new_data)
ma_conv=np.ma.masked_where(new_data<=2,new_data)
## Borders
tmp_data=ma_data.filled(0)
tmp_data[np.where(tmp_data!=0)]=255
bw_data=tmp_data.astype(np.uint8)
border = cv2.dilate(bw_data, None, iterations=5)
border = border - cv2.erode(border, None)
## Markers
tmp_conv=ma_conv.filled(0)
tmp_conv[np.where(tmp_conv!=0)]=255
bw_conv=tmp_conv.astype(np.uint8)
lbl, ncc = ndimage.label(bw_conv)
lbl = lbl * (255/ncc)
lbl[border == 255] = 255
lbl = lbl.astype(np.int32)
## Apply watershed
cv2.watershed(ma_data, lbl)
lbl[lbl == -1] = 0
lbl = lbl.astype(np.uint8)
result = 255 - lbl
I have the following error for the watershed in opencv-2.4.11/modules/imgproc/src/segmentation.cpp:
error: (-210) Only 8-bit, 3-channel input images are supported in function cvWatershed
For what I saw on the internet, this is due to the fact that the grayscale data is a 2D image and watershed needs a 3D image (from RGB). Indeed, I tried the script with a jpg image and I worked perfectly. This problem is mentionned here but the answer given was finally rejected. And I can't find any more recent link answering the question.
To try to solve this, I created a 3D array from the 2D new_data:
new_data = new_data[..., np.newaxis]
test=np.append(new_data, new_data, axis=2)
test=np.append(new_data, test, axis=2)
But, as expected, it didn't solve the problem (same error message).
I also tried to save the plot from matplotlib to get RGB data:
fig = mpl.figure()
fig.add_subplot(111)
fig.tight_layout(pad=0)
mpl.contourf(ma_data,levels=np.arange(0,255.1,0.1))
fig.canvas.draw()
test_data = np.fromstring(fig.canvas.tostring_rgb(), dtype=np.uint8, sep='')
test_data = test_data.reshape(fig.canvas.get_width_height()[::-1] + (3,))
But the size of the test_data created is different from ma_data (+ I can't get rid of the labels).
So, I am stuck here. Ideally, I want to apply the watershed on the 2D grayscale image directly and/or limit the number of operations as much as possible.
As yapws87 mentioned, there was indeed a problem with the format I was presenting to the watershed function.
Doing try_data=ma_data.astype(np.uint8) removed the error message.
Here is a minimal example that works now:
import os,sys
from netCDF4 import Dataset as nc
import cv2
import numpy as np
import scipy.ndimage as ndimage
from skimage.morphology import watershed
from scipy import ndimage
basename="/home/dcop696/Data/CMORPH/precip/CMORPH_V1.0/CRT/8km-30min/1999/"
filename=["Cmorph-1999_01_03.nc"]
fileslm=["/home/dcop696/Data/LSM/Cmorph_slm_8km.nc"]
nc_data=nc(basename+filename[0])
data=nc_data.variables["CMORPH"][23,0:250,250:750]
new_data=np.flipud(data)
ma_data=np.ma.masked_where(new_data<=0,new_data)
try_data=ma_data.astype(np.uint8)
## Building threshold
tmp_data=ma_data.filled(0)
tmp_data[np.where(tmp_data!=0)]=255
bw_data=tmp_data.astype(np.uint8)
## Building markers
ma_conv=np.ma.masked_where(new_data<=2,new_data)
tmp_conv=ma_conv.filled(0)
tmp_conv[np.where(tmp_conv!=0)]=255
bw_conv=tmp_conv.astype(np.uint8)
markers = ndimage.label(bw_conv)[0]
## Watershed
labels = watershed(-try_data, markers, mask=bw_data)