diff --git a/functions/sds.py b/functions/sds.py index e9b73fe..8bb8fa6 100644 --- a/functions/sds.py +++ b/functions/sds.py @@ -610,6 +610,7 @@ def convert_epsg(points, epsg_in, epsg_out): def classify_sand_unsupervised(im_ms_ps, im_pan, cloud_mask, wl_pix, buffer_size, min_beach_size, plot_bool): """ Classifies sand pixels using an unsupervised algorithm (Kmeans) + Set buffer size to False if you KV WRL 2018 @@ -623,8 +624,9 @@ def classify_sand_unsupervised(im_ms_ps, im_pan, cloud_mask, wl_pix, buffer_size 2D cloud mask with True where cloud pixels are wl_pix: list of np.ndarray list of arrays containig the pixel coordinates of the water line - buffer_size: int + buffer_size: int or False radius of the disk used to create a buffer around the water line + when False, the entire image is considered for kmeans min_beach_size: int minimum number of connected pixels belonging to a single beach plot_bool: boolean @@ -644,16 +646,19 @@ def classify_sand_unsupervised(im_ms_ps, im_pan, cloud_mask, wl_pix, buffer_size vec_features[:,[0,1,2,3]] = vec_ms_ps[:,[0,1,2,3]] vec_features[:,4] = vec_pan - # create binary image with ones where the detected water lines is - im_buffer = np.zeros((im_ms_ps.shape[0], im_ms_ps.shape[1])) - for i, contour in enumerate(wl_pix): - indices = [(int(_[0]), int(_[1])) for _ in list(np.round(contour))] - for j, idx in enumerate(indices): - im_buffer[idx] = 1 - # perform a dilation on the binary image - se = morphology.disk(buffer_size) - im_buffer = morphology.binary_dilation(im_buffer, se) - vec_buffer = (im_buffer == 1).reshape(im_ms_ps.shape[0] * im_ms_ps.shape[1]) + if buffer_size: + # create binary image with ones where the detected water lines is + im_buffer = np.zeros((im_ms_ps.shape[0], im_ms_ps.shape[1])) + for i, contour in enumerate(wl_pix): + indices = [(int(_[0]), int(_[1])) for _ in list(np.round(contour))] + for j, idx in enumerate(indices): + im_buffer[idx] = 1 + # perform a dilation on the binary image + se = morphology.disk(buffer_size) + im_buffer = morphology.binary_dilation(im_buffer, se) + vec_buffer = (im_buffer == 1).reshape(im_ms_ps.shape[0] * im_ms_ps.shape[1]) + else: + vec_buffer = np.ones((vec_pan.shape[0])) # add cloud mask to buffer vec_buffer= np.logical_and(vec_buffer, ~vec_mask) # perform kmeans (6 clusters) diff --git a/read_images.py b/read_images.py index 8960534..35ec388 100644 --- a/read_images.py +++ b/read_images.py @@ -40,6 +40,8 @@ plot_bool = False # if you want the plots prob_high = 99.9 # upper probability to clip and rescale pixel intensity min_contour_points = 100# minimum number of points contained in each water line output_epsg = 28356 # GDA94 / MGA Zone 56 +buffer_size = 10 # radius (in pixels) of disk for buffer (pixel classification) +min_beach_size = 50 # number of pixels in a beach (pixel classification) # load metadata (timestamps and epsg code) for the collection satname = 'L8' @@ -127,7 +129,7 @@ for i in range(N): idx_skipped.append(i) continue else: - del shorelines[idx_samedate] +# del shorelines[idx_samedate] del t[idx_samedate] del cloud_cover_ts[idx_samedate] del date_acquired_ts[idx_samedate] @@ -150,39 +152,43 @@ for i in range(N): # convert to output epsg spatial reference wl = sds.convert_epsg(wl_coords, input_epsg, output_epsg) - # plot a figure to select the correct water line and discard cloudy images - plt.figure() - cmap = cm.get_cmap('jet') - plt.subplot(121) - plt.imshow(im_ms_ps[:,:,[2,1,0]]) - for j,contour in enumerate(wl_pix): - colours = cmap(np.linspace(0, 1, num=len(wl_pix))) - plt.plot(contour[:, 1], contour[:, 0], linewidth=2, color=colours[j,:]) - plt.axis('image') - plt.title(file_names_pan[i]) - plt.subplot(122) - centroids = [] - for j,contour in enumerate(wl): - colours = cmap(np.linspace(0, 1, num=len(wl))) - centroids.append([np.mean(contour[:, 0]),np.mean(contour[:, 1])]) - plt.plot(contour[:, 0], contour[:, 1], linewidth=2, color=colours[j,:]) - plt.plot(np.mean(contour[:, 0]), np.mean(contour[:, 1]), 'o', color=colours[j,:]) - plt.plot(refpoints[:,0], refpoints[:,1], 'k.') - plt.axis('equal') - plt.title(file_names_pan[i]) - mng = plt.get_current_fig_manager() - mng.window.showMaximized() - plt.tight_layout() - plt.draw() - # click on the left image to discard, otherwise on the closest centroid in the right image - pt_in = np.array(ginput(n=1, timeout=1000)) - if pt_in[0][0] < 10000: - print('skip ' + str(i) + ' - manual') - idx_skipped.append(i) - continue - # get contour that was selected (clock closest to centroid) - dist_centroid = [np.linalg.norm(_ - pt_in) for _ in centroids] - shorelines.append(wl[np.argmin(dist_centroid)]) + # classify sand pixels + im_sand = sds.classify_sand_unsupervised(im_ms_ps, im_pan, cloud_mask, wl_pix, False, min_beach_size, True) + +# # plot a figure to select the correct water line and discard cloudy images +# plt.figure() +# cmap = cm.get_cmap('jet') +# plt.subplot(121) +# plt.imshow(im_ms_ps[:,:,[2,1,0]]) +# for j,contour in enumerate(wl_pix): +# colours = cmap(np.linspace(0, 1, num=len(wl_pix))) +# plt.plot(contour[:, 1], contour[:, 0], linewidth=2, color=colours[j,:]) +# plt.axis('image') +# plt.title(file_names_pan[i]) +# plt.subplot(122) +# centroids = [] +# for j,contour in enumerate(wl): +# colours = cmap(np.linspace(0, 1, num=len(wl))) +# centroids.append([np.mean(contour[:, 0]),np.mean(contour[:, 1])]) +# plt.plot(contour[:, 0], contour[:, 1], linewidth=2, color=colours[j,:]) +# plt.plot(np.mean(contour[:, 0]), np.mean(contour[:, 1]), 'o', color=colours[j,:]) +# plt.plot(refpoints[:,0], refpoints[:,1], 'k.') +# plt.axis('equal') +# plt.title(file_names_pan[i]) +# mng = plt.get_current_fig_manager() +# mng.window.showMaximized() +# plt.tight_layout() +# plt.draw() +# # click on the left image to discard, otherwise on the closest centroid in the right image +# pt_in = np.array(ginput(n=1, timeout=1000)) +# if pt_in[0][0] < 10000: +# print('skip ' + str(i) + ' - manual') +# idx_skipped.append(i) +# continue +# # get contour that was selected (click closest to centroid) +# dist_centroid = [np.linalg.norm(_ - pt_in) for _ in centroids] +# shorelines.append(wl[np.argmin(dist_centroid)]) + t.append(timestamps_sorted[i]) cloud_cover_ts.append(cloud_cover) acc_georef_ts.append(acc_georef_sorted[i])