Additional thresholding outputs for zoo workflow

[mlundine]

Was playing around with different Otsu thresholding outputs and am thinking we should have a few of these computed in the zoo workflow in addition to the deep learning model outputs. See little sample code below:

`
def compute_otsu_threshold(in_tiff, out_tiff):
"""
Otsu binary thresholding on a geotiff.
inputs:
in_tiff (str): path to the input geotiff
out_tiff (str): path to the output geotiff
outputs:
out_tiff (str): path to the output geotiff
"""
def rescale(arr):
arr_min = arr.min()
arr_max = arr.max()
return (arr - arr_min) / (arr_max - arr_min)

def average_img(images_list):
    # Alternative method using numpy mean function
    images = np.array(images_list)
    arr = np.array(np.mean(images, axis=(0)), dtype=np.uint8)
    return arr

with rasterio.open(in_tiff) as src:
    blue = src.read(1)
    green = src.read(2)
    red = src.read(3)
    nir = src.read(4)
    swir = src.read(5)
        
ndwi = (green - nir)/(green+nir)
mndwi = (green - swir)/(green+swir)
rgb = np.dstack([red,green,blue])
rgb_disp = 255.0*rescale(np.dstack([red,green,blue]))
stacked = np.dstack([red,green,blue,nir,swir,ndwi,mndwi])

# Compute Otsu's threshold
# Need to make nodata values zero or else the threshold will be just data vs. nodata    
swir[swir==src.meta['nodata']]=0
nir[nir==src.meta['nodata']]=0
blue[blue==src.meta['nodata']]=0
red[red==src.meta['nodata']]=0
green[green==src.meta['nodata']]=0

stacked[stacked==src.meta['nodata']]=0
rgb[rgb==src.meta['nodata']]=0

##compute thresholds
thresholds_swir = threshold_multiotsu(swir)
thresholds_nir = threshold_multiotsu(nir)

# Apply the threshold to create a binary image
binary_image_swir = swir > min(thresholds_swir)
binary_image_nir = nir > min(thresholds_nir)

# comparison plot, not much difference between swir and nir
plt.subplot(3,1,1)
plt.title('RGB')
plt.imshow(rgb_disp.astype('uint8'))
plt.xticks([],[])
plt.yticks([],[])
plt.subplot(3,1,2)
plt.title('Otsu Threshold SWIR')
plt.imshow(rgb_disp.astype('uint8'))
plt.imshow(binary_image_swir, alpha=0.25, cmap='jet')
plt.xticks([],[])
plt.yticks([],[])
plt.subplot(3,1,3)
plt.title('Otsu Threshold NIR')
plt.imshow(rgb_disp.astype('uint8'))
plt.imshow(binary_image_nir, alpha=0.25, cmap='jet')
plt.xticks([],[])
plt.yticks([],[])
plt.savefig(os.path.splitext(out_tiff)[0]+'_results.jpg',dpi=500)
plt.close('all')

# Define the metadata for the new geotiff
transform = from_origin(src.bounds.left, src.bounds.top, src.res[0], src.res[1])
new_meta = src.meta.copy()
new_meta.update({
    'dtype': 'uint8',
    'count': 1,
    'transform': transform,
    'nodata':None
})

# Save the binary image
with rasterio.open(out_tiff, 'w', **new_meta) as dst:
    dst.write(binary_image_nir.astype(np.uint8), 1)

return out_tiff

`

Basically, we could have zoo always output the NIR and SWIR Otsu threshold as a fallback/comparison method. So we could get at uncertainty bands as the average of the deep learning output, NIR Otsu, and SWIR Otsu.

[mlundine]

https://github.com/Doodleverse/SDStools/blob/main/src/sdstools/thresholding.py

[mlundine]

This was the toy script i was experimenting with. I am going to put together an averaging function on the shorelines for getting uncertainty bands.

[2320sharon]

This looks interesting for sure. So for the zoo workflow we would save these new shorelines in addition to the shorelines we save currently correct?

Thank you for including the functions in this issue as well as linking your script with the example at the bottom

Excerpt from https://github.com/Doodleverse/SDStools/blob/main/src/sdstools/thresholding.py

##merged_lines = os.path.join(out_folder, 'shorelines.geojson')
##batch_compute_otsu_threshold(images, out_folder_otsu)
##batch_binary_raster_to_vector(out_folder_otsu, out_folder_otsu_contours)
##batch_contours_to_lines(out_folder_otsu_contours,out_folder_otsu_lines)
##merge_lines(out_folder_otsu_lines,merged_lines)

[mlundine]

It might be better just as a function from SDSTools than as an output from coastseg right now, my thought is just to have a function like this: https://github.com/Doodleverse/SDStools/blob/main/src/sdstools/uncertainty_bands.py

[mlundine]

In this case, you'd have all of your outputs from coastseg (tidally_corrected_transect_time_series_merged.csv) and then just make a list of these and merge them to get the mean and conf interval outputs.

[mlundine]

the thresholding techniques could probably just be two additional zoo models. otsu(nir) and otsu(swir).

[dbuscombe-usgs]

Separate script in sds tools for sure while we research best implementation. That could take a while. Loads of potential pathways.

As for mean + confidence interval, that is one approach. Another would be ransac or similar. Robust methods.