Custom coastline shapefile from high resolution elevation data

Question

Struggling for days, I'm seeking help for creating a custom coastline shapefile based on a high-res digital elevation dataset in GeoTIFF format. I'm willing to do so because the default coastline is shifted and not accurate enough for my needs (climatology) ; see example screenshot.

Default coastline vs. DEM mask

The aim behind creating this shapefile is that I can re-use it while plotting climate data over the area.

For now, here's what I got :

import os
import rasterio
import numpy as np
import geopandas as gpd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from shapely.geometry import Point

def create_coastline_from_dem(dem_path, coastline_path):
    # Read DEM data
    with rasterio.open(dem_path) as dem:
        dem_data = dem.read(1)
        transform = dem.transform

    # Create coastline data
    coastline_data = np.zeros_like(dem_data, dtype=np.uint8)

    # Define a simple threshold to identify coastline
    threshold = 0  # Adjust this threshold based on your DEM characteristics

    coastline_data[dem_data <= threshold] = 1

    # Write coastline data to a new GeoTIFF file
    with rasterio.open(coastline_path, 'w', driver='GTiff', height=dem_data.shape[0],
                       width=dem_data.shape[1], count=1, dtype='uint8', crs=dem.crs,
                       transform=transform) as dst:
        dst.write(coastline_data, 1)

def create_coastline_geojson(coastline_path, output_geojson):
    # Read coastline data
    with rasterio.open(coastline_path) as coastline:
        transform = coastline.transform

        # Identify coastline pixels
        coastline_pixels = np.where(coastline.read(1) == 1)

        # Convert pixel coordinates to geographic coordinates
        lon, lat = rasterio.transform.xy(transform, coastline_pixels[0], coastline_pixels[1])

        # Create a GeoDataFrame with the coastline as points
        gdf = gpd.GeoDataFrame({'geometry': [Point(x, y) for x, y in zip(lon, lat)]}, crs=coastline.crs)

        # Save the GeoDataFrame to a GeoJSON file
        gdf.to_file(output_geojson, driver='GeoJSON')

def plot_coastline_geojson(geojson_path):
    # Read GeoDataFrame from GeoJSON
    gdf = gpd.read_file(geojson_path)

    # Calculate the central longitude and latitude
    central_lon = (gdf.bounds.minx + gdf.bounds.maxx) / 2
    central_lat = (gdf.bounds.miny + gdf.bounds.maxy) / 2

    # Create a plot with equal aspect ratio
    fig, ax = plt.subplots(subplot_kw={'projection': ccrs.PlateCarree(central_longitude=central_lon)})
    ax.set_aspect('equal', adjustable='datalim')

    # Plot coastline
    gdf.plot(ax=ax, markersize=0.1, color='black')

    # Add gridlines
    ax.gridlines(draw_labels=True)

    # Add coastlines
    ax.coastlines(resolution='10m')

    plt.show()

if __name__ == "__main__":
    # Get the path to the current directory
    current_directory = os.path.dirname(os.path.abspath(__file__))

    # Build the paths to the DEM and coastline files
    dem_file = os.path.join(current_directory, "tahi.tiff")
    coastline_file = os.path.join(current_directory, "coastline.tif")
    coastline_geojson = os.path.join(current_directory, "coastline.geojson")

    create_coastline_from_dem(dem_file, coastline_file)
    create_coastline_geojson(coastline_file, coastline_geojson)
    plot_coastline_geojson(coastline_geojson)

Edit ; that may be helpful :

gdalinfo tahi.tiff 
Driver: GTiff/GeoTIFF
Files: tahi.tiff
Size is 625, 601
Origin = (-150.261250000000217,-16.798750000000002)
Pixel Size = (0.002500000000000,-0.002500000000000)
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-150.2612500, -16.7987500) 
Lower Left  (-150.2612500, -18.3012500) 
Upper Right (-148.6987500, -16.7987500) 
Lower Right (-148.6987500, -18.3012500) 
Center      (-149.4800000, -17.5500000) 
Band 1 Block=625x3 Type=Float32, ColorInterp=Gray
  Min=0.000 Max=1980.548 
  Minimum=0.000, Maximum=1980.548, Mean=17.268, StdDev=106.653
  NoData Value=nan
  Metadata:
    STATISTICS_MAXIMUM=1980.5482177734
    STATISTICS_MEAN=17.268496548744
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=106.65328874517
    STATISTICS_VALID_PERCENT=100

2023/12/04 EDIT :

I did some more work and came closer, but am still facing some issues. For now, my script looks like :

import numpy as np
import matplotlib.pyplot as plt
from osgeo import gdal
import fiona
from shapely.geometry import Polygon, mapping
from shapely.ops import unary_union
from cartopy import crs as ccrs
from cartopy.feature import ShapelyFeature
from cartopy.io import shapereader

def read_geotiff(file_path):
    dataset = gdal.Open(file_path)
    if dataset is None:
        raise Exception("Error opening dataset!")

    band = dataset.GetRasterBand(1)
    data = band.ReadAsArray()
    geotransform = dataset.GetGeoTransform()

    return data, geotransform

def create_coastline_shapefile(data, geotransform, output_shapefile):
    contours = plt.contour(geotransform[0] + geotransform[1] * np.arange(data.shape[1]),
                           geotransform[3] + geotransform[5] * np.arange(data.shape[0]),
                           data, levels=[0])

    # Extract polygons from the contours
    polygons = [Polygon(line.vertices) for line in contours.collections[0].get_paths()]

    # Merge adjacent polygons
    merged_polygon = unary_union(polygons)

    # Write the merged polygon to a shapefile
    schema = {
        'geometry': 'Polygon',
        'properties': {'id': 'int'},
    }

    with fiona.open(output_shapefile, 'w', 'ESRI Shapefile', schema) as output:
        output.write({
            'geometry': mapping(merged_polygon),
            'properties': {'id': 1},
        })

def plot_custom_coastline(ax, shapefile):
    shp_feature = ShapelyFeature(shapereader.Reader(shapefile).geometries(),
                                 ccrs.PlateCarree(), edgecolor='black', facecolor='none')
    ax.add_feature(shp_feature, zorder=2)

def plot_default_coastline(ax):
    ax.coastlines(resolution='10m', color='black', linewidth=1, zorder=1)

def plot_comparison(data, geotransform, custom_shapefile):
    fig, ax = plt.subplots(subplot_kw={'projection': ccrs.PlateCarree()})
    ax.set_extent([geotransform[0], geotransform[0] + geotransform[1] * data.shape[1],
                   geotransform[3] + geotransform[5] * data.shape[0], geotransform[3]], crs=ccrs.PlateCarree())

    # Plot the default coastline
    plot_default_coastline(ax)

    # Plot the custom coastline
    plot_custom_coastline(ax, custom_shapefile)

    plt.show()

if __name__ == "__main__":
    tiff_file = "tahi.tiff"
    output_shapefile = "coastline_shapefile.shp"

    data, geotransform = read_geotiff(tiff_file)
    create_coastline_shapefile(data, geotransform, output_shapefile)
    plot_comparison(data, geotransform, output_shapefile)

This gets me these figures. Figure 1 looks pretty close to the shapefile I need. But when I plot the custom shapefile created, it shows some unwanted lines between lone pixels. Figure 1 and 2 screenshot