Bathymetry and Metadata

# Import Path to handle the raster path
from pathlib import Path

# Import NumPy for numerical array operations
import numpy as np

# Import Matplotlib for plotting
import matplotlib.pyplot as plt

# Help to create a nice colorbar
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib.ticker as mticker

# Formatter without scientific notation
formatter = mticker.FuncFormatter(lambda x, pos: f"{int(x):,}")

# Import RasterIO to read the raster data
import rasterio as rio

# Import plotting_extent to display the raster in map coordinates
from rasterio.plot import plotting_extent

import geopandas as gpd
from shapely import wkt
from shapely.geometry import Point,box
import folium

data_path = Path("../data/gifford_bathy.tif")

# rasterio as rio
dataset = rio.open(data_path)
dataset

# print metadata:
print("CRS:", dataset.crs)
print("Bounds:", dataset.bounds)
print("Width:", dataset.width)
print("Height:", dataset.height)
print("Resolution:", dataset.res)
print("Transform:", dataset.transform)
print("Number of bands:", dataset.count)
print("NoData value:", dataset.nodata)
print("Data type:", dataset.dtypes)

dict(dataset.profile)

band1 = dataset.read(1)
print("type :", type(band1).__name__)
print("shape:", band1.shape)
print("dtype:", band1.dtype)

bathy_data = dataset.read(1, masked=True)
bathy_data

bathy_nan = bathy_data.filled(np.nan)

print("shallowest:", np.nanmax(bathy_nan), "m")
print("deepest   :", np.nanmin(bathy_nan), "m")
print("mean depth:", np.nanmean(bathy_nan), "m")

# Pixel (row, col) -> map (x, y)
row, col = 200, 300
x, y = dataset.xy(row, col)
print(f"pixel (row={row}, col={col}) -> map (x={x:.1f}, y={y:.1f})")

# Map (x, y) -> pixel (row, col)
r, c = dataset.index(x, y)
print(f"map (x={x:.1f}, y={y:.1f}) -> pixel (row={r}, col={c})")

dataset.close()

bathy_path = '../data/gifford_bathy.tif'

# Raster reading
# --------------
with rio.open(bathy_path) as dataset:
    bathy = dataset.read(1, masked=True).astype(float).filled(np.nan)
    transform = dataset.transform
    crs = dataset.crs
    bounds = dataset.bounds   # <-- correct way to get raster bounds

    # print metadata:
    print("CRS:", dataset.crs)
    print("Bounds:", dataset.bounds)
    print("Width:", dataset.width)
    print("Height:", dataset.height)
    print("Resolution:", dataset.res)
    print("Transform:", dataset.transform)
    print("Number of bands:", dataset.count)
    print("NoData value:", dataset.nodata)
    print("Data type:", dataset.dtypes)

    
# Resolution
dx = abs(transform.a)
dy = abs(transform.e)

print("shape:", bathy.shape)
print("crs:", crs)
print("cell size (resolution):", (dx, dy))
print("depth range:", (np.nanmin(bathy), np.nanmax(bathy)))

crs

print(crs)

import random

class MysteriousObject:
    # Printable representation (more for developers)
    def __repr__(self):
        return f"MysteriousObject({random.randint(1, 100)})"

    # Printable String (more for users)
    def __str__(self):
        return f"I'm the mysterious object number {random.randint(1, 100)}"

# Create an instance
obj = MysteriousObject()

# Not using print (the __repr__ method is used)
obj

# Using print (the __str__ method is used)
print(obj)

print(crs.__repr__())

print(crs.to_wkt())

# Assuming o_data.crs.to_wkt() returns a valid WKT string
wkt_string = crs.to_wkt()

# Create a GeoDataFrame with a single geometry from the WKT
# A geometry is needed to define a GeoDataFrame, it is just a dummy value here
gdf = gpd.GeoDataFrame(geometry=[wkt.loads("POINT(0 0)")], crs=wkt_string)

# Display the CRS information in a nicely formatted way
print("Formatted CRS WKT:")
print(gdf.crs.to_wkt(pretty=True))

# Get the bounding box of the raster and compute its center
center_x = (bounds.left + bounds.right) / 2
center_y = (bounds.bottom + bounds.top) / 2

# Create a GeoDataFrame using the raster's native CRS
gdf = gpd.GeoDataFrame(
    geometry=[Point(center_x, center_y)],
    crs=crs
)

# Print the original CRS information
print("Original CRS WKT:")
print(gdf.crs.to_wkt(pretty=True))

# Reproject the GeoDataFrame to EPSG:4326 for proper display on Folium
gdf = gdf.to_crs(epsg=4326)
print("Reprojected CRS WKT (EPSG:4326):")
print(gdf.crs.to_wkt(pretty=True))

# Extract the reprojected coordinates (latitude, longitude)
center_point = gdf.geometry.iloc[0]
latitude, longitude = center_point.y, center_point.x

# Create a Folium map centered on the reprojected point
m = folium.Map(location=[latitude, longitude], zoom_start=10)
folium.Marker([latitude, longitude], popup="Raster Center").add_to(m)

# To save the map, uncomment the following line:
# m.save("raster_center_map.html")

# Display the map
m

