xkykai · Jul 25, 2026 · Jul 25, 2026 · Jul 25, 2026 · Jul 28, 2026 · Jul 28, 2026
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
@@ -39,6 +39,8 @@ env:
   COPERNICUSMARINE_SERVICE_PASSWORD: ${{ secrets.COPERNICUSMARINE_SERVICE_PASSWORD }}
   ECCO_USERNAME: ${{ secrets.ECCO_USERNAME }}
   ECCO_WEBDAV_PASSWORD: ${{ secrets.ECCO_WEBDAV_PASSWORD }}
+  EARTHDATA_USERNAME: ${{ secrets.EARTHDATA_USERNAME }}
+  EARTHDATA_PASSWORD: ${{ secrets.EARTHDATA_PASSWORD }}
   JULIA_CUDA_USE_COMPAT: false
 
 jobs:

diff --git a/Project.toml b/Project.toml
@@ -107,7 +107,7 @@ NCDatasets = "0.12, 0.13, 0.14"
 NaturalEarth = "0.1"
 NetworkOptions = "1"
 OceanBioME = "0.19.1, 0.20"
-Oceananigans = "0.111, 0.112, 0.113, 0.114"
+Oceananigans = "0.112, 0.113, 0.114"
 OffsetArrays = "1.14"
 PrecompileTools = "1"
 PythonCall = "0.9.28"

diff --git a/docs/make.jl b/docs/make.jl
@@ -53,6 +53,8 @@ examples = [
     # Example("Differentiable ERA5-forced slab land", "era5_forced_slab_land"; build_always=false, gpu=false),
     Example("ERA5-forced Terrarium land", "era5_forced_terrarium_land"; build_always=false, gpu=false),
     Example("ERA5 downscaling with Breeze", "breeze_downscaling_era5"; build_always=true, gpu=true),
+    Example("MODIS leaf area index", "modis_leaf_area_index"; build_always=false, gpu=false),
+    Example("MODIS land cover and gap filling", "modis_landcover_gap_fill"; build_always=false, gpu=false),
 ]
 
 # Developer examples from docs/src/developers/ directory

