diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index c2ba95283..1171dc75a 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -218,7 +218,7 @@ jobs:
               using Pkg;
               Pkg.test(; coverage=true,
                         julia_args=["--check-bounds=yes", "--compiled-modules=yes", "-O0"],
-                        test_args=["--verbose", "test_cds_downloading", "test_downloading"])
+                        test_args=["--verbose", "test_cds_downloading", "test_jra55_ecco_en4_etopo_downloading", "test_glorys_downloading"])
           '
       - uses: julia-actions/julia-processcoverage@v1
       - uses: codecov/codecov-action@v5
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index 849bff9f5..8d583b64b 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -25,6 +25,10 @@ env:
   NUMERICAL_EARTH_LABEL_BUILD_ALL_EXAMPLES: 'build all examples'
   CDSAPI_URL: "https://cds.climate.copernicus.eu/api"
   CDSAPI_KEY: ${{ secrets.CDSAPI_KEY }}
+  COPERNICUS_USERNAME: ${{ secrets.COPERNICUSMARINE_SERVICE_USERNAME }}
+  COPERNICUS_PASSWORD: ${{ secrets.COPERNICUSMARINE_SERVICE_PASSWORD }}
+  ECCO_USERNAME: ${{ secrets.ECCO_USERNAME }}
+  ECCO_WEBDAV_PASSWORD: ${{ secrets.ECCO_WEBDAV_PASSWORD }}
 
 jobs:
   build-docs:
@@ -63,10 +67,6 @@ jobs:
         run: julia --project=docs --color=yes -e 'using Pkg; Pkg.instantiate(verbose=true)'
       - name: Build documentation
         env:
-          ECCO_USERNAME: ${{ secrets.ECCO_USERNAME }}
-          ECCO_WEBDAV_PASSWORD: ${{ secrets.ECCO_WEBDAV_PASSWORD }}
-          COPERNICUS_USERNAME: ${{ secrets.COPERNICUS_SERVICE_USERNAME }}
-          COPERNICUS_PASSWORD: ${{ secrets.COPERNICUS_USERNAME_PASSWORD }}
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
           JULIA_DEBUG: Documenter
           JULIA_SSL_NO_VERIFY: "**"
diff --git a/.gitignore b/.gitignore
index edb030721..4acc5052a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -64,3 +64,6 @@ docs/src/literated/
 .CondaPkg
 CondaPkg.toml
 !docs/CondaPkg.toml
+
+# claude
+.claude
diff --git a/Project.toml b/Project.toml
index 3abcd1994..ca11979ef 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "NumericalEarth"
 uuid = "904d977b-046a-4731-8b86-9235c0d1ef02"
 license = "MIT"
-version = "0.2.3"
+version = "0.3.0"
 authors = ["NumericalEarth contributors"]
 
 [deps]
diff --git a/docs/CondaPkg.toml b/docs/CondaPkg.toml
index 3e90fc46d..5f3149eb6 100644
--- a/docs/CondaPkg.toml
+++ b/docs/CondaPkg.toml
@@ -1,4 +1,6 @@
-# Pin Python to < 3.14 to avoid the free-threaded build (cp314t)
-# which is ABI-incompatible with PythonCall.
+
+[deps.copernicusmarine]
+channel = "conda-forge"
+
 [deps.python]
 version = ">=3.10,<3.14"
diff --git a/docs/Project.toml b/docs/Project.toml
index c9257bca2..82347bf37 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -6,6 +6,7 @@ CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 ClimaSeaIce = "6ba0ff68-24e6-4315-936c-2e99227c95a4"
 CondaPkg = "992eb4ea-22a4-4c89-a5bb-47a3300528ab"
+CopernicusMarine = "cd43e856-93a3-40c8-bc9e-6146cdce14fa"
 DataDeps = "124859b0-ceae-595e-8997-d05f6a7a8dfe"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
diff --git a/docs/src/Metadata/metadata_overview.md b/docs/src/Metadata/metadata_overview.md
index 393aacea6..a2ed7a07f 100644
--- a/docs/src/Metadata/metadata_overview.md
+++ b/docs/src/Metadata/metadata_overview.md
@@ -64,9 +64,11 @@ The key ingredients stored in a [`Metadata`](@ref) or [`Metadatum`](@ref) object
 - the variable name (for example `:temperature` or `:u_velocity`);
 - the dataset (such as `EN4Monthly`, `ECCO2Daily`, or `GLORYSMonthly`);
 - the temporal coverage: either a single timestamp (`Metadatum`) or a range/vector of dates (`Metadata`);
-- an optional [`BoundingBox`](@ref NumericalEarth.DataWrangling.BoundingBox) describing regional subsets in
-  longitude, latitude, or depth;
-- the on-disk `dir`ectory where the dataset are be cached.
+- an optional `region` describing the spatial extent — either a
+  [`BoundingBox`](@ref NumericalEarth.DataWrangling.BoundingBox) for a rectangular sub-domain,
+  a [`Column`](@ref NumericalEarth.DataWrangling.Column) for a single horizontal location, or
+  `nothing` for the full global domain (see [Regions, locations, and FieldTimeSeries](@ref));
+- the on-disk `dir`ectory where the dataset files are cached.
 
 This bookkeeping lets downstream utilities (for example `set!` or `FieldTimeSeries`) request exactly the
 slices of data they need, and it keeps track of where those slices live so we do not redownload
diff --git a/examples/download_glorys_data.jl b/examples/download_glorys_data.jl
index e77228b72..b87052bb7 100644
--- a/examples/download_glorys_data.jl
+++ b/examples/download_glorys_data.jl
@@ -17,20 +17,20 @@ grid = LatitudeLongitudeGrid(arch;
                              latitude  = (35, 55),
                              longitude = (200, 220))
 
-bounding_box = NumericalEarth.DataWrangling.BoundingBox(longitude=(200, 220), latitude=(35, 55))
+region = NumericalEarth.DataWrangling.BoundingBox(longitude=(200, 220), latitude=(35, 55))
 
 # dataset = NumericalEarth.DataWrangling.Copernicus.GLORYSStatic()
-# static_meta = NumericalEarth.DataWrangling.Metadatum(:depth; dataset, bounding_box)
+# static_meta = NumericalEarth.DataWrangling.Metadatum(:depth; dataset, region)
 # coords_path = NumericalEarth.DataWrangling.download_dataset(static_meta)
 # @info "Downloaded coordinates data to $coords_path"
 
-# T_ecco = NumericalEarth.DataWrangling.ECCOMetadatum(:temperature; dataset, bounding_box)
+# T_ecco = NumericalEarth.DataWrangling.ECCOMetadatum(:temperature; dataset, region)
 # T_en4_meta = NumericalEarth.DataWrangling.EN4Metadatum(:temperature)
 # T_en4_path = NumericalEarth.DataWrangling.download_dataset(T_en4_meta)
 # T_en4 = Field(T_en4_meta)
 
 dataset = NumericalEarth.DataWrangling.Copernicus.GLORYSDaily()
-T_meta = NumericalEarth.DataWrangling.Metadatum(:temperature; dataset, bounding_box)
+T_meta = NumericalEarth.DataWrangling.Metadatum(:temperature; dataset, region)
 T_path = NumericalEarth.DataWrangling.download_dataset(T_meta)
 @info "Downloaded temperature data to $T_path"
 T = Field(T_meta, inpainting=nothing)
diff --git a/examples/single_column_os_papa_simulation.jl b/examples/single_column_os_papa_simulation.jl
index 39d8182ca..5c0439c12 100644
--- a/examples/single_column_os_papa_simulation.jl
+++ b/examples/single_column_os_papa_simulation.jl
@@ -14,6 +14,7 @@
 # pkg"add Oceananigans, NumericalEarth, CairoMakie"
 # ```
 
+using CopernicusMarine
 using NumericalEarth
 using Oceananigans
 using Oceananigans: prognostic_fields
@@ -29,7 +30,7 @@ using Printf
 
 # Ocean station papa location
 location_name = "ocean_station_papa"
-λ★, φ★ = 35.1, 50.1
+λ★, φ★ = -144.9, 50.1
 
 grid = RectilinearGrid(size = 200,
                        x = λ★,
@@ -48,10 +49,14 @@ ocean = ocean_simulation(grid; Δt=10minutes, coriolis=FPlane(latitude = φ★))
 
 ocean.model
 
-# We set initial conditions from ECCO4:
+# We set initial conditions from GLORYS, using a `Column` region so that
+# only the single water column at `(λ★, φ★)` is downloaded from the
+# Copernicus Marine Service.
 
-set!(ocean.model, T=Metadatum(:temperature, dataset=ECCO4Monthly()),
-                  S=Metadatum(:salinity, dataset=ECCO4Monthly()))
+col = Column(λ★, φ★; interpolation=Nearest())
+
+set!(ocean.model, T=Metadatum(:temperature, dataset=GLORYSMonthly(), region=col),
+                  S=Metadatum(:salinity,    dataset=GLORYSMonthly(), region=col))
 
 # # A prescribed atmosphere based on JRA55 re-analysis
 #
diff --git a/ext/NumericalEarthCDSAPIExt.jl b/ext/NumericalEarthCDSAPIExt.jl
index bac2aefe9..8232e6c77 100644
--- a/ext/NumericalEarthCDSAPIExt.jl
+++ b/ext/NumericalEarthCDSAPIExt.jl
@@ -77,7 +77,7 @@ function download_dataset(meta::ERA5Metadatum; skip_existing=true)
     )
 
     # Add area constraint from bounding box
-    area = build_era5_area(meta.bounding_box)
+    area = build_era5_area(meta.region)
     if !isnothing(area)
         request["area"] = area
     end
diff --git a/ext/NumericalEarthCopernicusMarineExt.jl b/ext/NumericalEarthCopernicusMarineExt.jl
index ec2013634..7105f8cb0 100644
--- a/ext/NumericalEarthCopernicusMarineExt.jl
+++ b/ext/NumericalEarthCopernicusMarineExt.jl
@@ -46,11 +46,12 @@ function download_dataset(meta::GLORYSMetadatum;
         (; start_datetime, end_datetime)
     end
 
-    lon_kw = longitude_bounds_kw(meta.bounding_box)
-    lat_kw = latitude_bounds_kw(meta.bounding_box)
-    z_kw = depth_bounds_kw(meta.bounding_box)
+    lon_kw = longitude_bounds_kw(meta.region)
+    lat_kw = latitude_bounds_kw(meta.region)
+    z_kw = depth_bounds_kw(meta.region)
+    selection_method = coordinates_selection_method(meta.region)
 
-    kw = (; coordinates_selection_method = "outside",
+    kw = (; coordinates_selection_method = selection_method,
           skip_existing,
           dataset_id,
           variables,
@@ -76,12 +77,38 @@ end
 longitude_bounds_kw(::Nothing) = NamedTuple()
 latitude_bounds_kw(::Nothing) = NamedTuple()
 depth_bounds_kw(::Nothing) = NamedTuple()
+coordinates_selection_method(::Nothing) = "outside"
 
 const BBOX = NumericalEarth.DataWrangling.BoundingBox
+const COL  = NumericalEarth.DataWrangling.Column
+const LIN  = NumericalEarth.DataWrangling.Linear
+const NR   = NumericalEarth.DataWrangling.Nearest
+
+longitude_bounds_kw(bbox::BBOX) = longitude_bounds_kw(bbox.longitude)
+latitude_bounds_kw(bbox::BBOX) = latitude_bounds_kw(bbox.latitude)
+depth_bounds_kw(bbox::BBOX) = depth_bounds_kw(bbox.z)
+coordinates_selection_method(::BBOX) = "outside"
+
+# Column with Nearest interpolation: download the single nearest point
+longitude_bounds_kw(col::COL{<:Any, <:Any, <:Any, NR}) = (; minimum_longitude = col.longitude, maximum_longitude = col.longitude)
+latitude_bounds_kw(col::COL{<:Any, <:Any, <:Any, NR})  = (; minimum_latitude = col.latitude, maximum_latitude = col.latitude)
+depth_bounds_kw(col::COL) = depth_bounds_kw(col.z)
+coordinates_selection_method(::COL{<:Any, <:Any, <:Any, NR}) = "nearest"
+
+# Column with Linear interpolation: expand by a small margin for interpolation
+longitude_bounds_kw(col::COL{<:Any, <:Any, <:Any, LIN}) = expand_longitude(col.longitude)
+latitude_bounds_kw(col::COL{<:Any, <:Any, <:Any, LIN})  = expand_latitude(col.latitude)
+coordinates_selection_method(::COL{<:Any, <:Any, <:Any, LIN}) = "outside"
+
+function expand_longitude(lon)
+    ε = 1/6  # slightly more than one GLORYS grid cell (1/12°)
+    return (; minimum_longitude = lon - ε, maximum_longitude = lon + ε)
+end
 
-longitude_bounds_kw(bounding_box::BBOX) = longitude_bounds_kw(bounding_box.longitude)
-latitude_bounds_kw(bounding_box::BBOX) = latitude_bounds_kw(bounding_box.latitude)
-depth_bounds_kw(bounding_box::BBOX) = depth_bounds_kw(bounding_box.z)
+function expand_latitude(lat)
+    ε = 1/6
+    return (; minimum_latitude = lat - ε, maximum_latitude = lat + ε)
+end
 
 function longitude_bounds_kw(longitude)
     minimum_longitude = longitude[1]
diff --git a/src/DataWrangling/DataWrangling.jl b/src/DataWrangling/DataWrangling.jl
index 0bd3c973b..92fa94cd7 100644
--- a/src/DataWrangling/DataWrangling.jl
+++ b/src/DataWrangling/DataWrangling.jl
@@ -5,11 +5,13 @@ restoring, or validation.
 module DataWrangling
 
 export Metadata, Metadatum, DatewiseFilename, ECCOMetadatum, EN4Metadatum, all_dates, first_date, last_date
+export BoundingBox, Column, Linear, Nearest
 export WOAClimatology, WOAAnnual, WOAMonthly
 export metadata_time_step, metadata_epoch
 export LinearlyTaperedPolarMask
 export DatasetRestoring, SurfaceFluxRestoring
 export ERA5Hourly, ERA5Monthly
+export native_grid
 
 using Oceananigans
 using Downloads
@@ -189,6 +191,7 @@ function z_interfaces end
 function longitude_interfaces end
 function latitude_interfaces end
 function reversed_vertical_axis end
+reversed_latitude_axis(dataset) = false
 function native_grid end
 function binary_data_grid end
 function binary_data_size end
diff --git a/src/DataWrangling/ECCO/ECCO.jl b/src/DataWrangling/ECCO/ECCO.jl
index 7eac889ec..576e81014 100644
--- a/src/DataWrangling/ECCO/ECCO.jl
+++ b/src/DataWrangling/ECCO/ECCO.jl
@@ -21,20 +21,25 @@ using NumericalEarth.DataWrangling:
     netrc_downloader,
     NearestNeighborInpainting,
     BoundingBox,
+    Column,
     metadata_path,
     GramPerKilogramMinus35,
     MicromolePerLiter,
     Metadata,
     Metadatum,
     download_progress,
-    InverseSign
+    InverseSign,
+    native_grid,
+    location,
+    compute_mask,
+    inpaint_mask!,
+    set_metadata_field!,
+    extract_column!
 
 using KernelAbstractions: @kernel, @index
 
 using Dates: year, month, day
 
-import Oceananigans: location
-
 import NumericalEarth.DataWrangling:
     default_download_directory,
     all_dates,
@@ -42,6 +47,7 @@ import NumericalEarth.DataWrangling:
     download_dataset,
     conversion_units,
     dataset_variable_name,
+    dataset_location,
     metaprefix,
     longitude_interfaces,
     latitude_interfaces,
@@ -247,7 +253,7 @@ end
 metaprefix(::ECCOMetadata) = "ECCOMetadata"
 
 # File name generation specific to each dataset
-function metadata_filename(::ECCO4Monthly, name, date, bounding_box)
+function metadata_filename(::ECCO4Monthly, name, date, region)
     shortname = ECCO4_dataset_variable_names[name]
     yearstr   = string(Dates.year(date))
     monthstr  = string(Dates.month(date), pad=2)
@@ -260,7 +266,7 @@ ecco2_is_three_dimensional(name) =
     name == :u_velocity ||
     name == :v_velocity
 
-function metadata_filename(dataset::Union{ECCO2Daily, ECCO2Monthly}, name, date, bounding_box)
+function metadata_filename(dataset::Union{ECCO2Daily, ECCO2Monthly}, name, date, region)
     shortname = ECCO2_dataset_variable_names[name]
     yearstr   = string(Dates.year(date))
     monthstr  = string(Dates.month(date), pad=2)
@@ -278,7 +284,7 @@ end
 dataset_variable_name(data::Metadata{<:ECCO2Daily})   = ECCO2_dataset_variable_names[data.name]
 dataset_variable_name(data::Metadata{<:ECCO2Monthly}) = ECCO2_dataset_variable_names[data.name]
 dataset_variable_name(data::Metadata{<:ECCO4Monthly}) = ECCO4_dataset_variable_names[data.name]
-location(data::ECCOMetadata) = ECCO_location[data.name]
+dataset_location(::ECCODataset, name) = ECCO_location[name]
 
 is_three_dimensional(data::ECCOMetadata) =
     data.name == :temperature ||
@@ -365,4 +371,37 @@ inpainted_metadata_path(metadata::ECCOMetadatum) = joinpath(metadata.dir, inpain
 
 include("ECCO_atmosphere.jl")
 
+#####
+##### Column Field for ECCO datasets (which always download globally)
+#####
+
+using Oceananigans.BoundaryConditions: fill_halo_regions!
+
+const ECCOColumnMetadatum = Metadatum{<:ECCODataset, <:Any, <:Column}
+
+function Oceananigans.Fields.Field(metadata::ECCOColumnMetadatum, arch=CPU();
+                                   inpainting = default_inpainting(metadata),
+                                   mask = nothing,
+                                   halo = (3, 3, 3),
+                                   cache_inpainted_data = true)
+
+    download_dataset(metadata)
+    column_grid = native_grid(metadata, arch; halo)
+
+    # Build a full-grid Field without a region to load the global data
+    global_metadatum = Metadatum(metadata.name;
+                                 dataset = metadata.dataset,
+                                 date = metadata.dates)
+
+    intermediate_field = Field(global_metadatum, arch; inpainting, mask, halo, cache_inpainted_data)
+    fill_halo_regions!(intermediate_field)
+
+    # Extract the column
+    _, _, LZ = location(metadata)
+    column_field = Field{Nothing, Nothing, LZ}(column_grid)
+    extract_column!(column_field, intermediate_field, metadata.region)
+
+    return column_field
+end
+
 end # Module
diff --git a/src/DataWrangling/ECCO/ECCO_darwin.jl b/src/DataWrangling/ECCO/ECCO_darwin.jl
index d672b82a3..4087731b6 100644
--- a/src/DataWrangling/ECCO/ECCO_darwin.jl
+++ b/src/DataWrangling/ECCO/ECCO_darwin.jl
@@ -26,14 +26,14 @@ all_dates(dataset::ECCO2DarwinMonthly, name) = metadata_epoch(dataset) : Month(1
 
 # File name generation specific to each Dataset dataset
 """
-    metadata_filename(dataset, name, date, bounding_box)
+    metadata_filename(dataset, name, date, region)
 
 Generate the filename for a given ECCO Darwin dataset and date.
 
 The filename is constructed using the dataset variable name, and the iteration number is calculated
 from the date and epoch.
 """
-function metadata_filename(dataset::Union{ECCO2DarwinMonthly, ECCO4DarwinMonthly}, name, date, bounding_box)
+function metadata_filename(dataset::Union{ECCO2DarwinMonthly, ECCO4DarwinMonthly}, name, date, region)
     shortname = ECCO_darwin_dataset_variable_names[name]
 
     reference_date = metadata_epoch(dataset)
@@ -52,8 +52,6 @@ default_mask_value(::ECCO2DarwinMonthly) = 0
 
 dataset_variable_name(data::Metadata{<:Union{ECCO2DarwinMonthly,ECCO4DarwinMonthly}}) = ECCO_darwin_dataset_variable_names[data.name]
 
-location(::Metadata{<:Union{ECCO2DarwinMonthly, ECCO4DarwinMonthly}}) = (Center, Center, Center)
-
 variable_is_three_dimensional(::Metadata{<:Union{ECCO2DarwinMonthly, ECCO4DarwinMonthly}}) = true
 
 ECCO_darwin_dataset_variable_names = Dict(
diff --git a/src/DataWrangling/EN4/EN4.jl b/src/DataWrangling/EN4/EN4.jl
index acc9ace30..851401439 100644
--- a/src/DataWrangling/EN4/EN4.jl
+++ b/src/DataWrangling/EN4/EN4.jl
@@ -50,7 +50,6 @@ import NumericalEarth.DataWrangling:
     inpainted_metadata_path,
     available_variables
 
-import Oceananigans.Fields: location
 
 download_EN4_cache::String = ""
 function __init__()
@@ -160,7 +159,7 @@ metaprefix(::EN4Metadatum) = "EN4Metadatum"
 
 # Note, EN4 files contain all variables, so the filenames do not
 # depend on name.
-function metadata_filename(::EN4Monthly, name, date, bounding_box)
+function metadata_filename(::EN4Monthly, name, date, region)
     yearstr  = string(Dates.year(date))
     monthstr = string(Dates.month(date), pad=2)
     return "EN.4.2.2.f.analysis.g10." * yearstr * lpad(string(monthstr), 2, '0') * ".nc"
@@ -168,7 +167,6 @@ end
 
 # Convenience functions
 dataset_variable_name(data::EN4Metadata) = EN4_dataset_variable_names[data.name]
-location(::EN4Metadata) = (Center, Center, Center)
 is_three_dimensional(::EN4Metadata) = true
 
 ## This function is explicitly for the downloader to check if the zip file/extracted file exists,
diff --git a/src/DataWrangling/ERA5/ERA5.jl b/src/DataWrangling/ERA5/ERA5.jl
index 09ba949c9..359a2d50e 100644
--- a/src/DataWrangling/ERA5/ERA5.jl
+++ b/src/DataWrangling/ERA5/ERA5.jl
@@ -11,11 +11,10 @@ using NumericalEarth.DataWrangling: Metadata, Metadatum, metadata_path
 using Dates
 using Dates: DateTime, Day, Month, Hour
 
-import Oceananigans.Fields: location
-
 import NumericalEarth.DataWrangling:
     all_dates,
     dataset_variable_name,
+    dataset_location,
     default_download_directory,
     longitude_interfaces,
     latitude_interfaces,
@@ -24,7 +23,8 @@ import NumericalEarth.DataWrangling:
     inpainted_metadata_path,
     available_variables,
     retrieve_data,
-    metadata_path
+    metadata_path,
+    reversed_latitude_axis
 
 import Base: eltype
 
@@ -48,6 +48,8 @@ struct ERA5Monthly <: ERA5Dataset end
 dataset_name(::ERA5Hourly) = "ERA5Hourly"
 dataset_name(::ERA5Monthly) = "ERA5Monthly"
 
+reversed_latitude_axis(::ERA5Dataset) = true
+
 # Wave variables are on a 0.5° grid (720×361), atmospheric variables on 0.25° (1440×721)
 const ERA5_wave_variables = Set([
     :eastward_stokes_drift, :northward_stokes_drift,
@@ -165,7 +167,9 @@ end
 function date_str(date::DateTime)
     y = Dates.year(date)
     m = lpad(Dates.month(date), 2, '0')
-    return "$(y)-$(m)"
+    d = lpad(Dates.day(date), 2, '0')
+    h = lpad(Dates.hour(date), 2, '0')
+    return "$(y)-$(m)-$(d)T$(h)"
 end
 
 start_date_str(date::DateTime) = date_str(date)
@@ -180,34 +184,39 @@ function bbox_strs(::Nothing)
     return "_nothing", "_nothing"
 end
 
+bbox_strs(c::Number) = @sprintf("_%.1f", c), @sprintf("_%.1f", c)
+
 function bbox_strs(c)
     first = @sprintf("_%.1f", c[1])
     second = @sprintf("_%.1f", c[2])
     return first, second
 end
 
-function metadata_prefix(dataset::ERA5Dataset, name, date, bounding_box)
+function region_suffix(::Nothing)
+    return ""
+end
+
+function region_suffix(region)
+    w, e = bbox_strs(region.longitude)
+    s, n = bbox_strs(region.latitude)
+    return string(w, e, s, n)
+end
+
+function metadata_prefix(dataset::ERA5Dataset, name, date, region)
     var = ERA5_dataset_variable_names[name]
     ds = dataset_name(dataset)
     start_date = start_date_str(date)
     end_date = end_date_str(date)
 
-    if !isnothing(bounding_box)
-        w, e = bbox_strs(bounding_box.longitude)
-        s, n = bbox_strs(bounding_box.latitude)
-        suffix = string(w, e, s, n)
-    else
-        suffix = ""
-    end
-
+    suffix = region_suffix(region)
     prefix = string(var, "_", ds, "_", start_date, "_", end_date, suffix)
     prefix = colon2dash(prefix)
     prefix = underscore_spaces(prefix)
     return prefix
 end
 
-function metadata_filename(dataset::ERA5Dataset, name, date, bounding_box)
-    prefix = metadata_prefix(dataset, name, date, bounding_box)
+function metadata_filename(dataset::ERA5Dataset, name, date, region)
+    prefix = metadata_prefix(dataset, name, date, region)
     return string(prefix, ".nc")
 end
 
@@ -222,7 +231,8 @@ inpainted_metadata_path(metadata::ERA5Metadatum) = joinpath(metadata.dir, inpain
 ##### Grid interfaces
 #####
 
-location(::ERA5Metadata) = (Center, Center, Center)
+# ERA5 is a 2D surface dataset — vertical location is Nothing
+dataset_location(::ERA5Dataset, name) = (Center, Center, Nothing)
 
 # ERA5 global coverage: 0-360 longitude, -90 to 90 latitude at 0.25 degree resolution
 longitude_interfaces(::ERA5Metadata) = (0, 360)
diff --git a/src/DataWrangling/ETOPO/ETOPO.jl b/src/DataWrangling/ETOPO/ETOPO.jl
index 56564d3a7..deca7a45f 100644
--- a/src/DataWrangling/ETOPO/ETOPO.jl
+++ b/src/DataWrangling/ETOPO/ETOPO.jl
@@ -53,7 +53,7 @@ const ETOPO_url = "https://www.dropbox.com/scl/fi/6pwalcuuzgtpanysn4h6f/" *
 
 z_interfaces(::ETOPOMetadatum) = (0, 1)
 metadata_url(::ETOPOMetadatum) = ETOPO_url
-metadata_filename(::ETOPO2022, name, date, bounding_box) = "ETOPO_2022_v1_60s_N90W180_surface.nc"
+metadata_filename(::ETOPO2022, name, date, region) = "ETOPO_2022_v1_60s_N90W180_surface.nc"
 
 function download_dataset(metadatum::ETOPOMetadatum)
     fileurl  = metadata_url(metadatum)
diff --git a/src/DataWrangling/GLORYS/GLORYS.jl b/src/DataWrangling/GLORYS/GLORYS.jl
index 0206dcecf..f7ec3953e 100644
--- a/src/DataWrangling/GLORYS/GLORYS.jl
+++ b/src/DataWrangling/GLORYS/GLORYS.jl
@@ -9,12 +9,10 @@ using Oceananigans.Fields: Center
 using NumericalEarth.DataWrangling: Metadata, Metadatum, metadata_path
 using Dates: DateTime, Day, Month
 
-import Oceananigans.Fields:
-    location
-
 import NumericalEarth.DataWrangling:
     all_dates,
     dataset_variable_name,
+    dataset_location,
     default_download_directory,
     longitude_interfaces,
     latitude_interfaces,
@@ -89,6 +87,7 @@ end_date_str(dates::AbstractVector) = last(dates) |> string
 dataset_variable_name(metadata::GLORYSMetadata) = GLORYS_dataset_variable_names[metadata.name]
 
 bbox_strs(::Nothing) = "_nothing", "_nothing"
+bbox_strs(c::Number) = @sprintf("_%.1f", c), @sprintf("_%.1f", c)
 
 function bbox_strs(c)
     first = @sprintf("_%.1f", c[1])
@@ -98,38 +97,42 @@ end
 
 colon2dash(s::String) = replace(s, ":" => "-")
 
-function metadata_prefix(dataset::GLORYSDataset, name, date, bounding_box)
+function region_suffix(::Nothing)
+    return ""
+end
+
+function region_suffix(region)
+    w, e = bbox_strs(region.longitude)
+    s, n = bbox_strs(region.latitude)
+    return string(w, e, s, n)
+end
+
+function metadata_prefix(dataset::GLORYSDataset, name, date, region)
     var = GLORYS_dataset_variable_names[name]
     ds = dataset_name(dataset)
     start_date = start_date_str(date)
     end_date = end_date_str(date)
-    if !isnothing(bounding_box)
-        w, e = bbox_strs(bounding_box.longitude)
-        s, n = bbox_strs(bounding_box.latitude)
-        suffix = string(w, e, s, n)
-    else
-        suffix = ""
-    end
+    suffix = region_suffix(region)
     return string(var, "_",
                   ds, "_",
                   start_date, "_",
                   end_date, suffix) |> colon2dash
 end
 
-function metadata_filename(dataset::GLORYSDataset, name, date, bounding_box)
-    prefix = metadata_prefix(dataset, name, date, bounding_box)
+function metadata_filename(dataset::GLORYSDataset, name, date, region)
+    prefix = metadata_prefix(dataset, name, date, region)
     return string(prefix, ".nc")
 end
 
 function inpainted_metadata_filename(metadata::GLORYSMetadatum)
-    prefix = metadata_prefix(metadata.dataset, metadata.name, metadata.dates, metadata.bounding_box)
+    prefix = metadata_prefix(metadata.dataset, metadata.name, metadata.dates, metadata.region)
     return string(prefix, "_inpainted.jld2")
 end
 
 inpainted_metadata_path(metadata::GLORYSMetadatum) = joinpath(metadata.dir, inpainted_metadata_filename(metadata))
 
-function location(metadata::GLORYSMetadata)
-    metadata.name == :free_surface && return (Center, Center, Nothing)
+function dataset_location(::GLORYSDataset, name)
+    name == :free_surface && return (Center, Center, Nothing)
     return (Center, Center, Center)
 end
 
@@ -140,7 +143,8 @@ latitude_interfaces(::GLORYSMetadata) = (-80, 90)
 
 function z_interfaces(metadata::GLORYSMetadata)
     metadata.name == :free_surface && return (-1.0, 0.0)
-    path = metadata_path(metadata)
+    paths = metadata_path(metadata)
+    path = paths isa AbstractVector ? first(paths) : paths
     ds = Dataset(path)
     zc = - reverse(ds["depth"][:])
     close(ds)
diff --git a/src/DataWrangling/JRA55/JRA55.jl b/src/DataWrangling/JRA55/JRA55.jl
index 4d703dd3f..1c53efbc5 100644
--- a/src/DataWrangling/JRA55/JRA55.jl
+++ b/src/DataWrangling/JRA55/JRA55.jl
@@ -25,6 +25,7 @@ using JLD2
 using Dates
 using Scratch
 
+using Oceananigans: location
 import Oceananigans.Fields: set!
 import Oceananigans.OutputReaders: new_backend, update_field_time_series!
 using Downloads: download
diff --git a/src/DataWrangling/JRA55/JRA55_metadata.jl b/src/DataWrangling/JRA55/JRA55_metadata.jl
index a32b5cc95..53f1fa58c 100644
--- a/src/DataWrangling/JRA55/JRA55_metadata.jl
+++ b/src/DataWrangling/JRA55/JRA55_metadata.jl
@@ -10,8 +10,6 @@ using NumericalEarth.DataWrangling: Metadata, metadata_path, download_progress,
 import Dates: year, month, day
 import Oceananigans.Fields: set!
 import Base
-
-import Oceananigans.Fields: set!, location
 import NumericalEarth.DataWrangling: all_dates, metadata_filename, build_filename, download_dataset, default_download_directory, available_variables
 
 struct MultiYearJRA55 end
@@ -59,13 +57,13 @@ end
 # File name generation specific to each Dataset dataset
 # Note that `RepeatYearJRA55` has only one file associated, so the filename
 # is independent of the date. Override the multi-date fallback to return a plain String.
-metadata_filename(::RepeatYearJRA55, name, date, bounding_box) =
+metadata_filename(::RepeatYearJRA55, name, date, region) =
     "RYF." * JRA55_dataset_variable_names[name] * ".1990_1991.nc"
 
-build_filename(::RepeatYearJRA55, name, dates::AbstractArray, bounding_box) =
+build_filename(::RepeatYearJRA55, name, dates::AbstractArray, region) =
     "RYF." * JRA55_dataset_variable_names[name] * ".1990_1991.nc"
 
-function metadata_filename(::MultiYearJRA55, name, date, bounding_box)
+function metadata_filename(::MultiYearJRA55, name, date, region)
     shortname = JRA55_dataset_variable_names[name]
     year      = Dates.year(date)
     suffix    = "_input4MIPs_atmosphericState_OMIP_MRI-JRA55-do-1-5-0_gr_"
@@ -86,8 +84,6 @@ end
 
 # Convenience functions
 dataset_variable_name(data::JRA55Metadata) = JRA55_dataset_variable_names[data.name]
-location(::JRA55Metadata) = (Center, Center, Center)
-
 available_variables(::MultiYearJRA55)  = JRA55_variable_names
 available_variables(::RepeatYearJRA55) = JRA55_variable_names
 
diff --git a/src/DataWrangling/ORCA/ORCA.jl b/src/DataWrangling/ORCA/ORCA.jl
index c824687ca..4a7f3f623 100644
--- a/src/DataWrangling/ORCA/ORCA.jl
+++ b/src/DataWrangling/ORCA/ORCA.jl
@@ -85,7 +85,7 @@ function metadata_url(metadatum::ORCA12Metadatum)
     end
 end
 
-function metadata_filename(::ORCA1, name, date, bounding_box)
+function metadata_filename(::ORCA1, name, date, region)
     if name == :mesh_mask
         return "eORCA1.2_mesh_mask.nc"
     elseif name == :bottom_height
@@ -95,7 +95,7 @@ function metadata_filename(::ORCA1, name, date, bounding_box)
     end
 end
 
-function metadata_filename(::ORCA12, name, date, bounding_box)
+function metadata_filename(::ORCA12, name, date, region)
     if name == :mesh_mask
         return "grid_mask_eORCA12-GO6.nc"
     elseif name == :bottom_height
diff --git a/src/DataWrangling/OSPapa/OSPapa_flux_observations.jl b/src/DataWrangling/OSPapa/OSPapa_flux_observations.jl
index 2fe93482d..75cad9f9d 100644
--- a/src/DataWrangling/OSPapa/OSPapa_flux_observations.jl
+++ b/src/DataWrangling/OSPapa/OSPapa_flux_observations.jl
@@ -89,9 +89,9 @@ end
 flux_uniform_filename(start_date, end_date) =
     "ocs_papa_flux_uniform_$(Dates.format(start_date, "yyyymmddTHHMMSS"))_$(Dates.format(end_date, "yyyymmddTHHMMSS")).nc"
 
-metadata_filename(::OSPapaFluxHourly, name, date, bounding_box) = flux_uniform_filename(date, date)
+metadata_filename(::OSPapaFluxHourly, name, date, region) = flux_uniform_filename(date, date)
 
-build_filename(::OSPapaFluxHourly, name, dates::AbstractArray, bounding_box) =
+build_filename(::OSPapaFluxHourly, name, dates::AbstractArray, region) =
     flux_uniform_filename(first(dates), last(dates))
 
 function download_dataset(md::OSPapaFluxMetadata)
diff --git a/src/DataWrangling/OSPapa/OSPapa_ocean_observations.jl b/src/DataWrangling/OSPapa/OSPapa_ocean_observations.jl
index 0a18054ce..7123d78f8 100644
--- a/src/DataWrangling/OSPapa/OSPapa_ocean_observations.jl
+++ b/src/DataWrangling/OSPapa/OSPapa_ocean_observations.jl
@@ -44,10 +44,10 @@ const OSPapa_depth_variable_names = Dict(
 )
 
 dataset_variable_name(data::OSPapaMetadata) = OSPapa_dataset_variable_names[data.name]
-
 location(::OSPapaMetadata) = (Center, Center, Center)
 is_three_dimensional(md::OSPapaMetadata) = md.name in (:temperature, :salinity, :eastward_velocity, :northward_velocity)
 reversed_vertical_axis(::OSPapaHourly) = true
+
 function conversion_units(metadatum::OSPapaMetadatum)
     name = metadatum.name
     name == :air_temperature    && return Celsius()
@@ -56,6 +56,7 @@ function conversion_units(metadatum::OSPapaMetadatum)
     name in (:eastward_velocity, :northward_velocity) && return CentimetersPerSecond()
     return nothing
 end
+
 default_inpainting(::OSPapaMetadata) = nothing
 
 #####
@@ -63,12 +64,9 @@ default_inpainting(::OSPapaMetadata) = nothing
 #####
 
 metadata_filename(::OSPapaMetadatum) = OSPAPA_FILENAME
-metadata_filename(::OSPapaHourly, name, date, bounding_box) = OSPAPA_FILENAME
+metadata_filename(::OSPapaHourly, name, date, region) = OSPAPA_FILENAME
 
-function download_dataset(metadata::OSPapaMetadata)
-    download_ospapa_file(metadata.dir)
-    return nothing
-end
+download_dataset(metadata::OSPapaMetadata) = download_ospapa_file(metadata.dir)
 
 function inpainted_metadata_path(metadata::OSPapaMetadata)
     filename = metadata_filename(first(metadata))
diff --git a/src/DataWrangling/WOA/WOA.jl b/src/DataWrangling/WOA/WOA.jl
index 55d2ee6a4..013b4371e 100644
--- a/src/DataWrangling/WOA/WOA.jl
+++ b/src/DataWrangling/WOA/WOA.jl
@@ -39,7 +39,6 @@ import NumericalEarth.DataWrangling:
     available_variables,
     retrieve_data
 
-import Oceananigans.Fields: location
 
 download_WOA_cache::String = ""
 function __init__()
@@ -148,12 +147,12 @@ metaprefix(::WOAMetadatum) = "WOAMetadatum"
 woa_period(::WOAAnnual, date) = 0
 woa_period(::WOAMonthly, date) = Dates.month(date)
 
-function metadata_filename(::WOAAnnual, name, date, bounding_box)
+function metadata_filename(::WOAAnnual, name, date, region)
     varname = WOA_variable_names[name]
     return "woa_$(varname)_annual.nc"
 end
 
-function metadata_filename(::WOAMonthly, name, date, bounding_box)
+function metadata_filename(::WOAMonthly, name, date, region)
     varname = WOA_variable_names[name]
     m = lpad(Dates.month(date), 2, '0')
     return "woa_$(varname)_monthly_$(m).nc"
@@ -162,7 +161,6 @@ end
 # WOA NetCDF variables are named "{tracer}_an" for the objectively analyzed field
 dataset_variable_name(data::WOAMetadata) = WOA_variable_names[data.name] * "_an"
 
-location(::WOAMetadata) = (Center, Center, Center)
 is_three_dimensional(::WOAMetadata) = true
 
 function inpainted_metadata_filename(metadata::WOAMetadatum)
diff --git a/src/DataWrangling/metadata.jl b/src/DataWrangling/metadata.jl
index 7b415742d..cb1153e16 100644
--- a/src/DataWrangling/metadata.jl
+++ b/src/DataWrangling/metadata.jl
@@ -25,6 +25,42 @@ Create a bounding box with `latitude`, `longitude`, and `z` bounds on the sphere
 BoundingBox(; longitude=nothing, latitude=nothing, z=nothing) =
     BoundingBox(longitude, latitude, z)
 
+#####
+##### Column region and interpolation types
+#####
+
+struct Linear end
+struct Nearest end
+
+"""
+    Column(longitude, latitude; z=nothing, interpolation=Linear())
+
+Create a column region at a single horizontal point `(longitude, latitude)`.
+When used as a `Metadata` region, `native_grid` returns a single-column
+`RectilinearGrid` and `location` reduces horizontal dimensions to `Nothing`.
+
+Keyword Arguments
+=================
+
+- `z`: depth range tuple `(z_bottom, z_top)` for restricting downloads
+  (used by CopernicusMarine/GLORYS). Default: `nothing` (full depth).
+- `interpolation`: method for extracting data from the surrounding grid
+  cells. `Linear()` (default) bilinearly interpolates to the exact point;
+  `Nearest()` selects the closest grid cell.
+"""
+struct Column{X, Y, Z, I}
+    longitude :: X
+    latitude :: Y
+    z :: Z
+    interpolation :: I
+end
+
+Column(longitude, latitude; z=nothing, interpolation=Linear()) =
+    Column(longitude, latitude, z, interpolation)
+
+Base.summary(col::Column) = string("Column(longitude=", prettysummary(col.longitude),
+                                   ", latitude=", prettysummary(col.latitude), ")")
+
 struct DatewiseFilename{A}
     filenames :: A
 end
@@ -37,16 +73,16 @@ getfilename(f::DatewiseFilename, i) = f.filenames[i]
 getfilename(f::String, i) = f
 getfilename(::Nothing, i) = nothing
 
-struct Metadata{V, D, B, S, F}
+struct Metadata{V, D, R, S, F}
     name :: S
     dataset :: V
     dates :: D
-    bounding_box :: B
+    region :: R
     dir :: String
     filename :: F
 end
 
-Metadata(name, dataset, dates, bbox, dir) = Metadata(name, dataset, dates, bbox, dir, nothing)
+Metadata(name, dataset, dates, region, dir) = Metadata(name, dataset, dates, region, dir, nothing)
 
 is_three_dimensional(::Metadata) = true
 z_interfaces(md::Metadata) = z_interfaces(md.dataset)
@@ -58,7 +94,7 @@ latitude_interfaces(md::Metadata) = latitude_interfaces(md.dataset)
              dataset,
              dates = all_dates(dataset, variable_name),
              dir = default_download_directory(dataset),
-             bounding_box = nothing,
+             region = nothing,
              filename = nothing,
              start_date = nothing,
              end_date = nothing)
@@ -88,7 +124,9 @@ Keyword Arguments
               (`Dates.AbstractDateTime` or `CFTime.AbstractCFDateTime`). If outside the
                 date range of the dataset, the last allowable date is chosen. Default: nothing.
 
-- `bounding_box`: Specifies the bounds of the dataset. See [`BoundingBox`](@ref).
+- `region`: Specifies the spatial region of the dataset. Can be a [`BoundingBox`](@ref)
+            for a rectangular region, a [`Column`](@ref) for a single horizontal location,
+            or `nothing` for the full domain.
 
 - `filename`: The filename(s) for the dataset. If `nothing`, the filename is computed from
               the dataset type. Can be a `String` (single file for all dates) or a
@@ -100,7 +138,7 @@ function Metadata(variable_name;
                   dataset,
                   dates = all_dates(dataset, variable_name),
                   dir = default_download_directory(dataset),
-                  bounding_box = nothing,
+                  region = nothing,
                   filename = nothing,
                   start_date = nothing,
                   end_date = nothing)
@@ -117,10 +155,10 @@ function Metadata(variable_name;
     end
 
     if isnothing(filename)
-        filename = build_filename(dataset, variable_name, dates, bounding_box)
+        filename = build_filename(dataset, variable_name, dates, region)
     end
 
-    return Metadata(variable_name, dataset, dates, bounding_box, dir, filename)
+    return Metadata(variable_name, dataset, dates, region, dir, filename)
 end
 
 const AnyDateTime  = Union{AbstractCFDateTime, Dates.AbstractDateTime}
@@ -140,7 +178,7 @@ end
 """
     Metadatum(variable_name;
               dataset,
-              bounding_box = nothing,
+              region = nothing,
               date = first_date(dataset, variable_name),
               filename = nothing,
               dir = default_download_directory(dataset))
@@ -149,7 +187,7 @@ A specialized constructor for a [`Metadata`](@ref) object with a single date, re
 """
 function Metadatum(variable_name;
                    dataset,
-                   bounding_box = nothing,
+                   region = nothing,
                    date = first_date(dataset, variable_name),
                    filename = nothing,
                    dir = default_download_directory(dataset))
@@ -164,10 +202,10 @@ function Metadatum(variable_name;
     end
 
     if isnothing(filename)
-        filename = metadata_filename(dataset, variable_name, date, bounding_box)
+        filename = metadata_filename(dataset, variable_name, date, region)
     end
 
-    return Metadata(variable_name, dataset, date, bounding_box, dir, filename)
+    return Metadata(variable_name, dataset, date, region, dir, filename)
 end
 
 datestr(md::Metadata) = string(first(md.dates), "--", last(md.dates))
@@ -192,9 +230,9 @@ function Base.show(io::IO, metadata::Metadata)
     "├── dataset: ", prettysummary(metadata.dataset), '\n',
     "├── dates: ", prettysummary(metadata.dates), '\n')
 
-    bbox = metadata.bounding_box
-    if !isnothing(bbox)
-        print(io, "├── bounding_box: ", summary(bbox), '\n')
+    rgn = metadata.region
+    if !isnothing(rgn)
+        print(io, "├── region: ", summary(rgn), '\n')
     end
 
     print(io, "├── filename: $(metadata.filename)", '\n')
@@ -213,17 +251,17 @@ Base.summary(md::Metadata) = string(metaprefix(md),
 Base.length(metadata::Metadatum) = 1
 
 @propagate_inbounds Base.getindex(m::Metadata, i::Int) =
-    Metadata(m.name, m.dataset, m.dates[i], m.bounding_box, m.dir, getfilename(m.filename, i))
+    Metadata(m.name, m.dataset, m.dates[i], m.region, m.dir, getfilename(m.filename, i))
 
 @propagate_inbounds Base.first(m::Metadata) =
-    Metadata(m.name, m.dataset, m.dates[1], m.bounding_box, m.dir, getfilename(m.filename, 1))
+    Metadata(m.name, m.dataset, m.dates[1], m.region, m.dir, getfilename(m.filename, 1))
 
 @propagate_inbounds Base.last(m::Metadata) =
-    Metadata(m.name, m.dataset, m.dates[end], m.bounding_box, m.dir, getfilename(m.filename, lastindex(m.dates)))
+    Metadata(m.name, m.dataset, m.dates[end], m.region, m.dir, getfilename(m.filename, lastindex(m.dates)))
 
 @inline function Base.iterate(m::Metadata, i=1)
    if (i % UInt) - 1 < length(m)
-        return Metadata(m.name, m.dataset, m.dates[i], m.bounding_box, m.dir, getfilename(m.filename, i)), i + 1
+        return Metadata(m.name, m.dataset, m.dates[i], m.region, m.dir, getfilename(m.filename, i)), i + 1
     else
         return nothing
     end
@@ -293,6 +331,15 @@ Return the name used for the variable `metadata.name` in its raw dataset file.
 """
 function dataset_variable_name end
 
+"""
+    dataset_location(dataset, variable_name)
+
+Return the native field location `(LX, LY, LZ)` for `variable_name` in
+`dataset`. Defaults to `(Center, Center, Center)`. Only datasets with
+staggered variables (e.g., ECCO velocity fields) need to extend this.
+"""
+dataset_location(dataset, variable_name) = (Center, Center, Center)
+
 # Note: all_dates needs to be extended for any new dataset.
 """
     all_dates(metadata)
@@ -323,7 +370,7 @@ Return the stored filename(s) of `metadata`.
 metadata_filename(metadata::Metadata) = metadata.filename
 
 """
-    metadata_filename(dataset, name, date, bounding_box)
+    metadata_filename(dataset, name, date, region)
 
 Compute the filename for a single date. Extended by each dataset module.
 """
@@ -331,12 +378,12 @@ function metadata_filename end
 
 # Internal: build filename for construction.
 # Single date: delegate to metadata_filename
-build_filename(dataset, name, date, bounding_box) =
-    metadata_filename(dataset, name, date, bounding_box)
+build_filename(dataset, name, date, region) =
+    metadata_filename(dataset, name, date, region)
 
 # Multi-date: one filename per date, wrapped in DatewiseFilename
-build_filename(dataset, name, dates::AbstractArray, bounding_box) =
-    DatewiseFilename([metadata_filename(dataset, name, d, bounding_box) for d in dates])
+build_filename(dataset, name, dates::AbstractArray, region) =
+    DatewiseFilename([metadata_filename(dataset, name, d, region) for d in dates])
 
 """
     available_variables(metadata)
diff --git a/src/DataWrangling/metadata_field.jl b/src/DataWrangling/metadata_field.jl
index 2f3bc1ba8..6a7c60b4d 100644
--- a/src/DataWrangling/metadata_field.jl
+++ b/src/DataWrangling/metadata_field.jl
@@ -2,45 +2,93 @@ using NCDatasets
 using JLD2
 using NumericalEarth.InitialConditions: interpolate!
 using Statistics: median
+using Oceananigans.Grids: λnodes, φnodes
+using Oceananigans.Architectures: on_architecture
+using Oceananigans.Fields: fractional_x_index, fractional_y_index
 
-import Oceananigans.Fields: set!, Field
+import Oceananigans.Fields: set!, Field, location
+
+#####
+##### Location with automatic restriction based on region
+#####
+
+location(metadata::Metadata) = restrict_location(dataset_location(metadata.dataset, metadata.name), metadata.region)
+
+restrict_location(loc, ::Nothing) = loc
+restrict_location(loc, ::BoundingBox) = loc
+restrict_location((LX, LY, LZ), ::Column) = (Nothing, Nothing, LZ)
+
+#####
+##### Native grid construction — dispatches on region type
+#####
 
 restrict(::Nothing, interfaces, N) = interfaces, N
 
 # TODO support stretched native grids
 function restrict(bbox_interfaces, interfaces, N)
-    Δ = interfaces[2] - interfaces[1]
+    extent = interfaces[end] - interfaces[1]
     rΔ = bbox_interfaces[2] - bbox_interfaces[1]
-    ϵ = rΔ / Δ
-    rN = ceil(Int, ϵ * N)  # Round up to ensure bounding box is covered
+    rN = round(Int, rΔ / extent * N)
+    rN = max(rN, 1)  # at least one cell
     return bbox_interfaces, rN
 end
 
 """
     native_grid(metadata::Metadata, arch=CPU(); halo = (3, 3, 3))
 
-Return a `LatitudeLongitudeGrid` on `arch` corresponding to the native grid of `metadata` with `halo` size.
+Return the native grid corresponding to `metadata` with `halo` size.
+Returns a `LatitudeLongitudeGrid` for global or `BoundingBox` regions,
+and a column `RectilinearGrid` for `Column` regions.
 """
-function native_grid(metadata::Metadata, arch=CPU(); halo = (3, 3, 3))
+native_grid(metadata::Metadata, arch=CPU(); halo=(3, 3, 3)) =
+    construct_native_grid(metadata, metadata.region, arch; halo)
+
+# Full global grid (no region restriction)
+function construct_native_grid(metadata, ::Nothing, arch; halo)
     Nx, Ny, Nz, _ = size(metadata)
     z = z_interfaces(metadata)
-
     FT = eltype(metadata)
+    longitude = longitude_interfaces(metadata)
+    latitude = latitude_interfaces(metadata)
+
+    grid = LatitudeLongitudeGrid(arch, FT; size = (Nx, Ny, Nz),
+                                 halo, longitude, latitude, z)
+    return grid
+end
 
+# BoundingBox-restricted LatitudeLongitudeGrid
+function construct_native_grid(metadata, bbox::BoundingBox, arch; halo)
+    Nx, Ny, Nz, _ = size(metadata)
+    z = z_interfaces(metadata)
+    FT = eltype(metadata)
     longitude = longitude_interfaces(metadata)
     latitude = latitude_interfaces(metadata)
 
-    # Restrict with BoundingBox
     # TODO: can we restrict in `z` as well?
-    bbox = metadata.bounding_box
-    if !isnothing(bbox)
-        longitude, Nx = restrict(bbox.longitude, longitude, Nx)
-        latitude, Ny = restrict(bbox.latitude, latitude, Ny)
-    end
+    longitude, Nx = restrict(bbox.longitude, longitude, Nx)
+    latitude, Ny = restrict(bbox.latitude, latitude, Ny)
+
+    # Clamp halo so it does not exceed grid size in any dimension
+    halo = min.(halo, (Nx, Ny, Nz))
 
     grid = LatitudeLongitudeGrid(arch, FT; size = (Nx, Ny, Nz),
                                  halo, longitude, latitude, z)
+    return grid
+end
+
+# Column RectilinearGrid
+function construct_native_grid(metadata, col::Column, arch; halo)
+    _, _, Nz, _ = size(metadata)
+    z = z_interfaces(metadata)
+    FT = eltype(metadata)
 
+    grid = RectilinearGrid(arch, FT;
+                           size = Nz,
+                           x = FT(col.longitude),
+                           y = FT(col.latitude),
+                           z,
+                           halo = halo[3],
+                           topology = (Flat, Flat, Bounded))
     return grid
 end
 
@@ -68,6 +116,13 @@ function retrieve_data(metadata::Metadatum)
     end
 
     close(ds)
+
+    # ERA5 (and some other datasets) store latitude north-to-south;
+    # flip to south-to-north to match the grid.
+    if reversed_latitude_axis(metadata.dataset)
+        data = reverse(data, dims=2)
+    end
+
     return data
 end
 
@@ -93,6 +148,14 @@ function Field(metadata::Metadatum, arch=CPU();
 
     download_dataset(metadata)
 
+    # Column regions need special handling: the downloaded file may contain
+    # more data than a single column (e.g. CopernicusMarine returns a small
+    # grid around the point). Load onto an intermediate grid from the file's
+    # actual dimensions, then extract the column.
+    if metadata.region isa Column
+        return column_field_from_file(metadata, arch; inpainting, mask, halo, cache_inpainted_data)
+    end
+
     grid = native_grid(metadata, arch; halo)
     LX, LY, LZ = location(metadata)
     field = Field{LX, LY, LZ}(grid)
@@ -178,9 +241,137 @@ function set!(target_field::Field, metadata::Metadatum; kw...)
     end
 
     interpolate!(target_field, meta_field)
+
     return target_field
 end
 
+#####
+##### Column field construction
+#####
+
+function column_field_from_file(metadata, arch; halo=(3, 3, 3), kw...)
+    column_grid = native_grid(metadata, arch; halo)
+
+    # Read the file's actual dimensions to build a matching intermediate grid
+    path = metadata_path(metadata)
+    ds = Dataset(path)
+    varname = dataset_variable_name(metadata)
+    var = ds[varname]
+    data_size = size(var)
+    Nx_file, Ny_file = data_size[1], data_size[2]
+
+    # Read coordinate arrays
+    lon_dimname = NCDatasets.dimnames(var)[1]
+    lat_dimname = NCDatasets.dimnames(var)[2]
+    λ = haskey(ds, lon_dimname) ? ds[lon_dimname][:] : ds["longitude"][:]
+    φ = haskey(ds, lat_dimname) ? ds[lat_dimname][:] : ds["latitude"][:]
+    close(ds)
+
+    if reversed_latitude_axis(metadata.dataset)
+        reverse!(φ)
+    end
+
+    _, _, Nz, _ = size(metadata)
+    z = z_interfaces(metadata)
+    FT = eltype(metadata)
+
+    # Build cell interfaces from centers
+    Δλ = Nx_file > 1 ? λ[2] - λ[1] : FT(1)
+    λf = range(λ[1] - Δλ/2, stop = λ[end] + Δλ/2, length = Nx_file + 1)
+
+    Δφ = Ny_file > 1 ? φ[2] - φ[1] : FT(1)
+    φf = range(φ[1] - Δφ/2, stop = φ[end] + Δφ/2, length = Ny_file + 1)
+
+    halo = min.(halo, (Nx_file, Ny_file, Nz))
+
+    intermediate_grid = LatitudeLongitudeGrid(arch, FT;
+                                              size = (Nx_file, Ny_file, Nz),
+                                              halo, longitude = λf, latitude = φf, z)
+
+    # Load data onto intermediate grid (no inpainting — columns have no horizontal neighbors)
+    LX, LY, LZ = dataset_location(metadata.dataset, metadata.name)
+    intermediate_field = Field{LX, LY, LZ}(intermediate_grid)
+
+    data = retrieve_data(metadata)
+    set_metadata_field!(intermediate_field, data, metadata)
+    fill_halo_regions!(intermediate_field)
+
+    # Extract column
+    _, _, LZ_col = location(metadata)
+    col_field = Field{Nothing, Nothing, LZ_col}(column_grid)
+    extract_column!(col_field, intermediate_field, metadata.region)
+
+    return col_field
+end
+
+#####
+##### Column extraction utilities
+#####
+
+# Dispatch extraction on interpolation method
+function extract_column!(column_field, intermediate_field, col::Column)
+    extract_column!(column_field, intermediate_field, col, col.interpolation)
+end
+
+function extract_column!(column_field, intermediate_field, col, ::Linear)
+    grid = intermediate_field.grid
+    arch = architecture(grid)
+    LX, LY, LZ = Oceananigans.Fields.location(intermediate_field)
+    locs = (LX(), LY(), LZ())
+
+    # Fractional indices (1-based, continuous)
+    fi = fractional_x_index(col.longitude, locs, grid)
+    fj = fractional_y_index(col.latitude,  locs, grid)
+
+    # Lower-left index and weights
+    i₁ = clamp(floor(Int, fi), 1, size(grid, 1))
+    j₁ = clamp(floor(Int, fj), 1, size(grid, 2))
+    i₂ = clamp(i₁ + 1, 1, size(grid, 1))
+    j₂ = clamp(j₁ + 1, 1, size(grid, 2))
+
+    wx = clamp(fi - floor(fi), 0, 1)
+    wy = clamp(fj - floor(fj), 0, 1)
+
+    launch!(arch, column_field.grid, :z, _bilinear_interpolate_column!,
+            column_field, intermediate_field, i₁, j₁, i₂, j₂, wx, wy)
+
+    return nothing
+end
+
+@kernel function _bilinear_interpolate_column!(column_field, source, i₁, j₁, i₂, j₂, wx, wy)
+    k = @index(Global, Linear)
+    @inbounds begin
+        v00 = source[i₁, j₁, k]
+        v10 = source[i₂, j₁, k]
+        v01 = source[i₁, j₂, k]
+        v11 = source[i₂, j₂, k]
+        column_field[1, 1, k] = (1 - wx) * (1 - wy) * v00 +
+                                     wx  * (1 - wy) * v10 +
+                                (1 - wx) *      wy  * v01 +
+                                     wx  *      wy  * v11
+    end
+end
+
+function extract_column!(column_field, intermediate_field, col, ::Nearest)
+    grid = intermediate_field.grid
+    arch = architecture(grid)
+    LX, LY, LZ = Oceananigans.Fields.location(intermediate_field)
+    locs = (LX(), LY(), LZ())  # fractional index functions expect instances, not types
+
+    # Use Oceananigans' fractional index machinery (handles cyclic longitude etc.)
+    i★ = round(Int, fractional_x_index(col.longitude, locs, grid))
+    j★ = round(Int, fractional_y_index(col.latitude,  locs, grid))
+
+    launch!(arch, column_field.grid, :z, copy_column!, column_field, intermediate_field, i★, j★)
+
+    return nothing
+end
+
+@kernel function copy_column!(column_field, source_field, i★, j★)
+    k = @index(Global, Linear)
+    @inbounds column_field[1, 1, k] = source_field[i★, j★, k]
+end
+
 # manglings
 struct ShiftSouth end
 struct AverageNorthSouth end
diff --git a/src/DataWrangling/metadata_field_time_series.jl b/src/DataWrangling/metadata_field_time_series.jl
index 0bfed83af..142eb6b07 100644
--- a/src/DataWrangling/metadata_field_time_series.jl
+++ b/src/DataWrangling/metadata_field_time_series.jl
@@ -116,7 +116,8 @@ function FieldTimeSeries(metadata::Metadata, grid::AbstractGrid;
     download_dataset(metadata)
 
     inpainting isa Int && (inpainting = NearestNeighborInpainting(inpainting))
-    backend = DatasetBackend(time_indices_in_memory, metadata; on_native_grid, inpainting, cache_inpainted_data)
+    is_native = grid == native_grid(metadata)
+    backend = DatasetBackend(time_indices_in_memory, metadata; on_native_grid=is_native, inpainting, cache_inpainted_data)
 
     times = native_times(metadata)
     loc = LX, LY, LZ = location(metadata)
@@ -131,11 +132,11 @@ function FieldTimeSeries(variable_name::Symbol;
                          dataset, dir,
                          architecture = CPU(),
                          start_date = first_date(dataset, variable_name),
-                         end_date = first_date(dataset, variable_name),
+                         end_date = last_date(dataset, variable_name),
                          kw...)
 
     native_dates = all_dates(dataset, variable_name)
     dates = compute_native_date_range(native_dates, start_date, end_date)
-    metadata = Metadata(variable_name, dataset, dates, dir)
+    metadata = Metadata(variable_name; dataset, dates, dir)
     return FieldTimeSeries(metadata, architecture; kw...)
 end
diff --git a/src/NumericalEarth.jl b/src/NumericalEarth.jl
index 21f61fc19..2a3ebb6aa 100644
--- a/src/NumericalEarth.jl
+++ b/src/NumericalEarth.jl
@@ -34,6 +34,8 @@ export
     regrid_bathymetry,
     Metadata,
     Metadatum,
+    BoundingBox,
+    Column, Linear, Nearest,
     ECCOMetadatum,
     EN4Metadatum,
     ETOPO2022,
@@ -60,9 +62,12 @@ export
     frazil_heat_flux, net_ocean_heat_flux, sea_ice_ocean_heat_flux, atmosphere_ocean_heat_flux,
     net_ocean_salinity_flux, sea_ice_ocean_salinity_flux, atmosphere_ocean_salinity_flux,
     net_ocean_freshwater_flux, sea_ice_ocean_freshwater_flux, atmosphere_ocean_freshwater_flux,
-    meridional_heat_transport
+    meridional_heat_transport,
+    location,
+    native_grid
 
 using Oceananigans
+import Oceananigans: location
 using Oceananigans.Operators: ℑxyᶠᶜᵃ, ℑxyᶜᶠᵃ
 using DataDeps
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 39e6985b0..3061c0c14 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -22,8 +22,9 @@ gpu_test = parse(Bool, get(ENV, "GPU_TEST", "false"))
 
 if filter_tests!(testsuite, args)
     # Always remove tests that are treated separately
-    delete!(testsuite, "test_downloading")
+    delete!(testsuite, "test_jra55_ecco_en4_etopo_downloading")
     delete!(testsuite, "test_cds_downloading")
+    delete!(testsuite, "test_glorys_downloading")
     delete!(testsuite, "test_distributed_utils")
     delete!(testsuite, "test_reactant")
 
diff --git a/test/test_cds_downloading.jl b/test/test_cds_downloading.jl
index 9b721da6d..9802c2be2 100644
--- a/test/test_cds_downloading.jl
+++ b/test/test_cds_downloading.jl
@@ -18,11 +18,11 @@ start_date = DateTime(2005, 2, 16, 12)
     dataset = ERA5Hourly()
 
     # Use a small bounding box to reduce download time
-    bounding_box = NumericalEarth.DataWrangling.BoundingBox(longitude=(0, 5), latitude=(40, 45))
+    region = NumericalEarth.DataWrangling.BoundingBox(longitude=(0, 5), latitude=(40, 45))
 
     @testset "Download ERA5 temperature data" begin
         variable = :temperature
-        metadatum = Metadatum(variable; dataset, bounding_box, date=start_date)
+        metadatum = Metadatum(variable; dataset, region, date=start_date)
 
         # Clean up any existing file
         filepath = metadata_path(metadatum)
@@ -84,13 +84,13 @@ start_date = DateTime(2005, 2, 16, 12)
 
     @testset "ERA5 metadata properties" begin
         variable = :temperature
-        metadatum = Metadatum(variable; dataset, bounding_box, date=start_date)
+        metadatum = Metadatum(variable; dataset, region, date=start_date)
 
         # Test metadata properties
         @test metadatum.name == :temperature
         @test metadatum.dataset isa ERA5Hourly
         @test metadatum.dates == start_date
-        @test metadatum.bounding_box == bounding_box
+        @test metadatum.region == region
 
         # Test size (should be global ERA5 size with 1 time step)
         Nx, Ny, Nz, Nt = size(metadatum)
@@ -141,7 +141,7 @@ start_date = DateTime(2005, 2, 16, 12)
 
         @testset "Field creation from ERA5 on $A" begin
             variable = :temperature
-            metadatum = Metadatum(variable; dataset, bounding_box, date=start_date)
+            metadatum = Metadatum(variable; dataset, region, date=start_date)
 
             # Download if not present (falls back to NumericalEarthArtifacts if CDS is unreachable)
             filepath = metadata_path(metadatum)
@@ -170,7 +170,7 @@ start_date = DateTime(2005, 2, 16, 12)
 
         @testset "Setting a field from ERA5 metadata on $A" begin
             variable = :temperature
-            metadatum = Metadatum(variable; dataset, bounding_box, date=start_date)
+            metadatum = Metadatum(variable; dataset, region, date=start_date)
 
             # Download if not present (falls back to NumericalEarthArtifacts if CDS is unreachable)
             filepath = metadata_path(metadatum)
diff --git a/test/test_column_field.jl b/test/test_column_field.jl
new file mode 100644
index 000000000..ea33bad6c
--- /dev/null
+++ b/test/test_column_field.jl
@@ -0,0 +1,233 @@
+include("runtests_setup.jl")
+
+using NumericalEarth.DataWrangling: Column, Linear, Nearest,
+                                    BoundingBox, native_grid,
+                                    restrict_location, dataset_location
+
+using NumericalEarth.DataWrangling: extract_column!
+
+using Oceananigans
+using Oceananigans.BoundaryConditions: fill_halo_regions!
+using Oceananigans.Fields: location
+using Oceananigans.Grids: λnodes, φnodes, topology, Flat, Bounded, Periodic
+
+# Test coordinates for end-to-end Column tests (ECCO4 ocean point)
+const test_longitude = 12.0
+const test_latitude = -50.0
+
+@testset "extract_column! with Nearest interpolation" begin
+    for arch in test_architectures
+        A = typeof(arch)
+        @testset "Nearest extraction on $A" begin
+            # Create a LatitudeLongitudeGrid with spatially varying data
+            intermediate_grid = LatitudeLongitudeGrid(arch;
+                size = (4, 4, 2),
+                longitude = (0, 4),
+                latitude = (0, 4),
+                z = (-20, 0))
+
+            intermediate_field = CenterField(intermediate_grid)
+
+            # Set distinct values at each horizontal point
+            for i in 1:4, j in 1:4, k in 1:2
+                @allowscalar intermediate_field[i, j, k] = 10 * i + j + 0.1 * k
+            end
+            fill_halo_regions!(intermediate_field)
+
+            # Column near grid point (3, 2) → lon≈2.5, lat≈1.5
+            col = Column(2.5, 1.5; interpolation=Nearest())
+            column_grid = RectilinearGrid(arch;
+                size = 2,
+                x = 2.5,
+                y = 1.5,
+                z = (-20, 0),
+                halo = 3,
+                topology = (Flat, Flat, Bounded))
+
+            column_field = Field{Nothing, Nothing, Center}(column_grid)
+
+            extract_column!(column_field, intermediate_field, col, Nearest())
+
+            # Find expected nearest indices
+            λnodes_arr = λnodes(intermediate_grid, Center(); with_halos=false)
+            φnodes_arr = φnodes(intermediate_grid, Center(); with_halos=false)
+            i★ = argmin(abs.(λnodes_arr .- 2.5))
+            j★ = argmin(abs.(φnodes_arr .- 1.5))
+
+            @allowscalar begin
+                for k in 1:2
+                    @test column_field[1, 1, k] == intermediate_field[i★, j★, k]
+                end
+            end
+        end
+
+        @testset "Nearest extraction preserves vertical profile on $A" begin
+            intermediate_grid = LatitudeLongitudeGrid(arch;
+                size = (3, 3, 5),
+                longitude = (10, 13),
+                latitude = (40, 43),
+                z = (-50, 0))
+
+            intermediate_field = CenterField(intermediate_grid)
+
+            # Set a vertical profile: value = depth level
+            for k in 1:5
+                interior(intermediate_field)[:, :, k] .= Float64(k)
+            end
+            fill_halo_regions!(intermediate_field)
+
+            col = Column(11.5, 41.5; interpolation=Nearest())
+            column_grid = RectilinearGrid(arch;
+                size = 5,
+                x = 11.5,
+                y = 41.5,
+                z = (-50, 0),
+                halo = 3,
+                topology = (Flat, Flat, Bounded))
+
+            column_field = Field{Nothing, Nothing, Center}(column_grid)
+            extract_column!(column_field, intermediate_field, col, Nearest())
+
+            @allowscalar begin
+                for k in 1:5
+                    @test column_field[1, 1, k] == k
+                end
+            end
+        end
+    end
+end
+
+@testset "extract_column! dispatch routes on interpolation type" begin
+    for arch in test_architectures
+        A = typeof(arch)
+        @testset "Dispatch on $A" begin
+            intermediate_grid = LatitudeLongitudeGrid(arch;
+                size = (4, 4, 2),
+                longitude = (0, 4),
+                latitude = (0, 4),
+                z = (-20, 0))
+
+            intermediate_field = CenterField(intermediate_grid)
+            interior(intermediate_field) .= 42.0
+            fill_halo_regions!(intermediate_field)
+
+            column_grid = RectilinearGrid(arch;
+                size = 2,
+                x = 2.0,
+                y = 2.0,
+                z = (-20, 0),
+                halo = 3,
+                topology = (Flat, Flat, Bounded))
+
+            # Column dispatch routes to the correct method
+            col_nearest = Column(2.0, 2.0; interpolation=Nearest())
+            cf = Field{Nothing, Nothing, Center}(column_grid)
+            extract_column!(cf, intermediate_field, col_nearest)
+
+            @allowscalar begin
+                @test cf[1, 1, 1] == 42.0
+                @test cf[1, 1, 2] == 42.0
+            end
+        end
+    end
+end
+
+@testset "End-to-end Column Field construction" begin
+    for arch in test_architectures
+        A = typeof(arch)
+
+        @testset "Column Field with Linear interpolation on $A" begin
+            col = Column(test_longitude, test_latitude; interpolation=Linear())
+            md = Metadatum(:temperature; dataset=ECCO4Monthly(), date=start_date, region=col)
+            field = Field(md, arch)
+
+            @test field.grid isa RectilinearGrid
+            @test topology(field.grid) == (Flat, Flat, Bounded)
+            @test location(field) == (Nothing, Nothing, Center)
+
+            # Field should have non-trivial data (not all zeros)
+            @allowscalar begin
+                @test any(!=(0), interior(field))
+            end
+        end
+
+        @testset "Column Field with Nearest interpolation on $A" begin
+            col = Column(test_longitude, test_latitude; interpolation=Nearest())
+            md = Metadatum(:temperature; dataset=ECCO4Monthly(), date=start_date, region=col)
+            field = Field(md, arch)
+
+            @test field.grid isa RectilinearGrid
+            @test location(field) == (Nothing, Nothing, Center)
+
+            @allowscalar begin
+                @test any(!=(0), interior(field))
+            end
+        end
+
+        @testset "set! with Column metadata on $A" begin
+            col = Column(test_longitude, test_latitude)
+            md = Metadatum(:temperature; dataset=ECCO4Monthly(), date=start_date, region=col)
+
+            # Build a target column field
+            column_grid = native_grid(md, arch)
+            target = Field{Nothing, Nothing, Center}(column_grid)
+
+            set!(target, md)
+
+            @allowscalar begin
+                @test any(!=(0), interior(target))
+            end
+        end
+
+        @testset "Column Linear vs Nearest give similar results on $A" begin
+            col_lin = Column(test_longitude, test_latitude; interpolation=Linear())
+            col_near = Column(test_longitude, test_latitude; interpolation=Nearest())
+
+            md_lin = Metadatum(:temperature; dataset=ECCO4Monthly(), date=start_date, region=col_lin)
+            md_near = Metadatum(:temperature; dataset=ECCO4Monthly(), date=start_date, region=col_near)
+
+            field_lin = Field(md_lin, arch)
+            field_near = Field(md_near, arch)
+
+            # Both should produce finite, non-zero vertical profiles
+            @allowscalar begin
+                @test all(isfinite, interior(field_lin))
+                @test all(isfinite, interior(field_near))
+            end
+        end
+    end
+end
+
+@testset "Column native_grid construction" begin
+    @testset "ECCO4 Column grid" begin
+        col = Column(35.1, 50.1)
+        md = Metadatum(:temperature; dataset=ECCO4Monthly(), region=col)
+        grid = native_grid(md)
+
+        @test grid isa RectilinearGrid
+        @test topology(grid) == (Flat, Flat, Bounded)
+        _, _, Nz, _ = size(md)
+        @test size(grid) == (1, 1, Nz)
+    end
+
+    @testset "ERA5 Column grid" begin
+        col = Column(200.0, 35.0)
+        md = Metadatum(:temperature; dataset=ERA5Hourly(),
+                       date=DateTime(2020, 1, 1), region=col)
+        grid = native_grid(md)
+
+        @test grid isa RectilinearGrid
+        @test topology(grid) == (Flat, Flat, Bounded)
+        # ERA5 has z = (0, 1), single level
+        @test size(grid) == (1, 1, 1)
+    end
+
+    @testset "Column grid uses Float32 for ECCO" begin
+        col = Column(123.4, -45.6)
+        md = Metadatum(:temperature; dataset=ECCO4Monthly(), region=col)
+        grid = native_grid(md)
+
+        # ECCO metadata has Float32 eltype
+        @test eltype(grid) == Float32
+    end
+end
diff --git a/test/test_distributed_utils.jl b/test/test_distributed_utils.jl
index a594bb4bb..bc366d87d 100644
--- a/test/test_distributed_utils.jl
+++ b/test/test_distributed_utils.jl
@@ -53,7 +53,7 @@ Base.size(::TrivalBathymetry, variable) = (Nλ, Nφ, 1)
 z_interfaces(::TrivalBathymetry) = (0, 1)
 longitude_interfaces(::TrivalBathymetry) = (-180, 180)
 latitude_interfaces(::TrivalBathymetry) = (0, 50)
-metadata_filename(::TrivalBathymetry, name, date, bounding_box) = "trivial_bathymetry.nc"
+metadata_filename(::TrivalBathymetry, name, date, region) = "trivial_bathymetry.nc"
 
 @testset "Distributed ECCO download" begin
     dates = DateTimeProlepticGregorian(1992, 1, 1) : Month(1) : DateTimeProlepticGregorian(1994, 4, 1)
diff --git a/test/test_glorys_downloading.jl b/test/test_glorys_downloading.jl
index 612be7fd3..d00b7c4f3 100644
--- a/test/test_glorys_downloading.jl
+++ b/test/test_glorys_downloading.jl
@@ -9,24 +9,22 @@ using Oceananigans.Fields: location
 
 @testset "Downloading GLORYS data" begin
     variables = (:temperature, :salinity, :u_velocity, :v_velocity, :free_surface)
-    bounding_box = BoundingBox(longitude=(200, 202), latitude=(35, 37))
+    region = BoundingBox(longitude=(200, 202), latitude=(35, 37))
     dataset = GLORYSDaily()
     for variable in variables
-        metadatum = Metadatum(variable; dataset, bounding_box)
+        metadatum = Metadatum(variable; dataset, region)
         filepath = NumericalEarth.DataWrangling.metadata_path(metadatum)
         isfile(filepath) && rm(filepath; force=true)
-        download_dataset_with_fallback(filepath; dataset_name="GLORYSDaily $variable") do
-            NumericalEarth.DataWrangling.download_dataset(metadatum)
-        end
+        NumericalEarth.DataWrangling.download_dataset(metadatum)
         @test isfile(filepath)
     end
 end
 
 @testset "Download and set GLORYS free_surface" begin
     for arch in test_architectures
-        bounding_box = BoundingBox(longitude=(200, 202), latitude=(35, 37))
+        region = BoundingBox(longitude=(200, 202), latitude=(35, 37))
         dataset = GLORYSDaily()
-        md = Metadatum(:free_surface; dataset, bounding_box)
+        md = Metadatum(:free_surface; dataset, region)
 
         @test !is_three_dimensional(md)
         @test location(md) === (Center, Center, Nothing)
diff --git a/test/test_downloading.jl b/test/test_jra55_ecco_en4_etopo_downloading.jl
similarity index 100%
rename from test/test_downloading.jl
rename to test/test_jra55_ecco_en4_etopo_downloading.jl
diff --git a/test/test_metadata.jl b/test/test_metadata.jl
new file mode 100644
index 000000000..e55fd2fdf
--- /dev/null
+++ b/test/test_metadata.jl
@@ -0,0 +1,143 @@
+include("runtests_setup.jl")
+
+using NumericalEarth.DataWrangling: Column, Linear, Nearest,
+                                    BoundingBox, dataset_location,
+                                    restrict_location, native_grid
+
+using Oceananigans: RectilinearGrid, LatitudeLongitudeGrid, location
+using Oceananigans.Grids: topology, Flat, Bounded, Periodic
+
+@testset "Column construction" begin
+    col = Column(35.1, 50.1)
+    @test col.longitude == 35.1
+    @test col.latitude == 50.1
+    @test col.z === nothing
+    @test col.interpolation isa Linear
+
+    col_nearest = Column(35.1, 50.1; interpolation=Nearest())
+    @test col_nearest.interpolation isa Nearest
+
+    col_z = Column(35.1, 50.1; z=(-400, 0))
+    @test col_z.z == (-400, 0)
+end
+
+@testset "Column isa checks" begin
+    @test Column(0, 0) isa Column
+    @test !(BoundingBox(longitude=(0, 10), latitude=(0, 10)) isa Column)
+    @test !(nothing isa Column)
+end
+
+@testset "restrict_location" begin
+    # Column reduces horizontal locations to Nothing
+    @test restrict_location((Center, Center, Center), Column(0, 0)) == (Nothing, Nothing, Center)
+    @test restrict_location((Face, Center, Center), Column(0, 0)) == (Nothing, Nothing, Center)
+    @test restrict_location((Center, Face, Center), Column(0, 0)) == (Nothing, Nothing, Center)
+    @test restrict_location((Center, Center, Nothing), Column(0, 0)) == (Nothing, Nothing, Nothing)
+
+    # BoundingBox and nothing leave location unchanged
+    bbox = BoundingBox(longitude=(0, 10), latitude=(0, 10))
+    @test restrict_location((Face, Center, Center), bbox) == (Face, Center, Center)
+    @test restrict_location((Center, Center, Center), nothing) == (Center, Center, Center)
+end
+
+@testset "dataset_location fallback" begin
+    # Default fallback returns (Center, Center, Center)
+    @test dataset_location(ECCO2Monthly(), :temperature) == (Center, Center, Center)
+    @test dataset_location(ECCO4Monthly(), :temperature) == (Center, Center, Center)
+
+    # ECCO staggered velocities
+    @test dataset_location(ECCO4Monthly(), :u_velocity) == (Face, Center, Center)
+    @test dataset_location(ECCO4Monthly(), :v_velocity) == (Center, Face, Center)
+
+    # ECCO 2D fields
+    @test dataset_location(ECCO4Monthly(), :free_surface) == (Center, Center, Nothing)
+
+    # Non-ECCO datasets use the generic fallback
+    @test dataset_location(JRA55.RepeatYearJRA55(), :temperature) == (Center, Center, Center)
+end
+
+@testset "location(metadata) with Column region" begin
+    # Column metadata: location is restricted
+    col = Column(35.1, 50.1)
+    md = Metadatum(:temperature; dataset=ECCO4Monthly(), region=col)
+    @test location(md) == (Nothing, Nothing, Center)
+
+    # ECCO velocity + Column: horizontal locations dropped
+    md_u = Metadatum(:u_velocity; dataset=ECCO4Monthly(), region=col)
+    @test location(md_u) == (Nothing, Nothing, Center)
+
+    # ECCO 2D field + Column
+    md_fs = Metadatum(:free_surface; dataset=ECCO4Monthly(), region=col)
+    @test location(md_fs) == (Nothing, Nothing, Nothing)
+
+    # No region: full dataset location
+    md_full = Metadatum(:u_velocity; dataset=ECCO4Monthly())
+    @test location(md_full) == (Face, Center, Center)
+
+    # BoundingBox: full dataset location
+    bbox = BoundingBox(longitude=(0, 10), latitude=(0, 10))
+    md_bbox = Metadatum(:u_velocity; dataset=ECCO4Monthly(), region=bbox)
+    @test location(md_bbox) == (Face, Center, Center)
+end
+
+@testset "native_grid with Column region" begin
+    col = Column(35.1, 50.1)
+    md = Metadatum(:temperature; dataset=ECCO4Monthly(), region=col)
+    grid = native_grid(md)
+
+    @test grid isa RectilinearGrid
+    @test topology(grid) == (Flat, Flat, Bounded)
+    _, _, Nz, _ = size(md)
+    @test size(grid) == (1, 1, Nz)
+end
+
+@testset "native_grid without region" begin
+    md = Metadatum(:temperature; dataset=ECCO4Monthly())
+    grid = native_grid(md)
+
+    @test grid isa LatitudeLongitudeGrid
+    Nx, Ny, Nz, _ = size(md)
+    @test size(grid) == (Nx, Ny, Nz)
+end
+
+@testset "native_grid with BoundingBox region" begin
+    bbox = BoundingBox(longitude=(0, 10), latitude=(0, 10))
+    md = Metadatum(:temperature; dataset=ECCO4Monthly(), region=bbox)
+    grid = native_grid(md)
+
+    @test grid isa LatitudeLongitudeGrid
+    # Grid should be smaller than the full global grid
+    Nx_full, Ny_full, _, _ = size(md)
+    Nx, Ny, Nz = size(grid)
+    @test Nx < Nx_full
+    @test Ny < Ny_full
+end
+
+@testset "Metadata region keyword" begin
+    # region keyword replaces BoundingBox
+    col = Column(35.1, 50.1)
+    md = Metadatum(:temperature; dataset=ECCO4Monthly(), region=col)
+    @test md.region isa Column
+    @test md.region.longitude == 35.1
+
+    bbox = BoundingBox(longitude=(0, 10), latitude=(0, 10))
+    md2 = Metadatum(:temperature; dataset=ECCO4Monthly(), region=bbox)
+    @test md2.region isa BoundingBox
+
+    md3 = Metadatum(:temperature; dataset=ECCO4Monthly())
+    @test md3.region === nothing
+end
+
+@testset "Metadata iteration propagates region" begin
+    col = Column(35.1, 50.1)
+    md = Metadata(:temperature; dataset=ECCO4Monthly(), region=col,
+                  start_date=DateTime(1992, 1, 1), end_date=DateTime(1992, 3, 1))
+
+    for sub_md in md
+        @test sub_md.region === col
+    end
+
+    @test first(md).region === col
+    @test last(md).region === col
+    @test md[1].region === col
+end