Shapefile ➜ Float32 GeoTIFF (DEPTH2D) — Interactive + Batch

import warnings, math
from pathlib import Path
import geopandas as gpd
from shapely import make_valid
import rasterio
from rasterio.transform import from_origin
from rasterio.features import rasterize
from rasterio.enums import MergeAlg
import matplotlib.pyplot as plt

# SINGLE TEST
test_shp = Path(r'D:/gis/single_test/one.shp')   # change

# BATCH
batch_root   = Path(r'D:/gis/batch_root')        # each sub‑folder = timestamp
batch_output = Path(r'D:/gis/tif_output')
batch_output.mkdir(parents=True, exist_ok=True)

print('Batch root :', batch_root.resolve())
print('Output root:', batch_output.resolve())

PIXEL_LIMIT     = 4096
PIXEL_SIZE_FINE = 1e-5    # deg pls defined by yourself
NODATA          = -9999.0
#MERGE_STRATEGY  = 'max'   # max/min/first
FRAME_FIELD     = 'FRAME' # attribute for frame number

def repair(gdf):
    def _f(geom):
        if geom is None or geom.is_empty:
            return None
        if not geom.is_valid:
            try:
                return make_valid(geom)
            except Exception:
                return geom.buffer(0)
        return geom
    gdf['geometry'] = gdf.geometry.apply(_f)
    return gdf[~gdf.geometry.is_empty & gdf.geometry.notnull()]

def choose_px(bounds):
    minx, miny, maxx, maxy = bounds
    dx, dy = maxx - minx, maxy - miny
    if dx/PIXEL_SIZE_FINE <= PIXEL_LIMIT and dy/PIXEL_SIZE_FINE <= PIXEL_LIMIT:
        return PIXEL_SIZE_FINE
    px = max(dx, dy)/PIXEL_LIMIT
    if px > PIXEL_SIZE_FINE:
        warnings.warn(f'Pixel size relaxed to {px:.6f}')
    return px

def rasterize_depth(gdf):
    minx, miny, maxx, maxy = gdf.total_bounds
    px = choose_px((minx, miny, maxx, maxy))
    width  = math.ceil((maxx - minx)/px)
    height = math.ceil((maxy - miny)/px)
    transform = from_origin(minx, maxy, px, px)

    shapes = ((row.geometry, float(row['DEPTH2D'])) for _, row in gdf.iterrows())
    arr = rasterize(shapes, out_shape=(height, width),
                    fill=NODATA, dtype='float32',
                    transform=transform)
    return arr, transform

if test_shp.exists():
    gdf = gpd.read_file(test_shp)
    if 'DEPTH2D' not in gdf.columns:
        raise KeyError('DEPTH2D missing')
    gdf = repair(gdf)
    arr, transform = rasterize_depth(gdf)

    minx, miny, maxx, maxy = gdf.total_bounds
    fig, ax = plt.subplots(figsize=(6,6))
    im = ax.imshow(arr, extent=[minx,maxx,miny,maxy], origin='upper', cmap='viridis', alpha=0.6)
    gdf.boundary.plot(ax=ax, edgecolor='white', linewidth=0.4)
    plt.colorbar(im, ax=ax, fraction=0.022, pad=0.01, label='DEPTH2D')
    ax.set_axis_off()
    plt.show()
else:
    print('Set `test_shp` path to an existing shapefile for preview.')

folders = [p for p in batch_root.iterdir() if p.is_dir()]
print(f'Folders found: {len(folders)}')

for ts_dir in folders:
    timestamp = ts_dir.name
    shp_paths = list(ts_dir.glob('*.shp'))
    if not shp_paths:
        print(f'[SKIP] {timestamp}: no shp')
        continue
    shp_path = shp_paths[0]

    try:
        gdf = gpd.read_file(shp_path)
        if 'DEPTH2D' not in gdf.columns:
            print(f'[SKIP] {shp_path.name}: DEPTH2D missing')
            continue
        gdf = repair(gdf)
        frame_val = (str(gdf[FRAME_FIELD].dropna().iloc[0])
                     if FRAME_FIELD in gdf.columns and not gdf[FRAME_FIELD].dropna().empty
                     else shp_path.stem)
        arr, transform = rasterize_depth(gdf)
        out_tif = batch_output / f'{timestamp}_{frame_val}.tif'

        meta = dict(driver='GTiff',
                    height=arr.shape[0], width=arr.shape[1],
                    count=1, dtype='float32', transform=transform,
                    crs=gdf.crs, nodata=NODATA)

        with rasterio.open(out_tif, 'w', **meta) as dst:
            dst.write(arr, 1)

        print(f'[OK] {timestamp}/{shp_path.name} → {out_tif.name}')
    except Exception as e:
        print(f'[ERR] {shp_path}: {e}')