diff --git a/Project.toml b/Project.toml index ab7f75323..b7b8532f1 100644 --- a/Project.toml +++ b/Project.toml @@ -33,6 +33,7 @@ Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" [weakdeps] +ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3" Breeze = "660aa2fb-d4c8-4359-a52c-9c057bc511da" CDSAPI = "8a7b9de3-9c00-473e-88b4-7eccd7ef2fea" CondaPkg = "992eb4ea-22a4-4c89-a5bb-47a3300528ab" @@ -44,6 +45,7 @@ WorldOceanAtlasTools = "04f20302-f1b9-11e8-29d9-7d841cb0a64a" XESMF = "2e0b0046-e7a1-486f-88de-807ee8ffabe5" [extensions] +NumericalEarthArchGDALExt = "ArchGDAL" NumericalEarthBreezeExt = "Breeze" NumericalEarthCDSAPIExt = "CDSAPI" NumericalEarthCopernicusMarineExt = "CopernicusMarine" @@ -54,6 +56,8 @@ NumericalEarthWOAExt = "WorldOceanAtlasTools" [compat] Adapt = "4" +ArchGDAL = "0.10" + Breeze = "0.4" CDSAPI = "2.2.1" CFTime = "0.1, 0.2" diff --git a/docs/src/Metadata/metadata_overview.md b/docs/src/Metadata/metadata_overview.md index a2ed7a07f..96fa06362 100644 --- a/docs/src/Metadata/metadata_overview.md +++ b/docs/src/Metadata/metadata_overview.md @@ -80,12 +80,18 @@ NumericalEarth currently ships connectors for the following data products: | Dataset | Supported Variables | Documentation Link | |--------------------|-----------------------------------------------------------|-----------------------------------------------------------------------------------------------------| +| **Bathymetry** | | | | `ETOPO2022` | [Supported variables](@ref dataset-etopo2022-vars) | [NOAA ETOPO 2022 overview](https://www.ncei.noaa.gov/products/etopo-global-relief-model) | +| `GEBCO2024` | [Supported variables](@ref dataset-gebco2024-vars) | [GEBCO 2024 overview](https://www.gebco.net/data_and_products/gridded_bathymetry_data/) | +| `IBCSOv2` | [Supported variables](@ref dataset-ibcsov2-vars) | [IBCSO overview](https://ibcso.org/ibcso-2024-annual-release/) | +| `IBCAOv5` | [Supported variables](@ref dataset-ibcaov5-vars) | [IBCAO overview](https://www.gebco.net/data_and_products/gridded_bathymetry_data/arctic_ocean/) | +| **Ocean reanalysis** | | | | `ECCO2Monthly` | [Supported variables](@ref dataset-ecco2monthly-vars) | [ECCO2 documentation](https://ecco.jpl.nasa.gov/products/all/) | | `ECCO2Daily` | [Supported variables](@ref dataset-ecco2daily-vars) | [ECCO2 documentation](https://ecco.jpl.nasa.gov/products/all/) | | `ECCO4Monthly` | [Supported variables](@ref dataset-ecco4monthly-vars) | [ECCO V4r4 product guide](https://ecco-group.org/products-ECCO-V4r4.htm) | | `EN4Monthly` | [Supported variables](@ref dataset-en4monthly-vars) | [Met Office EN4 overview](https://www.metoffice.gov.uk/hadobs/en4/) | | `GLORYSDaily` | [Supported variables](@ref dataset-glorysdaily-vars) | [Copernicus GLORYS product page](https://data.marine.copernicus.eu/product/GLOBAL_MULTIYEAR_PHY_001_030/description) | | `GLORYSMonthly` | [Supported variables](@ref dataset-glorysmonthly-vars) | [Copernicus GLORYS product page](https://data.marine.copernicus.eu/product/GLOBAL_MULTIYEAR_PHY_001_030/description) | +| **Atmospheric forcing** | | | | `RepeatYearJRA55` | [Supported variables](@ref dataset-repeatyearjra55-vars) | [JRA-55 Reanalysis](https://www.data.jma.go.jp/jra/html/JRA-55/index_en.html) | | `MultiYearJRA55` | [Supported variables](@ref dataset-multiyearjra55-vars) | [JRA-55 Reanalysis](https://www.data.jma.go.jp/jra/html/JRA-55/index_en.html) | diff --git a/docs/src/Metadata/supported_variables.md b/docs/src/Metadata/supported_variables.md index 7f7c2cb03..7e6cd7497 100644 --- a/docs/src/Metadata/supported_variables.md +++ b/docs/src/Metadata/supported_variables.md @@ -4,13 +4,19 @@ NumericalEarth currently ships connectors for the following data products: | Dataset | Supported Variables | Documentation Link | |--------------------|-----------------------------------------------------------|-----------------------------------------------------------------------------------------------------| +| **Bathymetry** | | | | `ETOPO2022` | [Supported variables](@ref dataset-etopo2022-vars) | [NOAA ETOPO 2022 overview](https://www.ncei.noaa.gov/products/etopo-global-relief-model) | +| `GEBCO2024` | [Supported variables](@ref dataset-gebco2024-vars) | [GEBCO 2024 overview](https://www.gebco.net/data_and_products/gridded_bathymetry_data/) | +| `IBCSOv2` | [Supported variables](@ref dataset-ibcsov2-vars) | [IBCSO overview](https://ibcso.org/ibcso-2024-annual-release/) | +| `IBCAOv5` | [Supported variables](@ref dataset-ibcaov5-vars) | [IBCAO overview](https://www.gebco.net/data_and_products/gridded_bathymetry_data/arctic_ocean/) | +| **Ocean reanalysis** | | | | `ECCO2Monthly` | [Supported variables](@ref dataset-ecco2monthly-vars) | [ECCO2 documentation](https://ecco.jpl.nasa.gov/products/all/) | | `ECCO2Daily` | [Supported variables](@ref dataset-ecco2daily-vars) | [ECCO2 documentation](https://ecco.jpl.nasa.gov/products/all/) | | `ECCO4Monthly` | [Supported variables](@ref dataset-ecco4monthly-vars) | [ECCO V4r4 product guide](https://ecco-group.org/products-ECCO-V4r4.htm) | | `EN4Monthly` | [Supported variables](@ref dataset-en4monthly-vars) | [Met Office EN4 overview](https://www.metoffice.gov.uk/hadobs/en4/) | | `GLORYSDaily` | [Supported variables](@ref dataset-glorysdaily-vars) | [Copernicus GLORYS product page](https://data.marine.copernicus.eu/product/GLOBAL_MULTIYEAR_PHY_001_030/description) | | `GLORYSMonthly` | [Supported variables](@ref dataset-glorysmonthly-vars) | [Copernicus GLORYS product page](https://data.marine.copernicus.eu/product/GLOBAL_MULTIYEAR_PHY_001_030/description) | +| **Atmospheric forcing** | | | | `RepeatYearJRA55` | [Supported variables](@ref dataset-repeatyearjra55-vars) | [JRA-55 Reanalysis](https://www.data.jma.go.jp/jra/html/JRA-55/index_en.html) | | `MultiYearJRA55` | [Supported variables](@ref dataset-multiyearjra55-vars) | [JRA-55 Reanalysis](https://www.data.jma.go.jp/jra/html/JRA-55/index_en.html) | @@ -114,3 +120,15 @@ NumericalEarth currently ships connectors for the following data products: - `:snow_freshwater_flux` - Precipitation flux from snow/ice (kg m⁻² s⁻¹). - `:river_freshwater_flux` - River discharge flux (kg m⁻² s⁻¹). - `:iceberg_freshwater_flux` - Iceberg calving flux (kg m⁻² s⁻¹). + +## [Supported variables for IBCSOv2](@id dataset-ibcsov2-vars) + +- `:bottom_height` - Southern Ocean bathymetry at 500 m resolution, south of 50°S (m). + +## [Supported variables for GEBCO2024](@id dataset-gebco2024-vars) + +- `:bottom_height` - Global bathymetry and topography at 15 arc-second resolution (m). + +## [Supported variables for IBCAOv5](@id dataset-ibcaov5-vars) + +- `:bottom_height` - Arctic Ocean bathymetry at 100 m resolution, north of 64°N, including Greenland ice sheet surface elevation (m). diff --git a/ext/NumericalEarthArchGDALExt.jl b/ext/NumericalEarthArchGDALExt.jl new file mode 100644 index 000000000..f8707a5e1 --- /dev/null +++ b/ext/NumericalEarthArchGDALExt.jl @@ -0,0 +1,50 @@ +module NumericalEarthArchGDALExt + +using NumericalEarth +using ArchGDAL +using NCDatasets + +import NumericalEarth.DataWrangling.IBCAO: reproject_ibcao_to_netcdf + +function reproject_ibcao_to_netcdf(tiff_path, nc_path) + ArchGDAL.read(tiff_path) do src + # Warp from EPSG:3996 (Polar Stereographic) to EPSG:4326 (WGS84) + # at 0.01° resolution, clipping to 64–90°N + ArchGDAL.gdalwarp([src], + ["-t_srs", "EPSG:4326", + "-te", "-180", "64", "180", "90", # xmin ymin xmax ymax + "-tr", "0.01", "0.01", # target resolution (degrees) + "-r", "bilinear", # resampling method + "-ot", "Float32"]) do warped + + # ArchGDAL returns data as (Nx, Ny) with y from north to south (GDAL convention) + data = Float32.(ArchGDAL.read(warped, 1)) + data = data[:, end:-1:1] # flip y: j=1 → 64°N (south), j=Ny → 90°N (north) + + Nx, Ny = size(data) # expected: (36000, 2600) + + NCDataset(nc_path, "c") do ds + defDim(ds, "lon", Nx) + defDim(ds, "lat", Ny) + + lon_var = defVar(ds, "lon", Float64, ("lon",); + attrib = ["units" => "degrees_east", + "long_name" => "longitude"]) + lat_var = defVar(ds, "lat", Float64, ("lat",); + attrib = ["units" => "degrees_north", + "long_name" => "latitude"]) + z_var = defVar(ds, "z", Float32, ("lon", "lat"); + attrib = ["long_name" => "elevation", + "units" => "m"]) + + lon_var[:] = range(-180 + 0.005, 180 - 0.005; length=Nx) + lat_var[:] = range(64 + 0.005, 90 - 0.005; length=Ny) + z_var[:, :] = data + end + end + end + + return nothing +end + +end # module NumericalEarthArchGDALExt diff --git a/src/Bathymetry/Bathymetry.jl b/src/Bathymetry/Bathymetry.jl index c290a7f44..9dabad723 100644 --- a/src/Bathymetry/Bathymetry.jl +++ b/src/Bathymetry/Bathymetry.jl @@ -19,7 +19,8 @@ using NCDatasets using Printf using Scratch -using ..DataWrangling: Metadatum, native_grid, metadata_path, download_dataset +using ..DataWrangling: Metadatum, native_grid, metadata_path, download_dataset, + dataset_variable_name, validate_dataset_coverage using ..DataWrangling.ETOPO: ETOPO2022 include("regrid_bathymetry.jl") diff --git a/src/Bathymetry/regrid_bathymetry.jl b/src/Bathymetry/regrid_bathymetry.jl index 430cad426..f3b7de2ab 100644 --- a/src/Bathymetry/regrid_bathymetry.jl +++ b/src/Bathymetry/regrid_bathymetry.jl @@ -186,6 +186,8 @@ function regrid_bathymetry(target_grid, metadata; major_basins = 1, cache = true) + validate_dataset_coverage(target_grid, metadata) + config = BathymetryRegridding(target_grid, metadata; height_above_water, minimum_depth, interpolation_passes, major_basins) @@ -240,7 +242,7 @@ function _regrid_bathymetry(target_grid, metadata; filepath = metadata_path(metadata) dataset = Dataset(filepath, "r") - z_data = convert(Array{FT}, dataset["z"][:, :]) + z_data = convert(Array{FT}, dataset[dataset_variable_name(metadata)][:, :]) close(dataset) if !isnothing(height_above_water) diff --git a/src/DataWrangling/DataWrangling.jl b/src/DataWrangling/DataWrangling.jl index 3741c6e8a..6db3f4853 100644 --- a/src/DataWrangling/DataWrangling.jl +++ b/src/DataWrangling/DataWrangling.jl @@ -5,6 +5,7 @@ restoring, or validation. module DataWrangling export Metadata, Metadatum, DatewiseFilename, ECCOMetadatum, EN4Metadatum, all_dates, first_date, last_date +export validate_dataset_coverage, metadata_filename export BoundingBox, Column, Linear, Nearest export WOAClimatology, WOAAnnual, WOAMonthly export metadata_time_step, metadata_epoch @@ -198,6 +199,29 @@ function binary_data_size end default_mask_value(dataset) = NaN +""" + AbstractStaticDataset + +Supertype for datasets without a time dimension. Provides default no-op implementations for the date-related interface +(`all_dates`, `first_date`, `last_date`). +""" +abstract type AbstractStaticDataset end + +all_dates(::AbstractStaticDataset, args...) = nothing +first_date(::AbstractStaticDataset, args...) = nothing +last_date(::AbstractStaticDataset, args...) = nothing + +""" + AbstractStaticBathymetry <: AbstractStaticDataset + +Supertype for static, two-dimensional bathymetry datasets (e.g. ETOPO, GEBCO, IBCSO, IBCAO). +Adds defaults for the degenerate vertical axis and a variable-agnostic `Base.size`. +""" +abstract type AbstractStaticBathymetry <: AbstractStaticDataset end + +z_interfaces(::AbstractStaticBathymetry) = (0, 1) +Base.size(dataset::AbstractStaticBathymetry, variable) = size(dataset) + # Fundamentals include("metadata.jl") include("metadata_field.jl") @@ -231,6 +255,9 @@ include("ORCA/ORCA.jl") include("WOA/WOA.jl") include("JRA55/JRA55.jl") include("OSPapa/OSPapa.jl") +include("IBCSO/IBCSO.jl") +include("GEBCO/GEBCO.jl") +include("IBCAO/IBCAO.jl") using .ETOPO using .ECCO @@ -241,6 +268,9 @@ using .ORCA using .WOA using .JRA55 using .OSPapa +using .IBCSO +using .GEBCO +using .IBCAO # Fallback: if no download extension is loaded, check that all files already exist function download_dataset(metadata::Metadata) diff --git a/src/DataWrangling/ETOPO/ETOPO.jl b/src/DataWrangling/ETOPO/ETOPO.jl index deca7a45f..f4a47e565 100644 --- a/src/DataWrangling/ETOPO/ETOPO.jl +++ b/src/DataWrangling/ETOPO/ETOPO.jl @@ -7,19 +7,15 @@ using Oceananigans using Oceananigans.DistributedComputations: @root using Scratch -using ..DataWrangling: download_progress, Metadatum, metadata_path +using ..DataWrangling: download_progress, Metadatum, metadata_path, AbstractStaticBathymetry import NumericalEarth.DataWrangling: metadata_filename, default_download_directory, - all_dates, - first_date, - last_date, dataset_variable_name, download_dataset, longitude_interfaces, latitude_interfaces, - z_interfaces, reversed_vertical_axis download_ETOPO_cache::String = "" @@ -31,18 +27,13 @@ ETOPO_bathymetry_variable_names = Dict( :bottom_height => "z", ) -struct ETOPO2022 end +struct ETOPO2022 <: AbstractStaticBathymetry end default_download_directory(::ETOPO2022) = download_ETOPO_cache reversed_vertical_axis(::ETOPO2022) = true longitude_interfaces(::ETOPO2022) = (-180, 180) latitude_interfaces(::ETOPO2022) = (-90, 90) Base.size(::ETOPO2022) = (21600, 10800, 1) -Base.size(dataset::ETOPO2022, variable) = size(dataset) - -all_dates(::ETOPO2022, args...) = nothing -first_date(::ETOPO2022, args...) = nothing -last_date(::ETOPO2022, args...) = nothing const ETOPOMetadatum = Metadatum{<:ETOPO2022} @@ -51,7 +42,6 @@ dataset_variable_name(data::ETOPOMetadatum) = ETOPO_bathymetry_variable_names[da const ETOPO_url = "https://www.dropbox.com/scl/fi/6pwalcuuzgtpanysn4h6f/" * "ETOPO_2022_v1_60s_N90W180_surface.nc?rlkey=2t7890ruyk4nd5t5eov5768lt&st=yfxsy1lu&dl=0" -z_interfaces(::ETOPOMetadatum) = (0, 1) metadata_url(::ETOPOMetadatum) = ETOPO_url metadata_filename(::ETOPO2022, name, date, region) = "ETOPO_2022_v1_60s_N90W180_surface.nc" diff --git a/src/DataWrangling/GEBCO/GEBCO.jl b/src/DataWrangling/GEBCO/GEBCO.jl new file mode 100644 index 000000000..58ba7e796 --- /dev/null +++ b/src/DataWrangling/GEBCO/GEBCO.jl @@ -0,0 +1,125 @@ +module GEBCO + +export GEBCO2024 + +using Downloads +using ZipFile +using Oceananigans +using Oceananigans.DistributedComputations: @root +using Scratch + +using ..DataWrangling: download_progress, Metadatum, metadata_path, AbstractStaticBathymetry + +import NumericalEarth.DataWrangling: + metadata_filename, + default_download_directory, + dataset_variable_name, + download_dataset, + longitude_interfaces, + latitude_interfaces, + reversed_vertical_axis + +download_GEBCO_cache::String = "" +function __init__() + global download_GEBCO_cache = @get_scratch!("GEBCO") +end + +GEBCO_bathymetry_variable_names = Dict( + :bottom_height => "elevation", # Variable name in NetCDF +) + +""" + GEBCO2024 + +General Bathymetric Chart of the Oceans 2024 release. +Global bathymetry and topography at 15 arc-second resolution. + +The GEBCO_2024 Grid is a global terrain model for ocean and land, +providing elevation data on a 15 arc-second interval grid. + +Reference: GEBCO Compilation Group (2024) GEBCO 2024 Grid +Data source: https://www.gebco.net/data_and_products/gridded_bathymetry_data/ +""" +struct GEBCO2024 <: AbstractStaticBathymetry end + +default_download_directory(::GEBCO2024) = download_GEBCO_cache +reversed_vertical_axis(::GEBCO2024) = false + +# GEBCO covers the entire globe +longitude_interfaces(::GEBCO2024) = (-180, 180) +latitude_interfaces(::GEBCO2024) = (-90, 90) + +# Grid size for 15 arc-second resolution +# 360° / (15/3600)° = 86400 points in longitude +# 180° / (15/3600)° = 43200 points in latitude +Base.size(::GEBCO2024) = (86400, 43200, 1) + +const GEBCOMetadatum = Metadatum{<:GEBCO2024} + +dataset_variable_name(data::GEBCOMetadatum) = GEBCO_bathymetry_variable_names[data.name] + +# GEBCO 2024 download URL from BODC +# Note: This is a large file (~8 GB zipped, ~22 GB unzipped) +const GEBCO_zip_url = "https://www.bodc.ac.uk/data/open_download/gebco/gebco_2024/zip/" + +# The expected NetCDF filename inside the ZIP +const GEBCO_nc_filename = "GEBCO_2024.nc" +metadata_filename(::GEBCO2024, name, date, bounding_box) = GEBCO_nc_filename + +function download_dataset(metadatum::GEBCOMetadatum) + filepath = metadata_path(metadatum) + download_dir = metadatum.dir + + @root if !isfile(filepath) + @info "Downloading GEBCO data: $(metadatum.name) to $download_dir..." + @info "Note: GEBCO is a large dataset (~8 GB download, ~22 GB uncompressed). This may take a while." + + # Download the ZIP file + zip_path = joinpath(download_dir, "GEBCO_2024.zip") + + try + @info "Downloading from BODC..." + Downloads.download(GEBCO_zip_url, zip_path; progress=download_progress) + + # Extract the NetCDF file from the ZIP using ZipFile.jl + @info "Extracting NetCDF from ZIP archive..." + zf = ZipFile.Reader(zip_path) + extracted = false + for f in zf.files + if endswith(f.name, GEBCO_nc_filename) + open(filepath, "w") do io + write(io, read(f)) + end + extracted = true + break + end + end + close(zf) + + if !extracted + error("Could not find $GEBCO_nc_filename in ZIP archive") + end + + if isfile(filepath) + @info "GEBCO data extracted successfully" + else + error("Failed to extract GEBCO NetCDF file") + end + + # Clean up ZIP file to save space + rm(zip_path; force=true) + + catch e + @warn "Failed to download GEBCO: $e" + + # Clean up any partial download + rm(zip_path; force=true) + + rethrow(e) + end + end + + return filepath +end + +end # module diff --git a/src/DataWrangling/IBCAO/IBCAO.jl b/src/DataWrangling/IBCAO/IBCAO.jl new file mode 100644 index 000000000..f26d1a4b2 --- /dev/null +++ b/src/DataWrangling/IBCAO/IBCAO.jl @@ -0,0 +1,102 @@ +module IBCAO + +export IBCAOv5 + +using Downloads +using Oceananigans +using Oceananigans.DistributedComputations: @root +using Scratch +using NCDatasets + +using ..DataWrangling: download_progress, Metadatum, metadata_path, AbstractStaticBathymetry + +import NumericalEarth.DataWrangling: + metadata_filename, + default_download_directory, + dataset_variable_name, + download_dataset, + longitude_interfaces, + latitude_interfaces, + reversed_vertical_axis, + validate_dataset_coverage + + +download_IBCAO_cache::String = "" +function __init__() + global download_IBCAO_cache = @get_scratch!("IBCAO") +end + +IBCAO_bathymetry_variable_names = Dict( + :bottom_height => "z", +) + +""" + IBCAOv5 + +International Bathymetric Chart of the Arctic Ocean Version 5.1 (2025). +100m resolution bathymetry for the Arctic Ocean (north of 64°N), including +Greenland ice sheet surface elevation. Data is provided in Polar Stereographic +projection (EPSG:3996) and reprojected to WGS84 geographic coordinates at 0.01° +resolution on first use. + +Reference: Jakobsson et al. (2024), https://doi.org/10.1038/s41597-024-04278-w +Data source: https://www.gebco.net/data-products/gridded-bathymetry-data/arctic-ocean/ +""" +struct IBCAOv5 <: AbstractStaticBathymetry end + +default_download_directory(::IBCAOv5) = download_IBCAO_cache +reversed_vertical_axis(::IBCAOv5) = false + +# Geographic bounds after reprojection to WGS84 +longitude_interfaces(::IBCAOv5) = (-180, 180) +latitude_interfaces(::IBCAOv5) = (64, 90) + +# 0.01° resolution: 36000 × 2600 +Base.size(::IBCAOv5) = (36000, 2600, 1) + +const IBCAOMetadatum = Metadatum{<:IBCAOv5} + +dataset_variable_name(data::IBCAOMetadatum) = IBCAO_bathymetry_variable_names[data.name] + +# CEDA BODC direct download — 100m, with Greenland ice sheet elevation (~25 GB) +const IBCAO_tiff_url = "https://dap.ceda.ac.uk/bodc/gebco/ibcao/ibcao_v5.1/" * + "greenland_ice_sheet_elevation_data/100mx100m_grid_cell_spacing/" * + "single_complete_bathymetric_grid/ibcao_5_1_2025_ice_100m.tiff?download=1" + +const IBCAO_tiff_filename = "ibcao_5_1_2025_ice_100m.tiff" +const IBCAO_nc_filename = "ibcao_v5_wgs84_0p01deg.nc" + +metadata_filename(::IBCAOv5, name, date, bounding_box) = IBCAO_nc_filename + +function validate_dataset_coverage(grid, ::IBCAOMetadatum) + φ_south, _ = Oceananigans.Grids.y_domain(grid) + if φ_south < 64 + error("IBCAOv5 only covers the Arctic Ocean (north of 64°N). " * + "The grid extends to $(round(φ_south, digits=1))°N. " * + "Use ETOPO2022() or GEBCO2024() for domains that extend south of 64°N.") + end +end + +function download_dataset(metadatum::IBCAOMetadatum) + nc_path = metadata_path(metadatum) + tiff_path = joinpath(metadatum.dir, IBCAO_tiff_filename) + + @root if !isfile(nc_path) + if !isfile(tiff_path) + @info "Downloading IBCAO V5.1 GeoTIFF (100m, with Greenland ice, ~25 GB)..." + Downloads.download(IBCAO_tiff_url, tiff_path; progress=download_progress) + end + + @info "Reprojecting IBCAO from Polar Stereographic (EPSG:3996) to WGS84 at 0.01°..." + reproject_ibcao_to_netcdf(tiff_path, nc_path) + @info "Reprojection complete. Removing raw GeoTIFF to save disk space..." + rm(tiff_path; force=true) + end + + return nc_path +end + +# Implemented in ext/NumericalEarthArchGDALExt.jl when ArchGDAL is loaded. +function reproject_ibcao_to_netcdf end + +end # module diff --git a/src/DataWrangling/IBCSO/IBCSO.jl b/src/DataWrangling/IBCSO/IBCSO.jl new file mode 100644 index 000000000..080855db8 --- /dev/null +++ b/src/DataWrangling/IBCSO/IBCSO.jl @@ -0,0 +1,87 @@ +module IBCSO + +export IBCSOv2 + +using Downloads +using Oceananigans +using Oceananigans.DistributedComputations: @root +using Scratch + +using ..DataWrangling: download_progress, Metadatum, metadata_path, AbstractStaticBathymetry + +import NumericalEarth.DataWrangling: + metadata_filename, + default_download_directory, + dataset_variable_name, + download_dataset, + longitude_interfaces, + latitude_interfaces, + reversed_vertical_axis, + validate_dataset_coverage + + +download_IBCSO_cache::String = "" +function __init__() + global download_IBCSO_cache = @get_scratch!("IBCSO") +end + +IBCSO_bathymetry_variable_names = Dict( + :bottom_height => "z", # Variable name in NetCDF +) + +""" + IBCSOv2 + +International Bathymetric Chart of the Southern Ocean Version 2 (2024 Annual Release). +High-resolution (500m) bathymetry for the Southern Ocean (south of 50°S). + +Reference: Dorschel et al. (2022), https://doi.org/10.1594/PANGAEA.937574 +Data source: https://ibcso.org/ibcso-2024-annual-release/ +""" +struct IBCSOv2 <: AbstractStaticBathymetry end + +default_download_directory(::IBCSOv2) = download_IBCSO_cache +reversed_vertical_axis(::IBCSOv2) = false + +# WGS84 version covers -180 to 180 longitude, -90 to -50 latitude +longitude_interfaces(::IBCSOv2) = (-180, 180) +latitude_interfaces(::IBCSOv2) = (-90, -50) + +# Grid size for WGS84 version (500m resolution) +# lon: 33812, lat: 3757 (from -180 to 180, -90 to -50) +Base.size(::IBCSOv2) = (33812, 3757, 1) + +const IBCSOMetadatum = Metadatum{<:IBCSOv2} + +dataset_variable_name(data::IBCSOMetadatum) = IBCSO_bathymetry_variable_names[data.name] + +const IBCSO_pangaea_url = "https://download.pangaea.de/dataset/937574/files/IBCSO_v2_bed_WGS84.nc" + +metadata_url(::IBCSOMetadatum) = IBCSO_pangaea_url + +# The expected NetCDF filename inside the ZIP or from PANGAEA +const IBCSO_nc_filename = "IBCSO_v2_bed_WGS84.nc" +metadata_filename(::IBCSOv2, name, date, bounding_box) = IBCSO_nc_filename + +function validate_dataset_coverage(grid, ::IBCSOMetadatum) + _, φ_north = Oceananigans.Grids.y_domain(grid) + if φ_north > -50 + error("IBCSOv2 only covers the Southern Ocean (south of 50°S). " * + "The grid extends to $(round(φ_north, digits=1))°. " * + "Use ETOPO2022() or GEBCO2024() for domains that include latitudes north of 50°S.") + end +end + +function download_dataset(metadatum::IBCSOMetadatum) + filepath = metadata_path(metadatum) + download_dir = metadatum.dir + + @root if !isfile(filepath) + @info "Downloading IBCSO data: $(metadatum.name) to $download_dir..." + Downloads.download(IBCSO_pangaea_url, filepath; progress=download_progress) + end + + return filepath +end + +end # module diff --git a/src/DataWrangling/metadata.jl b/src/DataWrangling/metadata.jl index a6f3fb092..4f973fd59 100644 --- a/src/DataWrangling/metadata.jl +++ b/src/DataWrangling/metadata.jl @@ -331,6 +331,16 @@ Return the name used for the variable `metadata.name` in its raw dataset file. """ function dataset_variable_name end +""" + validate_dataset_coverage(grid, metadata) + +Check that `grid` lies within the spatial coverage of `metadata`'s dataset. +Throws an error if the grid extends outside the dataset's domain. +The default implementation does nothing (all grids are accepted). +Dataset-specific methods can override this to enforce coverage constraints. +""" +validate_dataset_coverage(grid, metadata) = nothing + """ dataset_location(dataset, variable_name) diff --git a/test/Project.toml b/test/Project.toml index 54e377ccc..bfb90f138 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,4 +1,5 @@ [deps] +ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3" Breeze = "660aa2fb-d4c8-4359-a52c-9c057bc511da" CDSAPI = "8a7b9de3-9c00-473e-88b4-7eccd7ef2fea" CFTime = "179af706-886a-5703-950a-314cd64e0468" @@ -32,6 +33,7 @@ NumericalEarth = {path = ".."} Breeze = {url = "https://github.com/NumericalEarth/Breeze.jl/", rev = "main"} [compat] +ArchGDAL = "0.10" Breeze = "0.4" CDSAPI = "2.2.2" CFTime = "0.1, 0.2" diff --git a/test/test_polar_bathymetry.jl b/test/test_polar_bathymetry.jl new file mode 100644 index 000000000..f1b24b281 --- /dev/null +++ b/test/test_polar_bathymetry.jl @@ -0,0 +1,136 @@ +include("runtests_setup.jl") + +using NumericalEarth.DataWrangling.IBCSO +using NumericalEarth.DataWrangling.GEBCO +using NumericalEarth.DataWrangling.IBCAO +using NumericalEarth.DataWrangling: longitude_interfaces, latitude_interfaces, z_interfaces, + dataset_variable_name, validate_dataset_coverage, + metadata_filename +using NumericalEarth.Bathymetry: regrid_bathymetry + +@testset "Polar bathymetry metadata interfaces" begin + + @testset "IBCSOv2 metadata" begin + ds = IBCSOv2() + @test longitude_interfaces(ds) == (-180, 180) + @test latitude_interfaces(ds) == (-90, -50) + @test z_interfaces(ds) == (0, 1) + @test size(ds) == (33812, 3757, 1) + + meta = Metadatum(:bottom_height, dataset=ds) + @test dataset_variable_name(meta) == "z" + @test endswith(metadata_filename(ds, :bottom_height, nothing, nothing), ".nc") + end + + @testset "IBCSOv2 coverage validation" begin + meta = Metadatum(:bottom_height, dataset=IBCSOv2()) + + # Grid that extends north of -50°S should throw + grid_bad = LatitudeLongitudeGrid(CPU(); + size = (10, 10, 1), + longitude = (0, 360), + latitude = (-60, 0), + z = (-1, 0)) + @test_throws ErrorException validate_dataset_coverage(grid_bad, meta) + + # Grid entirely south of -50°S should pass + grid_ok = LatitudeLongitudeGrid(CPU(); + size = (10, 10, 1), + longitude = (0, 360), + latitude = (-90, -55), + z = (-1, 0)) + @test validate_dataset_coverage(grid_ok, meta) === nothing + end + + @testset "GEBCO2024 metadata" begin + ds = GEBCO2024() + @test longitude_interfaces(ds) == (-180, 180) + @test latitude_interfaces(ds) == (-90, 90) + @test z_interfaces(ds) == (0, 1) + Nx, Ny, Nz = size(ds) + @test Nx == 86400 # 360° at 15 arc-second + @test Ny == 43200 # 180° at 15 arc-second + @test Nz == 1 + + meta = Metadatum(:bottom_height, dataset=ds) + @test dataset_variable_name(meta) == "elevation" + @test endswith(metadata_filename(ds, :bottom_height, nothing, nothing), ".nc") + end + + @testset "IBCAOv5 metadata" begin + ds = IBCAOv5() + @test longitude_interfaces(ds) == (-180, 180) + @test latitude_interfaces(ds) == (64, 90) + @test z_interfaces(ds) == (0, 1) + @test size(ds) == (36000, 2600, 1) + + meta = Metadatum(:bottom_height, dataset=ds) + @test dataset_variable_name(meta) == "z" + @test endswith(metadata_filename(ds, :bottom_height, nothing, nothing), ".nc") + end + + @testset "IBCAOv5 coverage validation" begin + meta = Metadatum(:bottom_height, dataset=IBCAOv5()) + + # Grid that extends south of 64°N should throw + grid_bad = LatitudeLongitudeGrid(CPU(); + size = (10, 10, 1), + longitude = (-20, 20), + latitude = (50, 80), + z = (-1, 0)) + @test_throws ErrorException validate_dataset_coverage(grid_bad, meta) + + # Grid entirely north of 64°N should pass + grid_ok = LatitudeLongitudeGrid(CPU(); + size = (10, 10, 1), + longitude = (-20, 20), + latitude = (70, 85), + z = (-1, 0)) + @test validate_dataset_coverage(grid_ok, meta) === nothing + end + +end + +@testset "validate_dataset_coverage wired into regrid_bathymetry" begin + # Out-of-range grid should throw before any download occurs + meta = Metadatum(:bottom_height, dataset=IBCSOv2()) + grid_bad = LatitudeLongitudeGrid(CPU(); + size = (10, 10, 1), + longitude = (0, 360), + latitude = (-60, 0), + z = (-1, 0)) + @test_throws ErrorException regrid_bathymetry(grid_bad, meta) + + meta = Metadatum(:bottom_height, dataset=IBCAOv5()) + grid_bad = LatitudeLongitudeGrid(CPU(); + size = (10, 10, 1), + longitude = (-20, 20), + latitude = (50, 80), + z = (-1, 0)) + @test_throws ErrorException regrid_bathymetry(grid_bad, meta) +end + +@testset "IBCSO regridding" begin + @info "Testing IBCSO regridding (downloads ~1.5 GB on first run)..." + + # Drake Passage: open deep ocean well within IBCSO coverage + grid = LatitudeLongitudeGrid(CPU(); + size = (20, 20, 1), + longitude = (-70, -60), + latitude = (-60, -55), + z = (-6000, 0)) + + meta = Metadatum(:bottom_height, dataset=IBCSOv2()) + bathy = regrid_bathymetry(grid, meta; cache=false, height_above_water=0) + z = interior(bathy, :, :, 1) + + # All values should be finite (no NaN or Inf from interpolation gaps) + @test all(isfinite, z) + + # With height_above_water=0 all land cells are capped at 0 + @test maximum(z) ≤ 0 + + # Realistic ocean depths: deeper than 500 m, shallower than the deepest ocean + @test minimum(z) > -12000 + @test minimum(z) < -500 +end