diff --git a/docs/src/Metadata/metadata_overview.md b/docs/src/Metadata/metadata_overview.md
@@ -217,4 +217,6 @@ NumericalEarth currently ships connectors for the following data products:
 | `GHSBuiltS`        | `:built_up_fraction` — plan-area built-up fraction (10 m or 100 m) over a regional window; reprojected from GHSL World-Mollweide (needs `ArchGDAL`) | [GHSL GHS-BUILT-S](https://human-settlement.emergency.copernicus.eu) |
 | `GHSBuilt`         | `:built_up_fraction` and `:building_volume` (m³ per m²): the GHS-BUILT-S fraction and its product with the epoch-2018 ANBH height, paired on the built-up surface pixels over a regional window; reprojected from GHSL World-Mollweide (needs `ArchGDAL`) | [GHSL GHS-BUILT-S and GHS-BUILT-H](https://human-settlement.emergency.copernicus.eu) |
 | `ESAWorldCover`    | `:landcover_class`, `:vegetation_fraction`, and a per-class `:<class>_fraction` — 10 m land-cover classification aggregated block-wise over a regional window (needs `ArchGDAL`) | [ESA WorldCover](https://esa-worldcover.org/en/data-access) |
-
+| `MCD15A2H` | `:leaf_area_index`, `:fpar`, `:leaf_area_index_uncertainty`, `:landcover_code` — 8-day 500 m MODIS composites (needs `ArchGDAL`) | [MCD15A2H V061](https://www.earthdata.nasa.gov/data/catalog/lpcloud-mcd15a2h-061) |
+| `MODISLAIClimatology` | 46-period seasonal climatology of an [`MCD15A2H`](@ref) variable | [MCD15A2H V061](https://www.earthdata.nasa.gov/data/catalog/lpcloud-mcd15a2h-061) |
+| `MCD12Q1` | `:landcover_class`, `:quality_flag`, `:land_water_mask` — annual 500 m land cover on [`MCD15A2H`](@ref)'s lattice | [MCD12Q1 V061](https://www.earthdata.nasa.gov/data/catalog/lpcloud-mcd12q1-061) |
diff --git a/docs/src/appendix/notation.md b/docs/src/appendix/notation.md
@@ -220,6 +220,12 @@ Symbols introduced by [`VariablySaturatedHydrology`](@ref),
 | ``\alpha`` | `α` | albedo | Surface reflectivity (–) |
 | ``\epsilon`` | `ϵ` | emissivity | Surface emissivity (–) |
 
+## Vegetation properties
+
+| Math | Code | Property | Description |
+|:----:|:----:|:---------|:------------|
+| ``𝒜`` | `𝒜` | leaf area index | One-sided leaf area per unit ground area (m² m⁻²) |
+
 ## Similarity theory / surface layer
 
 | Math | Code | Property | Description |
@@ -366,6 +372,7 @@ Most symbols can be entered in the Julia REPL and in editors with Julia support
 | `𝒬` | `\scrQ` | Script Q (heat flux) |
 | `ℐ` | `\scrI` | Script I (radiative intensity) |
 | `ℒ` | `\scrL` | Script L (latent heat) |
+| `𝒜` | `\scrA` | Script A (area index) |
 | `𝓃` | `\scrn` | Script n (pore-size uniformity) |
 | `𝓂` | `\scrm` | Script m (van Genuchten `1 - 1/𝓃`) |
 | `τ` | `\tau` | Tau (kinematic stress) |

diff --git a/examples/modis_landcover_gap_fill.jl b/examples/modis_landcover_gap_fill.jl
diff --git a/examples/modis_leaf_area_index.jl b/examples/modis_leaf_area_index.jl
diff --git a/ext/NumericalEarthArchGDALExt/NumericalEarthArchGDALExt.jl b/ext/NumericalEarthArchGDALExt/NumericalEarthArchGDALExt.jl
@@ -32,6 +32,8 @@ using NumericalEarth.DataWrangling.GHSL: GHSBuilt, GHSBuiltH, GHSBuiltS, GHSLMet
                                          ghsl_regional_raster, built_surface_to_fraction,
                                          mask_building_height, ghsl_pixel_window, download_GHSL_cache,
                                          longitude_latitude_to_mollweide, MOLLWEIDE_RADIUS
+using NumericalEarth.DataWrangling.MODISLand: MODISLand, granule_urls, regional_lattice,
+                                              stored_granule_layers, write_lattice_file
 using NumericalEarth.DataWrangling.OpenLandMap: assemble_raster_window, raster_window_indices,
                                                 validate_wgs84_longitude_latitude, validate_geographic_northup
 using NumericalEarth.DataWrangling.WorldCover: ESAWorldCoverMetadatum, version_year, version_string,
@@ -48,6 +50,7 @@ include("ibcao.jl")
 include("asterged.jl")
 include("globfp3d.jl")
 include("ghsl.jl")
+include("modis_land.jl")
 include("openlandmap.jl")
 include("worldcover.jl")
 

diff --git a/ext/NumericalEarthArchGDALExt/modis_land.jl b/ext/NumericalEarthArchGDALExt/modis_land.jl
@@ -0,0 +1,73 @@
+#####
+##### MODIS land products — Earthdata granule fetch + sinusoidal HDF-EOS reprojection
+#####
+##### A region spans several sinusoidal tiles, so every granule covering it is fetched and
+##### the layers are mosaicked and warped to the region's latitude-longitude window in one
+##### `gdalwarp` call each. The digital numbers are written through unchanged; masking,
+##### quality screening, and scaling happen on read, where the product's decode rules live.
+#####
+
+# Locate the `HDF4_EOS:EOS_GRID:"…":<grid>:<layer>` subdataset without hardcoding the
+# product's grid-group name.
+function modis_subdataset(granule_path, layer)
+    ArchGDAL.read(granule_path) do ds
+        for entry in ArchGDAL.metadata(ds; domain = "SUBDATASETS")
+            occursin("_NAME=", entry) || continue
+            name = split(entry, "_NAME="; limit = 2)[2]
+            endswith(name, ":" * layer) && return name
+        end
+        error("Layer $(layer) not found among the HDF-EOS subdatasets of $(granule_path).")
+    end
+end
+
+# Mosaic one layer of `granule_paths` onto the regional lattice. Nearest-neighbor
+# resampling throughout: every layer is either a bit-packed quality byte or a digital number
+# whose out-of-range codes carry meaning, and averaging either would invent values. A lattice
+# cell is 463 m north-south but 463 cos φ m east-west, so away from the equator each 463 m
+# sinusoidal pixel fills about 1 / cos φ neighboring cells. Uncovered cells take the fill
+# code 255, which the read path rejects on both the digital numbers and the quality bytes.
+function warp_modis_layer(granule_paths, layer, lattice)
+    sources = ArchGDAL.IDataset[]
+
+    warped = try
+        for path in granule_paths
+            push!(sources, ArchGDAL.read(modis_subdataset(path, layer)))
+        end
+
+        ArchGDAL.gdalwarp(sources,
+            ["-t_srs", "EPSG:4326",
+             "-te", string(lattice.west), string(lattice.south),
+                    string(lattice.east), string(lattice.north),
+             "-ts", string(lattice.Nx), string(lattice.Ny),
+             "-r", "near",
+             "-ot", "Byte",
+             "-dstnodata", "255"]) do destination
+            ArchGDAL.read(destination, 1)
+        end
+    finally
+        foreach(ArchGDAL.destroy, sources)
+    end
+
+    # GDAL writes rows north→south; flip to the ascending latitude axis of the stored file.
+    return reverse(warped, dims = 2)
+end
+
+function MODISLand.modis_granules_to_netcdf(metadatum::MODISLand.MODISLandMetadatum, nc_path)
+    lattice = regional_lattice(metadatum)
+    urls = granule_urls(metadatum)
+
+    # A sinusoidal tile spans 10° and carries every layer of the product, so the same
+    # granules serve neighboring regions, sibling variables, and every year a climatology
+    # composites. Keeping them beside the warped files makes those warps a local operation.
+    granule_cache = mkpath(joinpath(dirname(nc_path), "granules"))
+    fetched = count(url -> !isfile(joinpath(granule_cache, basename(url))), urls)
+    fetched > 0 && @info string("Fetching ", fetched, " of ", length(urls),
+                                " MODIS granules for ", metadatum.dates, "...")
+
+    granule_paths = [earthdata_download_cached(url, granule_cache) for url in urls]
+    layers = Dict(layer => warp_modis_layer(granule_paths, layer, lattice)
+                  for layer in stored_granule_layers(metadatum.dataset))
+    write_lattice_file(nc_path, lattice, layers)
+
+    return nothing
+end
diff --git a/src/DataWrangling/CopernicusLandAlbedo/CopernicusLandAlbedo.jl b/src/DataWrangling/CopernicusLandAlbedo/CopernicusLandAlbedo.jl
@@ -2,15 +2,16 @@ module CopernicusLandAlbedo
 
 export CopernicusAlbedo, CopernicusAlbedoClimatology, build_monthly_climatology!
 
-using Dates: Dates, DateTime, Month, year, month, daysinmonth
+using Dates: Dates, DateTime, Day, Month, year, month, daysinmonth
 using DocStringExtensions: TYPEDSIGNATURES
 using Downloads: Downloads
 using NCDatasets: NCDataset, defVar, nomissing
 using Oceananigans: Center
 using Oceananigans.DistributedComputations: @root
 
 using ..DataWrangling: DataWrangling, Metadata, Metadatum, BoundingBox,
-                       metadata_path, default_download_directory, native_convention_longitude, unzip
+                       metadata_path, default_download_directory,
+                       native_convention_longitude, native_cell_range, write_atomically, unzip
 
 # TODO: move `group_by_calendar_month` and `is_zip` (product-agnostic download helpers) from
 # ERA5 into DataWrangling proper, so dataset modules don't reach across siblings.
@@ -59,6 +60,9 @@ series derived from SPOT/VGT and PROBA-V observations. Provides `:albedo`, the b
 blue-sky albedo blended from the black-sky (`AL_DH_BB`) and white-sky (`AL_BH_BB`)
 broadband albedos with diffuse fraction `diffuse_fraction`.
 
+Each sample uses a trailing 20-day composite ending on its nominal date. Successive
+composites overlap, and the retrieval also uses prior BRDF estimates.
+
 Files are global on a regular 1/112° latitude-longitude grid; build the `Metadata`
 with a lon/lat [`BoundingBox`](@ref) to window a region at read time. Downloads come
 from the C3S `satellite-albedo` catalog entry and require the CDSAPI backend:
@@ -171,6 +175,14 @@ DataWrangling.all_dates(::CopernicusAlbedo, variable) = copernicus_albedo_ten_da
 # 12 climatological months; the year is arbitrary, only the month matters.
 DataWrangling.all_dates(::CopernicusAlbedoClimatology, variable) = [DateTime(2018, m, 1) for m in 1:12]
 
+DataWrangling.is_seasonal_climatology(::CopernicusAlbedoClimatology) = true
+
+# C3S v2 PUGS §2.5: the stamp labels the last observation of a 20-day composite.
+DataWrangling.averaging_window(md::Metadatum{<:CopernicusAlbedo}) =
+    (DateTime(md.dates) - Day(19), DateTime(md.dates) + Day(1))
+
+DataWrangling.window_span(::Metadata{<:CopernicusAlbedo}) = nothing
+
 # Each climatological month averages every ten days falling in it, so it spans the calendar month.
 DataWrangling.averaging_window(md::Metadatum{<:CopernicusAlbedoClimatology}) =
     DataWrangling.calendar_month_window(md)
@@ -254,25 +266,24 @@ source are preserved.
 function repack_albedo_pair(blacksky, whitesky, destination_names, filepath, expected_size)
     blacksky_path, blacksky_name = blacksky
     whitesky_path, whitesky_name = whitesky
-    staging_path = filepath * ".tmp"
-
-    NCDataset(staging_path, "c") do destination
-        NCDataset(blacksky_path) do source
-            λ = nomissing(source[coordinate_name(source, ("lon", "longitude"))][:])
-            φ = nomissing(source[coordinate_name(source, ("lat", "latitude"))][:])
-            (length(λ), length(φ)) == expected_size ||
-                error("CGLS albedo file grid $((length(λ), length(φ))) does not match ",
-                      "the expected $expected_size; the dataset traits need updating.")
-            defVar(destination, "lon", λ, ("lon",))
-            defVar(destination, "lat", φ, ("lat",))
-            copy_packed_variable!(destination, source, blacksky_name, destination_names[1])
-        end
-        NCDataset(whitesky_path) do source
-            copy_packed_variable!(destination, source, whitesky_name, destination_names[2])
+
+    write_atomically(filepath) do staging_path
+        NCDataset(staging_path, "c") do destination
+            NCDataset(blacksky_path) do source
+                λ = nomissing(source[coordinate_name(source, ("lon", "longitude"))][:])
+                φ = nomissing(source[coordinate_name(source, ("lat", "latitude"))][:])
+                (length(λ), length(φ)) == expected_size ||
+                    error("CGLS albedo file grid $((length(λ), length(φ))) does not match ",
+                          "the expected $expected_size; the dataset traits need updating.")
+                defVar(destination, "lon", λ, ("lon",))
+                defVar(destination, "lat", φ, ("lat",))
+                copy_packed_variable!(destination, source, blacksky_name, destination_names[1])
+            end
+            NCDataset(whitesky_path) do source
+                copy_packed_variable!(destination, source, whitesky_name, destination_names[2])
+            end
         end
     end
-
-    mv(staging_path, filepath; force = true)
     return nothing
 end
 
@@ -286,25 +297,10 @@ end
 
 #####
 ##### Reading — blend the black-sky/white-sky pair into the blue-sky array. Regional
-##### reads hyperslab exactly the native-grid cell window off disk with no global
+##### reads take exactly the native-grid cell window off disk with no global
 ##### materialization, so `retrieve_data` and `read_file_coords` must window identically.
 #####
 
-# 1-based native cell range covered by `bbox` on the axis `(left, right)` split into `N`
-# cells — must match `restrict()` in metadata_field.jl (returned count = `i⁺ - i⁻`, so the
-# cell range `(i⁻+1):i⁺` has that same length, pinning the region offset to di = dj = 0).
-function albedo_cell_range(bbox, interfaces, N)
-    left, right = interfaces
-    Δ = (right - left) / N
-    i⁻ = clamp(floor(Int, (bbox[1] - left) / Δ - 1/2), 0, N)
-    i⁺ = clamp(ceil( Int, (bbox[2] - left) / Δ + 1/2), 0, N)
-    if i⁺ ≤ i⁻
-        i⁺ = min(i⁻ + 1, N)
-        i⁻ = max(i⁺ - 1, 0)
-    end
-    return (i⁻ + 1):i⁺
-end
-
 """
     albedo_read_window(metadatum)
 
@@ -326,8 +322,8 @@ function albedo_read_window(metadatum)
     left, right = native_longitude
     span = bbox_longitude[2] - bbox_longitude[1]
     (span == 360 || (bbox_longitude[1] ≥ left && bbox_longitude[2] > right)) && return nothing
-    icols = albedo_cell_range(bbox_longitude, native_longitude, Nx)
-    jrows = albedo_cell_range(region.latitude, DataWrangling.latitude_interfaces(metadatum), Ny)
+    icols = native_cell_range(bbox_longitude, native_longitude, Nx)
+    jrows = native_cell_range(region.latitude, DataWrangling.latitude_interfaces(metadatum), Ny)
 
     return icols, jrows
 end
@@ -357,7 +353,7 @@ function DataWrangling.retrieve_data(metadatum::CopernicusAlbedoMetadatum)
         icols, jrows = window
         _, Ny, _ = size(metadatum.dataset, metadatum.name)
         # The file stores latitude north→south, so the ascending cells `jrows` sit at file
-        # rows (Ny − last + 1):(Ny − first + 1). Only this hyperslab leaves disk.
+        # rows (Ny − last + 1):(Ny − first + 1). Only this window leaves disk.
         file_rows = (Ny - last(jrows) + 1):(Ny - first(jrows) + 1)
         NCDataset(path) do ds
             α_bs = copernicus_albedo_decode.(ds[blacksky_name][icols, file_rows])
@@ -435,35 +431,34 @@ function write_monthly_mean(filepath, source_paths, variable_names, latitude_chu
         nomissing(ds["lon"][:]), nomissing(ds["lat"][:])
     end
     Nx, Ny = length(λ), length(φ)
-    staging_path = filepath * ".tmp"
-
-    NCDataset(staging_path, "c") do destination
-        defVar(destination, "lon", λ, ("lon",))
-        defVar(destination, "lat", φ, ("lat",))
-
-        for variable_name in variable_names
-            monthly_mean = defVar(destination, variable_name, Float32, ("lon", "lat");
-                                  deflatelevel = 2, shuffle = true)
-
-            for j in 1:latitude_chunk:Ny
-                rows = j:min(j + latitude_chunk - 1, Ny)
-                Σα = zeros(Float32, Nx, length(rows))
-                n = zeros(Int32, Nx, length(rows))
-
-                for path in source_paths
-                    NCDataset(path) do ds
-                        α = copernicus_albedo_decode.(ds[variable_name][:, rows])
-                        @. Σα += ifelse(isnan(α), 0f0, α)
-                        @. n += !isnan(α)
+
+    write_atomically(filepath) do staging_path
+        NCDataset(staging_path, "c") do destination
+            defVar(destination, "lon", λ, ("lon",))
+            defVar(destination, "lat", φ, ("lat",))
+
+            for variable_name in variable_names
+                monthly_mean = defVar(destination, variable_name, Float32, ("lon", "lat");
+                                      deflatelevel = 2, shuffle = true)
+
+                for j in 1:latitude_chunk:Ny
+                    rows = j:min(j + latitude_chunk - 1, Ny)
+                    Σα = zeros(Float32, Nx, length(rows))
+                    n = zeros(Int32, Nx, length(rows))
+
+                    for path in source_paths
+                        NCDataset(path) do ds
+                            α = copernicus_albedo_decode.(ds[variable_name][:, rows])
+                            @. Σα += ifelse(isnan(α), 0f0, α)
+                            @. n += !isnan(α)
+                        end
                     end
-                end
 
-                monthly_mean[:, rows] = @. ifelse(n == 0, NaN32, Σα / n)
+                    monthly_mean[:, rows] = @. ifelse(n == 0, NaN32, Σα / n)
+                end
             end
         end
     end
-
-    mv(staging_path, filepath; force = true)
     return nothing
 end
 

diff --git a/src/DataWrangling/CopernicusLandAlbedo/README.md b/src/DataWrangling/CopernicusLandAlbedo/README.md
@@ -54,6 +54,11 @@ on rebuild.
 
 ## Notes
 
+- This adapter requests C3S version 2. Its [Product User Guide, section 2.5](https://dast.copernicus-climate.eu/documents/satellite-albedo/D3.3.4-v2.0_PUGS_CDR_SA_MULTI_SENSOR_v2.0_PRODUCTS_v1.2.pdf)
+  describes a trailing 20-day composite ending on the file's nominal date, with prior
+  BRDF estimates also contributing. Consecutive composites overlap, so `window_bounds`
+  rejects them. Monthly climatologies average the composite estimates assigned to each
+  nominal month; they are not exact calendar-month means of the underlying observations.
 - Ocean and unretrieved pixels are `NaN` (no inpainting by default — this is a land dataset).
 - Downloads are global and shared across regions; regional windowing happens at read
   time via the `BoundingBox` machinery.