Read ECCO-Darwin field from NASA Earthdata

Project.toml

@@ -17,6 +17,7 @@ ImageMorphology = "787d08f9-d448-5407-9aad-5290dd7ab264"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
+MeshArrays = "cb8c808f-1acf-59a3-9d2b-6e38d009f683"
 NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"

src/DataWrangling/ECCO/ECCO.jl

@@ -2,6 +2,7 @@ module ECCO
 
 export ECCOMetadata, ECCO_field, ECCO_mask, ECCO_immersed_grid, adjusted_ECCO_tracers, initialize!
 export ECCO2Monthly, ECCO4Monthly, ECCO2Daily
+export ECCO4DarwinMonthly
 export ECCORestoring, LinearlyTaperedPolarMask
 
 using ClimaOcean
@@ -32,6 +33,7 @@ end
 include("ECCO_metadata.jl")
 include("ECCO_mask.jl")
 include("ECCO_restoring.jl")
+include("ECCO_darwin.jl")
 
 # Vertical coordinate
 const ECCO_z = [
@@ -121,6 +123,26 @@ function empty_ECCO_field(metadata::ECCOMetadata;
     return Field{loc...}(grid)
 end
 
+"""
+    retrieve_data(metadata, path)
+
+Retrieve data from `path` according to `metadata`.
+"""
+function retrieve_data(metadata, path)
+    ds = Dataset(path)
+    shortname = short_name(metadata)
+
+    if variable_is_three_dimensional(metadata)
+        data = ds[shortname][:, :, :, 1]
+        data = reverse(data, dims=3)
+    else
+        data = ds[shortname][:, :, 1]
+    end        
+
+    close(ds)
+    return data
+end
+
 """
     ECCO_field(metadata::ECCOMetadata;
                architecture = CPU(),
@@ -162,17 +184,7 @@ function ECCO_field(metadata::ECCOMetadata;
 
     download_dataset(metadata)
     path = metadata_path(metadata)
-    ds = Dataset(path)
-    shortname = short_name(metadata)
-
-    if variable_is_three_dimensional(metadata)
-        data = ds[shortname][:, :, :, 1]
-        data = reverse(data, dims=3)
-    else
-        data = ds[shortname][:, :, 1]
-    end        
-
-    close(ds)
+    data = retrieve_data(metadata, path)
 
     # Convert data from Union(FT, missing} to FT
     FT = eltype(field)
@@ -188,7 +200,7 @@ function ECCO_field(metadata::ECCOMetadata;
     # ECCO4 data is on a -180, 180 longitude grid as opposed to ECCO2 data that
     # is on a 0, 360 longitude grid. To make the data consistent, we shift ECCO4
     # data by 180 degrees in longitude
-    if metadata.version isa ECCO4Monthly 
+    if metadata.version isa ECCO4Monthly || metadata.version isa ECCO4DarwinMonthly
         Nx = size(data, 1)
         if variable_is_three_dimensional(metadata)
             shift = (Nx ÷ 2, 0, 0)

src/DataWrangling/ECCO/ECCO_darwin.jl

@@ -0,0 +1,91 @@
+using MeshArrays
+
+Base.size(::ECCOMetadata{<:AbstractCFDateTime, <:ECCO4DarwinMonthly})   = (720, 360, 50, 1)
+Base.size(data::ECCOMetadata{<:Any, <:ECCO4DarwinMonthly}) = (720, 360, 50, length(data.dates))
+
+# File name generation specific to each Dataset version
+function metadata_filename(metadata::ECCOMetadata{<:AbstractCFDateTime, <:ECCO4DarwinMonthly})
+    shortname = short_name(metadata)
+
+    reference_date = DateTimeProlepticGregorian(1992, 1, 1, 12, 0, 0)
+    timestep_size  = 3600
+
+    iternum = Dates.value((metadata.dates - reference_date) / (timestep_size * 1e3))
+    iterstr = string(iternum, pad=10)
+
+    return shortname * "." * iterstr * ".data"
+end
+
+# Convenience functions
+short_name(data::ECCOMetadata{<:Any, <:ECCO4DarwinMonthly}) = ECCO_darwin_short_names[data.name]
+
+metadata_url(prefix, m::ECCOMetadata{<:Any, <:ECCO4DarwinMonthly}) = prefix * "/" * short_name(m) * "/" * metadata_filename(m)
+
+location(::ECCOMetadata{<:Any, <:ECCO4DarwinMonthly}) = (Center, Center, Center)
+
+variable_is_three_dimensional(::ECCOMetadata{<:Any, <:ECCO4DarwinMonthly}) = true
+
+ECCO_darwin_short_names = Dict(
+    :DIC => "DIC",
+    :ALK => "ALK",
+    :PO₄ => "PO4",
+    :NO₃ => "NO3",
+    :DOP => "DOP",
+    :POP => "POP",
+    :Fe  => "FeT",
+    :Siᵀ => "SiO2",
+)
+
+ECCO_darwin_scale_factor = Dict(
+    :DIC => 1e-3,
+    :ALK => 1e-3,
+    :PO₄ => 1e-3,
+    :NO₃ => 1e-3,
+    :DOP => 1e-3,
+    :POP => 1e-3,
+    :Fe  => 1e-3,
+    :Siᵀ => 1e-3,
+)
+
+# URLs for the ECCO datasets specific to each version
+urls(::ECCOMetadata{<:Any, <:ECCO4DarwinMonthly}) = "https://ecco.jpl.nasa.gov/drive/files/ECCO2/LLC90/ECCO-Darwin/monthly/"
+
+ECCO_darwin_native_grid(::ECCO4DarwinMonthly) = GridSpec("LatLonCap", MeshArrays.GRID_LLC90)
+ECCO_darwin_native_size(::ECCO4DarwinMonthly) = (90, 1170, 50)
+
+"""
+    ECCO_darwin_model_data(metafile)
+
+Read a ECCO4DarwinMonthly data file and regrid using MeshArrays on to regular lat-lon grid
+"""
+function ECCO_darwin_model_data(metadata, path)
+    native_size = ECCO_darwin_native_size(metadata.version)
+    native_grid = ECCO_darwin_native_grid(metadata.version)
+    native_data = zeros(Float32, prod(native_size)) # Native LLC90 grid at precision of the input binary file
+
+    read!(path, native_data)
+    native_data = bswap.(native_data)
+
+    meshed_data = read(reshape(native_data, native_size...), native_grid)
+    Nx, Ny, Nz, _ = size(metadata)
+    coeffs = interpolation_setup()
+    data = zeros(Float32, Nx, Ny, Nz) # Native LLC90 grid at precision of the input binary file
+
+    # Interpolate each layer
+    for k in 1:Nz
+        i, j, c = Interpolate(meshed_data[:, k], coeffs)
+        data[:, :, k] = c
+    end
+
+    # Reverse the z-axis
+    data = reverse(data, dims=3)
+
+    # Fill NaNs in Antarctica with zeros
+    data[isnan.(data)] .= 0.f0
+
+    scale_factor = ECCO_darwin_scale_factor[metadata.name]
+    # Scale data according to metadata.scale_factor
+    return data .* scale_factor
+end
+
+retrieve_data(metadata::ECCOMetadata{<:Any, <:ECCO4DarwinMonthly}, path) = ECCO_darwin_model_data(metadata, path)

src/DataWrangling/ECCO/ECCO_mask.jl

@@ -18,7 +18,7 @@ function ECCO_mask(metadata, architecture = CPU();
 
     # ECCO4 has zeros in place of the missing values, while
     # ECCO2 expresses missing values with values < -1e5
-    if metadata.version isa ECCO4Monthly 
+    if metadata.version isa ECCO4Monthly || metadata.version isa ECCO4DarwinMonthly
         _set_mask! = _set_ECCO4_mask!
     else
         _set_mask! = _set_ECCO2_mask!

src/DataWrangling/ECCO/ECCO_metadata.jl

@@ -15,6 +15,8 @@ struct ECCO2Monthly end
 struct ECCO2Daily end
 struct ECCO4Monthly end
 
+struct ECCO4DarwinMonthly end
+
 """
     struct ECCOMetadata{D, V}
 

test/test_ecco.jl

@@ -5,7 +5,7 @@ using Dates
 using ClimaOcean
 
 using ClimaOcean.ECCO
-using ClimaOcean.ECCO: ECCO_field, metadata_path, ECCO_times
+using ClimaOcean.ECCO: ECCO_field, metadata_path, ECCO_times, ECCO4DarwinMonthly
 using ClimaOcean.DataWrangling: NearestNeighborInpainting
 
 using Oceananigans.Grids: topology
@@ -138,6 +138,11 @@ end
             set!(field, ECCOMetadata(:salinity))
             true
         end
+
+        @test begin
+            set!(field, ECCOMetadata(:DIC; version=ECCO4DarwinMonthly()))
+            true
+        end
     end
 end