From fbaee80d15de833b61ead92680dd37044e892729 Mon Sep 17 00:00:00 2001 From: MicheleCeresoli Date: Fri, 5 Jan 2024 14:54:56 +0100 Subject: [PATCH] Reworked EOP parsing and loading --- src/IERSConventions.jl | 13 +- src/angles.jl | 2 - src/eop/data.jl | 204 +++++++++++++++++++++++++++++++ src/eop/loaders.jl | 129 +++++++++++++++++++ src/eop/parsers.jl | 204 +++++++++++++++++++++++++++++++ src/{eop.jl => eop/retrieval.jl} | 51 -------- 6 files changed, 546 insertions(+), 57 deletions(-) delete mode 100644 src/angles.jl create mode 100644 src/eop/data.jl create mode 100644 src/eop/loaders.jl create mode 100644 src/eop/parsers.jl rename src/{eop.jl => eop/retrieval.jl} (70%) diff --git a/src/IERSConventions.jl b/src/IERSConventions.jl index 9edd002..ab3a808 100644 --- a/src/IERSConventions.jl +++ b/src/IERSConventions.jl @@ -1,14 +1,16 @@ module IERSConventions - using Tempo - + using DelimitedFiles using PrecompileTools using ReferenceFrameRotations using StaticArrays + using Tempo + using JSMDInterfaces.Math: AbstractInterpolationMethod, interpolate + using JSMDUtils.Math: InterpAkima, arcsec2rad + # Basic definitions include("models.jl") - include("angles.jl") include("poisson.jl") # Fundamental arguments @@ -32,7 +34,10 @@ module IERSConventions include("polar.jl") # EOP functions - include("eop.jl") + include("eop/data.jl") + include("eop/parsers.jl") + include("eop/loaders.jl") + include("eop/retrieval.jl") # Rotation functions include("rotations.jl") diff --git a/src/angles.jl b/src/angles.jl deleted file mode 100644 index 7c5920a..0000000 --- a/src/angles.jl +++ /dev/null @@ -1,2 +0,0 @@ - -arcsec2rad(x) = x/648000*π \ No newline at end of file diff --git a/src/eop/data.jl b/src/eop/data.jl new file mode 100644 index 0000000..12832ce --- /dev/null +++ b/src/eop/data.jl @@ -0,0 +1,204 @@ + +export eop_filename + +# EOP Data +# ============================ + +""" + NutationCorrections + +Container to hold the nutation corrections and celestial pole offsets associated to a +given IERS model (e.g., 2000A) + +### Fields +- `δX, δY`: Celestial pole offsets referred to the model IAU2000A (rad). +- `δΔψ, δΔϵ`: Nutation corrections in longitude and obliquity (rad) + +### See also +See also [`EOPData`](@ref). +""" +struct NutationCorrections{T} + δX::Vector{T} + δY::Vector{T} + δΔψ::Vector{T} + δΔϵ::Vector{T} +end + +function NutationCorrections(::Type{T}=Float64) where T + return NutationCorrections(T[], T[], T[], T[]) +end + + +""" + EOPData{T} + +EOP Data container. Data is parameterised by time expressed in Terrestrial Time (TT) days. + +### Fields +- `filename` : File where the EOP data are stored. +- `days_UTC`: UTC days since J2000. +- `days_TT` : TT days since J2000. +- `UT1_TT`: UT1 minus TT offset, in seconds. +- `LOD`: Length of day offset (s). +- `xp, yp`: Polar motion with respect to the crust (rad). +- `nut1996, nut2003, nut2006`: nutation and CIP corrections. +""" +mutable struct EOPData{T} + + filename::String + + days_UTC::Vector{T} + days_TT::Vector{T} + + UT1_TT::Vector{T} + LOD::Vector{T} + + xp::Vector{T} + yp::Vector{T} + + nut1996::NutationCorrections{T} + nut2003::NutationCorrections{T} + nut2010::NutationCorrections{T} + +end + +function EOPData(::Type{T}=Float64) where T + return EOPData( + "", T[], T[], T[], T[], T[], T[], + NutationCorrections(T), NutationCorrections(T), NutationCorrections(T) + ) +end + +function Base.show(io::IO, eop::EOPData) + + if isempty(eop.days_UTC) + println(io, "EOPData()") + else + println( + io, + "EOPData(filename=\"$(eop.filename), \"beg=\"$(eop.days_UTC[1]) (UTC)\", "* + "end=\"$(eop.days_UTC[end]) (UTC)\")" + ) + end + +end + +""" + set_nutation_corr!(eop::EOPData, m::IERSModel, nc::NutationCorrections) + +Set the nutation and cip corrections associated to the IERS model `m` to `nc`. +""" +set_nutation_corr!(eop::EOPData, ::IERS1996, nc::NutationCorrections) = eop.nut1996 = nc +set_nutation_corr!(eop::EOPData, ::IERS2003, nc::NutationCorrections) = eop.nut2003 = nc +set_nutation_corr!(eop::EOPData, ::IERS2010, nc::NutationCorrections) = eop.nut2010 = nc + +""" + eop_filename(eop::EOPData) + +Get the loaded Earth Orientation Parameters (EOP) filename. +""" +function eop_filename(eop::EOPData) + + if isempty(eop.filename) + throw(ErrorException("Unable to retrieve filename, no EOP data has been loaded.")) + end + + return eop.filename +end + + +# Interpolators +# ============================ + +""" + NutCorrectionsInterpolator + +Container to store the interpolators for the nutation corrections and celestial pole offsets +associated to a given IERS model. + +### See also +See also [`EOPInterpolator`](@ref) +""" +struct NutCorrectionsInterpolator{T <: AbstractInterpolationMethod} + δX::T + δY::T + δΔψ::T + δΔϵ::T +end + +function NutCorrectionsInterpolator(::Type{T}) where T + NutCorrectionsInterpolator( + _akima_init(T), _akima_init(T), _akima_init(T), _akima_init(T) + ) +end + +function NutCorrectionsInterpolator(days_tt, nc::NutationCorrections) + NutCorrectionsInterpolator( + InterpAkima(days_tt, nc.δX), + InterpAkima(days_tt, nc.δY), + InterpAkima(days_tt, nc.δΔψ), + InterpAkima(days_tt, nc.δΔϵ) + ) +end + +""" + EOPInterpolator + +Container to store the interpolators for the loaded EOP data. +""" +mutable struct EOPInterpolator{T <: AbstractInterpolationMethod} + + init::Bool + + xp::T + yp::T + + ut1_tt::T + lod::T + + nut1996::NutCorrectionsInterpolator{T} + nut2003::NutCorrectionsInterpolator{T} + nut2010::NutCorrectionsInterpolator{T} + +end + +function EOPInterpolator(::Type{T}=Float64) where T + EOPInterpolator( + false, + _akima_init(T), _akima_init(T), _akima_init(T), _akima_init(T), + NutCorrectionsInterpolator(T), + NutCorrectionsInterpolator(T), + NutCorrectionsInterpolator(T) + ) +end + +function Base.show(io::IO, i::EOPInterpolator) + println(io, "EOPInterpolator(init=$(i.init))") +end + +# Function to initialise a dummy Akima interpolator +_akima_init(::Type{T}) where T = InterpAkima([1, 2, 3, 4, 5, 6], zeros(T, 6)) + + +# IERS Constants +# ============================ + +""" + IERS_EOP_DATA + +Earth Orientation Parameters Data. + +### See also +See also [`EOPData`](@ref). +""" +const IERS_EOP_DATA = EOPData(); + +""" + IERS_EOP + +Earth Orientation Parameters (EOP) interpolators. + +### See also +See also [`load_eop!`](@ref) and [`EOPInterpolator`](@ref). +""" +const IERS_EOP = EOPInterpolator(); \ No newline at end of file diff --git a/src/eop/loaders.jl b/src/eop/loaders.jl new file mode 100644 index 0000000..5410eee --- /dev/null +++ b/src/eop/loaders.jl @@ -0,0 +1,129 @@ + +export eop_load_data! + +""" + eop_load_data!(filename, m::IERS2010) + +Initialise the Earth Orientation Parameters (EOP) from a dedicated JSMD `.eop.dat` file. + +!!! note + Currently, only EOP data associated to the IAU2006/2000A precession-nutation model + is supported. +""" +function eop_load_data!(filename, m::IERS2010) + + # Set the EOP data + eop_set_data!(filename, m) + + # Initialise the interpolators + IERS_EOP.init = true + + # Set polar motion + IERS_EOP.xp = InterpAkima(IERS_EOP_DATA.days_TT, IERS_EOP_DATA.xp) + IERS_EOP.yp = InterpAkima(IERS_EOP_DATA.days_TT, IERS_EOP_DATA.yp) + + IERS_EOP.lod = InterpAkima(IERS_EOP_DATA.days_TT, IERS_EOP_DATA.LOD) + IERS_EOP.ut1_tt = InterpAkima(IERS_EOP_DATA.days_TT, IERS_EOP_DATA.UT1_TT) + + # Set nutation correction interpolators + # TODO: ensure these are set! + # nci1996 = NutCorrectionsInterpolator(IERS_EOP_DATA.days_TT, IERS_EOP_DATA.nut1996) + # nci2003 = NutCorrectionsInterpolator(IERS_EOP_DATA.days_TT, IERS_EOP_DATA.nut2003) + nci2010 = NutCorrectionsInterpolator(IERS_EOP_DATA.days_TT, IERS_EOP_DATA.nut2010) + + # IERS_EOP.nut1996 = nci1996 + # IERS_EOP.nut2003 = nci2003 + IERS_EOP.nut2010 = nci2010 + + @info "EOP initialised from file `$(filename)`." + nothing + +end + + +""" + eop_set_data!(filename, m::IERSModel) + +Set Earth Orientation Parameters (EOP) to be used for frames transformations from a JSMD +`.eop.dat` file, given a reference model `m`. +""" +function eop_set_data!(filename, m::IERSModel) + + oldfile = IERS_EOP_DATA.filename + days_utc, days_tt, xp, yp, ut1_tt, lod, δX, δY, δΔψ, δΔϵ = eop_read_data(filename) + + if (!isempty(IERS_EOP_DATA.days_UTC)) + @warn "Existing EOP data from \'$(oldfile)\' will be overwritten by \'$(filename)\'." + end + + # Set data time stamps + IERS_EOP_DATA.filename = filename + IERS_EOP_DATA.days_UTC = days_utc + IERS_EOP_DATA.days_TT = days_tt + + # Set TT-to-UT1 offset and LOD + IERS_EOP_DATA.UT1_TT = ut1_tt + IERS_EOP_DATA.LOD = lod + + # Set polar motion coordinates + IERS_EOP_DATA.xp = xp + IERS_EOP_DATA.yp = yp + + # Set the nutation corrections + nc = NutationCorrections(δX, δY, δΔψ, δΔϵ) + set_nutation_corr!(IERS_EOP_DATA, m, nc) + + # TODO: convert these corrections to those for the other models! + nothing + +end + + +""" + eop_read_data(filename) + +Read Earth Orientation Parameters (EOP) from a dedicated JSMD `.eop.dat` file. + +!!! note + JSMD's `.eop.dat` EOP files are meant to have a fixed structure to ease the + retrieval of the relevant EOP data. +""" +function eop_read_data(filename) + + if !endswith(filename, "eop.dat") + throw( + ArgumentError( + "EOP reader support only '.eop.dat' files! Please prepare " * + "the data with \'eop_parse_csv\' and retry." + ) + ) + end + + # Load the file + data = readdlm(filename; header=false) + + # Retrieve UT1-UTC + Δut1 = @view(data[:, 4]) + + # Convert UTC days to TT seconds + days_utc = @view(data[:, 1]) + days_tai = map(t->Tempo.utc2tai(Tempo.DJ2000, t)[2], days_utc) + days_tt = days_tai .+ Tempo.OFFSET_TAI_TT/Tempo.DAY2SEC + + # Compute UT1-TT, in seconds + sec_ut1 = days_utc/Tempo.DAY2SEC + Δut1 + ut1_tt = sec_ut1 - days_tt*Tempo.DAY2SEC + + # Retrieve other quantities + xp = @view(data[:, 2]) + yp = @view(data[:, 3]) + lod = @view(data[:, 5]) + + δX = data[:, 6] + δY = data[:, 7] + δΔψ = data[:, 6] + δΔϵ = data[:, 7] + + return days_utc, days_tt, xp, yp, ut1_tt, lod, δX, δY, δΔψ, δΔϵ + +end \ No newline at end of file diff --git a/src/eop/parsers.jl b/src/eop/parsers.jl new file mode 100644 index 0000000..2d54699 --- /dev/null +++ b/src/eop/parsers.jl @@ -0,0 +1,204 @@ + +export eop_parse_csv + +""" + eop_parse_csv(m::IERSModel, inputfile, outputfile) + + +Parse CSV files containing IERS EOP data and extracts the relevant information to a dedicated +JSMD `.eop.dat` file. Supported formats are the EOP C04 series and the Rapid Data +prediction (finals). + +!!! note + The `outputfile` name should not include the file extension, which is automatically + added by this function. + +!!! note + Depending on the type of file, either the CIP or the nutation corrections may be present. + This function applies a conversion algorithm to automatically retrieve the missing data. + +!!! warning + The rapid data prediction files store the LOD, CIP and nutation corrections in + milliseconds and milliarcseconds, respectively, whereas the EOP C04 files do not. This + routine automatically detects the correct unit of measure by analysing the format of the + input filename. Thus, the filename should be kept equal to the one released by the IERS. + +### References +# TODO: add websites +""" +function eop_parse_csv(m::IERSModel, inputfile, outputfile) + + # Check whether the file is a Combined Series (C04) or the Rapid Data Series. + isfinal = startswith(splitdir(inputfile)[2], "finals") + + # Multiplicative factor to bring all quantities to seconds/arcseconds, since the in the + # rapid data series the LOD is given in milliseconds, and the dX, dY, dPsi and dEps + # corrections are in milliarcseconds. + fct = isfinal ? 1e-3 : 1.0 + + # Read the data in the file + data, labels = readdlm(inputfile, ';'; header=true) + headers = collect(String.(labels[:])) + + # Instead of relying on fixed column indexes, we retrieve the column ID associated + # to a specific quantity by checking the header strings. If the required string is + # unavailable, then the CSV file is not supported! + + cols = ["MJD", "LOD", "UT1-UTC", "x_pole", "y_pole", "dX", "dY", "dPsi", "dEpsilon"] + + idxs = zeros(Int, length(cols)) + for (j, label) in enumerate(cols) + # Check whether a valid index has been found! + idx = findfirst(headers .== label) + if !isnothing(idx) + idxs[j] = idx + else + throw( + ErrorException( + "Unsupported or invalid EOP file. Could not find the '$(label)' column." + ) + ) + end + end + + # Retrieve the last row in which there is available data for ΔUT1, for predictions, + # the other parameters usually end before. + last_row = findlast(!isempty, @view(data[:, idxs[3]])) + + mjd = convert(Vector{Float64}, @view(data[1:last_row, idxs[1]])) + Δut1 = convert(Vector{Float64}, @view(data[1:last_row, idxs[3]])) + + # Convert UTC days to TT centuries (for the corrections conversion) + days_utc = mjd .- Tempo.DMJD + days_tai = map(t->Tempo.utc2tai(Tempo.DJ2000, t)[2], days_utc) + cent_tt = (days_tai .+ Tempo.OFFSET_TAI_TT/Tempo.DAY2SEC)/Tempo.CENTURY2DAY + + # Retrieve pole coordinates + xp = convert(Vector{Float64}, @view(data[1:last_row, idxs[4]])) + yp = convert(Vector{Float64}, @view(data[1:last_row, idxs[5]])) + + # Retrieve columns with optional missing data + raw_lod = @view(data[1:last_row, idxs[2]]) + raw_δX = @view(data[1:last_row, idxs[6]]) + raw_δY = @view(data[1:last_row, idxs[7]]) + raw_δΔψ = @view(data[1:last_row, idxs[8]]) + raw_δΔϵ = @view(data[1:last_row, idxs[9]]) + + # Fill missing LOD elements + lod = fct*fill_eop_data(raw_lod) + + k = π/648000 + + # Retrieve the latest rows with valid data + lrow_nut = findlast(!isempty, raw_δΔψ) + if isnothing(lrow_nut) + # Need to convert CIP into nutation corrections + δX = fct*fill_eop_data(raw_δX) + δY = fct*fill_eop_data(raw_δY) + + corr = map((t, dx, dy)->δcip_to_δnut(m, t, dx*k, dy*k), cent_tt, δX, δY) + δΔψ, δΔϵ = map(x->x[1]/k, corr), map(x->x[2]/k, corr) + + else + # Need to convert nutation into CIP corrections + δΔψ = fct*fill_eop_data(raw_δΔψ) + δΔϵ = fct*fill_eop_data(raw_δΔϵ) + + corr = map((t, dp, de)->δnut_to_δcip(m, t, dp*k, de*k), cent_tt, δΔψ, δΔϵ) + δΔψ, δΔϵ = map(x->x[1]/k, corr), map(x->x[2]/k, corr) + end + + # Write the EOP data to the desired file + eop_write_data( + hcat( + days_utc, xp, yp, Δut1, + round.(lod; digits=7), round.(δX; digits=7), round.(δY; digits=7), + round.(δΔψ; digits=7), round.(δΔϵ; digits=7) + ), + outputfile + ) + +end + + +""" + eop_write_data(data, output_filename) + +Write the EOP data stored in the matrix `data` to a dedicated JSMD `.eop-dat` file. + +!!! note + The `output_filename` should not include the file extension, which is automatically + added by this function. + +### See also +See also [`eop_parse_csv`](@ref). +""" +function eop_write_data(data, output_filename) + writedlm(output_filename * ".eop.dat", data) + @info "IERS EOP file '$(filename)' converted to '$(output_filename).eop.dat'." + nothing +end + + +function fill_eop_data(raw_data) + + # Fill all the raws with empty data with the latest valid EOP element and + # return the output as a Float64 vector. + + data = zeros(size(raw_data)) + + # Find index of the latest raw with valid data + lrow = findlast(!isempty, raw_data) + isnothing(lrow) && return data + + # Fill the missing elements with the latest available data + data[1:lrow] = raw_data[1:lrow] + data[lrow+1:end] .= raw_data[lrow] + + return data + +end + + +# EOP Conversion functions +# ========================================= + +# Function to convert nutation corrections to CIP corrections and viceversa +function δnut_to_δcip(m::IERSModel, t::Number, δΔψ::Number, δΔϵ::Number) + + # Compute the precession angles + ϵ₀, ψₐ, _, χₐ = precession_angles_rot4(m, t) + + # Compute sine\cosine of mean obliquity + se = sin(iers_obliquity(m, t)) + ce = cos(ϵ₀) + + c = ψₐ*ce - χₐ + + # Convert nutation corrections + δx = δΔψ*se + c*δΔϵ + δy = δΔϵ - c*se*δΔψ + + return δx, δy + +end + +function δcip_to_δnut(m::IERSModel, t::Number, δx::Number, δy::Number) + + # Compute the precession angles + ϵ₀, ψₐ, _, χₐ = precession_angles_rot4(m, t) + + # Compute sine\cosine of mean obliquity + se = sin(iers_obliquity(m, t)) + ce = cos(ϵ₀) + + c = ψₐ*ce - χₐ + d = 1 + c^2 + + # Convert CIP corrections + δΔψ = (δx - c*δy)/se/d + δΔϵ = (δy + c*δx)/d + + return δΔψ, δΔϵ + +end diff --git a/src/eop.jl b/src/eop/retrieval.jl similarity index 70% rename from src/eop.jl rename to src/eop/retrieval.jl index a683afa..d8ed0d0 100644 --- a/src/eop.jl +++ b/src/eop/retrieval.jl @@ -1,13 +1,4 @@ -# EOP Data structures -# ========================================= - -# TODO: write me - - -# EOP Retrival functions -# ========================================= - """ eop_δΔψ(m::IERSModel, t::Number) @@ -136,45 +127,3 @@ function offset_tt2ut1(seconds) end -# EOP Conversion functions -# ========================================= - -# Function to convert nutation corrections to CIP corrections and viceversa -function δnut_to_δcip(m::IERSModel, t::Number, δΔψ::Number, δΔϵ::Number) - - # Compute the precession angles - ϵ₀, ψₐ, _, χₐ = precession_angles_rot4(m, t) - - # Compute sine\cosine of mean obliquity - se = sin(orient_obliquity(m, t)) - ce = cos(ϵ₀) - - c = ψₐ*ce - χₐ - - # Convert nutation corrections - δx = δΔψ*se + c*δΔϵ - δy = δΔϵ - c*se*δΔψ - - return δx, δy - -end - -function δcip_to_δnut(m::IERSModel, t::Number, δx::Number, δy::Number) - - # Compute the precession angles - ϵ₀, ψₐ, _, χₐ = precession_angles_rot4(m, t) - - # Compute sine\cosine of mean obliquity - se = sin(orient_obliquity(m, t)) - ce = cos(ϵ₀) - - c = ψₐ*ce - χₐ - d = 1 + c^2 - - # Convert CIP corrections - δΔψ = (δx - c*δy)/se/d - δΔϵ = (δy + c*δx)/d - - return δΔψ, δΔϵ - -end