|
| 1 | +using MeshArrays |
| 2 | +using JLD2 |
| 3 | +using Glob |
| 4 | + |
| 5 | +struct ECCO2DarwinMonthly <:ECCODataset end |
| 6 | +struct ECCO4DarwinMonthly <:ECCODataset end |
| 7 | + |
| 8 | +# URLs for the ECCO datasets specific to each version |
| 9 | +const ECCO4Darwin_url = "https://ecco.jpl.nasa.gov/drive/files/ECCO2/LLC90/ECCO-Darwin/" |
| 10 | +const ECCO2Darwin_url = "https://ecco.jpl.nasa.gov/drive/files/ECCO2/LLC270/ECCO-Darwin_extension/" |
| 11 | + |
| 12 | +Base.size(data::Metadata{<:ECCO4DarwinMonthly}) = (720, 360, 50, length(data.dates)) |
| 13 | +Base.size(::Metadatum{<:ECCO4DarwinMonthly}) = (720, 360, 50, 1) |
| 14 | +Base.size(data::Metadata{<:ECCO2DarwinMonthly}) = (1440, 720, 50, length(data.dates)) |
| 15 | +Base.size(::Metadatum{<:ECCO2DarwinMonthly}) = (1440, 720, 50, 1) |
| 16 | + |
| 17 | +metadata_time_step(::Metadatum{<:ECCO4DarwinMonthly}) = 3600 |
| 18 | +metadata_epoch(::Metadatum{<:ECCO4DarwinMonthly}) = DateTime(1992, 1, 1, 12, 0, 0) |
| 19 | + |
| 20 | +metadata_time_step(::Metadatum{<:ECCO2DarwinMonthly}) = 1200 |
| 21 | +metadata_epoch(::Metadatum{<:ECCO2DarwinMonthly}) = DateTime(1992, 1, 1, 0, 0, 0) |
| 22 | + |
| 23 | +# The whole range of dates in the different dataset datasets |
| 24 | +all_dates(dataset::ECCO4DarwinMonthly, name) = metadata_epoch(dataset) : Month(1) : DateTime(2023, 3, 1) |
| 25 | +all_dates(dataset::ECCO2DarwinMonthly, name) = metadata_epoch(dataset) : Month(1) : DateTime(2025, 5, 1) |
| 26 | + |
| 27 | +# File name generation specific to each Dataset dataset |
| 28 | +""" |
| 29 | + metadata_filename(metadata::Metadatum{<:Union{ECCO2DarwinMonthly, ECCO4DarwinMonthly}}) |
| 30 | +
|
| 31 | +Generate the filename for a given ECCO Darwin dataset and date. |
| 32 | +
|
| 33 | +The filename is constructed using the dataset variable name, and the iteration number is calculated |
| 34 | +from the date and epoch. |
| 35 | +""" |
| 36 | +function metadata_filename(metadata::Metadatum{<:Union{ECCO2DarwinMonthly, ECCO4DarwinMonthly}}) |
| 37 | + shortname = dataset_variable_name(metadata) |
| 38 | + |
| 39 | + reference_date = metadata_epoch(metadata) |
| 40 | + timestep_size = metadata_time_step(metadata) |
| 41 | + |
| 42 | + # Explicitly convert to Int to avoid return of a float |
| 43 | + iternum = Int(Dates.value((metadata.dates - reference_date) / (timestep_size * 1e3))) |
| 44 | + iterstr = string(iternum, pad=10) |
| 45 | + |
| 46 | + return shortname * "." * iterstr * ".data" |
| 47 | +end |
| 48 | + |
| 49 | +# Convenience functions |
| 50 | +default_mask_value(::ECCO4DarwinMonthly) = 0 |
| 51 | +default_mask_value(::ECCO2DarwinMonthly) = 0 |
| 52 | + |
| 53 | +dataset_variable_name(data::Metadata{<:Union{ECCO2DarwinMonthly,ECCO4DarwinMonthly}}) = ECCO_darwin_dataset_variable_names[data.name] |
| 54 | + |
| 55 | +location(::Metadata{<:Union{ECCO2DarwinMonthly, ECCO4DarwinMonthly}}) = (Center, Center, Center) |
| 56 | + |
| 57 | +variable_is_three_dimensional(::Metadata{<:Union{ECCO2DarwinMonthly, ECCO4DarwinMonthly}}) = true |
| 58 | + |
| 59 | +ECCO_darwin_dataset_variable_names = Dict( |
| 60 | + :temperature => "THETA", |
| 61 | + :salinity => "SALTanom", |
| 62 | + :dissolved_inorganic_carbon => "DIC", |
| 63 | + :alkalinity => "ALK", |
| 64 | + :phosphate => "PO4", |
| 65 | + :nitrate => "NO3", |
| 66 | + :dissolved_organic_phosphorus => "DOP", |
| 67 | + :particulate_organic_phosphorus => "POP", |
| 68 | + :dissolved_iron => "FeT", |
| 69 | + :dissolved_silicate => "SiO2", |
| 70 | + :dissolved_oxygen => "O2", |
| 71 | +) |
| 72 | + |
| 73 | +""" |
| 74 | + concentration_units(metadatum::Metadatum{<:Union{ECCO2DarwinMonthly, ECCO4DarwinMonthly}}) |
| 75 | +
|
| 76 | +Set up conversion from the ECCODarwin output data to standard units |
| 77 | + - salinity = SALTanom + 35 |
| 78 | + - biogeochemical tracer concentrations are in uL => umol/L in the output files from Darwin |
| 79 | +""" |
| 80 | +function concentration_units(metadatum::Metadatum{<:Union{ECCO2DarwinMonthly, ECCO4DarwinMonthly}}) |
| 81 | + if dataset_variable_name(metadatum) == "SALTanom" |
| 82 | + return GramPerKilogramMinus35() |
| 83 | + elseif dataset_variable_name(metadatum) != "THETA" |
| 84 | + return MicromolePerLiter() |
| 85 | + else |
| 86 | + return nothing |
| 87 | + end |
| 88 | +end |
| 89 | + |
| 90 | +function default_download_directory(::ECCO4DarwinMonthly) |
| 91 | + path = joinpath(download_ECCO_cache, "v4_darwin", "monthly") |
| 92 | + return mkpath(path) |
| 93 | +end |
| 94 | + |
| 95 | +function default_download_directory(::ECCO2DarwinMonthly) |
| 96 | + path = joinpath(download_ECCO_cache, "v2_darwin", "monthly") |
| 97 | + return mkpath(path) |
| 98 | +end |
| 99 | + |
| 100 | +metadata_url(m::Metadata{<:ECCO4DarwinMonthly}) = ECCO4Darwin_url * "monthly/" * dataset_variable_name(m) * "/" * metadata_filename(m) |
| 101 | +metadata_url(m::Metadata{<:ECCO2DarwinMonthly}) = ECCO2Darwin_url * "monthly/" * dataset_variable_name(m) * "/" * metadata_filename(m) |
| 102 | + |
| 103 | +# Functions for reading the ECCO binary files using MeshArrays |
| 104 | +binary_data_grid(::ECCO4DarwinMonthly) = GridSpec(ID=:LLC90) |
| 105 | +binary_data_size(::ECCO4DarwinMonthly) = (90, 1170, 50) |
| 106 | +binary_data_grid(::ECCO2DarwinMonthly) = GridSpec(ID=:LLC270) |
| 107 | +binary_data_size(::ECCO2DarwinMonthly) = (270, 3510, 50) |
| 108 | + |
| 109 | +longitude_interfaces(::ECCO4DarwinMonthly) = (-180, 180) |
| 110 | + |
| 111 | +""" |
| 112 | + retrieve_data(metadata::Metadatum{<:ECCO4DarwinMonthly}) |
| 113 | +
|
| 114 | +Read a ECCO4DarwinMonthly data file and regrid using MeshArrays on to regular lat-lon grid |
| 115 | +""" |
| 116 | +function retrieve_data(metadata::Metadatum{<:Union{ECCO4DarwinMonthly, ECCO2DarwinMonthly}}) |
| 117 | + native_size = binary_data_size(metadata.dataset) |
| 118 | + native_grid = binary_data_grid(metadata.dataset) |
| 119 | + native_data = zeros(Float32, prod(native_size)) # Native LLC grid at precision of the input binary file |
| 120 | + |
| 121 | + read!(metadata_path(metadata), native_data) |
| 122 | + native_data = bswap.(native_data) |
| 123 | + |
| 124 | + meshed_data = read(reshape(native_data, native_size...), native_grid) |
| 125 | + Nx, Ny, Nz, _ = size(metadata) |
| 126 | + data = zeros(Float32, Nx, Ny, Nz) # Native LLC grid at precision of the input binary file |
| 127 | + mask = zeros(Float32, Nx, Ny, Nz) |
| 128 | + |
| 129 | + # Download the native grid data from MeshArrays repo (only if not in already in datadeps) |
| 130 | + native_grid_coords = GridLoad(native_grid; option="full") |
| 131 | + |
| 132 | + # Check if the interpolation coefficients are already calculated |
| 133 | + interp_file = joinpath(dirname(metadata_path(metadata)),"native_interp_coeffs.jld2") |
| 134 | + if !isfile(interp_file) |
| 135 | + # Calculate coefficients to interpolate from native grid to regular lat-lon grid (as in the ECCO netcdf files) |
| 136 | + resolution_X = 360/Nx |
| 137 | + resolution_Y = 180/Ny |
| 138 | + |
| 139 | + # Regular lat-lon grid |
| 140 | + longitudes = longitude_interfaces(metadata.dataset) |
| 141 | + latitudes = latitude_interfaces(metadata.dataset) |
| 142 | + lon = [i for i = longitudes[1]+resolution_X/2:resolution_X:longitudes[2]-resolution_X/2, |
| 143 | + j = latitudes[1]+resolution_Y/2:resolution_Y:latitudes[2]-resolution_Y/2] |
| 144 | + lat = [j for i = longitudes[1]+resolution_X/2:resolution_X:longitudes[2]-resolution_X/2, |
| 145 | + j = latitudes[1]+resolution_Y/2:resolution_Y:latitudes[2]-resolution_Y/2] |
| 146 | + |
| 147 | + # Interpolation factors for the native grid (writes out to a file "interp_file" for later use) |
| 148 | + coeffs = interpolation_setup(; Γ=native_grid_coords, lat, lon, filename=interp_file) |
| 149 | + else |
| 150 | + # Read the coefficients from the file that was previously calculated |
| 151 | + coeffs = interpolation_setup(interp_file) |
| 152 | + end |
| 153 | + |
| 154 | + # Read continental mask on the native model grid |
| 155 | + native_grid_fac_center = GridLoadVar("hFacC", native_grid) |
| 156 | + |
| 157 | + # Interpolate each masked layer on to the native lat lon grid |
| 158 | + for k in 1:Nz |
| 159 | + i, j, c = MeshArrays.Interpolate( |
| 160 | + meshed_data[:, k], |
| 161 | + coeffs, |
| 162 | + ) |
| 163 | + data[:, :, k] = c |
| 164 | + i, j, c = MeshArrays.Interpolate( |
| 165 | + land_mask(native_grid_fac_center[:, k]), |
| 166 | + coeffs, |
| 167 | + ) |
| 168 | + mask[:, :, k] = c |
| 169 | + end |
| 170 | + |
| 171 | + # Reverse the z-axis |
| 172 | + data = reverse(data, dims=3) |
| 173 | + mask = reverse(mask, dims=3) |
| 174 | + |
| 175 | + # Fill NaNs in Antarctica with zeros |
| 176 | + data[isnan.(data)] .= default_mask_value(metadata.dataset) |
| 177 | + |
| 178 | + return data .* mask |
| 179 | +end |
0 commit comments