# Create polygon from raster bounds in native CRS
bounds_polygon = box(bounds.left, bounds.bottom, bounds.right, bounds.top)

bounds_gdf = gpd.GeoDataFrame(
    geometry=[bounds_polygon],
    crs=crs
)

# Reproject to WGS84 for Folium
bounds_gdf_4326 = bounds_gdf.to_crs(epsg=4326)

# Reproject center point too, if not already done
gdf_4326 = gdf.to_crs(epsg=4326)

center_point = gdf_4326.geometry.iloc[0]
latitude, longitude = center_point.y, center_point.x

m = folium.Map(location=[latitude, longitude], zoom_start=5)

folium.GeoJson(
    bounds_gdf_4326,
    name="Raster Bounds",
    tooltip="Raster footprint"
).add_to(m)

popup_text = f"""
<b>Raster Center</b><br>
CRS: {crs}<br>
Bounds: {bounds}<br>
Resolution: {transform.a:.2f}, {abs(transform.e):.2f}<br>
Shape: {bathy.shape}
"""

folium.Marker(
    [latitude, longitude],
    popup=folium.Popup(popup_text, max_width=400)
).add_to(m)

folium.LayerControl().add_to(m)

m

DATASETS = {
    "Gifford Canyon": Path("../data/gifford_bathy.tif"),
    "Oceanic Shoals": Path("../data/os_bathy.tif"),
    "Point Cloates": Path("../data/pc_bathy.tif"),
}

def get_raster_metadata(path):
    """Read raster metadata and return footprint in EPSG:4326 plus metadata."""
    with rio.open(path) as src:
        data = src.read(1, masked=True)

        bounds_polygon = box(
            src.bounds.left,
            src.bounds.bottom,
            src.bounds.right,
            src.bounds.top,
        )

        footprint_native = gpd.GeoDataFrame(
            geometry=[bounds_polygon],
            crs=src.crs,
        )

        footprint_4326 = footprint_native.to_crs(epsg=4326)

        metadata = {
            "file": path.name,
            "crs": str(src.crs),
            "bounds": src.bounds,
            "resolution_x": src.transform.a,
            "resolution_y": abs(src.transform.e),
            "shape": src.shape,
            "width": src.width,
            "height": src.height,
            "band_count": src.count,
            "nodata": src.nodata,
            "dtype": src.dtypes[0],
            "depth_min": float(data.min()),
            "depth_max": float(data.max()),
        }

    return footprint_4326, metadata

site_footprints = {}

for site_name, path in DATASETS.items():
    footprint_4326, metadata = get_raster_metadata(path)
    site_footprints[site_name] = {
        "footprint": footprint_4326,
        "metadata": metadata,
    }

def raster_bounds_to_wgs84(path):
    """Create a WGS84 footprint polygon from raster metadata."""
    with rio.open(path) as src:
        bounds = src.bounds
        crs = src.crs
        transform = src.transform

        bounds_polygon = box(
            bounds.left,
            bounds.bottom,
            bounds.right,
            bounds.top,
        )

        bounds_gdf = gpd.GeoDataFrame(
            geometry=[bounds_polygon],
            crs=crs,
        )

        bounds_gdf_4326 = bounds_gdf.to_crs(epsg=4326)

        metadata = {
            "crs": crs,
            "bounds": bounds,
            "resolution": (transform.a, abs(transform.e)),
            "shape": src.shape,
            "width": src.width,
            "height": src.height,
            "nodata": src.nodata,
            "dtype": src.dtypes[0],
            "band_count": src.count,
        }

    return bounds_gdf_4326, metadata

# Use all site centroids to center the map
centroids = [
    item["footprint"].geometry.iloc[0].centroid
    for item in site_footprints.values()
]

center_latitude = sum(point.y for point in centroids) / len(centroids)
center_longitude = sum(point.x for point in centroids) / len(centroids)


m = folium.Map(
    location=[center_latitude, center_longitude],
    zoom_start=3,
    tiles="CartoDB positron",
)

for site_name, item in site_footprints.items():
    footprint = item["footprint"]
    metadata = item["metadata"]

    center_point = footprint.geometry.iloc[0].centroid
    latitude, longitude = center_point.y, center_point.x

    popup_text = f"""
    <b>{site_name}</b><br>
    CRS: {metadata["crs"]}<br>
    Bounds: {metadata["bounds"]}<br>
    Resolution: {metadata["resolution_x"]:.2f} m, {metadata["resolution_y"]:.2f} m<br>
    Shape: {metadata["shape"]}<br>
    Width × height: {metadata["width"]} × {metadata["height"]}<br>
    Bands: {metadata["band_count"]}<br>
    Data type: {metadata["dtype"]}<br>
    NoData: {metadata["nodata"]}
    """

    folium.GeoJson(
        footprint,
        name=site_name,
        tooltip=f"{site_name} raster footprint",
        style_function=lambda feature: {
            "fillOpacity": 0.2,
            "weight": 2,
        },
    ).add_to(m)

    folium.Marker(
        location=[latitude, longitude],
        popup=folium.Popup(popup_text, max_width=450),
        tooltip=site_name,
    ).add_to(m)


folium.LayerControl().add_to(m)

m

# Save to HTML file
m.save("./images/map.html")