Evaluate HARP2 BRDF¶
load libraries and define helper functions¶
In [63]:
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()
In [98]:
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")
Get HARP2 L1C and FastMAPOL L2 data¶
Example scene from rapid response page:
Feb 3 at ocean color website:
- https://oceandata.sci.gsfc.nasa.gov/directdataaccess/Level-2/PACE-HARP2/2026/03-Feb-2026/
- https://oceandata.sci.gsfc.nasa.gov/directdataaccess/Level-1C/PACE-HARP2/2026/03-Feb-2026/
direct data link:
- http://oceandata.sci.gsfc.nasa.gov/getfile/PACE_HARP2.20260203T053050.L1C.MAPOL_OCEAN.V3.5km.nc
- http://oceandata.sci.gsfc.nasa.gov/getfile/PACE_HARP2.20260203T053050.L2.MAPOL_OCEAN.V3_0.NRT.nc
tutorial on visualize fastmapol products:
In [68]:
#download data l1c data
OB_DAAC_PROVISIONAL = "https://oceandata.sci.gsfc.nasa.gov/cgi/getfile/"
HARP2_L2_MAPOL_FILENAME = "PACE_HARP2.20260203T053050.L1C.V3.5km.nc"
fs.get(f"{OB_DAAC_PROVISIONAL}/{HARP2_L2_MAPOL_FILENAME}", "data/")
path1 = list(Path("data").glob("*L1C*.nc"))
path1
Out[68]:
[PosixPath('data/PACE_HARP2.20260203T053050.L1C.V3.5km.nc')]
In [69]:
#download data l2 data
OB_DAAC_PROVISIONAL = "https://oceandata.sci.gsfc.nasa.gov/cgi/getfile/"
HARP2_L2_MAPOL_FILENAME = "PACE_HARP2.20260203T053050.L2.MAPOL_OCEAN.V3_0.NRT.nc"
fs.get(f"{OB_DAAC_PROVISIONAL}/{HARP2_L2_MAPOL_FILENAME}", "data/")
path2 = list(Path("data").glob("*L2*.nc"))
path2
Out[69]:
[PosixPath('data/PACE_HARP2.20260203T053050.L2.MAPOL_OCEAN.V3_0.NRT.nc')]
explore l2 data¶
In [81]:
datatree2 = xr.open_datatree(path2[0])
dataset2 = xr.merge(datatree2.to_dict().values())
lat2 = dataset2["latitude"].values
lon2 = dataset2["longitude"].values
plot_range2 = [lon2.min(), lon2.max(), lat2.min(), lat2.max()]
#per angle wavelength
wavelength2 = np.concat([np.repeat(550,10), np.repeat(670,60), np.repeat(870,10), np.repeat(440,10)])
angle2= dataset2['sensor_view_angle'].values
print("angles per view: ", angle2)
print("wavelegnth per view: ", wavelength2)
wavelength = dataset2["wavelength"].values
print("list of unique wavelength: ", wavelength)
#load Rrs
Rrs1 = dataset2['Rrs1'].values
Rrs2 = dataset2['Rrs2'].values
print("Rrs dimension: ", Rrs1.shape)
angles per view: [ 56.26 43.32 32.57 20.11 5.99 -8.25 -22.02 -34.22 -44.52 -53.32 55.61 54.35 53.05 51.63 50.21 48.77 47.24 45.72 44.19 42.55 40.8 39.05 37.3 35.43 33.56 31.57 29.58 27.58 25.47 23.47 21.23 18.99 16.75 14.22 11.86 9.6 7.12 4.75 2.38 -0.09 -2.38 -4.75 -7.12 -9.49 -11.86 -14.22 -16.41 -18.66 -20.9 -23.02 -25.14 -27.25 -29.36 -31.35 -33.23 -35.1 -36.97 -38.72 -40.48 -42.12 -43.75 -45.28 -46.81 -48.33 -49.86 -51.25 -52.66 -53.97 -55.28 -56.53 55.01 46.48 30.58 17.87 3.62 -10.62 -17.54 -30.36 -46.05 -54.63 53.75 44.96 34.44 22.35 8.36 -5.88 -19.78 -32.23 -42.88 -55.88] wavelegnth per view: [550 550 550 550 550 550 550 550 550 550 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 670 870 870 870 870 870 870 870 870 870 870 440 440 440 440 440 440 440 440 440 440] list of unique wavelength: [440. 550. 670. 870.] Rrs dimension: (396, 519, 90, 1)
define angle index¶
In [108]:
angle_index = 1
wv_angle_str = f"index:{angle_index}, wv:{wavelength2[angle_index]:.4g}nm, angle:{angle2[angle_index]:.4g}"
print(wv_angle_str)
index:1, wv:550nm, angle:43.32
check chla¶
In [85]:
angle_index = 1
key = 'chla'
data = dataset2[key].values
label = key
title = f"{label}"
plot_l2_product(
data, plot_range=plot_range, label=label, title=title, cmap="jet", log_scale=True
)
check Rrs¶
In [109]:
key = 'Rrs1'
data = dataset2[key].values[:, :, angle_index, 0]
label = key
title = f"{label}, {wv_angle_str}"
plot_l2_product(
data, plot_range=plot_range, label=label, title=title, cmap="jet", log_scale=False
)
In [88]:
key = 'Rrs2'
data = dataset2[key].values[:, :, angle_index, 0]
label = key
title = f"{label}, {wv_angle_str}"
plot_l2_product(
data, plot_range=plot_range, label=label, title=title, cmap="jet", log_scale=False
)
get brdf correction factor: Rrs2/Rrs1¶
In [110]:
key1, key2 = 'Rrs1', 'Rrs2'
data = dataset2[key2].values[:, :, angle_index, 0]/dataset2[key1].values[:, :, angle_index, 0]
label = f"{key2}/{key1}"
title = f"{label}, {wv_angle_str}"
plot_l2_product(
data, plot_range=plot_range, label=label, title=title, cmap="jet", log_scale=False
)
load L1C data¶
In [73]:
datatree1 = xr.open_datatree(path1[0],decode_timedelta=False)
dataset1 = xr.merge(datatree1.to_dict().values())
In [82]:
lat1 = dataset2["latitude"].values
lon1 = dataset2["longitude"].values
plot_range1 = [lon1.min(), lon1.max(), lat1.min(), lat1.max()]
plot_range1, plot_range2
Out[82]:
([np.float32(99.83995), np.float32(130.17134), np.float32(-34.135452), np.float32(-11.234323)], [np.float32(99.83995), np.float32(130.17134), np.float32(-34.135452), np.float32(-11.234323)])
In [99]:
#L1 and L2 share the same lat and lon
ensure_same_latlon(lat1, lon1, lat2, lon2)
YES
check geometry:
- solar_zenith_angle
- solar_azimuth_angle
- sensor_zenith_angle
- sensor_azimuth_angle
In [101]:
dataset1['sensor_zenith_angle'].shape
Out[101]:
(396, 519, 90)
In [103]:
angle_index = 1
wv_angle_str = f"index:{angle_index}, wv:{wavelength2[angle_index]:.4g}nm, angle:{angle2[angle_index]:.4g}"
print(wv_angle_str)
index:1, wv:550nm, angle:43.32
In [112]:
key1 = 'sensor_zenith_angle'
for key1 in ['sensor_zenith_angle', 'sensor_azimuth_angle', 'solar_zenith_angle', 'solar_azimuth_angle']:
data = dataset1[key1].values[:, :, angle_index]
label = f"{key1}"
title = f"{label},{wv_angle_str}"
plot_l2_product(
data, plot_range=plot_range, label=label, title=title, cmap="jet", log_scale=False
)
In [ ]:


