Evaluate HARP2 BRDF¶

from pathlib import Path

import cartopy.crs as ccrs
import earthaccess
import matplotlib.pyplot as plt
import numpy as np
import requests
import xarray as xr
from matplotlib.colors import LogNorm


auth = earthaccess.login(persist=True)
fs = earthaccess.get_fsspec_https_session()

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from matplotlib.colors import LogNorm


def plot_l2_product(
    data, plot_range, label, title,
    vmin=None, vmax=None,
    figsize=(12, 4),
    cmap="viridis",
    log_scale=False,
    land_color="#d9d9d9",
    ocean_color="#dbe9f6",
):
    """Make map + histogram with optional log color scaling and land/ocean background.

    Notes:
      - Assumes lon/lat are 2D or 1D arrays available in the outer scope
        (or change signature to pass them in).
      - For log_scale=True, only positive values are used for autoscaling and histogram.
    """

    # ------------------
    # Determine vmin / vmax if not given
    # ------------------
    valid = data[np.isfinite(data)]
    if valid.size == 0:
        raise ValueError("No finite values in `data`.")

    if log_scale:
        valid = valid[valid > 0]
        if valid.size == 0:
            raise ValueError("log_scale=True but `data` has no positive finite values.")

        if vmin is None:
            vmin = np.percentile(valid, 2)
        if vmax is None:
            vmax = np.percentile(valid, 98)

        # Safety: avoid invalid/degenerate bounds
        if (vmin is None) or (vmax is None) or (vmin <= 0) or (vmin >= vmax):
            vmin = float(np.min(valid))
            vmax = float(np.max(valid))

    else:
        if vmin is None:
            vmin = np.percentile(valid, 2)
        if vmax is None:
            vmax = np.percentile(valid, 98)

        if (vmin is None) or (vmax is None) or (vmin >= vmax):
            vmin = float(np.min(valid))
            vmax = float(np.max(valid))

    # ------------------
    # Figure layout
    # ------------------
    fig = plt.figure(figsize=figsize)
    gs = fig.add_gridspec(1, 2, width_ratios=[3, 1], wspace=0.3)

    # ------------------
    # Map subplot
    # ------------------
    ax_map = fig.add_subplot(gs[0], projection=ccrs.PlateCarree())
    ax_map.set_extent(plot_range, crs=ccrs.PlateCarree())

    # Land / ocean background (behind data)
    ax_map.add_feature(cfeature.OCEAN, facecolor=ocean_color, zorder=0)
    ax_map.add_feature(cfeature.LAND, facecolor=land_color, zorder=1)

    ax_map.coastlines(resolution="110m", color="black", linewidth=0.8)
    ax_map.gridlines(draw_labels=True)

    norm = LogNorm(vmin=vmin, vmax=vmax) if log_scale else None

    pm = ax_map.pcolormesh(
        lon, lat, data,
        norm=norm,
        vmin=None if log_scale else vmin,
        vmax=None if log_scale else vmax,
        transform=ccrs.PlateCarree(),
        cmap=cmap,
        zorder=2
    )

    cbar = plt.colorbar(pm, ax=ax_map, orientation="vertical", pad=0.1)
    cbar.set_label(label)

    ax_map.set_title(title, fontsize=12)

    # ------------------
    # Histogram subplot
    # ------------------
    ax_hist = fig.add_subplot(gs[1])

    hist_data = data[np.isfinite(data)]
    if log_scale:
        hist_data = hist_data[hist_data > 0]
        bins = np.logspace(np.log10(vmin), np.log10(vmax), 40)
        ax_hist.hist(hist_data, bins=bins, color="gray", edgecolor="black")
        ax_hist.set_xscale("log")
    else:
        ax_hist.hist(
            hist_data, bins=40, range=[vmin, vmax],
            color="gray", edgecolor="black"
        )

    ax_hist.set_xlabel(label)
    ax_hist.set_ylabel("Count")
    ax_hist.set_title(f"Histogram: N={hist_data.size}")

    plt.show()

def ensure_same_latlon(lat1, lon1, lat2, lon2, rtol=1e-6, atol=1e-8):
    assert lat1.shape == lat2.shape, "Latitude shape mismatch"
    assert lon1.shape == lon2.shape, "Longitude shape mismatch"

    assert np.allclose(lat1, lat2, rtol=rtol, atol=atol), (
        f"Latitude values differ (max diff = {np.nanmax(np.abs(lat1 - lat2))})"
    )
    assert np.allclose(lon1, lon2, rtol=rtol, atol=atol), (
        f"Longitude values differ (max diff = {np.nanmax(np.abs(lon1 - lon2))})"
    )
    print("YES")

datatree1 = xr.open_datatree(path1[0],decode_timedelta=False)
dataset1 = xr.merge(datatree1.to_dict().values())