From bc842c31add37e9014d23ed17dcd7c42956a20e1 Mon Sep 17 00:00:00 2001 From: Steve Marshall Date: Wed, 10 Dec 2025 13:25:06 -0500 Subject: [PATCH] Addition of the RefNLPEP2D GNSSRO forward operator This is a new 2D GNSSRO forward operator based on refractivity that simulates a path-integrated variable called nonLocalPseudoExcessPhase. Summary of changes: * New code for the RefNLPEP2D GNSSRO forward operator * Additional GNSSRO utility classes related to ray path generation, GeoVaLs access in C++, and refractivity computation. * Two new variable transforms: NonLocalPseudoExcessPhase and GnssroRefractivityGradient. * Extension of GeoVaLs code to allow reading and writing GeoVaLs where the number of profiles differs from the number of locations. * Extension of the NCEP refractivity error model in ROobserror to handle the nonLocalPseudoExcessPhase variable. * Units tests for the forward operator and variable transforms. --- src/ufo/operators/gnssro/CMakeLists.txt | 5 +- .../gnssro/QC/ufo_roobserror_mod.F90 | 63 +- .../gnssro/RefNLPEP2D/CMakeLists.txt | 22 + .../gnssro/RefNLPEP2D/ObsGnssroRefNLPEP2D.cc | 338 +++++++++++ .../gnssro/RefNLPEP2D/ObsGnssroRefNLPEP2D.h | 87 +++ .../ObsGnssroRefNLPEP2DParameters.h | 80 +++ .../RefNLPEP2D/ObsGnssroRefNLPEP2DTLAD.cc | 538 ++++++++++++++++++ .../RefNLPEP2D/ObsGnssroRefNLPEP2DTLAD.h | 84 +++ src/ufo/operators/gnssro/utils/CMakeLists.txt | 10 + .../operators/gnssro/utils/GnssroGeoVaLs.cc | 124 ++++ .../operators/gnssro/utils/GnssroGeoVaLs.h | 133 +++++ .../gnssro/utils/GnssroProfileExtractor.cc | 78 +++ .../gnssro/utils/GnssroProfileExtractor.h | 56 ++ .../gnssro/utils/GnssroProfileRayPaths.cc | 164 ++++++ .../gnssro/utils/GnssroProfileRayPaths.h | 90 +++ .../gnssro/utils/GnssroProfileSlice.h | 77 +++ .../operators/gnssro/utils/GnssroRayPath.cc | 144 +++++ .../operators/gnssro/utils/GnssroRayPath.h | 142 +++++ .../gnssro/utils/GnssroRayPathGenerator.cc | 42 ++ .../gnssro/utils/GnssroRayPathGenerator.h | 64 +++ .../gnssro/utils/GnssroRayPathOrchestrator.cc | 389 +++++++++++++ .../gnssro/utils/GnssroRayPathOrchestrator.h | 165 ++++++ .../gnssro/utils/GnssroRayPathParameters.cc | 122 ++++ .../gnssro/utils/GnssroRayPathParameters.h | 70 +++ .../gnssro/utils/GnssroRayTrajectory.cc | 89 +++ .../gnssro/utils/GnssroRayTrajectory.h | 78 +++ .../gnssro/utils/RefractivityCalculator.cc | 151 +++++ .../gnssro/utils/RefractivityCalculator.h | 69 +++ .../utils/StraightLineRayPathGenerator.cc | 302 ++++++++++ .../utils/StraightLineRayPathGenerator.h | 113 ++++ src/ufo/ufo_geovals_mod.F90 | 256 +++++++-- src/ufo/variabletransforms/CMakeLists.txt | 2 + .../Cal_GnssroRefractivityGradient.cc | 184 ++++++ .../Cal_GnssroRefractivityGradient.h | 100 ++++ .../Cal_NonLocalPseudoExcessPhase.cc | 165 ++++++ .../Cal_NonLocalPseudoExcessPhase.h | 92 +++ .../unit_tests/operators/CMakeLists.txt | 21 +- .../operators/gnssrorefnlpep2d.yaml | 35 ++ .../variabletransforms/CMakeLists.txt | 22 + ...transforms_gnssrorefractivitygradient.yaml | 81 +++ ...etransforms_nonlocalpseudoexcessphase.yaml | 38 ++ 41 files changed, 4830 insertions(+), 55 deletions(-) create mode 100644 src/ufo/operators/gnssro/RefNLPEP2D/CMakeLists.txt create mode 100644 src/ufo/operators/gnssro/RefNLPEP2D/ObsGnssroRefNLPEP2D.cc create mode 100644 src/ufo/operators/gnssro/RefNLPEP2D/ObsGnssroRefNLPEP2D.h create mode 100644 src/ufo/operators/gnssro/RefNLPEP2D/ObsGnssroRefNLPEP2DParameters.h create mode 100644 src/ufo/operators/gnssro/RefNLPEP2D/ObsGnssroRefNLPEP2DTLAD.cc create mode 100644 src/ufo/operators/gnssro/RefNLPEP2D/ObsGnssroRefNLPEP2DTLAD.h create mode 100644 src/ufo/operators/gnssro/utils/GnssroGeoVaLs.cc create mode 100644 src/ufo/operators/gnssro/utils/GnssroGeoVaLs.h create mode 100644 src/ufo/operators/gnssro/utils/GnssroProfileExtractor.cc create mode 100644 src/ufo/operators/gnssro/utils/GnssroProfileExtractor.h create mode 100644 src/ufo/operators/gnssro/utils/GnssroProfileRayPaths.cc create mode 100644 src/ufo/operators/gnssro/utils/GnssroProfileRayPaths.h create mode 100644 src/ufo/operators/gnssro/utils/GnssroProfileSlice.h create mode 100644 src/ufo/operators/gnssro/utils/GnssroRayPath.cc create mode 100644 src/ufo/operators/gnssro/utils/GnssroRayPath.h create mode 100644 src/ufo/operators/gnssro/utils/GnssroRayPathGenerator.cc create mode 100644 src/ufo/operators/gnssro/utils/GnssroRayPathGenerator.h create mode 100644 src/ufo/operators/gnssro/utils/GnssroRayPathOrchestrator.cc create mode 100644 src/ufo/operators/gnssro/utils/GnssroRayPathOrchestrator.h create mode 100644 src/ufo/operators/gnssro/utils/GnssroRayPathParameters.cc create mode 100644 src/ufo/operators/gnssro/utils/GnssroRayPathParameters.h create mode 100644 src/ufo/operators/gnssro/utils/GnssroRayTrajectory.cc create mode 100644 src/ufo/operators/gnssro/utils/GnssroRayTrajectory.h create mode 100644 src/ufo/operators/gnssro/utils/RefractivityCalculator.cc create mode 100644 src/ufo/operators/gnssro/utils/RefractivityCalculator.h create mode 100644 src/ufo/operators/gnssro/utils/StraightLineRayPathGenerator.cc create mode 100644 src/ufo/operators/gnssro/utils/StraightLineRayPathGenerator.h create mode 100644 src/ufo/variabletransforms/Cal_GnssroRefractivityGradient.cc create mode 100644 src/ufo/variabletransforms/Cal_GnssroRefractivityGradient.h create mode 100644 src/ufo/variabletransforms/Cal_NonLocalPseudoExcessPhase.cc create mode 100644 src/ufo/variabletransforms/Cal_NonLocalPseudoExcessPhase.h create mode 100644 test/testinput/unit_tests/operators/gnssrorefnlpep2d.yaml create mode 100644 test/testinput/unit_tests/variabletransforms/variabletransforms_gnssrorefractivitygradient.yaml create mode 100644 test/testinput/unit_tests/variabletransforms/variabletransforms_nonlocalpseudoexcessphase.yaml diff --git a/src/ufo/operators/gnssro/CMakeLists.txt b/src/ufo/operators/gnssro/CMakeLists.txt index 2a0c19b44..f5cfda7cc 100644 --- a/src/ufo/operators/gnssro/CMakeLists.txt +++ b/src/ufo/operators/gnssro/CMakeLists.txt @@ -7,6 +7,7 @@ add_subdirectory( BndROPP2D ) add_subdirectory( utils ) add_subdirectory( QC ) add_subdirectory( BendMetOffice ) +add_subdirectory( RefNLPEP2D ) PREPEND( _p_refncep_files "operators/gnssro/RefNCEP" ${refncep_src_files} ) PREPEND( _p_refmetoffice_files "operators/gnssro/RefMetOffice" ${refmetoffice_src_files} ) @@ -16,6 +17,7 @@ PREPEND( _p_bndropp2d_files "operators/gnssro/BndROPP2D" ${bndropp2d_src_ PREPEND( _p_utils_files "operators/gnssro/utils" ${utils_src_files} ) PREPEND( _p_qc_files "operators/gnssro/QC" ${qc_src_files} ) PREPEND( _p_metoffice_files "operators/gnssro/BendMetOffice" ${bendmetoffice_src_files} ) +PREPEND( _p_refnlpep2d_files "operators/gnssro/RefNLPEP2D" ${refnlpep2d_files} ) set ( gnssro_src_files ${_p_refncep_files} @@ -26,5 +28,6 @@ set ( gnssro_src_files ${_p_utils_files} ${_p_qc_files} ${_p_metoffice_files} - PARENT_SCOPE + ${_p_refnlpep2d_files} + PARENT_SCOPE ) diff --git a/src/ufo/operators/gnssro/QC/ufo_roobserror_mod.F90 b/src/ufo/operators/gnssro/QC/ufo_roobserror_mod.F90 index d044a1bf7..7ad28645c 100644 --- a/src/ufo/operators/gnssro/QC/ufo_roobserror_mod.F90 +++ b/src/ufo/operators/gnssro/QC/ufo_roobserror_mod.F90 @@ -147,6 +147,9 @@ subroutine ufo_roobserror_prior(self) real(kind_real), allocatable :: obsLocR(:) ! The local radius of curvature at the observation location real(kind_real), allocatable :: obsValue(:) real(kind_real), allocatable :: obsErr(:) +real(kind_real), allocatable :: obsRefr(:) ! observed atmosphericRefractivity (N-units) +real(kind_real), allocatable :: obsNLPEP(:) ! observed nonLocalPseudoExcessPhase (m) +real(kind_real), allocatable :: errScaleFactor(:) ! Factor to map one kind of error to another real(kind_real), allocatable :: averageTemp(:) ! The average model temperature integer(c_int), allocatable :: obsSaid(:) integer(c_int), allocatable :: QCflags(:) @@ -376,8 +379,66 @@ subroutine ufo_roobserror_prior(self) call abor1_ftn(err_msg) end select +case ("nonLocalPseudoExcessPhase") + + ! Get scale factors to convert refractivity errors into non-local pseudo excess phase errors. + ! This part is common to any technique that repurposes a refractivity error model. + allocate(errScaleFactor(nobs)) + if (obsspace_has(self%obsdb, "ObsValue", "atmosphericRefractivity") .AND. & + & obsspace_has(self%obsdb, "ObsValue", "nonLocalPseudoExcessPhase")) then + allocate(obsRefr(nobs)) + allocate(obsNLPEP(nobs)) + call obsspace_get_db(self%obsdb, "ObsValue", "atmosphericRefractivity", obsRefr) + call obsspace_get_db(self%obsdb, "ObsValue", "nonLocalPseudoExcessPhase", obsNLPEP) + + do iob = 1, nobs + errScaleFactor(iob) = obsNLPEP(iob) / obsRefr(iob) + end do + deallocate(obsRefr) + deallocate(obsNLPEP) + else + write(err_msg,*) "ufo_roobserror_mod: could not access observations of & + &atmosphericRefractivity and nonLocalPseudoExcessPhase while computing & + &errors for nonLocalPseudoExcessPhase from refractivity; error cannot & + &be properly scaled" + call abor1_ftn(err_msg) + end if + + select case (trim(self%errmodel)) + + case ("NCEP") + + allocate(obsZ(nobs)) + allocate(obsLat(nobs)) + call obsspace_get_db(self%obsdb, "MetaData", "height", obsZ) + call obsspace_get_db(self%obsdb, "MetaData", "latitude", obsLat) + call refractivity_obserr_NCEP(obsLat, obsZ, nobs, obsErr, QCflags, missing) + do iob = 1, nobs + if (QCflags(iob) .eq. 0) then + ! Scale refractivity error to nonLocalPseudoExcessPhase error. + obsErr(iob) = obsErr(iob) * errScaleFactor(iob) + end if + end do + + write(err_msg,*) 'ufo_roobserror_mod: setting up nonLocalPseudoExcessPhase & + &obs error with NCEP method' + call fckit_log%debug(err_msg) + deallocate(obsZ) + deallocate(obsLat) + ! up date obs error + call obsspace_put_db(self%obsdb, "FortranERR", trim(self%variable), obsErr) + + case default + write(err_msg,*) 'ufo_roobserror_mod: only NCEP error model is available & + &for variable nonLocalPseudoExcessPhase' + call abor1_ftn(err_msg) + end select + + deallocate(errScaleFactor) + case default - call abor1_ftn("ufo_roobserror_prior: variable has to be bending_angle or refractivity") + call abor1_ftn('ufo_roobserror_prior: variable has to be bendingAngle, & + &atmosphericRefractivity, or nonLocalPseudoExcessPhase') end select deallocate(QCflags) diff --git a/src/ufo/operators/gnssro/RefNLPEP2D/CMakeLists.txt b/src/ufo/operators/gnssro/RefNLPEP2D/CMakeLists.txt new file mode 100644 index 000000000..fba608409 --- /dev/null +++ b/src/ufo/operators/gnssro/RefNLPEP2D/CMakeLists.txt @@ -0,0 +1,22 @@ +# (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +set ( refnlpep2d_files + ObsGnssroRefNLPEP2D.h + ObsGnssroRefNLPEP2D.cc + ObsGnssroRefNLPEP2DTLAD.h + ObsGnssroRefNLPEP2DTLAD.cc + ObsGnssroRefNLPEP2DParameters.h +PARENT_SCOPE +) diff --git a/src/ufo/operators/gnssro/RefNLPEP2D/ObsGnssroRefNLPEP2D.cc b/src/ufo/operators/gnssro/RefNLPEP2D/ObsGnssroRefNLPEP2D.cc new file mode 100644 index 000000000..aed28f803 --- /dev/null +++ b/src/ufo/operators/gnssro/RefNLPEP2D/ObsGnssroRefNLPEP2D.cc @@ -0,0 +1,338 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#include +#include +#include +#include + +#include "eckit/exception/Exceptions.h" +#include "ioda/ObsVector.h" + +#include "oops/base/Locations.h" +#include "oops/base/Variables.h" +#include "oops/interface/SampledLocations.h" +#include "oops/util/Logger.h" +#include "oops/util/missingValues.h" + +#include "ufo/GeoVaLs.h" +#include "ufo/ObsDiagnostics.h" +#include "ufo/ObsTraits.h" +#include "ufo/operators/gnssro/RefNLPEP2D/ObsGnssroRefNLPEP2D.h" +#include "ufo/operators/gnssro/utils/GnssroGeoVaLs.h" +#include "ufo/operators/gnssro/utils/GnssroProfileExtractor.h" +#include "ufo/operators/gnssro/utils/GnssroRayPathGenerator.h" +#include "ufo/operators/gnssro/utils/RefractivityCalculator.h" +#include "ufo/SampledLocations.h" + +namespace ufo { + +// ----------------------------------------------------------------------------- +static ObsOperatorMaker makerGnssroRefNLPEP2D_( + "GnssroRefNLPEP2D"); +// ----------------------------------------------------------------------------- + +ObsGnssroRefNLPEP2D::ObsGnssroRefNLPEP2D(const ioda::ObsSpace & odb, + const Parameters_ & parameters) + : ObsOperatorBase(odb) + , odb_(odb) + , varin_() + , opts_(parameters.options.value()) + , rpParams_(opts_.rayPathGenType, opts_.approxRayLength, + opts_.top2D, opts_.res, opts_.nHoriz) + , orchestrator_(odb_, rpParams_) + , refrCalc_(RefractivityCalculator::create(opts_.refrAlgorithm, opts_.useCompress)) +{ + oops::Log::trace() << "ObsGnssroRefNLPEP2D - c'tor body starting" << std::endl; + const std::vector vv{GnssroGeoVaLs::VAR_TEMP, GnssroGeoVaLs::VAR_SPHUM, + GnssroGeoVaLs::VAR_PRES, GnssroGeoVaLs::VAR_GPH}; + varin_.reset(new oops::Variables(vv)); + + // Provide aggregated information about the ray creation process + oops::Log::info() << "ObsGnssroRefNLPEP2D ctor summary: " << orchestrator_ << std::endl; + + oops::Log::trace() << "ObsGnssroRefNLPEP2DParameters: " << parameters.toConfiguration() + << std::endl; + oops::Log::trace() << "ObsGnssroRefNLPEP2D created." << std::endl; +} + +// ----------------------------------------------------------------------------- + +ObsGnssroRefNLPEP2D::~ObsGnssroRefNLPEP2D() { + oops::Log::trace() << "ObsGnssroRefNLPEP2D destroyed" << std::endl; +} + +// ----------------------------------------------------------------------------- + +void ObsGnssroRefNLPEP2D::simulateObs(const GeoVaLs & gv, ioda::ObsVector & ovec, + ObsDiagnostics & d, const QCFlags_t & qc_flags) const { + const char funcName[] = "ObsGnssroRefNLPEP2D::simulateObs"; + oops::Log::trace() << funcName << " started" << std::endl; + + oops::Log::debug() << funcName << ": GeoVaLs has " << gv.nlocs() << " locations and " + << gv.getVars().size() << " variables:" << std::endl; + for (std::size_t varIdx = 0; varIdx < gv.getVars().size(); ++varIdx) + { + const oops::Variable& var = (gv.getVars())[varIdx]; + oops::Log::debug() << var.name() << ": " << gv.nlevs(var) << " levels, " + << gv.nprofiles(var) << " profiles/paths" << std::endl; + } + oops::Log::debug() << funcName << ": ovec has " << ovec.nlocs() << " locations, " + << ovec.size() << " array elements, " << ovec.nvars() << " variables" << std::endl; + + // Confirm our sizing for number of obs is what we expect. + std::size_t totalNumRays = orchestrator_.totalNumRays(); + if (totalNumRays != ovec.nlocs()) + { + oops::Log::error() << funcName << ": size mismatch between ObsVector (" + << ovec.nlocs() << ") and Orchestrator's total number of rays (" + << totalNumRays << ")" << std::endl; + throw eckit::BadValue("size mismatch between ObsVector and Orchestrator in " + "ObsGnssroRefNLPEP2D", Here()); + } + + // Return early if there are no geoVals to operate on. + if (gv.nlocs() == 0) + { + return; + } + + // Extract the model data from the GeoVaLs. + GnssroGeoVaLs ggv(gv); + + // Get geometric height and latitude from the observations. + std::vector latitudes(odb_.nlocs()); + std::vector longitudes(odb_.nlocs()); + std::vector geomHeights(odb_.nlocs()); + odb_.get_db("MetaData", "latitude", latitudes); + odb_.get_db("MetaData", "longitude", longitudes); + odb_.get_db("MetaData", "height", geomHeights); + + // Iterate over RO profiles. + double dmiss = util::missingValue(); + std::size_t ovecIdx = 0; + std::size_t profileIdx = 0; + const GnssroRayPathOrchestrator::Profiles_& profiles = orchestrator_.profiles(); + for (const auto& roProfile : profiles) + { + int seqNum = roProfile->seqNum(); + + // Determine range of observation indices in this profile. + std::size_t startOvecIdx = ovecIdx; + std::size_t endOvecIdx = ovecIdx + roProfile->numRays(); + ovecIdx = endOvecIdx; + + // Get range of model profile indices in this RO profile + std::size_t startProfileIdx = profileIdx; + std::size_t endProfileIdx = startProfileIdx + roProfile->totalSampledLocations(); + profileIdx = endProfileIdx; + + oops::Log::debug() << funcName << ": seqNum=" << seqNum << ": processing profile " + << "with numRays=" << roProfile->numRays() << " ovecIdx range=[" + << startOvecIdx << ":" << endOvecIdx << "), profile-wide profileIdx range=[" + << startProfileIdx << "," << endProfileIdx << "), starting at ovecIdx=" + << ovecIdx - 1 << " and proceeding backwards" << std::endl; + + // Iterate over rays within this profile from top to bottom (reverse order) + for (GnssroProfileRayPaths::const_reverse_iterator rayItr = roProfile->crbegin(); + rayItr != roProfile->crend(); ++rayItr) + { + // Get the integer rayIdx for use in diagnostic messages. + // Get a forward iterator pointing to one position after the reverse iterator. + // Then convert the forward iterator into an integer position. + // This prevents us from even computing a negative rayIdx, which would happen + // if we used an integer index as a reverse iterator. + --ovecIdx; + auto fwdItr = rayItr.base(); + std::size_t rayIdx = std::distance(roProfile->cbegin(), fwdItr) - 1; + const GnssroProfileRayPaths::Ray_& ray = *rayItr; + + // Get the model profile range for this ray + std::size_t rayEndProfileIdx = profileIdx; + std::size_t rayStartProfileIdx = rayEndProfileIdx - ray->numLocations(); + profileIdx = rayStartProfileIdx; + oops::Log::debug() << funcName << ": seqNum=" << seqNum << ", rayIdx= " + << rayIdx << ", numSampLocs=" << ray->numLocations() + << ", ray profileIdx range=[" << rayStartProfileIdx + << "," << rayEndProfileIdx << "), starting at profileIdx=" + << profileIdx << std::endl; + + std::size_t tpNodeIdx = ray->tpNodeIdx(); + + oops::Log::debug() << funcName << ": seqNum=" << seqNum << ", rayIdx=" << rayIdx + << ", ovecIdx=" << ovecIdx << ", tpNodeIdx=" << tpNodeIdx + << ", computing NLPEP over numNodes=" << ray->numNodes() + << ", totalLength=" << ray->totalLength() + << ", tpHgt=" << geomHeights[ovecIdx] << std::endl; + + // Iterate over nodes in the ray. We accumulate: + // fwdNLPEP = Model 2D Non-local pseudo excess phase. This is the forward part + // from integrating N(x,y) over the ray path. + double fwdNLPEP = 0.0; + for (std::size_t nodeIdx = 0; nodeIdx < ray->numNodes(); ++nodeIdx) + { + // Unnormalized segment length for this node of the ray in meters + double segmentLength = ray->segLen(nodeIdx); + + // Get the geometric height of the node and observed latitude of tangent point. + // Compute corresponding geopotential height of this observation sample. + // + // Use the average of the edge heights at all locations except the tangent point. + // The edge average will not find the lowest height at the tangent point. + // Instead, it will use two edges that are both higher than the tangent point. + // So we use the observed geometric height for the tangent point instead. + // + double avgGeomHgt; + if (nodeIdx == tpNodeIdx) { + avgGeomHgt = geomHeights[ovecIdx]; + } else { + float thisGeomHgt = ray->geomHgt(nodeIdx); + float nextGeomHgt = ray->geomHgt(nodeIdx + 1); + avgGeomHgt = static_cast(0.5 * (thisGeomHgt + nextGeomHgt)); + } + double refr = GnssroRayPathOrchestrator::modelRefractivity( + avgGeomHgt, latitudes[ovecIdx], ggv, profileIdx, refrCalc_); + + // Integrate refractivity by the segment length of this ray node + // to get nonLocalPseudoExcessPhase contribution from the horizontally + // varying model data. + double fwdNLPEPincr = (refr * segmentLength); + fwdNLPEP += fwdNLPEPincr; + oops::Log::debug() << funcName << ": seqNum=" << seqNum << ", rayIdx=" << rayIdx + << ", nodeIdx=" << nodeIdx << ", profileIdx=" << profileIdx + << ", geomHgt=" << avgGeomHgt + << ", increment fwdNLPEP=" << fwdNLPEP << ", refr=" << refr + << ", segmentLength=" << segmentLength + << ", NLPEPincr=" << fwdNLPEPincr << std::endl; + + oops::Log::debug() << funcName << ": seqNum=" << seqNum << ", rayIdx=" << rayIdx + << ", nodeIdx=" << nodeIdx << ", ovecIdx=" << ovecIdx + << ": Node complete: refr=" << refr << ", segLen=" << segmentLength + << ", fwdNLPEP=" << fwdNLPEP << std::endl; + + // Only change location for next node if we are not assuming + // spherical symmetry, which means all nodes have the same location. + if (!ray->assumeSphericalSymmetry()) { + ++profileIdx; + } + } // end iteration over nodes in a ray + + // Scale from refractivity N to excess index of refraction (n-1). + // This is a unit conversion from N * m to m. + ovec[ovecIdx] = fwdNLPEP * RefractivityCalculator::N_TO_EXCESS_IOR; + + // Prepare for the next lower ray by setting the profileIdx to point to + // the start of the model profiles for the current ray, which is also the + // end of the range for the next lower ray within the profile. + profileIdx = rayStartProfileIdx; + } // end iteration over rays in a single RO profile + + // Update ob vector index to skip to the next profile. + oops::Log::debug() << funcName << ": seqNum=" << seqNum << " DONE: " + << "ovecIdx=" << ovecIdx << std::endl; + ovecIdx = endOvecIdx; + if (profileIdx != startProfileIdx) { + oops::Log::error() << funcName << ": seqNum=" << seqNum + << ": model profile indexing error: after RO profile processed," + << " profileIdx=" << profileIdx << "; should be same as " + << "startProfileIdx=" << startProfileIdx << std::endl; + } + profileIdx = endProfileIdx; + oops::Log::debug() << funcName << ": seqNum=" << seqNum << " DONE: " + << "skipping to next profile ovecIdx=" << ovecIdx + << ", profileIdx=" << profileIdx << std::endl; + } // end iteration over RO profiles. + + oops::Log::trace() << funcName << " complete" << std::endl; +} +// ----------------------------------------------------------------------------- +ObsGnssroRefNLPEP2D::Locations_ +ObsGnssroRefNLPEP2D::locations() const { + typedef oops::SampledLocations SampledLocations_; + oops::Log::trace() << "ObsGnssroRefNLPEP2D::locations() enter" << std::endl; + + // Get the arrays of latitudes, longitudes, and datetimes used to create the + // sampled locations. This method will allocate the vectors properly. + std::vector lons; + std::vector lats; + std::vector times; + std::size_t numSampledLocs = orchestrator_.getSampledLatLonTimes(lats, lons, times); + + // Set range of sampled locations used by each observed ray in each profile. + std::vector> pathsGroupedByLocation; + orchestrator_.fillPathsGroupedByLocation(odb_.nlocs(), pathsGroupedByLocation); + + SampledLocations_ sampledLocations( + std::make_unique(lons, lats, times, odb_.distribution(), + std::move(pathsGroupedByLocation))); + oops::Log::debug() << "ObsGnssroRefNLPEP2D::locations: creating SampledLocations for " + << odb_.nlocs() << " ob samples to make " << numSampledLocs + << " model paths; SampledLocs/ob=" + << static_cast(numSampledLocs)/static_cast(odb_.nlocs()) + << ", locationsSampledOnceAndInOrder=" + << sampledLocations.sampledLocations().areLocationsSampledOnceAndInOrder() + << std::endl; + + const std::vector>& pathRanges = + sampledLocations.sampledLocations().pathsGroupedByLocation(); + if (pathRanges.size() != odb_.nlocs()) + { + oops::Log::error() << "ObsGnssroRefNLPEP2D::locations: SampledLocations has " + << pathRanges.size() << " entries, but must have one range for each of " + << odb_.nlocs() << " observations" << std::endl; + } + + std::vector seqNums(odb_.nlocs()); + odb_.get_db("MetaData", "sequenceNumber", seqNums); + + for (std::size_t obIdx = 0; obIdx < odb_.nlocs(); ++obIdx) + { + const util::Range & pathRange = pathRanges[obIdx]; + oops::Log::debug() << "SampledLoc for obIdx " << obIdx << " in profile seqNum=" + << seqNums[obIdx] << " has " << pathRange.end - pathRange.begin + << " locations from " << pathRange.begin << " to " << pathRange.end << "; " + << "lat/lon[" << pathRange.begin << "]=(" + << (sampledLocations.latitudes())[pathRange.begin] + << "," << (sampledLocations.longitudes())[pathRange.begin] << "), " + << "lat/lon[" << pathRange.end - 1 << "]=(" + << (sampledLocations.latitudes())[pathRange.end-1] + << "," << (sampledLocations.longitudes())[pathRange.end-1] << ")" << std::endl; + } + return sampledLocations; +} + +// ----------------------------------------------------------------------------- + +void ObsGnssroRefNLPEP2D::computeReducedVars(const oops::Variables & reducedVars, + GeoVaLs & /*geovals*/) const { + // No method for reducing the set of nhoriz_ profiles sampling each location into a single + // profile has been implemented so far, so when this obs operator is in use, neither it nor + // any obs filters or bias predictors can request variables in the reduced format. + if (reducedVars.size() != 0) + throw eckit::NotImplemented("ObsGnssroRefNLPEP2D is unable to compute the reduced " + "representation of GeoVaLs", Here()); +} + +// ----------------------------------------------------------------------------- + +void ObsGnssroRefNLPEP2D::print(std::ostream & os) const { + os << "ObsGnssroRefNLPEP2D::print not implemented"; +} + +// ----------------------------------------------------------------------------- + +} // namespace ufo diff --git a/src/ufo/operators/gnssro/RefNLPEP2D/ObsGnssroRefNLPEP2D.h b/src/ufo/operators/gnssro/RefNLPEP2D/ObsGnssroRefNLPEP2D.h new file mode 100644 index 000000000..b5c18034e --- /dev/null +++ b/src/ufo/operators/gnssro/RefNLPEP2D/ObsGnssroRefNLPEP2D.h @@ -0,0 +1,87 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#ifndef UFO_OPERATORS_GNSSRO_REFNLPEP2D_OBSGNSSROREFNLPEP2D_H_ +#define UFO_OPERATORS_GNSSRO_REFNLPEP2D_OBSGNSSROREFNLPEP2D_H_ + +#include +#include +#include + +#include "ioda/ObsDataVector.h" + +#include "oops/base/Variables.h" + +#include "ufo/ObsOperatorBase.h" +#include "ufo/operators/gnssro/RefNLPEP2D/ObsGnssroRefNLPEP2DParameters.h" +#include "ufo/operators/gnssro/utils/GnssroProfileRayPaths.h" +#include "ufo/operators/gnssro/utils/GnssroRayPathOrchestrator.h" +#include "ufo/operators/gnssro/utils/GnssroRayPathParameters.h" +#include "ufo/operators/gnssro/utils/RefractivityCalculator.h" + +/// Forward declarations +namespace ioda { + class ObsSpace; + class ObsVector; +} + +namespace ufo { + class GeoVaLs; + class ObsDiagnostics; + class GnssroRayPathGenerator; + +// ----------------------------------------------------------------------------- +/// GnssroRefNLPEP2D observation operator class +class ObsGnssroRefNLPEP2D : public ObsOperatorBase { + public: + /// The type of parameters accepted by the constructor of this operator. + /// This typedef is used by the ObsOperatorFactory. + typedef ObsGnssroRefNLPEP2DParameters Parameters_; + typedef ObsGnssroRefNLPEP2DOptions ObsOptions_; + typedef ioda::ObsDataVector QCFlags_t; + + static const std::string classname() {return "ufo::ObsGnssroRefNLPEP2D";} + + ObsGnssroRefNLPEP2D(const ioda::ObsSpace &, const Parameters_ &); + ~ObsGnssroRefNLPEP2D() override; + + // Obs Operator + void simulateObs(const GeoVaLs &, ioda::ObsVector &, ObsDiagnostics &, + const QCFlags_t &) const override; + + // Other + const oops::Variables & requiredVars() const override {return *varin_;} + + void computeReducedVars(const oops::Variables & reducedVars, GeoVaLs & geovals) const override; + Locations_ locations() const override; + + private: + void print(std::ostream &) const override; + + // Member data + const ioda::ObsSpace& odb_; + std::unique_ptr varin_; + ObsOptions_ opts_; + GnssroRayPathParameters rpParams_; + GnssroRayPathOrchestrator orchestrator_; + std::unique_ptr refrCalc_; +}; + +// ----------------------------------------------------------------------------- + +} // namespace ufo +#endif // UFO_OPERATORS_GNSSRO_REFNLPEP2D_OBSGNSSROREFNLPEP2D_H_ diff --git a/src/ufo/operators/gnssro/RefNLPEP2D/ObsGnssroRefNLPEP2DParameters.h b/src/ufo/operators/gnssro/RefNLPEP2D/ObsGnssroRefNLPEP2DParameters.h new file mode 100644 index 000000000..2da4c7a0f --- /dev/null +++ b/src/ufo/operators/gnssro/RefNLPEP2D/ObsGnssroRefNLPEP2DParameters.h @@ -0,0 +1,80 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#ifndef UFO_OPERATORS_GNSSRO_REFNLPEP2D_OBSGNSSROREFNLPEP2DPARAMETERS_H_ +#define UFO_OPERATORS_GNSSRO_REFNLPEP2D_OBSGNSSROREFNLPEP2DPARAMETERS_H_ + +#include + +#include "oops/util/parameters/OptionalParameter.h" +#include "oops/util/parameters/Parameter.h" +#include "oops/util/parameters/RequiredParameter.h" +#include "ufo/ObsOperatorParametersBase.h" +#include "ufo/operators/gnssro/utils/GnssroRayPathParameters.h" +#include "ufo/utils/parameters/ParameterTraitsVariable.h" + +namespace ufo { + +/// Configuration options recognized by the GnssroRefNLPEP2D operator. +class ObsGnssroRefNLPEP2DOptions : public oops::Parameters { + OOPS_CONCRETE_PARAMETERS(ObsGnssroRefNLPEP2DOptions, Parameters) + + public: + /// \brief GnssroRayPathGenerator to use to create ray path. + oops::Parameter rayPathGenType{"ray_path_gen_type", + GnssroRayPathParameters::DEFAULT_RAY_PATH_GEN_TYPE, this}; + + /// \brief Approximate ray length in km. + /// Ray length and res will be used to override the n_horiz value if ray_length + /// is set to a value > 0. The effective ray length may be adjusted slightly lower + /// or higher to meet the requirement of some ray path generators, such as having + /// an odd number of nodes. + oops::Parameter approxRayLength{"ray_length", + GnssroRayPathParameters::DEFAULT_APPROX_RAY_LENGTH_KM, this}; + + /// \brief horizontal resolution of nodes in ray path in km. + /// This should be close to the native resolution of the NWP model. + oops::Parameter res{"res", + GnssroRayPathParameters::DEFAULT_HORIZONTAL_RES_KM, this}; + + /// \brief The highest geometric height in km that will use 2D rays. + /// Levels above this height will be treated as 1D rays (one node at tangent point, + /// assuming the atmosphere is spherically symmetric around tangent point). + oops::Parameter top2D{"top_2d", + GnssroRayPathParameters::DEFAULT_TOP2D_KM, this}; + + /// \brief number of nodes in each 2D ray. Must be odd. + oops::Parameter nHoriz{"n_horiz", + GnssroRayPathParameters::DEFAULT_NHORIZ, this}; + + /// \brief The algorithm to use for computing refractivity and its derivatives. + oops::Parameter refrAlgorithm{"refr_algo", "RuegerBevis", this}; + + /// \brief true if the refractivity algorithm should assume a compressible atmosphere. + /// Set to false for an uncompressible atmosphere. + oops::Parameter useCompress{"use_compress", true, this}; +}; + +class ObsGnssroRefNLPEP2DParameters: public ObsOperatorParametersBase { + OOPS_CONCRETE_PARAMETERS(ObsGnssroRefNLPEP2DParameters, ObsOperatorParametersBase) + public: + oops::Parameter options{"obs options", {}, this}; +}; + + +} // namespace ufo +#endif // UFO_OPERATORS_GNSSRO_REFNLPEP2D_OBSGNSSROREFNLPEP2DPARAMETERS_H_ diff --git a/src/ufo/operators/gnssro/RefNLPEP2D/ObsGnssroRefNLPEP2DTLAD.cc b/src/ufo/operators/gnssro/RefNLPEP2D/ObsGnssroRefNLPEP2DTLAD.cc new file mode 100644 index 000000000..61686c411 --- /dev/null +++ b/src/ufo/operators/gnssro/RefNLPEP2D/ObsGnssroRefNLPEP2DTLAD.cc @@ -0,0 +1,538 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#include +#include +#include +#include +#include "eckit/exception/Exceptions.h" +#include "ufo/operators/gnssro/RefNLPEP2D/ObsGnssroRefNLPEP2DTLAD.h" +#include "ufo/operators/gnssro/utils/GnssroGeoVaLs.h" +#include "ufo/operators/gnssro/utils/GnssroProfileExtractor.h" +#include "ufo/operators/gnssro/utils/GnssroRayPathGenerator.h" +#include "ufo/utils/Constants.h" +#include "ufo/utils/VertInterp.interface.h" +#include "ufo/variabletransforms/Formulas.h" + + +#include "ioda/ObsSpace.h" +#include "ioda/ObsVector.h" + +#include "oops/util/Logger.h" + +#include "ufo/GeoVaLs.h" +#include "ufo/ObsDiagnostics.h" + +namespace ufo { + +// ----------------------------------------------------------------------------- +static LinearObsOperatorMaker makerGnssroRefNLPEP2D_TL_( + "GnssroRefNLPEP2D"); + +// ----------------------------------------------------------------------------- + +ObsGnssroRefNLPEP2DTLAD::ObsGnssroRefNLPEP2DTLAD(const ioda::ObsSpace & odb, + const Parameters_ & parameters) + : LinearObsOperatorBase(odb) + , varin_() + , opts_(parameters.options.value()) + , rpParams_(opts_.rayPathGenType, opts_.approxRayLength, + opts_.top2D, opts_.res, opts_.nHoriz) + , orchestrator_(obsspace(), rpParams_) + , refrCalc_(RefractivityCalculator::create(opts_.refrAlgorithm, opts_.useCompress)) + , trajectorySet_(false) +{ + const std::vector vv{GnssroGeoVaLs::VAR_TEMP, GnssroGeoVaLs::VAR_SPHUM, + GnssroGeoVaLs::VAR_PRES, GnssroGeoVaLs::VAR_GPH}; + varin_.reset(new oops::Variables(vv)); + + oops::Log::trace() << "ObsGnssroRefNLPEP2DTLAD created" << std::endl; +} + +// ----------------------------------------------------------------------------- + +ObsGnssroRefNLPEP2DTLAD::~ObsGnssroRefNLPEP2DTLAD() { + oops::Log::trace() << "ObsGnssroRefNLPEP2DTLAD destroyed" << std::endl; +} + +// ----------------------------------------------------------------------------- + +void ObsGnssroRefNLPEP2DTLAD::setTrajectory( + const GeoVaLs & gv, ObsDiagnostics &, const QCFlags_t &) { + const char funcName[] = "ObsGnssroRefNLPEP2DTLAD::setTrajectory"; + oops::Log::trace() << funcName << ": started" << std::endl; + + // Return early if there are no geoVals to operate on. + if (gv.nlocs() == 0) + { + oops::Log::trace() << funcName << ": trajectory not set, numGeoVaLs=" + << gv.nlocs() << std::endl; + return; + } + + // Extract the model data from the GeoVaLs. + GnssroGeoVaLs ggv(gv); + + // Get geometric height and location from the observations. + std::size_t nlocs = obsspace().nlocs(); + std::vector latitudes(nlocs); + std::vector geomHeights(nlocs); + obsspace().get_db("MetaData", "latitude", latitudes); + obsspace().get_db("MetaData", "height", geomHeights); + + std::size_t profileIdx = 0; + std::size_t ovecIdx = 0; // Index into locations in ObsSpace (nlocs) + double missing = util::missingValue(); + + // Iterate over RO profiles. + const GnssroRayPathOrchestrator::Profiles_& profiles = orchestrator_.profiles(); + for (const auto& roProfile : profiles) + { + int seqNum = roProfile->seqNum(); + oops::Log::debug() << funcName << ": processing profile seqNum=" << seqNum << std::endl; + + // Iterate over rays within this profile. + for (std::size_t rayIdx = 0; rayIdx < roProfile->numRays(); ++rayIdx, ++ovecIdx) + { + const GnssroRayPathOrchestrator::Ray_& ray = roProfile->ray(rayIdx); + + std::size_t tpNodeIdx = ray->tpNodeIdx(); + + // Initialize the trajectories for this ray. + ray->traj().initialize(ray->numNodes()); + + for (std::size_t nodeIdx = 0; nodeIdx < ray->numNodes(); ++nodeIdx) + { + // Include segment length and unit conversion in factor for refractivity Jacobian. + double segmentLength = ray->segLen(nodeIdx) * RefractivityCalculator::N_TO_EXCESS_IOR; + + // Get geometric height of this node + double geomHgt; + if (nodeIdx == tpNodeIdx) { + geomHgt = geomHeights[ovecIdx]; + } else { // Average of edges of segment. + double thisGeomHgt = ray->geomHgt(nodeIdx); + double nextGeomHgt = ray->geomHgt(nodeIdx + 1); + geomHgt = 0.5 * (thisGeomHgt + nextGeomHgt); + } + + // Get trajectory information from the model profile. + GnssroRayPathOrchestrator::TrajTuple trajTuple = + GnssroRayPathOrchestrator::modelTrajectory( + geomHgt, latitudes[ovecIdx], ggv, profileIdx, refrCalc_); + + // Save the weighted Jacobian values. + ray->traj().jacobianT(nodeIdx) = segmentLength * trajTuple.dNdT(); + ray->traj().jacobianP(nodeIdx) = segmentLength * trajTuple.dNdP(); + ray->traj().jacobianQ(nodeIdx) = segmentLength * trajTuple.dNdQ(); + + // Save the vertical interpolation weights in the trajectory. + ray->traj().wf(nodeIdx) = trajTuple.wf(); + ray->traj().wi(nodeIdx) = trajTuple.wi(); + oops::Log::debug() << funcName << ": seqNum=" << seqNum << ", rayIdx=" << rayIdx + << ", nodeIdx=" << nodeIdx << " (" << tpNodeIdx << "): 2D Traj" + << ", jacT=" << ray->traj().jacobianT(nodeIdx) + << ", jacQ=" << ray->traj().jacobianQ(nodeIdx) + << ", jacP=" << ray->traj().jacobianP(nodeIdx) + << ", wf=" << ray->traj().wf(nodeIdx) + << ", wi=" << ray->traj().wi(nodeIdx) << std::endl; + + // If using spherical symmetry, all nodes use the same model column. + if (!ray->assumeSphericalSymmetry()) { + ++profileIdx; + } + } // end of iteration of nodes in a single ray. + + // Mark trajectory for this ray as complete and go to next observed location. + ray->traj().finalize(); + } // end of iteration of rays in a single RO profile + + oops::Log::debug() << funcName << ": seqNum=" << seqNum << ", preparing for next profile" + << ", ovecIdx=" << ovecIdx << std::endl; + } // end of iteration over RO profiles + + trajectorySet_ = true; + oops::Log::trace() << funcName << ": trajectory set" << std::endl; +} + +// ----------------------------------------------------------------------------- + +void ObsGnssroRefNLPEP2DTLAD::simulateObsTL( + const GeoVaLs & gv, ioda::ObsVector & ovec, const QCFlags_t & qc_flags) const +{ + const char funcName[] = "ObsGnssroRefNLPEP2DTLAD::simulateObsTL"; + oops::Log::trace() << funcName << ": started" << std::endl; + + // Return early if there are no geoVals to operate on. + if (gv.nlocs() == 0) + { + oops::Log::trace() << funcName << ": not run because numGeoVaLs=" + << gv.nlocs() << std::endl; + return; + } + + // Validate our state. + if (!trajectorySet_) + { + oops::Log::error() << funcName << ": trajectory was not set, nlocs=" << gv.nlocs() + << std::endl; + throw eckit::SeriousBug("ObsGnssroRefNLPEP2DTLAD::simulateObsTL: trajectory was not set", + Here()); + } + + if (gv.nlocs() != ovec.size()) + { + oops::Log::error() << funcName << ": GeoVaLs:nlocs=" << gv.nlocs() << " not equal to " + << "size of ObsVector=" << ovec.size() << std::endl; + throw eckit::SeriousBug("ObsGnssroRefNLPEP2DTLAD::simulateObsTL: nlocs in geovals " + "not equal to size of ObsVector", Here()); + } + + // Extract the model data from the GeoVaLs. + GnssroGeoVaLs ggv(gv); + + const double missing = util::missingValue(); + std::size_t profileIdx = 0; + std::size_t ovecIdx = 0; // Index into locations in ObsVector (nlocs) + + // Iterate over RO profiles. + const GnssroRayPathOrchestrator::Profiles_& profiles = orchestrator_.profiles(); + for (const auto& roProfile : profiles) + { + int seqNum = roProfile->seqNum(); + oops::Log::debug() << funcName << ": processing profile seqNum=" << seqNum << std::endl; + + // Determine range of observation indices in this profile. + std::size_t startOvecIdx = ovecIdx; + std::size_t endOvecIdx = ovecIdx + roProfile->numRays(); + ovecIdx = endOvecIdx; + + // Get range of model profile indices in this RO profile + std::size_t startProfileIdx = profileIdx; + std::size_t endProfileIdx = startProfileIdx + roProfile->totalSampledLocations(); + profileIdx = endProfileIdx; + + oops::Log::debug() << funcName << ": seqNum=" << seqNum << ": processing profile " + << "with numRays=" << roProfile->numRays() << " ovecIdx range=[" + << startOvecIdx << ":" << endOvecIdx << "), profile-wide profileIdx range=[" + << startProfileIdx << "," << endProfileIdx << "), starting at ovecIdx=" + << ovecIdx - 1 << " and proceeding backwards" << std::endl; + + // Iterate over rays within this profile from top to bottom (reverse order) + for (GnssroProfileRayPaths::const_reverse_iterator rayItr = roProfile->crbegin(); + rayItr != roProfile->crend(); ++rayItr) + { + // Get the integer rayIdx for use in diagnostic messages. + // Get a forward iterator pointing to one position after the reverse iterator. + // Then convert the forward iterator into an integer position. + // This prevents us from even computing a negative rayIdx, which would happen + // if we used an integer index as a reverse iterator. + --ovecIdx; + auto fwdItr = rayItr.base() - 1; + std::size_t rayIdx = std::distance(roProfile->cbegin(), fwdItr); + const GnssroProfileRayPaths::Ray_& ray = *rayItr; + + // Get the model profile range for this ray + std::size_t rayEndProfileIdx = profileIdx; + std::size_t rayStartProfileIdx = rayEndProfileIdx - ray->numLocations(); + profileIdx = rayStartProfileIdx; + oops::Log::debug() << funcName << ": seqNum=" << seqNum << ", rayIdx=" + << rayIdx << ", numNodes=" << ray->numNodes() << ", numSampLocs=" + << ray->numLocations() << ", ray profileIdx range=[" << rayStartProfileIdx + << "," << rayEndProfileIdx << "), starting at profileIdx=" + << profileIdx << std::endl; + + std::size_t tpNodeIdx = ray->tpNodeIdx(); + + double fwdTL = 0.0; + ovec[ovecIdx] = 0.0; // Accumulate contribution to TL hofx from each node. + for (std::size_t nodeIdx = 0; nodeIdx < ray->numNodes(); ++nodeIdx) + { + // Get references to model vertical profiles (columns) for this profileIdx. + const GnssroGeoVaLs::Profile_& tempProf = ggv.tempProfile(profileIdx); + const GnssroGeoVaLs::Profile_& sphumProf = ggv.sphumProfile(profileIdx); + const GnssroGeoVaLs::Profile_& presProf = ggv.presProfile(profileIdx); + + // Vertically interpolate using previously computed weights for this node. + int wi = ray->traj().wi(nodeIdx); + double wf = ray->traj().wf(nodeIdx); + double temp; // Interpolated temperature perturbation in K + double sphum; // Interpolated specific humidity perturbation in kg/kg + double pres; // Interpolated pressure perturbation in Pa + vert_interp_apply_f90(ggv.nlevs(), tempProf.data(), temp, wi, wf); + vert_interp_apply_f90(ggv.nlevs(), sphumProf.data(), sphum, wi, wf); + vert_interp_apply_f90(ggv.nlevs(), presProf.data(), pres, wi, wf); + + // Compute the forward part + // Note: Jacobian at nodes is already node-weighted. + double fwdIncr = ( + ray->traj().jacobianT(nodeIdx) * temp + + ray->traj().jacobianP(nodeIdx) * pres + + ray->traj().jacobianQ(nodeIdx) * sphum); + fwdTL += fwdIncr; + + // If using spherical symmetry, all nodes use the same model column. + if (!ray->assumeSphericalSymmetry()) { + ++profileIdx; + } + } // end of loop over nodes + + ovec[ovecIdx] = fwdTL; + + // Prepare for the next lower ray by setting the profileIdx to point to + // the start of the model profiles for the current ray, which is also the + // end of the range for the next lower ray within the profile. + profileIdx = rayStartProfileIdx; + } // end of iteration of rays in a single RO profile + + // Prepare for iteration into next RO profile. + oops::Log::debug() << funcName << ": seqNum=" << seqNum << " DONE: " + << "ovecIdx=" << ovecIdx << std::endl; + ovecIdx = endOvecIdx; + if (profileIdx != startProfileIdx) { + oops::Log::error() << funcName << ": seqNum=" << seqNum + << ": model profile indexing error: after RO profile processed," + << " profileIdx=" << profileIdx << "; should be same as " + << "startProfileIdx=" << startProfileIdx << std::endl; + } + profileIdx = endProfileIdx; + oops::Log::debug() << funcName << ": seqNum=" << seqNum << " DONE: " + << "skipping to next profile ovecIdx=" << ovecIdx + << ", profileIdx=" << profileIdx << std::endl; + } // end of iteration over RO profiles + + oops::Log::trace() << funcName << ": complete" << std::endl; +} + +//---------------------------------------------------------------------------- + +void ObsGnssroRefNLPEP2DTLAD::simulateObsAD(GeoVaLs & gv, const ioda::ObsVector & ovec, + const QCFlags_t & qc_flags) const +{ + const char funcName[] = "ObsGnssroRefNLPEP2DTLAD::simulateObsAD"; + oops::Log::trace() << funcName << ": started" << std::endl; + + // Return early if there are no geoVals to operate on. + if (gv.nlocs() == 0) + { + oops::Log::trace() << funcName << ": not run because numGeoVaLs=" + << gv.nlocs() << std::endl; + return; + } + + // Validate our state. + if (!trajectorySet_) + { + oops::Log::error() << funcName << ": trajectory was not set, nlocs=" << gv.nlocs() + << std::endl; + throw eckit::SeriousBug("ObsGnssroRefNLPEP2DTLAD::simulateObsTL: trajectory was not set", + Here()); + } + + if (gv.nlocs() != ovec.size()) + { + oops::Log::error() << funcName << ": GeoVaLs:nlocs=" << gv.nlocs() << " not equal to " + << "size of ObsVector=" << ovec.size() << std::endl; + throw eckit::SeriousBug("ObsGnssroRefNLPEP2DTLAD::simulateObsTL: nlocs in geovals " + "not equal to size of ObsVector", Here()); + } + + // Extract the model data from the GeoVaLs. + GnssroGeoVaLs ggv(gv); + + const double missing = util::missingValue(); + std::size_t profileIdx = 0; + std::size_t ovecIdx = 0; // Index into locations in ObsVector (nlocs) + // Iterate over RO profiles. + const GnssroRayPathOrchestrator::Profiles_& profiles = orchestrator_.profiles(); + for (const auto& roProfile : profiles) + { + int seqNum = roProfile->seqNum(); + oops::Log::debug() << funcName << ": processing profile seqNum=" << seqNum << std::endl; + + // Determine range of observation indices in this profile. + std::size_t startOvecIdx = ovecIdx; + std::size_t endOvecIdx = ovecIdx + roProfile->numRays(); + ovecIdx = endOvecIdx; + + // Get range of model profile indices in this RO profile + std::size_t startProfileIdx = profileIdx; + std::size_t endProfileIdx = startProfileIdx + roProfile->totalSampledLocations(); + profileIdx = endProfileIdx; + + oops::Log::debug() << funcName << ": seqNum=" << seqNum << ": processing profile " + << "with numRays=" << roProfile->numRays() << " ovecIdx range=[" + << startOvecIdx << ":" << endOvecIdx << "), profile-wide profileIdx range=[" + << startProfileIdx << "," << endProfileIdx << "), starting at ovecIdx=" + << ovecIdx - 1 << " and proceeding backwards" << std::endl; + + // Iterate over rays within this profile from top to bottom (reverse order) + for (GnssroProfileRayPaths::const_reverse_iterator rayItr = roProfile->crbegin(); + rayItr != roProfile->crend(); ++rayItr) + { + // Get the integer rayIdx for use in diagnostic messages. + // Get a forward iterator pointing to one position after the reverse iterator. + // Then convert the forward iterator into an integer position. + // This prevents us from even computing a negative rayIdx, which would happen + // if we used an integer index as a reverse iterator. + --ovecIdx; + auto fwdItr = rayItr.base() - 1; // fwdItr points to same loc as reverse itr + std::size_t rayIdx = std::distance(roProfile->cbegin(), fwdItr); + const GnssroProfileRayPaths::Ray_& ray = *rayItr; + + // Get the model profile range for this ray + std::size_t rayEndProfileIdx = profileIdx; + std::size_t rayStartProfileIdx = rayEndProfileIdx - ray->numLocations(); + profileIdx = rayStartProfileIdx; + oops::Log::debug() << funcName << ": seqNum=" << seqNum << ", rayIdx=" + << rayIdx << ", ovecIdx=" << ovecIdx << ", ovec.size=" << ovec.size() + << ", numNodes=" << ray->numNodes() << ", numSampLocs=" + << ray->numLocations() << ", ray profileIdx range=[" << rayStartProfileIdx + << "," << rayEndProfileIdx << "), starting at profileIdx=" + << profileIdx << std::endl; + + std::size_t tpNodeIdx = ray->tpNodeIdx(); + + for (std::size_t nodeIdx = 0; nodeIdx < ray->numNodes(); ++nodeIdx) + { + // Compute the perturbations in each variable for the forward part + double dTemp = ovec[ovecIdx] * ray->traj().jacobianT(nodeIdx); + double dSphum = ovec[ovecIdx] * ray->traj().jacobianQ(nodeIdx); + double dPres = ovec[ovecIdx] * ray->traj().jacobianP(nodeIdx); + + // These perturbations should be divided over one or two vertical + // levels in the profile (wi and wi+1), with relative amounts + // determined by the wf value. This logic follows the Fortran + // function vert_interp_apply_ad. + // NOTE: when weighting the contributions, we want to include both + // the horizontal weight (segmentLength) and the vertical weight (wf). + // However, segmentLength is already part of the value returned + // by the node-specific Jacobians. + int wi = ray->traj().wi(nodeIdx); + int cwi = wi - 1; // 0-based index equivalent to wi. + double wf = ray->traj().wf(nodeIdx); + if (cwi >= 0) // Do not go below lowest vertical level in model + { + std::size_t levelIdx = cwi; + double weight = wf; + oops::Log::debug() << funcName << ": seqNum=" << seqNum + << ", rayIdx=" << rayIdx << ", nodeIdx=" << nodeIdx + << ", prof=" << profileIdx << ", lev=" << levelIdx + << ", incr temp=" << ggv.temp(profileIdx, levelIdx) + << " by " << (weight * dTemp) << ", incr sphum=" + << ggv.sphum(profileIdx, levelIdx) << " by " + << (weight * dSphum) << ", incr pres=" + << ggv.pres(profileIdx, levelIdx) << " by " + << (weight * dPres) << std::endl; + double & geoval_temp = ggv.temp(profileIdx, levelIdx); + if (geoval_temp != missing) { + geoval_temp += (weight * dTemp); + } + double & geoval_sphum = ggv.sphum(profileIdx, levelIdx); + if (geoval_sphum != missing) { + geoval_sphum += (weight * dSphum); + } + double & geoval_pres = ggv.pres(profileIdx, levelIdx); + if (geoval_pres != missing) { + geoval_pres += (weight * dPres); + } + oops::Log::debug() << funcName << ": seqNum=" << seqNum + << ", rayIdx=" << rayIdx << ", nodeIdx=" << nodeIdx + << ", prof=" << profileIdx << ", lev=" << levelIdx + << ", set temp=" << ggv.temp(profileIdx, levelIdx) + << ", set sphum=" << ggv.sphum(profileIdx, levelIdx) + << ", set pres=" << ggv.pres(profileIdx, levelIdx) + << std::endl; + } + + if (cwi + 1 < ggv.nlevs()) // Do not go above highest level in model + { + std::size_t levelIdx = cwi + 1; + double weight = (1.0 - wf); + oops::Log::debug() << funcName << ": seqNum=" << seqNum + << ", rayIdx=" << rayIdx << ", nodeIdx=" << nodeIdx + << ", prof=" << profileIdx << ", lev=" << levelIdx + << ", incr temp=" << ggv.temp(profileIdx, levelIdx) + << " by " << (weight * dTemp) << ", incr sphum=" + << ggv.sphum(profileIdx, levelIdx) << " by " + << (weight * dSphum) << ", incr pres=" + << ggv.pres(profileIdx, levelIdx) << " by " + << (weight * dPres) << std::endl; + double & geoval_temp = ggv.temp(profileIdx, levelIdx); + if (geoval_temp != missing) { + geoval_temp += (weight * dTemp); + } + double & geoval_sphum = ggv.sphum(profileIdx, levelIdx); + if (geoval_sphum != missing) { + geoval_sphum += (weight * dSphum); + } + double & geoval_pres = ggv.pres(profileIdx, levelIdx); + if (geoval_pres != missing) { + geoval_pres += (weight * dPres); + } + oops::Log::debug() << funcName << ": seqNum=" << seqNum + << ", rayIdx=" << rayIdx << ", nodeIdx=" << nodeIdx + << ", prof=" << profileIdx << ", lev=" << levelIdx + << ", set temp=" << ggv.temp(profileIdx, levelIdx) + << ", set sphum=" << ggv.sphum(profileIdx, levelIdx) + << ", set pres=" << ggv.pres(profileIdx, levelIdx) + << std::endl; + } + + // If using spherical symmetry, all nodes use the same model column. + if (!ray->assumeSphericalSymmetry()) { + ++profileIdx; + } + } // end of iteration over nodes for adjoint. + + // Prepare for the next lower ray by setting the profileIdx to point to + // the start of the model profiles for the current ray, which is also the + // end of the range for the next lower ray within the profile. + profileIdx = rayStartProfileIdx; + } // end of iteration of rays in a single RO profile + + // Prepare for iteration into next RO profile. + oops::Log::debug() << funcName << ": seqNum=" << seqNum << " DONE: " + << "ovecIdx=" << ovecIdx << std::endl; + ovecIdx = endOvecIdx; + if (profileIdx != startProfileIdx) { + oops::Log::error() << funcName << ": seqNum=" << seqNum + << ": model profile indexing error: after RO profile processed," + << " profileIdx=" << profileIdx << "; should be same as " + << "startProfileIdx=" << startProfileIdx << std::endl; + } + profileIdx = endProfileIdx; + oops::Log::debug() << funcName << ": seqNum=" << seqNum << " DONE: " + << "skipping to next profile ovecIdx=" << ovecIdx + << ", profileIdx=" << profileIdx << std::endl; + } // end of iteration over RO profiles + + // Save the data from the GnssroGeoVaLs to the GeoVaLs passed to us. + ggv.saveTo(gv); + + oops::Log::trace() << funcName << ": complete" << std::endl; +} + +// ----------------------------------------------------------------------------- + +void ObsGnssroRefNLPEP2DTLAD::print(std::ostream & os) const { + os << "ObsGnssroRefNLPEP2DTLAD::print not implemented" << std::endl; +} + +// ----------------------------------------------------------------------------- + +} // namespace ufo diff --git a/src/ufo/operators/gnssro/RefNLPEP2D/ObsGnssroRefNLPEP2DTLAD.h b/src/ufo/operators/gnssro/RefNLPEP2D/ObsGnssroRefNLPEP2DTLAD.h new file mode 100644 index 000000000..dd696332a --- /dev/null +++ b/src/ufo/operators/gnssro/RefNLPEP2D/ObsGnssroRefNLPEP2DTLAD.h @@ -0,0 +1,84 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#ifndef UFO_OPERATORS_GNSSRO_REFNLPEP2D_OBSGNSSROREFNLPEP2DTLAD_H_ +#define UFO_OPERATORS_GNSSRO_REFNLPEP2D_OBSGNSSROREFNLPEP2DTLAD_H_ + +#include +#include +#include + +#include "ioda/ObsDataVector.h" + +#include "oops/base/Variables.h" + +#include "ufo/LinearObsOperatorBase.h" +#include "ufo/operators/gnssro/RefNLPEP2D/ObsGnssroRefNLPEP2DParameters.h" +#include "ufo/operators/gnssro/utils/GnssroProfileRayPaths.h" +#include "ufo/operators/gnssro/utils/GnssroRayPathOrchestrator.h" +#include "ufo/operators/gnssro/utils/GnssroRayPathParameters.h" +#include "ufo/operators/gnssro/utils/RefractivityCalculator.h" + +// Forward declarations +namespace ioda { + class ObsSpace; + class ObsVector; +} + +namespace ufo { + class GeoVaLs; + +// ----------------------------------------------------------------------------- +/// GnssroRefNLPEP2D TL/AD observation operator class +class ObsGnssroRefNLPEP2DTLAD : public LinearObsOperatorBase { + public: + /// The type of parameters accepted by the constructor of this operator. + /// This typedef is used by the LinearObsOperatorFactory. + typedef ObsGnssroRefNLPEP2DParameters Parameters_; + typedef ObsGnssroRefNLPEP2DOptions ObsOptions_; + typedef ioda::ObsDataVector QCFlags_t; + + static const std::string classname() {return "ufo::ObsGnssroRefNLPEP2DTLAD";} + + + ObsGnssroRefNLPEP2DTLAD(const ioda::ObsSpace &, const Parameters_ &); + ~ObsGnssroRefNLPEP2DTLAD() override; + + // Obs Operators + void setTrajectory(const GeoVaLs &, ObsDiagnostics &, const QCFlags_t &) override; + void simulateObsTL(const GeoVaLs &, ioda::ObsVector &, const QCFlags_t &) const override; + void simulateObsAD(GeoVaLs &, const ioda::ObsVector &, const QCFlags_t &) const override; + + // Other + const oops::Variables & requiredVars() const override {return *varin_;} + + private: + void print(std::ostream &) const override; + + // Member data + std::unique_ptr varin_; + ObsOptions_ opts_; + GnssroRayPathParameters rpParams_; + GnssroRayPathOrchestrator orchestrator_; + std::unique_ptr refrCalc_; + bool trajectorySet_; +}; + +// ----------------------------------------------------------------------------- + +} // namespace ufo +#endif // UFO_OPERATORS_GNSSRO_REFNLPEP2D_OBSGNSSROREFNLPEP2DTLAD_H_ diff --git a/src/ufo/operators/gnssro/utils/CMakeLists.txt b/src/ufo/operators/gnssro/utils/CMakeLists.txt index 30142bad5..8ae301b3c 100644 --- a/src/ufo/operators/gnssro/utils/CMakeLists.txt +++ b/src/ufo/operators/gnssro/utils/CMakeLists.txt @@ -10,6 +10,16 @@ set ( utils_src_files gnssro_mod_constants.F90 gnssro_mod_transform.F90 gnssro_mod_obserror.F90 + GnssroGeoVaLs.cc + GnssroProfileExtractor.cc + GnssroProfileRayPaths.cc + GnssroRayPath.cc + GnssroRayPathGenerator.cc + GnssroRayPathOrchestrator.cc + GnssroRayPathParameters.cc + GnssroRayTrajectory.cc + RefractivityCalculator.cc + StraightLineRayPathGenerator.cc PARENT_SCOPE ) diff --git a/src/ufo/operators/gnssro/utils/GnssroGeoVaLs.cc b/src/ufo/operators/gnssro/utils/GnssroGeoVaLs.cc new file mode 100644 index 000000000..fb17d5646 --- /dev/null +++ b/src/ufo/operators/gnssro/utils/GnssroGeoVaLs.cc @@ -0,0 +1,124 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#include +#include "oops/base/Variable.h" +#include "oops/util/Logger.h" +#include "ufo/operators/gnssro/utils/GnssroGeoVaLs.h" + +namespace ufo { + +// oops variables used in GeoVaLs interface. +static const oops::Variable sTempVar(GnssroGeoVaLs::VAR_TEMP); +static const oops::Variable sSphumVar(GnssroGeoVaLs::VAR_SPHUM); +static const oops::Variable sPresVar(GnssroGeoVaLs::VAR_PRES); +static const oops::Variable sGphVar(GnssroGeoVaLs::VAR_GPH); + +// ----------------------------------------------------------------------------- +GnssroGeoVaLs::GnssroGeoVaLs(const GeoVaLs & cgv) + : nprofiles_(validateProfiles(cgv)) + , nlevs_(validateLevels(cgv)) + , temp_() + , sphum_() + , pres_() + , gph_() +{ + temp_.reserve(nprofiles_); + sphum_.reserve(nprofiles_); + pres_.reserve(nprofiles_); + gph_.reserve(nprofiles_); + + for (std::size_t profIdx = 0; profIdx < nprofiles_; ++profIdx) + { + getProfile(sTempVar, profIdx, temp_, cgv); + getProfile(sSphumVar, profIdx, sphum_, cgv); + getProfile(sPresVar, profIdx, pres_, cgv); + getProfile(sGphVar, profIdx, gph_, cgv); + } +} +// ----------------------------------------------------------------------------- + +GnssroGeoVaLs::~GnssroGeoVaLs() +{ +} +// ----------------------------------------------------------------------------- + +void GnssroGeoVaLs::saveTo(GeoVaLs & gv) +{ + putProfiles(sTempVar, temp_, gv); + putProfiles(sSphumVar, sphum_, gv); + putProfiles(sPresVar, pres_, gv); + putProfiles(sGphVar, gph_, gv); +} +// ----------------------------------------------------------------------------- + +std::size_t GnssroGeoVaLs::validateProfiles(const GeoVaLs& cgv) +{ + std::size_t nprofTemp = cgv.nprofiles(sTempVar); + std::size_t nprofSphum = cgv.nprofiles(sSphumVar); + std::size_t nprofPres = cgv.nprofiles(sPresVar); + std::size_t nprofGph = cgv.nprofiles(sGphVar); + if (nprofTemp != nprofSphum || nprofTemp != nprofPres || nprofTemp != nprofGph) + { + throw eckit::BadValue("GnssroGeoVaLs: number of profiles for " + sTempVar.name() + + ", " + sSphumVar.name() + ", " + sPresVar.name() + ", and " + sGphVar.name() + + " must be the same", Here()); + } + return nprofTemp; +} +// ----------------------------------------------------------------------------- + +std::size_t GnssroGeoVaLs::validateLevels(const GeoVaLs& cgv) +{ + std::size_t nlevsTemp = cgv.nlevs(sTempVar); + std::size_t nlevsSphum = cgv.nlevs(sSphumVar); + std::size_t nlevsPres = cgv.nlevs(sPresVar); + std::size_t nlevsGph = cgv.nlevs(sGphVar); + if (nlevsTemp != nlevsSphum || nlevsTemp != nlevsPres || nlevsTemp != nlevsGph) + { + throw eckit::BadValue("GnssroGeoVaLs: number of vertical levels for " + sTempVar.name() + + ", " + sSphumVar.name() + ", " + sPresVar.name() + ", and " + sGphVar.name() + + " must be the same", Here()); + } + return nlevsTemp; +} +// ----------------------------------------------------------------------------- + +void GnssroGeoVaLs::getProfile( + const oops::Variable& var, + std::size_t profIdx, + Profiles_ & profiles, + const GeoVaLs & cgv) +{ + ProfilePtr_ profPtr = std::make_unique(nlevs_, 0.0); + cgv.getProfile(*profPtr, var, profIdx); + profiles.push_back(std::move(profPtr)); +} +// ----------------------------------------------------------------------------- + +void GnssroGeoVaLs::putProfiles(const oops::Variable& var, const Profiles_& profiles, GeoVaLs& gv) +{ + for (std::size_t profIdx = 0; profIdx < nprofiles_; ++profIdx) + { + const ProfilePtr_& profPtr = profiles.at(profIdx); + gv.putProfile(*profPtr, var, profIdx); + } +} +// ----------------------------------------------------------------------------- + +} // namespace ufo + diff --git a/src/ufo/operators/gnssro/utils/GnssroGeoVaLs.h b/src/ufo/operators/gnssro/utils/GnssroGeoVaLs.h new file mode 100644 index 000000000..97c26da80 --- /dev/null +++ b/src/ufo/operators/gnssro/utils/GnssroGeoVaLs.h @@ -0,0 +1,133 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#ifndef UFO_OPERATORS_GNSSRO_UTILS_GNSSROGEOVALS_H_ +#define UFO_OPERATORS_GNSSRO_UTILS_GNSSROGEOVALS_H_ + +#include +#include +#include +#include "ufo/GeoVaLs.h" + +namespace ufo { + +// ----------------------------------------------------------------------------- +/// GnssroGeoVaLs - +/// C++ Interface to GeoVaLs variables needed by GNSSRO forward operators. +/// NOTE: this does not account for FO's that need some variables on +/// different vertical levels, e.g. air_pressure_levels variable. +/// +class GnssroGeoVaLs { + public: + // Type definitions + typedef std::vector Profile_; + typedef std::unique_ptr ProfilePtr_; + typedef std::vector Profiles_; + + // String constants + static constexpr const char* CLASSNAME = "GnssroGeoVaLs"; + static constexpr const char* VAR_TEMP = "air_temperature"; + static constexpr const char* VAR_SPHUM = "water_vapor_mixing_ratio_wrt_moist_air"; + static constexpr const char* VAR_PRES = "air_pressure"; + static constexpr const char* VAR_GPH = "geopotential_height"; + static constexpr const char* VAR_SFC_ALT = "height_above_mean_sea_level_at_surface"; + static constexpr std::array PROFILES_VARS = { + VAR_TEMP, VAR_SPHUM, VAR_PRES, VAR_GPH}; + + static std::string name() {return CLASSNAME;} + + // Constructor/Destructor + explicit GnssroGeoVaLs(const GeoVaLs & gv); + ~GnssroGeoVaLs(); + // Disallow default, copy and move construction and assignment. + GnssroGeoVaLs() = delete; + GnssroGeoVaLs(const GnssroGeoVaLs & other) = delete; + GnssroGeoVaLs(GnssroGeoVaLs && other) noexcept = delete; + GnssroGeoVaLs & operator=(const GnssroGeoVaLs & other) = delete; + GnssroGeoVaLs & operator=(GnssroGeoVaLs && other) noexcept = delete; + + // Accessors for dimensions + std::size_t nprofiles() const {return nprofiles_;} + std::size_t nlevs() const {return nlevs_;} + + // Profile Read Accessors + const Profile_ & tempProfile(std::size_t profileIdx) const {return *(temp_[profileIdx]);} + const Profile_ & sphumProfile(std::size_t profileIdx) const {return *(sphum_[profileIdx]);} + const Profile_ & presProfile(std::size_t profileIdx) const {return *(pres_[profileIdx]);} + const Profile_ & gphProfile(std::size_t profileIdx) const {return *(gph_[profileIdx]);} + + // Profile Read-write Accessors + Profile_ & tempProfile(std::size_t profileIdx) {return *(temp_[profileIdx]);} + Profile_ & sphumProfile(std::size_t profileIdx) {return *(sphum_[profileIdx]);} + Profile_ & presProfile(std::size_t profileIdx) {return *(pres_[profileIdx]);} + Profile_ & gphProfile(std::size_t profileIdx) {return *(gph_[profileIdx]);} + + // Scalar element read accessors + double temp(std::size_t profileIdx, std::size_t levIdx) const { + return (*(temp_[profileIdx]))[levIdx]; + } + double sphum(std::size_t profileIdx, std::size_t levIdx) const { + return (*(sphum_[profileIdx]))[levIdx]; + } + double pres(std::size_t profileIdx, std::size_t levIdx) const { + return (*(pres_[profileIdx]))[levIdx]; + } + double gph(std::size_t profileIdx, std::size_t levIdx) const { + return (*(gph_[profileIdx]))[levIdx]; + } + + // Scalar element read-write accessors + double & temp(std::size_t profileIdx, std::size_t levIdx) { + return (*(temp_[profileIdx]))[levIdx]; + } + double & sphum(std::size_t profileIdx, std::size_t levIdx) { + return (*(sphum_[profileIdx]))[levIdx]; + } + double & pres(std::size_t profileIdx, std::size_t levIdx) { + return (*(pres_[profileIdx]))[levIdx]; + } + double & gph(std::size_t profileIdx, std::size_t levIdx) { + return (*(gph_[profileIdx]))[levIdx]; + } + + // Save values to GeoVaLs + void saveTo(GeoVaLs & gv); + + private: + // Implementation methods + std::size_t validateProfiles(const GeoVaLs & cgv); + std::size_t validateLevels(const GeoVaLs & cgv); + void getProfile(const oops::Variable& var, + std::size_t profIdx, + Profiles_ & profiles, + const GeoVaLs& cgv); + void putProfiles(const oops::Variable& var, + const Profiles_& profiles, + GeoVaLs& gv); + + // Member data + std::size_t nprofiles_; // Number of profiles (sampled locations) + std::size_t nlevs_; // Number of levels in each profile + Profiles_ temp_; // Air temperature profiles. + Profiles_ sphum_; // Specific humidity profiles + Profiles_ pres_; // Air pressure profiles + Profiles_ gph_; // Geopotential height profiles. +}; + +} // namespace ufo + +#endif // UFO_OPERATORS_GNSSRO_UTILS_GNSSROGEOVALS_H_ diff --git a/src/ufo/operators/gnssro/utils/GnssroProfileExtractor.cc b/src/ufo/operators/gnssro/utils/GnssroProfileExtractor.cc new file mode 100644 index 000000000..f1af07300 --- /dev/null +++ b/src/ufo/operators/gnssro/utils/GnssroProfileExtractor.cc @@ -0,0 +1,78 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#include "oops/util/Logger.h" +#include "ufo/operators/gnssro/utils/GnssroProfileExtractor.h" + +namespace ufo { + +// ----------------------------------------------------------------------------- + +GnssroProfileExtractor::GnssroProfileExtractor(const ioda::ObsSpace & odb) + : slices_() +{ + oops::Log::trace() << "GnssroProfileExtractor() created" << std::endl; + std::size_t sumSliceSizes = 0; + if (odb.nlocs() > 0) + { + // Get sequenceNumber from the observations. + std::vector seqNums(odb.nlocs()); + std::vector times(odb.nlocs()); + odb.get_db("MetaData", "sequenceNumber", seqNums); + odb.get_db("MetaData", "dateTime", times); + + // If the data is empty, there is nothing to do. + if (seqNums.empty()) { + oops::Log::trace() << "GnssroProfileExtractor(): constructor complete: " + << "found no sequenceNumber values" << std::endl; + return; + } + oops::Log::debug() << "GnssroProfileExtractor(): found " << seqNums.size() + << " sequenceNumber values" << std::endl; + + // Find the range associated with each unique sequence number. + // NOTE: due to filtering and distribution of work between ranks, + // there may be gaps in the integer values and the sequence numbers + // will not be in order. + size_t start = 0; + int seqNum = seqNums[0]; + util::DateTime profileTime = times[0]; + for (std::size_t iloc = 1; iloc < seqNums.size(); ++iloc) { + if (seqNums[iloc] != seqNum) { + slices_.push_back(GnssroProfileSlice(seqNum, profileTime, start, iloc)); + start = iloc; + seqNum = seqNums[iloc]; + profileTime = times[iloc]; + } + } + slices_.push_back(GnssroProfileSlice(seqNum, profileTime, start, seqNums.size())); + sumSliceSizes += seqNums.size(); + } + + oops::Log::trace() << "GnssroProfileExtractor() constructor complete: " << slices_.size() + << " profiles found with total of " << sumSliceSizes << " obs; nlocs=" + << odb.nlocs() << " expected" << std::endl; +} +// ----------------------------------------------------------------------------- + +GnssroProfileExtractor::~GnssroProfileExtractor() +{ +} +// ----------------------------------------------------------------------------- + +} // namespace ufo + diff --git a/src/ufo/operators/gnssro/utils/GnssroProfileExtractor.h b/src/ufo/operators/gnssro/utils/GnssroProfileExtractor.h new file mode 100644 index 000000000..3d0f17845 --- /dev/null +++ b/src/ufo/operators/gnssro/utils/GnssroProfileExtractor.h @@ -0,0 +1,56 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#ifndef UFO_OPERATORS_GNSSRO_UTILS_GNSSROPROFILEEXTRACTOR_H_ +#define UFO_OPERATORS_GNSSRO_UTILS_GNSSROPROFILEEXTRACTOR_H_ + +#include // For std::size_t +#include +#include "ioda/ObsSpace.h" +#include "ufo/operators/gnssro/utils/GnssroProfileSlice.h" + +namespace ufo { + +// ----------------------------------------------------------------------------- +/// Class to extract the GnssroProfileSlices from a set of GNSSRO observation samples. +/// Each slice references a distinct RO profile. +class GnssroProfileExtractor { + public: + typedef std::vector Slices; + typedef Slices::const_iterator const_iterator; + + explicit GnssroProfileExtractor(const ioda::ObsSpace & odb); + ~GnssroProfileExtractor(); + + const Slices & slices() const {return slices_;} + const_iterator cbegin() const {return slices_.cbegin();} + const_iterator cend() const {return slices_.cend();} + bool is_empty() const {return slices_.empty();} + + private: + GnssroProfileExtractor() = delete; + GnssroProfileExtractor(const GnssroProfileExtractor &) = delete; + GnssroProfileExtractor(GnssroProfileExtractor &&) noexcept = delete; + GnssroProfileExtractor& operator=(const GnssroProfileExtractor &) = delete; + GnssroProfileExtractor& operator=(GnssroProfileExtractor &&) noexcept = delete; + + Slices slices_; +}; + +} // namespace ufo + +#endif // UFO_OPERATORS_GNSSRO_UTILS_GNSSROPROFILEEXTRACTOR_H_ diff --git a/src/ufo/operators/gnssro/utils/GnssroProfileRayPaths.cc b/src/ufo/operators/gnssro/utils/GnssroProfileRayPaths.cc new file mode 100644 index 000000000..8e547b746 --- /dev/null +++ b/src/ufo/operators/gnssro/utils/GnssroProfileRayPaths.cc @@ -0,0 +1,164 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#include // for std::copy, std::nth_element +#include +#include // for NaN +#include +#include "oops/util/Logger.h" +#include "ufo/operators/gnssro/utils/GnssroProfileRayPaths.h" +#include "ufo/operators/gnssro/utils/GnssroRayPath.h" + +namespace ufo { + +// Helper function to approximate a median calculation. +// This sorts the input vector just enough to determine the central value. +// However, for vectors with an even number of elements, it does +// not take the mean of the two values closest to the center. +static float approx_median( + const std::vector & v, + std::size_t start_offset, + std::size_t end_offset, + std::vector & tmp_v) +{ + // Copy the subset of the source vector to the temporary vector. + std::copy(v.begin() + start_offset, v.begin() + end_offset, tmp_v.begin()); + + // Determine approximate location of the median in sorted order. + // If the size is even, the median should be the average of the + // two center values. However, we don't need to be that precise. + std::size_t median_loc = tmp_v.size() / 2; + + // Use nth_element to sort the data just enough to place + // the value for the median location in the right place. + // This is more efficient than sorting the entire vector. + std::nth_element(tmp_v.begin(), tmp_v.begin() + median_loc, tmp_v.end()); + return tmp_v[median_loc]; +} + +// ----------------------------------------------------------------------------- +/// Class to represent represent the path of a single GNSSRO ray. +/// Typically, each sample in a set of RO observations becomes a ray. +GnssroProfileRayPaths::GnssroProfileRayPaths() + : seqNum_(-1) + , profileLat_(std::numeric_limits::quiet_NaN()) + , profileLon_(std::numeric_limits::quiet_NaN()) + , profileAzimuth_(std::numeric_limits::quiet_NaN()) + , rays_() +{ + oops::Log::trace() << "GnssroProfileRayPaths default c'tor" << std::endl; +} +// ----------------------------------------------------------------------------- + +GnssroProfileRayPaths::GnssroProfileRayPaths( + const GnssroProfileSlice & slice, + const std::vector & lats, + const std::vector & lons, + const std::vector & azimuths) + : seqNum_(slice.seqNum()) + , profileLat_(std::numeric_limits::quiet_NaN()) + , profileLon_(std::numeric_limits::quiet_NaN()) + , profileAzimuth_(std::numeric_limits::quiet_NaN()) + , rays_() +{ + oops::Log::trace() << "GnssroProfileRayPaths parameterized c'tor" << std::endl; + setProfileValues(slice, lats, lons, azimuths); + + // The size of the profile is the maximum number of rays we will have, + // assuming no observations are filtered out. + rays_.reserve(slice.size()); +} +// ----------------------------------------------------------------------------- + +GnssroProfileRayPaths::~GnssroProfileRayPaths() +{ +} +// ----------------------------------------------------------------------------- + +void GnssroProfileRayPaths::setProfileValues( + const GnssroProfileSlice & slice, + const std::vector & lats, + const std::vector & lons, + const std::vector & azimuths) +{ + oops::Log::trace() << "GnssroProfileRayPaths::setProfileValues - enter" << std::endl; + if (lats.size() == 0) { + return; + } + + // Define a scratch vector, then compute an approximate median for each array. + std::vector tmp_v(slice.size()); + profileLat_ = approx_median(lats, slice.start(), slice.end(), tmp_v); + profileLon_ = approx_median(lons, slice.start(), slice.end(), tmp_v); + profileAzimuth_ = approx_median(azimuths, slice.start(), slice.end(), tmp_v); + + // Lats, lons, and azimuths should all have the same size (number of samples + // in all the obs passed to this rank). The slice size is the number of samples + // in a single profile, which must be <= this size. + assert(lons.size() == lats.size()); + assert(azimuths.size() == lats.size()); + assert(slice.size() <= lats.size()); + oops::Log::debug() << "GnssroProfileRayPaths seqNum=" << slice.seqNum() << ": lat " + << *(std::min_element(lats.cbegin() + slice.start(), lats.cbegin() + slice.end())) << " to " + << *(std::max_element(lats.cbegin() + slice.start(), lats.cbegin() + slice.end())) << " (" + << profileLat_ << "), lon " + << *(std::min_element(lons.cbegin() + slice.start(), lons.cbegin() + slice.end())) << " to " + << *(std::max_element(lons.cbegin() + slice.start(), lons.cbegin() + slice.end())) << " (" + << profileLon_ << "), azimuth " + << *(std::min_element(azimuths.cbegin() + slice.start(), azimuths.cbegin() + slice.end())) + << " to " + << *(std::max_element(azimuths.cbegin() + slice.start(), azimuths.cbegin() + slice.end())) + << " (" << profileAzimuth_ << ")" << std::endl; + oops::Log::trace() << "GnssroProfileRayPaths::setProfileValues - exit" << std::endl; +} +// ----------------------------------------------------------------------------- + +void GnssroProfileRayPaths::addRay(std::unique_ptr && ray) noexcept +{ + rays_.push_back(std::move(ray)); +} +// ----------------------------------------------------------------------------- + +// Count all nodes in all rays in the profile. +std::size_t GnssroProfileRayPaths::totalNodes() const +{ + std::size_t nodeCnt = 0; + for (const_iterator rayItr = cbegin(); rayItr != cend(); ++rayItr) + { + const GnssroProfileRayPaths::Ray_& ray = *rayItr; + nodeCnt += ray->numNodes(); + } + return nodeCnt; +} +// ----------------------------------------------------------------------------- + +// Count all sampled locations in all rays in the profile. +std::size_t GnssroProfileRayPaths::totalSampledLocations() const +{ + std::size_t sampLocs = 0; + for (const_iterator rayItr = cbegin(); rayItr != cend(); ++rayItr) + { + const GnssroProfileRayPaths::Ray_& ray = *rayItr; + sampLocs += ray->numLocations(); + } + return sampLocs; +} +// ----------------------------------------------------------------------------- + + +} // namespace ufo + diff --git a/src/ufo/operators/gnssro/utils/GnssroProfileRayPaths.h b/src/ufo/operators/gnssro/utils/GnssroProfileRayPaths.h new file mode 100644 index 000000000..0e840b538 --- /dev/null +++ b/src/ufo/operators/gnssro/utils/GnssroProfileRayPaths.h @@ -0,0 +1,90 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#ifndef UFO_OPERATORS_GNSSRO_UTILS_GNSSROPROFILERAYPATHS_H_ +#define UFO_OPERATORS_GNSSRO_UTILS_GNSSROPROFILERAYPATHS_H_ + +#include // For std::size_t +#include +#include +#include "ufo/operators/gnssro/utils/GnssroProfileSlice.h" +#include "ufo/operators/gnssro/utils/GnssroRayPath.h" + +namespace ufo { + +// ----------------------------------------------------------------------------- +/// Class to represent represent the ray paths associated with a single GNSSRO +/// profile. +class GnssroProfileRayPaths +{ + public: + typedef std::unique_ptr Ray_; + typedef std::vector Rays_; + typedef std::vector::const_iterator const_iterator; + typedef std::vector::const_reverse_iterator const_reverse_iterator; + + GnssroProfileRayPaths(); + GnssroProfileRayPaths(const GnssroProfileSlice & slice, + const std::vector & lats, + const std::vector & lons, + const std::vector & azimuths); + ~GnssroProfileRayPaths(); + GnssroProfileRayPaths(const GnssroProfileRayPaths & other) = delete; + GnssroProfileRayPaths(GnssroProfileRayPaths && other) noexcept = delete; + GnssroProfileRayPaths & operator=(const GnssroProfileRayPaths & other) = delete; + GnssroProfileRayPaths & operator=(GnssroProfileRayPaths && other) noexcept = delete; + + void setProfileValues(const GnssroProfileSlice & slice, + const std::vector & lats, + const std::vector & lons, + const std::vector & azimuths); + + // Ray accessors. + void addRay(std::unique_ptr && ray) noexcept; + std::size_t numRays() const {return rays_.size();} + const Ray_ & ray(std::size_t rayIdx) const {return rays_[rayIdx];} + Ray_ & ray(std::size_t rayIdx) {return rays_[rayIdx];} + + // Ray iterators + const Rays_& rays() const {return rays_;} + const_iterator cbegin() const {return rays_.cbegin();} + const_iterator cend() const {return rays_.cend();} + const_reverse_iterator crbegin() const {return rays_.crbegin();} + const_reverse_iterator crend() const {return rays_.crend();} + + // Profile metadata accessors. + int seqNum() const {return seqNum_;} + float profileLat() const {return profileLat_;} + float profileLon() const {return profileLon_;} + float profileAzimuth() const {return profileAzimuth_;} + + // Aggregate quantities over all rays in RO profile + std::size_t totalNodes() const; + std::size_t totalSampledLocations() const; + + private: + // These are representative values for the entire profile (all in degrees). + int seqNum_; // Profile identifier: unique within an input obs file. + float profileLat_; + float profileLon_; + float profileAzimuth_; + Rays_ rays_; +}; + +} // namespace ufo + +#endif // UFO_OPERATORS_GNSSRO_UTILS_GNSSROPROFILERAYPATHS_H_ diff --git a/src/ufo/operators/gnssro/utils/GnssroProfileSlice.h b/src/ufo/operators/gnssro/utils/GnssroProfileSlice.h new file mode 100644 index 000000000..c2a3c823e --- /dev/null +++ b/src/ufo/operators/gnssro/utils/GnssroProfileSlice.h @@ -0,0 +1,77 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#ifndef UFO_OPERATORS_GNSSRO_UTILS_GNSSROPROFILESLICE_H_ +#define UFO_OPERATORS_GNSSRO_UTILS_GNSSROPROFILESLICE_H_ + +#include // For std::size_t +#include "oops/util/DateTime.h" + +namespace ufo { + +// ----------------------------------------------------------------------------- +/// Class to represent represent the boundaries in the samples of observed +/// GNSSRO data between separate, largely vertical profiles. +/// These slices are created by the GnssroProfileExtractor class. +/// +class GnssroProfileSlice { + public: + GnssroProfileSlice() + : seqNum_(-1) + , time_(19000101, 0) + , start_(0) + , end_(0) {} + GnssroProfileSlice(int seqNum, const util::DateTime& time, size_t start, size_t end) + : seqNum_(seqNum) + , time_(time) + , start_(start) + , end_(end) {} + GnssroProfileSlice(const GnssroProfileSlice & other) + : seqNum_(other.seqNum_) + , time_(other.time_) + , start_(other.start_) + , end_(other.end_) {} + ~GnssroProfileSlice() {} + GnssroProfileSlice & operator=(const GnssroProfileSlice & other) { + if (this != &other) { + seqNum_ = other.seqNum_; + time_ = other.time_; + start_ = other.start_; + end_ = other.end_; + } + return *this; + } + + // Accessors. + int seqNum() const {return seqNum_;} + std::size_t start() const {return start_;} + util::DateTime time() const {return time_;} + std::size_t end() const {return end_;} + std::size_t size() const {return end_ - start_;} + + private: + // A slice includes all indices start_ <= idx < end_, + // i.e. inclusive of start and exclusive of end. + int seqNum_; + util::DateTime time_; + std::size_t start_; + std::size_t end_; +}; + +} // namespace ufo + +#endif // UFO_OPERATORS_GNSSRO_UTILS_GNSSROPROFILESLICE_H_ diff --git a/src/ufo/operators/gnssro/utils/GnssroRayPath.cc b/src/ufo/operators/gnssro/utils/GnssroRayPath.cc new file mode 100644 index 000000000..27003429c --- /dev/null +++ b/src/ufo/operators/gnssro/utils/GnssroRayPath.cc @@ -0,0 +1,144 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +// Define compilation symbol DEBUG_GNSSRO_RAY_PATH (e.g. with -D compiler option) +// to enable debug messages in this file. +#include +#include +#include +#include "oops/util/Logger.h" +#include "ufo/operators/gnssro/utils/GnssroRayPath.h" + +namespace ufo { + +template +static void eraseVectorFromBothEnds(std::vector & v, std::uint32_t numNodesToRemove) +{ + // Erase from end first, since this does not require moving values. + v.erase(v.end() - numNodesToRemove, v.end()); + // Erase from beginning second. + v.erase(v.begin(), v.begin() + numNodesToRemove); +} + +// ----------------------------------------------------------------------------- +/// Class to represent represent the path of a single GNSSRO ray. +/// Typically, each sample in a set of RO observations becomes a ray. +GnssroRayPath::GnssroRayPath( + std::uint32_t numNodes, + float tpLat, + float tpLon, + float azimuth, + float tpAlt, + bool assumeSphericalSymmetry) + : tpLat_(tpLat) + , tpLon_(tpLon) + , azimuth_(azimuth) + , tpAlt_(tpAlt) + , totalLength_(0.0) + , tpNodeIdx_(std::numeric_limits::max()) + , assumeSphericalSymmetry_(assumeSphericalSymmetry) + , geomHeights_(numNodes + 1) // Dimensioned by edges + , latitudes_() // Dimensioned by nodes (not used w/ SphericalSym) + , longitudes_() // Dimensioned by nodes (not used w/ SphericalSym) + , segmentLengths_(numNodes) // Dimensioned by nodes (center points) + , beginLocIdx_(0) + , endLocIdx_(0) + , traj_() +{ + // Only allocate lat-lon vectors if we are not assuming spherical symmetry. + // Otherwise, we assume all nodes are located at the tangent point. + if (!assumeSphericalSymmetry_) { + latitudes_.assign(numNodes, 0.0); + longitudes_.assign(numNodes, 0.0); + } +} + +// ----------------------------------------------------------------------------- + +GnssroRayPath::~GnssroRayPath() +{ +} +// ----------------------------------------------------------------------------- + +void GnssroRayPath::removeUnusedNodes(std::uint32_t edgeOffset) +{ + // Compute nodes to remove from the front and the back of each vector. + std::uint32_t tpIdx = centerNodeIdx(); + std::uint32_t numNodesToRemove = tpIdx - edgeOffset; + + if (numNodesToRemove == 0) + return; + + eraseVectorFromBothEnds(geomHeights_, numNodesToRemove); + eraseVectorFromBothEnds(segmentLengths_, numNodesToRemove); + if (!assumeSphericalSymmetry_) { + eraseVectorFromBothEnds(latitudes_, numNodesToRemove); + eraseVectorFromBothEnds(longitudes_, numNodesToRemove); + + // Check for off-by-one error loading these arrays. + if (geomHeights_.size() > 0 && latitudes_[0] == 0.0 && longitudes_[0] == 0.0 + && geomHeights_[0] == 0.0) + { + oops::Log::error() << "GnssroRayPath: first lat/lon/hgt == 0.0 after erasing numNodes=" + << numNodesToRemove << " from each side of tp to shrink ray to " << latitudes_.size() + << " nodes" << std::endl; + } + } + + // Reset the tangent point node index. Since we removed nodes symmetrically + // from both ends, the tangent point will still be the center node. + setTpNodeIdx(centerNodeIdx()); +} +// ----------------------------------------------------------------------------- + +void GnssroRayPath::shrinkNodesTo(std::uint32_t newNumNodes) +{ + if (newNumNodes < numNodes()) { + uint32_t numToRemove = numNodes() - newNumNodes; + geomHeights_.erase(geomHeights_.end() - numToRemove, geomHeights_.end()); + segmentLengths_.erase(segmentLengths_.end() - numToRemove, segmentLengths_.end()); + if (!assumeSphericalSymmetry_) { + latitudes_.erase(latitudes_.end() - numToRemove, latitudes_.end()); + longitudes_.erase(longitudes_.end() - numToRemove, longitudes_.end()); + } + } else if (newNumNodes > numNodes()) { + oops::Log::error() << "GnssroRayPath::shrinkNodesTo: newNumNodes=" << newNumNodes + << " > oldNumNodes=" << numNodes() << "; not changing ray size" << std::endl; + } + // Compute the total length of the ray. We only removed nodes after the last used node, + // so the tangent point node will not have changed. + updateTotalLength(); + return; +} +// ----------------------------------------------------------------------------- + +void GnssroRayPath::updateTotalLength() +{ + double newLength = std::accumulate(segmentLengths_.begin(), segmentLengths_.end(), 0.0); + totalLength_ = newLength; +} +// ----------------------------------------------------------------------------- + +void GnssroRayPath::setLocIdxRange(std::uint32_t beginLocIdx, std::uint32_t endLocIdx) +{ + beginLocIdx_ = beginLocIdx; + endLocIdx_ = endLocIdx; +} +// ----------------------------------------------------------------------------- + +} // namespace ufo + diff --git a/src/ufo/operators/gnssro/utils/GnssroRayPath.h b/src/ufo/operators/gnssro/utils/GnssroRayPath.h new file mode 100644 index 000000000..c9036e493 --- /dev/null +++ b/src/ufo/operators/gnssro/utils/GnssroRayPath.h @@ -0,0 +1,142 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#ifndef UFO_OPERATORS_GNSSRO_UTILS_GNSSRORAYPATH_H_ +#define UFO_OPERATORS_GNSSRO_UTILS_GNSSRORAYPATH_H_ + +#include +#include // For std::size_t +#include // For uint32_t +#include +#include "ufo/operators/gnssro/utils/GnssroRayTrajectory.h" + +namespace ufo { + +// ----------------------------------------------------------------------------- +/// Class to represent represent the path of a single GNSSRO ray. +/// Typically, each sample in a set of RO observations becomes a ray. +class GnssroRayPath +{ + public: + typedef std::vector DoubleArray; + typedef std::vector RealArray; + typedef std::vector UIntArray; + + GnssroRayPath(uint32_t num_nodes, float tpLat, float tpLon, float azimuth, + float tpAlt, bool assumeSphericalSymmetry = false); + GnssroRayPath() = default; + GnssroRayPath(const GnssroRayPath & other) = delete; + GnssroRayPath(GnssroRayPath && other) noexcept = delete; + GnssroRayPath & operator=(const GnssroRayPath & other) = delete; + GnssroRayPath & operator=(GnssroRayPath && other) noexcept = delete; + ~GnssroRayPath(); + + // Accessors. + std::uint32_t centerNodeIdx() const {return segmentLengths_.size() / 2;} + std::uint32_t numEdges() const {return geomHeights_.size();} + std::uint32_t numNodes() const {return segmentLengths_.size();} + std::uint32_t numLocations() const { + return (assumeSphericalSymmetry_ ? 1 : numNodes()); + } + std::uint32_t edgesInHalfRay() const {return numEdges() / 2;} + bool assumeSphericalSymmetry() const {return assumeSphericalSymmetry_;} + const GnssroRayTrajectory& traj() const {return traj_;} + GnssroRayTrajectory& traj() {return traj_;} + + // Read-only accessors to array data. + float geomHgt(uint32_t edgeIdx) const {return geomHeights_[edgeIdx];} + float lat(uint32_t nodeIdx) const { + return (assumeSphericalSymmetry_ ? tpLat_ : latitudes_[nodeIdx]); + } + float lon(uint32_t nodeIdx) const { + return (assumeSphericalSymmetry_ ? tpLon_ : longitudes_[nodeIdx]); + } + double segLen(uint32_t nodeIdx) const {return segmentLengths_[nodeIdx];} + + // Read only accessors for the vectors + const RealArray& lats() const {return latitudes_;} + const RealArray& lons() const {return longitudes_;} + + // Read-write accessors to array data. + float & geomHgt(uint32_t edgeIdx) { + assert(edgeIdx < geomHeights_.size()); + return geomHeights_[edgeIdx]; + } + float & lat(uint32_t nodeIdx) { + if (assumeSphericalSymmetry_) { + return tpLat_; + } else { + assert(nodeIdx < latitudes_.size()); + return latitudes_[nodeIdx]; + } + } + float & lon(uint32_t nodeIdx) { + if (assumeSphericalSymmetry_) { + return tpLon_; + } else { + assert(nodeIdx < longitudes_.size()); + return longitudes_[nodeIdx]; + } + } + double & segLen(uint32_t nodeIdx) { + assert(nodeIdx < segmentLengths_.size()); + return segmentLengths_[nodeIdx]; + } + + // Remove elements from front and back of vector so there are edgeOffset points + // around the center. + void removeUnusedNodes(std::uint32_t edgeOffset); + + // Resize to specified number of nodes by removing elements from the end, + // if necessary. + void shrinkNodesTo(std::uint32_t newNumNodes); + + // Method to compute the total length from the segment lengths. + void updateTotalLength(); + + float tpLat() const {return tpLat_;} + float tpLon() const {return tpLon_;} + float azimuth() const {return azimuth_;} + float tpAlt() const {return tpAlt_;} + float totalLength() const {return totalLength_;} + std::uint32_t tpNodeIdx() const {return tpNodeIdx_;} + void setTpNodeIdx(std::uint32_t tpNodeIdx) {tpNodeIdx_ = tpNodeIdx;} + std::uint32_t beginLocIdx() const {return beginLocIdx_;} + std::uint32_t endLocIdx() const {return endLocIdx_;} + void setLocIdxRange(std::uint32_t beginLocIdx, std::uint32_t endLocIdx); + + private: + // Member data. + float tpLat_; // Latitude of tangent point defining this ray. + float tpLon_; // Longitude of tangent point defining this ray. + float azimuth_; // Azimuth angle defining this ray. + float tpAlt_; // Altitude of tangent point (Geometric height above MSL in m) + double totalLength_; // sum of all segment lengths. + uint32_t tpNodeIdx_; // Node index of the tangent point. + bool assumeSphericalSymmetry_; // Each node with use the tangent point lat, lon + RealArray geomHeights_; // For bottom edges of nodes (m). Last edge is top of highest node + RealArray latitudes_; // For edges (degrees) + RealArray longitudes_; // For edges (degrees) + DoubleArray segmentLengths_; // Length of segments centered on node (m) + std::uint32_t beginLocIdx_; // begin index for pathsGroupedByLocation + std::uint32_t endLocIdx_; // end index for pathsGroupedByLocation + GnssroRayTrajectory traj_; // 2D Trajectory used in TLAD +}; + +} // namespace ufo + +#endif // UFO_OPERATORS_GNSSRO_UTILS_GNSSRORAYPATH_H_ diff --git a/src/ufo/operators/gnssro/utils/GnssroRayPathGenerator.cc b/src/ufo/operators/gnssro/utils/GnssroRayPathGenerator.cc new file mode 100644 index 000000000..bebbe747d --- /dev/null +++ b/src/ufo/operators/gnssro/utils/GnssroRayPathGenerator.cc @@ -0,0 +1,42 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#include +#include +#include "eckit/exception/Exceptions.h" +#include "ufo/operators/gnssro/utils/GnssroRayPathGenerator.h" +#include "ufo/operators/gnssro/utils/StraightLineRayPathGenerator.h" + +namespace ufo { + +// Factory method for GnssroRayPathGenerator objects. +std::unique_ptr +GnssroRayPathGenerator::create(const ioda::ObsSpace & odb, + const GnssroRayPathParameters & params) +{ + if (params.rayPathGenType() == StraightLineRayPathGenerator::name()) { + return std::make_unique(odb, params); + } else { + throw eckit::BadValue( + "Unsupported GnssroRayPathGenerator subclass: " + params.rayPathGenType() + + "; the only supported type is " + StraightLineRayPathGenerator::name(), + Here()); + } +} + +} // namespace ufo + diff --git a/src/ufo/operators/gnssro/utils/GnssroRayPathGenerator.h b/src/ufo/operators/gnssro/utils/GnssroRayPathGenerator.h new file mode 100644 index 000000000..eefc9b9de --- /dev/null +++ b/src/ufo/operators/gnssro/utils/GnssroRayPathGenerator.h @@ -0,0 +1,64 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#ifndef UFO_OPERATORS_GNSSRO_UTILS_GNSSRORAYPATHGENERATOR_H_ +#define UFO_OPERATORS_GNSSRO_UTILS_GNSSRORAYPATHGENERATOR_H_ + +#include +#include + +#include "ioda/ObsSpace.h" +#include "ufo/operators/gnssro/utils/GnssroProfileRayPaths.h" +#include "ufo/operators/gnssro/utils/GnssroProfileSlice.h" +#include "ufo/operators/gnssro/utils/GnssroRayPathParameters.h" + +namespace ufo { + +// ----------------------------------------------------------------------------- +/// Base class for GnssroRayPathGenerator +/// +/// The developer should define a subclass of GnssroRayPathGenerator that +/// implements the makeRayPath method. +/// +class GnssroRayPathGenerator { + public: + static std::string name() {return "GnssroRayPathGenerator";} + + // Factory method. + static std::unique_ptr create( + const ioda::ObsSpace & odb, const GnssroRayPathParameters & params); + + // Constructor/Destructor + explicit GnssroRayPathGenerator(const ioda::ObsSpace & odb, + const GnssroRayPathParameters & params) + : odb_(odb) + , params_(params) {} + virtual ~GnssroRayPathGenerator() {} + + /// \brief Return a GnssroRayPath object based on observation data read for the specified + /// GNSSRO profile from the ObsSpace. + virtual std::unique_ptr makeProfileRayPaths( + const GnssroProfileSlice & profileSlice) = 0; + + protected: + const ioda::ObsSpace& odb_; + GnssroRayPathParameters params_; +}; + +} // namespace ufo + +#endif // UFO_OPERATORS_GNSSRO_UTILS_GNSSRORAYPATHGENERATOR_H_ diff --git a/src/ufo/operators/gnssro/utils/GnssroRayPathOrchestrator.cc b/src/ufo/operators/gnssro/utils/GnssroRayPathOrchestrator.cc new file mode 100644 index 000000000..156a9cbae --- /dev/null +++ b/src/ufo/operators/gnssro/utils/GnssroRayPathOrchestrator.cc @@ -0,0 +1,389 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#include +#include "ioda/ObsSpace.h" +#include "ufo/operators/gnssro/utils/GnssroGeoVaLs.h" +#include "ufo/operators/gnssro/utils/GnssroRayPathOrchestrator.h" +#include "ufo/operators/gnssro/utils/RefractivityCalculator.h" +#include "ufo/utils/Constants.h" +#include "ufo/utils/VertInterp.interface.h" +#include "ufo/variabletransforms/Formulas.h" + +namespace ufo { + +// ----------------------------------------------------------------------------- +static const char className[] = "GnssroRayPathOrchestrator"; + +GnssroRayPathOrchestrator::GnssroRayPathOrchestrator( + const ioda::ObsSpace & odb, + const GnssroRayPathParameters & params) + : odb_(odb) + , params_(params) + , extractor_(odb_) + , generator_() + , profiles_() +{ + oops::Log::trace() << className << " ctor enter" << std::endl; + + oops::Log::debug() << className << ": " << extractor_.slices().size() + << " profiles found" << std::endl; + + if (extractor_.is_empty()) { + return; // No observations to process. + } + + // Construct the ray path generator. + generator_ = GnssroRayPathGenerator::create(odb, params_); + + // Create a set of ray paths for each profile. + std::size_t numRays = 0; + for (GnssroProfileExtractor::const_iterator itr = extractor_.cbegin(); + itr != extractor_.cend(); ++itr) + { + oops::Log::debug() << className << ": Profile slice " + << itr->seqNum() << " from " << itr->start() + << " to " << itr->end() << std::endl; + Profile_ profileRayPaths = generator_->makeProfileRayPaths(*itr); + numRays += profileRayPaths->numRays(); + profiles_.push_back(std::move(profileRayPaths)); + } + + oops::Log::debug() << className << ": generated " << numRays << " rays for " + << profiles_.size() << " profiles" << std::endl; + oops::Log::trace() << className << " ctor exit" << std::endl; +} +// ----------------------------------------------------------------------------- + +GnssroRayPathOrchestrator::~GnssroRayPathOrchestrator() +{ + oops::Log::trace() << className << " dtor called" << std::endl; +} +// ----------------------------------------------------------------------------- + +std::size_t GnssroRayPathOrchestrator::numSampledLocs() const +{ + // We will have one sampled location for each node. + // Iterate over RO profiles + std::size_t numLocs = 0; + for (Profiles_::const_iterator p = profiles_.cbegin(); p != profiles_.cend(); ++p) + { + // Iterate over rays in a profile + const Profile_& roProfile = *p; + for (std::size_t rayIdx = 0; rayIdx < roProfile->numRays(); ++rayIdx) + { + numLocs += roProfile->ray(rayIdx)->numLocations(); + } + } + return numLocs; +} +// ----------------------------------------------------------------------------- + +std::size_t GnssroRayPathOrchestrator::getSampledLatLonTimes( + std::vector& lats, + std::vector& lons, + std::vector& times) const +{ + const char funcName[] = "GnssroRayPathOrchestrator::getSampledLatLonTimes"; + lats.clear(); + lons.clear(); + times.clear(); + std::size_t numLocs = numSampledLocs(); + lats.reserve(numLocs); + lons.reserve(numLocs); + times.reserve(numLocs); + + oops::Log::debug() << funcName << ": Number of profiles in extractor=" + << extractor_.slices().size() << ", and number of profiles in profiles= " + << profiles_.size() << std::endl; + + GnssroProfileExtractor::const_iterator extItr = extractor_.cbegin(); + for (Profiles_::const_iterator itr = profiles_.cbegin(); itr != profiles_.cend(); ++itr) + { + const Profile_& roProfile = *itr; + for (std::size_t rayIdx = 0; rayIdx < roProfile->numRays(); ++rayIdx) + { + const Ray_& ray = roProfile->ray(rayIdx); + + // Insert the lat-lon values, as well as the characteristic time of the profile + // for each sampled location associated with the profile. + if (ray->assumeSphericalSymmetry()) { + // All nodes share one location and time. + lats.insert(lats.end(), ray->tpLat()); + lons.insert(lons.end(), ray->tpLon()); + times.insert(times.end(), extItr->time()); + } else { + // Each node has a separate location, but the same time. + lats.insert(lats.end(), ray->lats().cbegin(), ray->lats().cend()); + lons.insert(lons.end(), ray->lons().cbegin(), ray->lons().cend()); + times.insert(times.end(), ray->lats().size(), extItr->time()); + } + } + ++extItr; + } + + if (lats.size() != times.size()) + { + oops::Log::error() << funcName << "Number of locations and times are inconsistent" + << ": lats (" << lats.size() << "), times (" << times.size() + << "); they must be the same" << std::endl; + assert(lats.size() == times.size()); + } + + if (lats.size() != numLocs) + { + oops::Log::error() << funcName << "Number of locations and computed numLocs are " + << "inconsistent: lats (" << lats.size() << "), numLocs (" << numLocs + << "); they must be the same" << std::endl; + assert(lats.size() == times.size()); + } + + return numLocs; +} +// ----------------------------------------------------------------------------- + +void GnssroRayPathOrchestrator::fillPathsGroupedByLocation( + std::size_t numObs, + std::vector> & pathsGroupedByLocation) const +{ + const char funcName[] = "GnssroRayPathOrchestrator::fillPathsGroupedByLocation"; + pathsGroupedByLocation.resize(numObs); + + std::size_t obIdx = 0; + std::size_t locIdx = 0; + for (Profiles_::const_iterator itr = profiles_.cbegin(); itr != profiles_.cend(); ++itr) + { + const Profile_& roProfile = *itr; + for (std::size_t rayIdx = 0; rayIdx < roProfile->numRays(); ++rayIdx) + { + const Ray_& ray = roProfile->ray(rayIdx); + + // Do not go out of bounds of the array, but keep a count of the + // number of rays in all the profiles. + if (obIdx >= numObs) + { + ++obIdx; + continue; + } + + // Set the range of sampled location indices used by this ray. + // It will have one sampled location for each lat-lon location in the ray. + // This is either 1 location per node, or a single node for the entire + // ray if the spherical symmetry assumption is in place. + std::uint32_t rayNumLocs = ray->numLocations(); + pathsGroupedByLocation[obIdx].begin = locIdx; + pathsGroupedByLocation[obIdx].end = locIdx + rayNumLocs; + oops::Log::debug() << funcName << ": obIdx=" << obIdx << " of " + << numObs << ", add " << rayNumLocs << " sampled locs " + << "[" << pathsGroupedByLocation[obIdx].begin + << "," << pathsGroupedByLocation[obIdx].end << ")" << std::endl; + locIdx += rayNumLocs; + ++obIdx; + } + } + + if (obIdx != numObs) + { + oops::Log::error() << funcName << ": found " << obIdx + << " ray paths to specify in pathsGroupedByLocation, but there are " + << numObs << " observations to specify; these must be the same." + << std::endl; + } + return; +} +// ----------------------------------------------------------------------------- + +std::size_t GnssroRayPathOrchestrator::totalNumRays() const +{ + std::size_t numRays = 0; + for (Profiles_::const_iterator p = profiles_.cbegin(); p != profiles_.cend(); ++p) + { + const Profile_& profile = *p; + numRays += profile->numRays(); + } + return numRays; +} +// ----------------------------------------------------------------------------- + +double GnssroRayPathOrchestrator::pressureInterpApply( + int nlevs, + const double * pres, + const double * temp, + int wi, + double wf) +{ + // Interpolate pressure in a way that ensures both the bottom + // and top of the layer are used. By contrast, the hydrostatic formulation + // used in RefNCEP only uses pressure at the bottom level. Using pressure + // at both bottom and top ensures better agreement between non-linear and + // linear models if pressure is only perturbed at one of those levels. + // + // Ideally, we use a form of the hypsometric equation with a + // non-zero lapse rate, but we fall back to isothermal interpolation + // if temperature is effectively constant over the layer. + // This formulation is similar to pressure interpolation used in ROPP. + const double NEAR_ZERO = 1e-10; + int cwi = wi - 1; // convert 1-based index wi to 0-based cwi + double interpPres = util::missingValue(); + if (std::abs(temp[wi] - temp[cwi]) < NEAR_ZERO) { + // layer is isothermal: use isothermal interpolation + if (pres[cwi] <= 0.0 || pres[wi] <= 0.0) { + // Min pressure <= 0.0: use linear interpolation + interpPres = wf * pres[cwi] + (1.0 - wf) * pres[wi]; + } else { + // Pressure is reasonable: use exponential interpolation + interpPres = pres[cwi] * std::exp(wf * std::log(pres[wi]/pres[cwi])); + } + } else { + // Use adapted hypsometric equation, scaled to temperature and pressure + // values at the layer's bottom and top. + double gamma = std::log(temp[wi]/temp[cwi]) + / std::log(pres[wi]/pres[cwi]); + double interpTemp = wf * temp[cwi] + (1.0 - wf) * temp[wi]; + interpPres = pres[cwi] * std::pow(interpTemp / temp[cwi], 1.0 / gamma); + } + return interpPres; +} + +// ----------------------------------------------------------------------------- +// Class method for computing refractivity from interpolated T, P, and Q. +// +double GnssroRayPathOrchestrator::modelRefractivity( + float geomHgt, float lat, const GnssroGeoVaLs & ggv, std::size_t profileIdx, + const std::unique_ptr & refrCalc) +{ + // Get references to model vertical profiles (columns) for this profileIdx. + const GnssroGeoVaLs::Profile_& gphProf = ggv.gphProfile(profileIdx); + const GnssroGeoVaLs::Profile_& tempProf = ggv.tempProfile(profileIdx); + const GnssroGeoVaLs::Profile_& sphumProf = ggv.sphumProfile(profileIdx); + const GnssroGeoVaLs::Profile_& presProf = ggv.presProfile(profileIdx); + + // Convert geometric height to geopotential height + float gph = formulas::Geometric_to_Geopotential_Height(lat, geomHgt); + + // Compute vertical interpolation weights within a model profile for the + // target geopotential height + int wi; // Base vertical index for vertical interpolation (1-based). + double wf; // Weight for base index (rest of weight goes to idx wi + 1). + vert_interp_weights_f90(ggv.nlevs(), gph, gphProf.data(), wi, wf); + + // Linearly interpolate temp and sphum to target gph + double temp; // Interpolated air temperature in K + double sphum; // Interpolated specific humidity in kg/kg + vert_interp_apply_f90(ggv.nlevs(), tempProf.data(), temp, wi, wf); + vert_interp_apply_f90(ggv.nlevs(), sphumProf.data(), sphum, wi, wf); + double pres = pressureInterpApply( + ggv.nlevs(), presProf.data(), tempProf.data(), wi, wf); + + // Compute model refractivity from the interpolated temp,sphum,pres. + double refr = refrCalc->N(temp, sphum, pres); + return refr; +} + +// ----------------------------------------------------------------------------- +// Class method for computing derivatives of refractivity from interpolated T, P, and Q. +GnssroRayPathOrchestrator::TrajTuple +GnssroRayPathOrchestrator::modelTrajectory( + float geomHgt, float lat, const GnssroGeoVaLs & ggv, std::size_t profileIdx, + const std::unique_ptr & refrCalc) +{ + // Get references to model vertical profiles (columns) for this profileIdx. + const GnssroGeoVaLs::Profile_& gphProf = ggv.gphProfile(profileIdx); + const GnssroGeoVaLs::Profile_& tempProf = ggv.tempProfile(profileIdx); + const GnssroGeoVaLs::Profile_& sphumProf = ggv.sphumProfile(profileIdx); + const GnssroGeoVaLs::Profile_& presProf = ggv.presProfile(profileIdx); + + // Convert geometric height to geopotential height + float gph = formulas::Geometric_to_Geopotential_Height(lat, geomHgt); + + // Compute vertical interpolation weights within a model profile for the + // target geopotential height + int wi; // Base vertical index for vertical interpolation (1-based). + double wf; // Weight for base index (rest of weight goes to idx wi + 1). + vert_interp_weights_f90(ggv.nlevs(), gph, gphProf.data(), wi, wf); + + // Linearly interpolate temp and sphum to target gph + double temp; // Interpolated air temperature in K + double sphum; // Interpolated specific humidity in kg/kg + vert_interp_apply_f90(ggv.nlevs(), tempProf.data(), temp, wi, wf); + vert_interp_apply_f90(ggv.nlevs(), sphumProf.data(), sphum, wi, wf); + double pres = pressureInterpApply( + ggv.nlevs(), presProf.data(), tempProf.data(), wi, wf); + + // Compute derivatives of model refractivity from the interpolated temp,sphum,pres. + double dNdT = refrCalc->dNdT(temp, sphum, pres); + double dNdP = refrCalc->dNdP(temp, sphum, pres); + double dNdQ = refrCalc->dNdQ(temp, sphum, pres); + + // Return the trajectory elements as a tuple. + TrajTuple traj(dNdT, dNdQ, dNdP, wf, wi); + return traj; +} + +// ----------------------------------------------------------------------------- +// Associated ostream operator +std::ostream& operator<<(std::ostream& os, const GnssroRayPathOrchestrator& orc) +{ + // Accumulate stats over all the profiles managed by the orchestrator + const GnssroRayPathOrchestrator::Profiles_& profiles = orc.profiles(); + std::size_t numProfiles = profiles.size(); // Total RO profiles/seqNums + std::size_t numRays = 0; // Total observation samples turned into rays. + std::size_t numNodes = 0; + std::size_t numSampledLocations = 0; + std::size_t maxNodesForRay = 0; + int maxNodeSeqNum = 0; + std::size_t maxNodeRayIdx = 0; + + for (GnssroRayPathOrchestrator::Profiles_::const_iterator itr = profiles.cbegin(); + itr != profiles.cend(); ++itr) + { + const GnssroRayPathOrchestrator::Profile_& roProfile = *itr; + numRays += roProfile->numRays(); + for (std::size_t rayIdx = 0; rayIdx < roProfile->numRays(); ++rayIdx) + { + const GnssroRayPathOrchestrator::Ray_& ray = roProfile->ray(rayIdx); + + numSampledLocations += ray->numLocations(); + numNodes += ray->numNodes(); + if (ray->numNodes() > maxNodesForRay) { + maxNodesForRay = ray->numNodes(); + maxNodeSeqNum = roProfile->seqNum(); + maxNodeRayIdx = rayIdx; + } + } + } + + // Ratio of sampledLocations (columns from the model data) to the number of + // tangent point observations. A 1D operator would return 1.0, so this is + // a metric of how efficiently we are sampling the model data. + float sampledLocationsPerRay = static_cast(numSampledLocations) + / static_cast(numRays); + + os << "ROprofiles=" << numProfiles + << ", rays=" << numRays + << ", numSampLocs=" << numSampledLocations + << ", sampLocs/ray=" << sampledLocationsPerRay + << ", numRayNodes=" << numNodes + << ", maxNodesForRay=" << maxNodesForRay + << " (seqNum=" << maxNodeSeqNum + << ", rayIdx=" << maxNodeRayIdx << ")"; + return os; +} +// ----------------------------------------------------------------------------- + +} // namespace ufo + diff --git a/src/ufo/operators/gnssro/utils/GnssroRayPathOrchestrator.h b/src/ufo/operators/gnssro/utils/GnssroRayPathOrchestrator.h new file mode 100644 index 000000000..d62e18381 --- /dev/null +++ b/src/ufo/operators/gnssro/utils/GnssroRayPathOrchestrator.h @@ -0,0 +1,165 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#ifndef UFO_OPERATORS_GNSSRO_UTILS_GNSSRORAYPATHORCHESTRATOR_H_ +#define UFO_OPERATORS_GNSSRO_UTILS_GNSSRORAYPATHORCHESTRATOR_H_ + +#include +#include +#include +#include "ioda/ObsSpace.h" +#include "oops/util/DateTime.h" +#include "oops/util/Range.h" +#include "ufo/operators/gnssro/utils/GnssroProfileExtractor.h" +#include "ufo/operators/gnssro/utils/GnssroProfileRayPaths.h" +#include "ufo/operators/gnssro/utils/GnssroProfileSlice.h" +#include "ufo/operators/gnssro/utils/GnssroRayPathGenerator.h" +#include "ufo/operators/gnssro/utils/GnssroRayPathParameters.h" + +namespace ufo { +// Forward class references +class RefractivityCalculator; +class GnssroGeoVaLs; + +// ----------------------------------------------------------------------------- +/// GnssroRayPathOrchestrator +/// +class GnssroRayPathOrchestrator { + public: + // Type definitions + typedef std::unique_ptr Generator_; + typedef std::unique_ptr Profile_; + typedef std::vector Profiles_; + typedef GnssroProfileRayPaths::Ray_ Ray_; + + // Nested class for holding elements of a trajectory for one + // node of a 2D ray path. + class TrajTuple + { + public: + TrajTuple() + : dNdT_(0.0) + , dNdQ_(0.0) + , dNdP_(0.0) + , wf_(0.0) + , wi_(0) + {} + + TrajTuple(double dNdT, double dNdQ, double dNdP, double wf, int wi) + : dNdT_(dNdT) + , dNdQ_(dNdQ) + , dNdP_(dNdP) + , wf_(wf) + , wi_(wi) + {} + + // Read-only accessors. + double dNdT() const {return dNdT_;} + double dNdQ() const {return dNdQ_;} + double dNdP() const {return dNdP_;} + double wf() const {return wf_;} + int wi() const {return wi_;} + + private: + // Member data. + double dNdT_; + double dNdQ_; + double dNdP_; + double wf_; + int wi_; + }; // end of nested class definition for TrajTuple + + /// \brief Class method to apply pressure interpolation. + /// This is consistent with interpolation weights and indices + /// produced by vert_interp_f90. Returns interpolated value + /// from pres using base index wi (1-based) and weight wf + /// as well as information from the temperature profile. + static double pressureInterpApply( + int nlevs, // Number of vertical levels + const double * pres, // Column of pressure values + const double * temp, // Column of temperature values + int wi, // Vertical interp base index (1-based) + double wf); // Linear vertical interp factor + + /// \brief Class method for computing model refractivity. + static double modelRefractivity( + float geomHgt, float lat, const GnssroGeoVaLs & ggv, std::size_t profileIdx, + const std::unique_ptr & refrCalc); + + /// \brief Class method for computing model trajectory. + static TrajTuple modelTrajectory( + float geomHgt, float lat, const GnssroGeoVaLs & ggv, std::size_t profileIdx, + const std::unique_ptr & refrCalc); + + // Constructor/Destructor + explicit GnssroRayPathOrchestrator( + const ioda::ObsSpace & odb, + const GnssroRayPathParameters & params); + ~GnssroRayPathOrchestrator(); + GnssroRayPathOrchestrator() = delete; + GnssroRayPathOrchestrator(const GnssroRayPathOrchestrator &) = delete; + GnssroRayPathOrchestrator(GnssroRayPathOrchestrator &&) noexcept = delete; + GnssroRayPathOrchestrator& operator=(const GnssroRayPathOrchestrator &) = delete; + GnssroRayPathOrchestrator& operator=(GnssroRayPathOrchestrator &&) noexcept = delete; + + // Read-only Accessors + const GnssroRayPathParameters & parameters() const {return params_;} + const GnssroProfileExtractor & extractor() const {return extractor_;} + const Generator_ & generator() const {return generator_;} + const Profiles_ & profiles() const {return profiles_;} + + // Number of paths (vertical columns) requested from the model. + std::size_t numSampledLocs() const; + + // Method to allocate and fill the location and time arrays for the SampledLocations. + // It clears and reallocates the vectors to the necessary number of sampled locations. + // It returns the common allocation size used for all three vectors. + std::size_t getSampledLatLonTimes( + std::vector& lats, + std::vector& lons, + std::vector& times) const; + + /// Sets the range of SampledLocations used by each ray path. + void fillPathsGroupedByLocation( + std::size_t numSampledLocs, + std::vector> & pathsGroupedByLocation) const; + + // Number of rays computed within profile ray paths. + std::size_t totalNumRays() const; + + // Read-write Accessors + Generator_ & generator() {return generator_;} + Profiles_ & profiles() {return profiles_;} + + // ostream operator can access private members. + friend std::ostream& operator<<(std::ostream& os, const GnssroRayPathOrchestrator& orc); + + private: + // Member data + const ioda::ObsSpace& odb_; + GnssroRayPathParameters params_; + GnssroProfileExtractor extractor_; + Generator_ generator_; + Profiles_ profiles_; + std::vector sampledLats_; + std::vector sampledLons_; + std::vector sampledTimes_; +}; + +} // namespace ufo + +#endif // UFO_OPERATORS_GNSSRO_UTILS_GNSSRORAYPATHORCHESTRATOR_H_ diff --git a/src/ufo/operators/gnssro/utils/GnssroRayPathParameters.cc b/src/ufo/operators/gnssro/utils/GnssroRayPathParameters.cc new file mode 100644 index 000000000..4d9149a8f --- /dev/null +++ b/src/ufo/operators/gnssro/utils/GnssroRayPathParameters.cc @@ -0,0 +1,122 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#include +#include +#include "eckit/exception/Exceptions.h" +#include "oops/util/Logger.h" +#include "ufo/operators/gnssro/utils/GnssroRayPathParameters.h" + +namespace ufo { + +static constexpr double KM_TO_METERS = 1000.0; + +// ----------------------------------------------------------------------------- +GnssroRayPathParameters::GnssroRayPathParameters() + : rayPathGenType_(DEFAULT_RAY_PATH_GEN_TYPE) + , approxRayLength_(DEFAULT_APPROX_RAY_LENGTH_KM * KM_TO_METERS) + , expectedRayLength_(0.0) + , top2D_(DEFAULT_TOP2D_KM * KM_TO_METERS) + , horizontalRes_(DEFAULT_HORIZONTAL_RES_KM * KM_TO_METERS) + , nHoriz_(DEFAULT_NHORIZ) +{ + computeDerivedParameters(); +} +// ----------------------------------------------------------------------------- + +GnssroRayPathParameters::GnssroRayPathParameters( + const std::string & rayPathGenType, + double approxRayLengthKm, + double top2DKm, + double resKm, + int nHoriz +) : rayPathGenType_(rayPathGenType) + , approxRayLength_(approxRayLengthKm * KM_TO_METERS) + , expectedRayLength_(0.0) + , top2D_(top2DKm * KM_TO_METERS) + , horizontalRes_(resKm * KM_TO_METERS) + , nHoriz_(nHoriz) +{ + computeDerivedParameters(); +} +// ----------------------------------------------------------------------------- + +void GnssroRayPathParameters::computeDerivedParameters() +{ + if (approxRayLength_ > 0.0) { + // Compute nHoriz_ if max_ray_length was specified. + int originalNHoriz = nHoriz_; + + // Ensure nHoriz_ is an odd number and within reasonable bounds. + int nHorizFloor = static_cast(std::floor(approxRayLength_ / horizontalRes_)); + int nHorizCeil = static_cast(std::ceil(approxRayLength_ / horizontalRes_)); + if (nHorizFloor % 2 != 0) { // Floor is odd + nHoriz_ = nHorizFloor; + } else if (nHorizCeil % 2 != 0) { // Ceiling is odd + nHoriz_ = nHorizCeil; + } else { // Both ceiling and floor are even (original number was even int) + nHoriz_ = nHorizFloor + 1; + } + + // Emit warning if we changed the configured n_horiz. + if (originalNHoriz != DEFAULT_NHORIZ && nHoriz_ != DEFAULT_NHORIZ) { + oops::Log::warning() << "GnssroRayPathParameters: overriding configured nHoriz=" + << originalNHoriz << " with computed value=" << nHoriz_ + << " from approxRayLength=" << approxRayLength_ + << " and horizontalRes=" << horizontalRes_ << std::endl; + } + + // Ensure nHoriz_ is an odd number and within reasonable bounds. + if (nHoriz_ <= 0) { + oops::Log::error() << "GnssroRayPathParameters: adjusting nHoriz=" + << nHoriz_ << " up to 1; calculated from approxRayLength=" + << approxRayLength_ << " and horizontalRes=" << horizontalRes_ + << std::endl; + nHoriz_ = 1; + } + if (nHoriz_ > MAX_NHORIZ) { + oops::Log::error() << "GnssroRayPathParameters: adjusting nHoriz=" + << nHoriz_ << " down to " << MAX_NHORIZ + << "; calculated from approxRayLength=" + << approxRayLength_ << " and horizontalRes=" << horizontalRes_ + << std::endl; + nHoriz_ = MAX_NHORIZ; + } + + oops::Log::debug() << "GnssroRayPathParameters: nHoriz=" << nHoriz_ + << " from approxRayLength=" << approxRayLength_ << " and horizontalRes=" + << horizontalRes_ << std::endl; + } else { + // ray length not configured; use nHoriz + // Ensure nHoriz_ is an odd number and within reasonable bounds. + if (nHoriz_ <= 0 || nHoriz_ > MAX_NHORIZ || nHoriz_ % 2 == 0) { + throw eckit::BadValue( + "GnssroRayPathParameters: Bad value for nhoriz: " + + std::to_string(nHoriz_) + + "; must be an odd number between 1 and " + + std::to_string(MAX_NHORIZ) + " (inclusive)", Here()); + } + } + + // Compute expected ray length + expectedRayLength_ = nHoriz_ * horizontalRes_; + return; +} +// ----------------------------------------------------------------------------- + +} // namespace ufo + diff --git a/src/ufo/operators/gnssro/utils/GnssroRayPathParameters.h b/src/ufo/operators/gnssro/utils/GnssroRayPathParameters.h new file mode 100644 index 000000000..33674c6f8 --- /dev/null +++ b/src/ufo/operators/gnssro/utils/GnssroRayPathParameters.h @@ -0,0 +1,70 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#ifndef UFO_OPERATORS_GNSSRO_UTILS_GNSSRORAYPATHPARAMETERS_H_ +#define UFO_OPERATORS_GNSSRO_UTILS_GNSSRORAYPATHPARAMETERS_H_ + +#include + +namespace ufo { + +// ----------------------------------------------------------------------------- +class GnssroRayPathParameters { + public: + // Class to hold configurable parameters for GnssroRayPathGenerator and + // GnssroRayPathOrchestrator classes. + // This should hold the super set of parameters needed by all subclasses. + static constexpr const char DEFAULT_RAY_PATH_GEN_TYPE[] = "StraightLine"; + static constexpr double DEFAULT_APPROX_RAY_LENGTH_KM = -1.0; + static constexpr double DEFAULT_TOP2D_KM = 90.0; + static constexpr double DEFAULT_HORIZONTAL_RES_KM = 11.0; + static constexpr int DEFAULT_NHORIZ = 11; + static constexpr int MAX_NHORIZ = 67; + + GnssroRayPathParameters(); + GnssroRayPathParameters( + const std::string & rayPathGenType, + double approxRayLengthKm, + double top2DKm, + double resKm, + int nHoriz); + ~GnssroRayPathParameters() {} + GnssroRayPathParameters(const GnssroRayPathParameters & other) = default; + GnssroRayPathParameters & operator=(const GnssroRayPathParameters & other) = default; + + void computeDerivedParameters(); + + const std::string & rayPathGenType() const {return rayPathGenType_;} + double approxRayLength() const {return approxRayLength_;} + double expectedRayLength() const {return expectedRayLength_;} + double top2D() const {return top2D_;} + double horizontalRes() const {return horizontalRes_;} + int nHoriz() const {return nHoriz_;} + + private: + // Member data + std::string rayPathGenType_; // Name for GnssroRayPathGenerator subclass + double approxRayLength_; // Configured ray length (m) + double expectedRayLength_; // ray length (m) resolved from n_horiz and res + double top2D_; // Highest geom height (m) for 2D rays + double horizontalRes_; // Horizontal resolution (m) + int nHoriz_; // Number of nodes in each 2D ray +}; + +} // namespace ufo + +#endif // UFO_OPERATORS_GNSSRO_UTILS_GNSSRORAYPATHPARAMETERS_H_ diff --git a/src/ufo/operators/gnssro/utils/GnssroRayTrajectory.cc b/src/ufo/operators/gnssro/utils/GnssroRayTrajectory.cc new file mode 100644 index 000000000..7a95ce3cd --- /dev/null +++ b/src/ufo/operators/gnssro/utils/GnssroRayTrajectory.cc @@ -0,0 +1,89 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#include +#include +#include +#include "oops/util/Logger.h" +#include "oops/util/missingValues.h" +#include "ufo/operators/gnssro/utils/GnssroRayTrajectory.h" + +namespace ufo { + +// ----------------------------------------------------------------------------- +GnssroRayTrajectory::GnssroRayTrajectory() + : jacT_() + , jacP_() + , jacQ_() + , wf_() + , wi_() + , isSet_(false) +{ +} + +// ----------------------------------------------------------------------------- +GnssroRayTrajectory::GnssroRayTrajectory(const GnssroRayTrajectory& other) + : jacT_(other.jacT_) + , jacP_(other.jacP_) + , jacQ_(other.jacQ_) + , wf_(other.wf_) + , wi_(other.wi_) + , isSet_(other.isSet_) +{ +} + +// ----------------------------------------------------------------------------- +GnssroRayTrajectory & GnssroRayTrajectory::operator=(const GnssroRayTrajectory& other) +{ + if (this != &other) + { + jacT_ = other.jacT_; + jacP_ = other.jacP_; + jacQ_ = other.jacQ_; + wf_ = other.wf_; + wi_ = other.wi_; + isSet_ = other.isSet_; + } + return *this; +} + +// ----------------------------------------------------------------------------- +GnssroRayTrajectory::~GnssroRayTrajectory() +{ +} + +// ----------------------------------------------------------------------------- +void GnssroRayTrajectory::initialize(std::size_t numNodes) +{ + isSet_ = false; + jacT_.assign(numNodes, util::missingValue()); + jacP_.assign(numNodes, util::missingValue()); + jacQ_.assign(numNodes, util::missingValue()); + wf_.assign(numNodes, util::missingValue()); + wi_.assign(numNodes, 0); +} + +// ----------------------------------------------------------------------------- +void GnssroRayTrajectory::finalize() +{ + isSet_ = true; + return; +} +// ----------------------------------------------------------------------------- + +} // namespace ufo + diff --git a/src/ufo/operators/gnssro/utils/GnssroRayTrajectory.h b/src/ufo/operators/gnssro/utils/GnssroRayTrajectory.h new file mode 100644 index 000000000..b82c04308 --- /dev/null +++ b/src/ufo/operators/gnssro/utils/GnssroRayTrajectory.h @@ -0,0 +1,78 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#ifndef UFO_OPERATORS_GNSSRO_UTILS_GNSSRORAYTRAJECTORY_H_ +#define UFO_OPERATORS_GNSSRO_UTILS_GNSSRORAYTRAJECTORY_H_ + +#include +#include // For std::size_t +#include +#include + +namespace ufo { + +// ----------------------------------------------------------------------------- +/// Class to represent represent the trajectory of a single GNSSRO ray, +/// used to implement the Tangent Linear and Adjoint models. +/// Each tangent-point sample in a set of RO observations becomes a ray. +class GnssroRayTrajectory +{ + public: + typedef std::vector DoubleArray_; + typedef std::vector IntArray_; + + GnssroRayTrajectory(); + GnssroRayTrajectory(const GnssroRayTrajectory & other); + GnssroRayTrajectory & operator=(const GnssroRayTrajectory & other); + ~GnssroRayTrajectory(); + + // Accessors + bool isSet() const {return isSet_;} + std::size_t numNodes() {return wf_.size();} + + // Read-only accessors to array data. + double jacobianT(uint32_t nodeIdx) const {return jacT_[nodeIdx];} + double jacobianP(uint32_t nodeIdx) const {return jacP_[nodeIdx];} + double jacobianQ(uint32_t nodeIdx) const {return jacQ_[nodeIdx];} + double wf(uint32_t nodeIdx) const {return wf_[nodeIdx];} + int wi(uint32_t nodeIdx) const {return wi_[nodeIdx];} + + // Read-write accessors to array data. + double & jacobianT(uint32_t nodeIdx) {assert(nodeIdx < jacT_.size()); return jacT_[nodeIdx];} + double & jacobianP(uint32_t nodeIdx) {assert(nodeIdx < jacP_.size()); return jacP_[nodeIdx];} + double & jacobianQ(uint32_t nodeIdx) {assert(nodeIdx < jacQ_.size()); return jacQ_[nodeIdx];} + double & wf(uint32_t nodeIdx) {assert(nodeIdx < wf_.size()); return wf_[nodeIdx];} + int & wi(uint32_t nodeIdx) {assert(nodeIdx < wi_.size()); return wi_[nodeIdx];} + + // Methods to call before and after populating the trajectory. + // Trajectory values should be set using the read-write accessors. + void initialize(std::size_t numNodes); + void finalize(); + + private: + // Member data. + DoubleArray_ jacT_; // Derivative of sim var wrt temperature for each node + DoubleArray_ jacP_; // Derivative of sim var wrt pressure for each node + DoubleArray_ jacQ_; // Derivative of sim var wrt specific humidity for each node + DoubleArray_ wf_; // vertical interpolation weight for each node + IntArray_ wi_; // vertical interpolation base index (1-based) for each node + bool isSet_; // True if trajectory variables that follow are populated. +}; + +} // namespace ufo + +#endif // UFO_OPERATORS_GNSSRO_UTILS_GNSSRORAYTRAJECTORY_H_ diff --git a/src/ufo/operators/gnssro/utils/RefractivityCalculator.cc b/src/ufo/operators/gnssro/utils/RefractivityCalculator.cc new file mode 100644 index 000000000..86c5836bb --- /dev/null +++ b/src/ufo/operators/gnssro/utils/RefractivityCalculator.cc @@ -0,0 +1,151 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#include +#include +#include +#include "eckit/exception/Exceptions.h" +#include "ufo/operators/gnssro/utils/RefractivityCalculator.h" +#include "ufo/utils/Constants.h" + +namespace ufo { + +// ----------------------------------------------------------------------------- +// Specific refractivity algorithms are all declared and implemented in this +// file. The only intended interface if through the RefractivityCalculator +// base class, constructed through the create() factory method. +// ----------------------------------------------------------------------------- +// ----------------------------------------------------------------------------- +// Rueger, 2002 (compressible) or Bevis, 1994 (non-compressible) +// ----------------------------------------------------------------------------- +class RuegerBevisRefractivityCalculator : public RefractivityCalculator +{ + public: + static constexpr const char ALGORITHM_NAME[] = "RuegerBevis"; + explicit RuegerBevisRefractivityCalculator(bool useCompress = true); + ~RuegerBevisRefractivityCalculator() override; + const char * algorithmName() const override; + double N(double T, double Q, double P) const override; + double dNdT(double T, double Q, double P) const override; + double dNdQ(double T, double Q, double P) const override; + double dNdP(double T, double Q, double P) const override; + private: + double n_a_; + double n_b_; + double n_c_; + double oneMinusRdOverRv_; +}; +// Implementation of RuegerBevisRefractivityCalculator + +RuegerBevisRefractivityCalculator::RuegerBevisRefractivityCalculator(bool useCompress) + : RefractivityCalculator() + , n_a_(0.776890) // Default compressible, Rueger 2002 + , n_b_(3.75463e3) + , n_c_(0.712952) + , oneMinusRdOverRv_(1.0 - Constants::rd_over_rv) +{ + // This combination of refractivity algorithm used in cucurull 2010, Healy 2011 + // Default is for compressible atmosphere, per Rueger 2002. + if (!useCompress) { + // Non-compressible atmosphere, per Bevis et al 1994 + n_a_ = 0.7760; + n_b_ = 3.739e3; + n_c_ = 0.704; + } + n_c_ -= n_a_; +} + +RuegerBevisRefractivityCalculator::~RuegerBevisRefractivityCalculator() +{ +} + +const char * RuegerBevisRefractivityCalculator::algorithmName() const +{ + return ALGORITHM_NAME; +} + +double RuegerBevisRefractivityCalculator::N(double T, double Q, double P) const +{ + double tfact = oneMinusRdOverRv_ * Q + Constants::rd_over_rv; + double refr1 = n_a_ * P / T; + double refr2 = n_b_ * Q * P / (T * T * tfact); + double refr3 = n_c_ * Q * P / (T * tfact); + double refr = refr1 + refr2 + refr3; + return refr; +} + +double RuegerBevisRefractivityCalculator::dNdT(double T, double Q, double P) const +{ + // Algorithm taken from Jacobian calculation in ufo_gnssro_refncep_tlad_mod.F90 + // self%jac_t(iobs) = - n_a*gesP/gesT**2 & + // - n_b*two*gesP*gesQ/( ((1-rd_over_rv)*gesQ+rd_over_rv)*gesT**3 ) & + // - n_c*gesP*gesQ/( ((1-rd_over_rv)*gesQ+rd_over_rv)*gesT**2 ) + double T_sq = T * T; + double PxQ = P * Q; + double tfact = oneMinusRdOverRv_ * Q + Constants::rd_over_rv; + double jacT = - n_a_ * P / T_sq + - n_b_ * 2.0 * PxQ / (tfact * T_sq * T) + - n_c_ * PxQ / (tfact * T_sq); + return jacT; +} + +double RuegerBevisRefractivityCalculator::dNdQ(double T, double Q, double P) const +{ + // Algorithm taken from Jacobian calculation in ufo_gnssro_refncep_tlad_mod.F90 + // self%jac_q(iobs) = n_b*gesP/( gesT**2*( (1-rd_over_rv)*gesQ+rd_over_rv)**2 ) & + // * rd_over_rv & + // + n_c*gesP/( gesT *( (1-rd_over_rv)*gesQ+rd_over_rv)**2 ) & + // * rd_over_rv + double T_sq = T * T; + double tfact = oneMinusRdOverRv_ * Q + Constants::rd_over_rv; + double tfact_sq = tfact * tfact; + double jacQ = n_b_ * P / (T_sq * tfact_sq) * Constants::rd_over_rv + + n_c_ * P / (T * tfact_sq) * Constants::rd_over_rv; + return jacQ; +} + +double RuegerBevisRefractivityCalculator::dNdP(double T, double Q, double P) const +{ + // Algorithm taken from Jacobian calculation in ufo_gnssro_refncep_tlad_mod.F90 + // self%jac_prs(iobs) = n_a/gesT & + // + n_b*gesQ/ ( ((1-rd_over_rv)*gesQ+rd_over_rv)*gesT**2 ) & + // + n_c*gesQ/ ( ((1-rd_over_rv)*gesQ+rd_over_rv)*gesT ) + double tfact = oneMinusRdOverRv_ * Q + Constants::rd_over_rv; + double jacP = n_a_ / T + + n_b_ * Q / (tfact * T * T) + + n_c_ * Q / (tfact * T); + return jacP; +} + +// ----------------------------------------------------------------------------- +// Base class implementation. +// ----------------------------------------------------------------------------- +// Factory method. +std::unique_ptr RefractivityCalculator::create( + const std::string& algorithmName, bool useCompress) +{ + if (algorithmName == RuegerBevisRefractivityCalculator::ALGORITHM_NAME) { + return std::make_unique(useCompress); + } else { + throw eckit::BadValue( + "Unsupported algorithm name for RefractivityCalculator: " + algorithmName, + Here()); + } +} + +} // namespace ufo + diff --git a/src/ufo/operators/gnssro/utils/RefractivityCalculator.h b/src/ufo/operators/gnssro/utils/RefractivityCalculator.h new file mode 100644 index 000000000..1beaf4bcd --- /dev/null +++ b/src/ufo/operators/gnssro/utils/RefractivityCalculator.h @@ -0,0 +1,69 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#ifndef UFO_OPERATORS_GNSSRO_UTILS_REFRACTIVITYCALCULATOR_H_ +#define UFO_OPERATORS_GNSSRO_UTILS_REFRACTIVITYCALCULATOR_H_ + +#include +#include + +namespace ufo { + +// ----------------------------------------------------------------------------- +/// Base class for RefractivityCalculator. +/// +/// Refractivity algorithms are identified by a textual name. Each algorithm +/// is implemented as a separate subclass of RefractivityCalculator. +/// In addition, the same algorithm can be identified as using a compressible +/// or non-compressible atmosphere using the useCompress argument. +/// +/// The RefractivityCalculator can compute refractivity from air temperature +/// (T) in Kelvin, specific humidity (Q) in kg/kg and pressure (P) in Pa. +/// It can also compute the derivatives of refractivity with respect to each +/// of these variables. +/// +class RefractivityCalculator { + public: + // Refractivity units support: + // Conversions between excess index of refraction (n-1) to Refractivity N. + static constexpr double EXCESS_IOR_TO_N = 1e6; + static constexpr double N_TO_EXCESS_IOR = 1e-6; + + // Factory method. + static std::unique_ptr create( + const std::string& algorithmName, bool useCompress); + + // Constructor/Destructor + RefractivityCalculator() = default; + virtual ~RefractivityCalculator() = default; + virtual const char * algorithmName() const = 0; + + /// \brief Return refractivity as a function of temperature, specific + /// humidity, and pressure. + virtual double N(double T, double Q, double P) const = 0; + + /// \brief Return derivative of refractivity with respect to temperature + virtual double dNdT(double T, double Q, double P) const = 0; + /// \brief Return derivative of refractivity with respect to specific humidity + virtual double dNdQ(double T, double Q, double P) const = 0; + /// \brief Return derivative of refractivity with respect to pressure + virtual double dNdP(double T, double Q, double P) const = 0; +}; + +} // namespace ufo + +#endif // UFO_OPERATORS_GNSSRO_UTILS_REFRACTIVITYCALCULATOR_H_ diff --git a/src/ufo/operators/gnssro/utils/StraightLineRayPathGenerator.cc b/src/ufo/operators/gnssro/utils/StraightLineRayPathGenerator.cc new file mode 100644 index 000000000..01ea41576 --- /dev/null +++ b/src/ufo/operators/gnssro/utils/StraightLineRayPathGenerator.cc @@ -0,0 +1,302 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#include +#include +#include // for std::size_t +#include +#include +#include +#include "oops/util/Logger.h" +#include "oops/util/missingValues.h" +#include "ufo/operators/gnssro/utils/GnssroRayPath.h" +#include "ufo/operators/gnssro/utils/StraightLineRayPathGenerator.h" + +namespace ufo { + +static const double RAD2DEG = 180.0 / M_PI; +static const double DEG2RAD = M_PI / 180.0; + +template +static T normalizeLongitude_180(T lon) +{ + // Normalize longitude in degrees to -180 to 180 range. + T angle = std::fmod(lon + static_cast(180), static_cast(360)); + if (angle < static_cast(0)) + angle += static_cast(360); + angle -= static_cast(180); + return angle; +} + +template +static T normalizeLongitude_360(T lon) +{ + static const T three_sixty = static_cast(360); + return std::fmod(std::fmod(lon, three_sixty) + three_sixty, three_sixty); +} + +/* + * Geometry for straightline calculations. + * d = rayDist + * ----- + * \ | + * edgeRadius = R \ | R0 = radius0 + * \ | + * \| + * Angle between R and R0 is baseTheta. + */ +StraightLineRayPathGenerator::NodeCalculator::NodeCalculator( + double tpGeomHeight, + double localEarthRadius, + float azimuth, + float tpLat, + float tpLon, + const GnssroRayPathParameters & params) + : tpGeomHeight_(tpGeomHeight) + , localEarthRadius_(localEarthRadius) + , radius0_(tpGeomHeight_ + localEarthRadius_) + , radius0_sq_(radius0_ * radius0_) + , azimuth_((azimuth > 180.0) ? azimuth - 360.0 : azimuth) + , tpLat_(tpLat) + , tpLon_(tpLon) + , params_(params) + , azimuthRad_(azimuth_ * DEG2RAD) + , tpLatRad_(tpLat_ * DEG2RAD) + , tpLonRad_(tpLon_ * DEG2RAD) + , cosAzimuth_(std::cos(azimuthRad_)) + , signAzimuth_((azimuth_ >= 0.0) ? 1 : -1) + , cosTpLat_(std::cos(tpLatRad_)) + , cosTpLon_(std::cos(tpLonRad_)) + , sinTpLat_(std::sin(tpLatRad_)) + , sinTpLon_(std::sin(tpLonRad_)) + , edgeIncrement_(0.5 * params_.horizontalRes()) // increment for rayDist + , edgeDist_(0.0) // Length of half ray (from TP) + , edgeRadius_(radius0_) // Radius to current edge + , baseTheta_(0.0) // Angle to current node (not edge) +{ +} + +StraightLineRayPathGenerator::NodeCalculator::~NodeCalculator() +{ +} + +// Compute increment ray distance and recompute quantities that depend on it. +void StraightLineRayPathGenerator::NodeCalculator::extendRay() +{ + // Compute the distance to the next edge between nodes. + edgeDist_ += edgeIncrement_; + + // After the first ray distance computation, we increment by full resolution. + edgeIncrement_ = params_.horizontalRes(); + + // Compute the edge radius, used to determine geometric height of an edge + edgeRadius_ = std::sqrt(edgeDist_ * edgeDist_ + radius0_sq_); + + // Compute the angle to the current node, which is half of a segment length + // closer to the tangent point than the current edge. + // This angle is 0 at the tangent point. + double nodeDist = edgeDist_ - 0.5 * edgeIncrement_; + baseTheta_ = std::atan2(nodeDist, radius0_); + return; +} + +std::pair StraightLineRayPathGenerator::NodeCalculator::lonLatAlongRay(double theta) +{ + // This algorithm is a copy of the Fortran logic from ROPP, except we use a + // variable theta, rather than a constant value based on distances at the equator. + double sinTheta = std::sin(theta); + double cosTheta = std::cos(theta); + + double sinLat = sinTheta * cosAzimuth_ * cosTpLat_ + sinTpLat_ * cosTheta; + double latRad = std::asin(sinLat); + + double cosDelta = 1.0; // For North/South poles + if (std::abs(latRad) != M_PI / 2.0) // Not North/South pole + { + cosDelta = (cosTheta - (sinTpLat_ * sinLat)) / (cosTpLat_ * std::cos(latRad)); + if (cosDelta > 1.0) + cosDelta = 1.0; + if (cosDelta < -1.0) + cosDelta = -1.0; + } + + int signTheta = (theta >= 0.0) ? 1 : -1; + double lonRad; + if (signTheta == signAzimuth_) { + lonRad = tpLonRad_ + std::acos(cosDelta); + } else { + lonRad = tpLonRad_ - std::acos(cosDelta); + } + + // Convert lat, lon to degrees; normalize lon to -pi to pi range. + double latDeg = latRad * RAD2DEG; + double lonDeg = lonRad * RAD2DEG; + lonDeg = normalizeLongitude_180(lonDeg); + return std::pair(lonDeg, latDeg); +} + +////////////////////////////////////////////////// +// Implementation of StraightLineRayPathGenerator +////////////////////////////////////////////////// +StraightLineRayPathGenerator::StraightLineRayPathGenerator( + const ioda::ObsSpace & odb, + const GnssroRayPathParameters & params) + : GnssroRayPathGenerator(odb, params) + , obsLat_(odb.nlocs()) + , obsLon_(odb.nlocs()) + , obsAlt_(odb.nlocs()) + , obsRoc_(odb.nlocs()) + , obsUndul_(odb.nlocs()) + , obsAzimuth_(odb.nlocs()) +{ + oops::Log::trace() << name() << " entering constructor" << std::endl; + + // Get observation data. This may contain multiple profiles. + odb.get_db("MetaData", "latitude", obsLat_); + odb.get_db("MetaData", "longitude", obsLon_); + odb.get_db("MetaData", "height", obsAlt_); + odb.get_db("MetaData", "earthRadiusCurvature", obsRoc_); + odb.get_db("MetaData", "geoidUndulation", obsUndul_); + odb.get_db("MetaData", "sensorAzimuthAngle", obsAzimuth_); + + oops::Log::trace() << name() << " exiting constructor" << std::endl; +} + +StraightLineRayPathGenerator::~StraightLineRayPathGenerator() +{ +} + +std::unique_ptr +StraightLineRayPathGenerator::makeProfileRayPaths(const GnssroProfileSlice & profileSlice) +{ + oops::Log::trace() << name() << "::makeProfileRayPaths - enter" << std::endl; + auto profileRayPaths = std::make_unique( + profileSlice, obsLat_, obsLon_, obsAzimuth_); + + oops::Log::trace() << name() << "::makeProfileRayPaths - created empty profileRayPaths" + << std::endl; + std::size_t lastObBelowTop = profileSlice.end() - 1; + std::size_t obIdx = profileSlice.start(); + for (; obIdx <= lastObBelowTop; ++obIdx) { + // Number of nodes will be constant for all levels. + // However, lat-lon for all nodes of a ray will be set + // to the tangent point location for rays above top2D. + double alt = obsAlt_[obIdx]; + bool assumeSphericalSymmetry = (alt > params_.top2D()); + size_t numNodes = params_.nHoriz(); + oops::Log::debug() << name() << "seqNum=" << profileSlice.seqNum() << ", obIdx=" << obIdx + << ": checking geometric height=" << alt << " against top2D=" << params_.top2D() + << ", numNodes=" << numNodes + << ", assumeSphericalSymmetry=" << assumeSphericalSymmetry + << std::endl; + + double localEarthRadius = obsRoc_[obIdx] + obsUndul_[obIdx]; + NodeCalculator nodeCalculator(alt, + localEarthRadius, + obsAzimuth_[obIdx], + obsLat_[obIdx], + obsLon_[obIdx], + params_); + oops::Log::debug() << name() << ": obIdx " << obIdx << ": created NodeCalculator" + << " with tpLat/Lon " << nodeCalculator.tpLat() << "," << nodeCalculator.tpLon() + << ", azimuth " << nodeCalculator.azimuth() << std::endl; + + auto rayPath = std::make_unique(numNodes, + nodeCalculator.tpLat(), + nodeCalculator.tpLon(), + nodeCalculator.azimuth(), + alt, + assumeSphericalSymmetry); + + // Set the node index for the tangent point. + size_t centerNodeIdx = rayPath->centerNodeIdx(); + rayPath->setTpNodeIdx(centerNodeIdx); + + oops::Log::debug() << name() << ": obIdx " << obIdx << ": created GnssroRayPath with " + << rayPath->numNodes() << " nodes, edgesInHalfRay=" << rayPath->edgesInHalfRay() + << std::endl; + + for (std::size_t edgeOffset = 0; edgeOffset < rayPath->edgesInHalfRay(); ++edgeOffset) + { + nodeCalculator.extendRay(); + double edgeGeomHeight = nodeCalculator.edgeGeomHeight(); + + // Set edge-oriented geometric heights, which are symmetric + // around the tangent point. + std::size_t rightEdgeIdx = centerNodeIdx + edgeOffset + 1; + std::size_t leftEdgeIdx = centerNodeIdx - edgeOffset; + rayPath->geomHgt(rightEdgeIdx) = edgeGeomHeight; + rayPath->geomHgt(leftEdgeIdx) = edgeGeomHeight; + + // Set node-oriented values. + if (edgeOffset == 0) { + // Tangent point + assert(rayPath->numNodes() > centerNodeIdx); + rayPath->segLen(centerNodeIdx) = params_.horizontalRes(); + if (!assumeSphericalSymmetry) { + rayPath->lat(centerNodeIdx) = nodeCalculator.tpLat(); + rayPath->lon(centerNodeIdx) = nodeCalculator.tpLon(); + } + } else { + // Off-tangent points to right and left of tangent point + std::size_t rightNodeIdx = centerNodeIdx + edgeOffset; + std::size_t leftNodeIdx = centerNodeIdx - edgeOffset; + assert(rayPath->numNodes() > rightNodeIdx); + assert(rayPath->numNodes() > leftNodeIdx); + rayPath->segLen(rightNodeIdx) = params_.horizontalRes(); + rayPath->segLen(leftNodeIdx) = params_.horizontalRes(); + + // Set lat-lons if we are not assuming spherical symmetry. + if (!assumeSphericalSymmetry) { + std::pair rightLonLat = nodeCalculator.lonLatAlongRay( + nodeCalculator.rightTheta()); + std::pair leftLonLat = nodeCalculator.lonLatAlongRay( + nodeCalculator.leftTheta()); + if (rightLonLat.first == 0.0 && rightLonLat.second == 0.0) { + oops::Log::warning() << name() << "seqNum=" << profileSlice.seqNum() + << ", obIdx=" << obIdx << ", edgeIdx=" << rightEdgeIdx << " of " + << rayPath->numEdges() << ", edgeOffset=" << edgeOffset + << ": rightLonLat values are both 0.0" << std::endl; + } + if (leftLonLat.first == 0.0 && leftLonLat.second == 0.0) { + oops::Log::warning() << name() << "seqNum=" << profileSlice.seqNum() + << ", obIdx=" << obIdx << ", edgeIdx=" << leftEdgeIdx << " of " + << rayPath->numEdges() << ", edgeOffset=" << edgeOffset + << ": leftLonLat values are both 0.0 " << std::endl; + } + + rayPath->lon(rightNodeIdx) = rightLonLat.first; + rayPath->lat(rightNodeIdx) = rightLonLat.second; + rayPath->lon(leftNodeIdx) = leftLonLat.first; + rayPath->lat(leftNodeIdx) = leftLonLat.second; + } // end of setting lat-lons + } + } // end loop over edges in the ray path. + + rayPath->updateTotalLength(); + oops::Log::debug() << name() << ": obIdx " << obIdx << ": loaded GnssroRayPath with " + << "numNodes=" << rayPath->numNodes() << ", totalLength=" << rayPath->totalLength() + << std::endl; + profileRayPaths->addRay(std::move(rayPath)); + } // end of loop over observation samples in the profile (except the topmost obs). + + oops::Log::trace() << name() << "::makeProfileRayPaths - exit " + << profileRayPaths->numRays() << " ray paths created" << std::endl; + return profileRayPaths; +} + +} // namespace ufo diff --git a/src/ufo/operators/gnssro/utils/StraightLineRayPathGenerator.h b/src/ufo/operators/gnssro/utils/StraightLineRayPathGenerator.h new file mode 100644 index 000000000..6cfddf8c9 --- /dev/null +++ b/src/ufo/operators/gnssro/utils/StraightLineRayPathGenerator.h @@ -0,0 +1,113 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#ifndef UFO_OPERATORS_GNSSRO_UTILS_STRAIGHTLINERAYPATHGENERATOR_H_ +#define UFO_OPERATORS_GNSSRO_UTILS_STRAIGHTLINERAYPATHGENERATOR_H_ + +#include +#include +#include +#include +#include "ufo/operators/gnssro/utils/GnssroRayPathGenerator.h" +#include "ufo/operators/gnssro/utils/GnssroRayPathParameters.h" + +namespace ufo { + +// ----------------------------------------------------------------------------- +/// Subclass of GnssroRayPathGenerator for producing straight-line ray paths. +/// +// +class StraightLineRayPathGenerator: public GnssroRayPathGenerator { + public: + class NodeCalculator { + public: + NodeCalculator(double tpGeomHeight, + double localEarthRadius, + float azimuth, + float tpLat, + float tpLon, + const GnssroRayPathParameters & params); + ~NodeCalculator(); + + // Accessors + double tpLat() const {return tpLat_;} + double tpLon() const {return tpLon_;} + double azimuth() const {return azimuth_;} + double edgeDist() const {return edgeDist_;} + double edgeRadius() const {return edgeRadius_;} + double edgeGeomHeight() const {return edgeRadius_ - localEarthRadius_;} + double rightTheta() const {return baseTheta_;} + double leftTheta() const {return -1.0 * baseTheta_;} + + /// \brief Method to extend ray properties by one more node in + // each direction from the tangent point. + void extendRay(); + + /// Returns lon, lat in degrees that is arc-length theta from the + // tangent point along the ray path defined by the azimuth angle. + std::pair lonLatAlongRay(double theta); + + private: + double tpGeomHeight_; + double localEarthRadius_; + double radius0_; + double radius0_sq_; + double azimuth_; + double tpLat_; + double tpLon_; + GnssroRayPathParameters params_; + double azimuthRad_; + double tpLatRad_; + double tpLonRad_; + double cosAzimuth_; + int signAzimuth_; + double cosTpLat_; + double cosTpLon_; + double sinTpLat_; + double sinTpLon_; + double edgeIncrement_; + double edgeDist_; + double edgeRadius_; + double baseTheta_; + }; // end NodeCalculator nested class + + static constexpr const char GENTYPE[] = "StraightLine"; + static const char * name() {return GENTYPE;} + + // Constructor/Destructor + explicit StraightLineRayPathGenerator( + const ioda::ObsSpace & odb, + const GnssroRayPathParameters & params); + ~StraightLineRayPathGenerator() override; + + /// \brief Return a GnssroRayPath object based on observation data read for + /// the specified GNSSRO profile from the ObsSpace. + std::unique_ptr makeProfileRayPaths( + const GnssroProfileSlice & profileSlice) override; + + private: + std::vector obsLat_; // observed latitude (degrees) + std::vector obsLon_; // observed longitude (degrees) + std::vector obsAlt_; // observed geometric height (m) + std::vector obsRoc_; // observed earth radius of curvature (m) + std::vector obsUndul_; // observed geoid undulation (m) + std::vector obsAzimuth_; // observed azimuth angle/bearing (deg clockwise from north) +}; + +} // namespace ufo + +#endif // UFO_OPERATORS_GNSSRO_UTILS_STRAIGHTLINERAYPATHGENERATOR_H_ diff --git a/src/ufo/ufo_geovals_mod.F90 b/src/ufo/ufo_geovals_mod.F90 index 300f7ef84..2b500f57a 100644 --- a/src/ufo/ufo_geovals_mod.F90 +++ b/src/ufo/ufo_geovals_mod.F90 @@ -174,6 +174,7 @@ end subroutine ufo_geovals_update_linit !> @param npaths_by_method !> An array of length `nsampling_methods` whose ith element is the number of paths in the set !> of interpolation paths obtained with the ith sampling method. +!if (self%nsampling_methods .gt. 0) then !> @param sampling_method_by_var !> An array of length `nvars` whose ith element is the index of the observation location !> sampling method producing the set of paths along which the ith variable will be @@ -1675,7 +1676,7 @@ subroutine ufo_geovals_read_netcdf(self, filename, loc_multiplier, levels_are_to type(c_ptr), intent(in) :: c_obspace type(oops_variables), intent(in) :: vars -integer :: global_npaths, var_global_npaths +integer :: global_npaths, var_global_npaths, global_nlocs integer :: nval integer :: obs_nlocs integer :: obs_all_nlocs @@ -1683,7 +1684,8 @@ subroutine ufo_geovals_read_netcdf(self, filename, loc_multiplier, levels_are_to character(len=MAXVARLEN), pointer :: varname type(oops_variables) :: reduced_vars -integer :: ncid, dimid, varid, vartype, ndims +logical :: use_nprofiles = .false. +integer :: ncid, dimid, varid, vartype, ndims, dimid_nlocs integer, dimension(3) :: dimids integer :: ivar integer :: ierr @@ -1695,27 +1697,43 @@ subroutine ufo_geovals_read_netcdf(self, filename, loc_multiplier, levels_are_to integer(c_size_t), allocatable, dimension(:) :: global_loc_by_my_loc integer(c_size_t), allocatable, dimension(:) :: global_path_by_my_path type(ufo_index_range), allocatable :: path_ranges_by_loc(:) +integer, allocatable, dimension(:) :: pblb_indices +integer, allocatable, dimension(:) :: pble_indices real, allocatable :: field2d(:,:), field1d(:) -! At the time being, the NetCDF GeoVaLs file format requires all GeoVaLs to be interpolated +! For the most part, the NetCDF GeoVaLs file format requires all GeoVaLs to be interpolated ! along the same set of paths composed of as many paths as there are observation locations. ! The format will be made more general in future. +! There are two exception to this limitation: +! * when loc_multiplier is used to specify multiple geoval profiles per observation. +! * when the nprofiles dimension exists in the geovals file to specify the +! dimension of all the geovals variables, distinct from nlocs. integer, parameter :: nsampling_methods = 1 integer(c_size_t), allocatable :: sampling_method_by_var(:) logical(c_bool) :: is_sampling_method_trivial(1) ! open netcdf file call check('nf90_open', nf90_open(trim(filename),nf90_nowrite,ncid)) - -! find how many locs are in the file -ierr = nf90_inq_dimid(ncid, "nlocs", dimid) +ierr = nf90_inq_dimid(ncid, "nlocs", dimid_nlocs) if(ierr /= nf90_noerr) then write(err_msg,*) "Error: Dimension nlocs not found in ", trim(filename) call abor1_ftn(err_msg) -endif +end if -call check('nf90_inq_dimid', nf90_inq_dimid(ncid, "nlocs", dimid)) -call check('nf90_inquire_dimension', nf90_inquire_dimension(ncid, dimid, len = global_npaths)) +call check('nf90_inq_dimid', nf90_inq_dimid(ncid, "nlocs", dimid_nlocs)) +call check('nf90_inquire_dimension', nf90_inquire_dimension(ncid, dimid_nlocs, len = global_nlocs)) + +! Determine if the geovals uses the nprofiles or nlocs dimension for variables. +ierr = nf90_inq_dimid(ncid, "nprofiles", dimid) +if(ierr == nf90_noerr) then + ! The geovals file has an nprofiles dimension, distinct from nlocs + use_nprofiles = .true. + call check('nf90_inq_dimid', nf90_inq_dimid(ncid, "nprofiles", dimid)) + call check('nf90_inquire_dimension', nf90_inquire_dimension(ncid, dimid, len = global_npaths)) +else + dimid = dimid_nlocs + global_npaths = global_nlocs +end if !> round-robin distribute the observations to PEs !> Calculate how many obs. on each PE @@ -1728,48 +1746,86 @@ subroutine ufo_geovals_read_netcdf(self, filename, loc_multiplier, levels_are_to ! single location in the obs file. There needs to be at least ! loc_multiplier * obs_all_nlocs locations in the geovals file. -if (global_npaths .lt. (loc_multiplier * obs_all_nlocs)) then - write(obs_nlocs_str, *) loc_multiplier * obs_all_nlocs - write(geo_nlocs_str, *) global_npaths - write(err_msg,'(7a)') & - "Error: Number of locations in the geovals file (", & - trim(adjustl(geo_nlocs_str)), ") must be greater than or equal to ", & - "the product of loc_multiplier and number of locations in the ", & - "obs file (", trim(adjustl(obs_nlocs_str)), ")" - call abor1_ftn(err_msg) +if (.not. use_nprofiles) then + if (global_npaths .lt. (loc_multiplier * obs_all_nlocs)) then + write(obs_nlocs_str, *) loc_multiplier * obs_all_nlocs + write(geo_nlocs_str, *) global_npaths + write(err_msg,'(7a)') & + "Error: Number of locations in the geovals file (", & + trim(adjustl(geo_nlocs_str)), ") must be greater than or equal to ", & + "the product of loc_multiplier and number of locations in the ", & + "obs file (", trim(adjustl(obs_nlocs_str)), ")" + call abor1_ftn(err_msg) + endif endif -! We have enough locations in the geovals file to cover the span of the -! number of locations in the obs file. Generate the global_path_by_my_path map according -! to the loc_multiplier and obs_nlocs values. allocate(path_ranges_by_loc(obs_nlocs)) -if (loc_multiplier > 0) then - my_npaths = loc_multiplier * obs_nlocs +if (use_nprofiles) then + ! Read path_ranges_by_loc from the geovals file. The values in the file + ! should already be global values, so the adjustment from local to global + ! is only needed for my_loc, not for the begin and end offsets. + my_npaths = global_npaths + allocate(pblb_indices(obs_all_nlocs)) + allocate(pble_indices(obs_all_nlocs)) allocate(global_path_by_my_path(my_npaths)) - my_path = 1 + ierr = nf90_inq_varid(ncid, "paths_by_loc_begin", varid) + if(ierr /= nf90_noerr) then + write(err_msg,*) "Error: Variable paths_by_loc_begin not found in ", trim(filename) + call abor1_ftn(err_msg) + endif + call check('nf90_get_var_pblb', nf90_get_var(ncid, varid, pblb_indices)) + + ierr = nf90_inq_varid(ncid, "paths_by_loc_end", varid) + if(ierr /= nf90_noerr) then + write(err_msg,*) "Error: Variable paths_by_loc_end not found in ", trim(filename) + call abor1_ftn(err_msg) + endif + call check('nf90_get_var_pble', nf90_get_var(ncid, varid, pble_indices)) + do my_loc = 1,obs_nlocs - global_path_start = ((global_loc_by_my_loc(my_loc) - 1) * loc_multiplier) + 1 - global_path_end = global_loc_by_my_loc(my_loc) * loc_multiplier - path_ranges_by_loc(my_loc)%begin = my_path - do global_path = global_path_start, global_path_end - global_path_by_my_path(my_path) = global_path - my_path = my_path + 1 - enddo - path_ranges_by_loc(my_loc)%end = my_path + ! global_loc = global_loc_by_my_loc(my_loc) + !write(*,*) "map my_loc=", my_loc, " to global_loc=", global_loc + path_ranges_by_loc(my_loc)%begin = pblb_indices(my_loc) + path_ranges_by_loc(my_loc)%end = pble_indices(my_loc) enddo -else - ! Negative loc_multipliers are no longer supported. They used to map each location to a set - ! of paths with non-consecutive indices, which cannot be represented with the current data - ! structures. - write(err_msg, *) "Error: loc_multiplier must be positive" - call abor1_ftn(err_msg) -end if + do my_path = 1, my_npaths + global_path_by_my_path(my_path) = my_path + enddo + + deallocate(pblb_indices) + deallocate(pble_indices) +else + ! We have enough locations in the geovals file to cover the span of the + ! number of locations in the obs file. Generate the global_path_by_my_path map according + ! to the loc_multiplier and obs_nlocs values. + if (loc_multiplier > 0) then + my_npaths = loc_multiplier * obs_nlocs + allocate(global_path_by_my_path(my_npaths)) + my_path = 1 + do my_loc = 1,obs_nlocs + global_path_start = ((global_loc_by_my_loc(my_loc) - 1) * loc_multiplier) + 1 + global_path_end = global_loc_by_my_loc(my_loc) * loc_multiplier + path_ranges_by_loc(my_loc)%begin = my_path + do global_path = global_path_start, global_path_end + global_path_by_my_path(my_path) = global_path + my_path = my_path + 1 + enddo + path_ranges_by_loc(my_loc)%end = my_path + enddo + else + ! Negative loc_multipliers are no longer supported. They used to map each location to a set + ! of paths with non-consecutive indices, which cannot be represented with the current data + ! structures. + write(err_msg, *) "Error: loc_multiplier must be positive" + call abor1_ftn(err_msg) + endif +endif ! define a common sampling method for all variables allocate(sampling_method_by_var(vars%nvars())) sampling_method_by_var(:) = 1 -is_sampling_method_trivial(1) = (loc_multiplier == 1) +is_sampling_method_trivial(1) = (loc_multiplier == 1 .and. .not. use_nprofiles) ! Currently NetCDF files store variables only in the sampled format. But if loc_multiplier is ! equal to 1, there is no difference between the sampled and reduced format, so we can set up a @@ -1777,7 +1833,7 @@ subroutine ufo_geovals_read_netcdf(self, filename, loc_multiplier, levels_are_to ! VaLs object storing variables only in the sampled format. ! ! In future we may want to extend the NetCDF file layout to accommodate both formats at once. -if (loc_multiplier == 1) then +if (loc_multiplier == 1 .and. .not. use_nprofiles) then reduced_vars = oops_variables(vars) else reduced_vars = oops_variables() @@ -1877,6 +1933,7 @@ subroutine ufo_geovals_read_netcdf(self, filename, loc_multiplier, levels_are_to end subroutine ufo_geovals_read_netcdf + ! ------------------------------------------------------------------------------ subroutine ufo_geovals_write_netcdf(self, filename) use netcdf @@ -1884,45 +1941,140 @@ subroutine ufo_geovals_write_netcdf(self, filename) type(ufo_geovals), intent(inout) :: self character(max_string), intent(in) :: filename -integer :: i -integer :: ncid, dimid_nlocs, dimid_nval, dims(2) +integer :: i, nlocs +integer :: ncid, dimid_nlocs, dimid_nval, dimid_nprofile, dims(2) integer, allocatable :: ncid_var(:) +integer :: ncid_pbl_begin, ncid_pbl_end +integer, allocatable :: pbl_indices(:) +integer :: nprofiles +logical :: nprofiles_and_nlocs_differ +nlocs = self%nlocs allocate(ncid_var(self%nvar)) call check('nf90_create', nf90_create(trim(filename),nf90_hdf5,ncid)) -! TODO(wsmigaj): define a new format with nlocs replaced by nprofiles and stored -! separately for each variable -!------------------------------------------------------------------------------ -!call check('nf90_def_dim', nf90_def_dim(ncid,'nlocs',self%nlocs, dimid_nlocs)) ! This modification accounts for two scenarios: ! the number of observations (nlocs) for the trivial sampling method, where self%geovals(1)%nprofiles -! is equal to nlocs, thus not affecting previous cases; and the number of extended geovals required for non-trivial cases. -! Further refinement is needed for storing each variable separately. -!------------------------------------------------------------------------------ +! is equal to nlocs, thus not affecting previous cases; and the number of extended geovals +! required for non-trivial cases. Further refinement is needed for storing each variable with its +! own separate dimensions. call check('nf90_def_dim', nf90_def_dim(ncid,'nlocs',self%geovals(1)%nprofiles, dimid_nlocs)) dims(2) = dimid_nlocs - +nprofiles = -1 +nprofiles_and_nlocs_differ = .false. do i = 1, self%nvar + ! Check if number of profiles and observed locs are the same + ! If not, we will create an nprofiles dimension for use with the geoval vars. + ! This code requires the nprofile dim of all variables to be the same. + if (self%geovals(i)%nprofiles /= self%nlocs) then + if (nprofiles == -1) then + nprofiles = self%geovals(i)%nprofiles + nprofiles_and_nlocs_differ = .true. + ! create a profile dimension for this var to use in place of nlocs. + call check('nf90_def_dim', & + nf90_def_dim(ncid,"nprofiles",self%geovals(i)%nprofiles, dimid_nprofile)) + else if (self%geovals(i)%nprofiles /= nprofiles) then + call abor1_ftn("ufo_geovals_write_netcdf: all variables must share the same"& + // " nprofiles dimension") + end if + dims(2) = dimid_nprofile + else + dims(2) = dimid_nlocs + end if + + ! Create a separate nval dimension for each variable call check('nf90_def_dim', & nf90_def_dim(ncid,trim(self%variables(i))//"_nval",self%geovals(i)%nval, dimid_nval)) dims(1) = dimid_nval - call check('nf90_def_var', & + + call check('nf90_def_var', & nf90_def_var(ncid,trim(self%variables(i)),nf90_float,dims,ncid_var(i))) enddo +! Define variables for paths_by_loc mapping if the relationship +! between paths and observed locations is not 1-to-1. +if (nprofiles_and_nlocs_differ) then + ! For now, we just define the paths_by_loc for the first sampling method. If + ! the geovals use sampled and reduced format, there could be more than one. + call check('nf90_def_pblb_var', & + nf90_def_var(ncid,"paths_by_loc_begin",nf90_int,dimid_nlocs,ncid_pbl_begin)) + call check('nf90_def_pble_var', & + nf90_def_var(ncid,"paths_by_loc_end",nf90_int,dimid_nlocs,ncid_pbl_end)) +end if + call check('nf90_enddef', nf90_enddef(ncid)) do i = 1, self%nvar call check('nf90_put_var', nf90_put_var(ncid,ncid_var(i),self%geovals(i)%vals(:,:))) enddo +! Write the paths_by_loc mapping if the relationship between +! paths and observed locations is not 1-to-1. +if (nprofiles_and_nlocs_differ) then + ! For now, we write just the paths_by_loc for the first sampling method. If + ! the geovals use sampled and reduced format, there could be more than one. + ! TODO(smarshall) - if the geovals are written from multiple PEs, the + ! pbl_indices will be relative to only the paths and locs in that PE. + ! Ultimately, these indices should be global to all locations across all PEs. + if (self%nsampling_methods .eq. 0) then + call abor1_ftn("ufo_geovals_write_netcdf: num sampling method must be " & + // "> 0 when nprofiles and nlocs differ") + end if + allocate(pbl_indices(self%nlocs)) + do i = 1, self%nlocs + pbl_indices(i) = self%sampling_methods(1)%paths_by_loc(i)%begin + enddo + call check('nf90_put_var', nf90_put_var(ncid,ncid_pbl_begin,pbl_indices(:))) + do i = 1, self%nlocs + pbl_indices(i) = self%sampling_methods(1)%paths_by_loc(i)%end + enddo + call check('nf90_put_var', nf90_put_var(ncid,ncid_pbl_end,pbl_indices(:))) + deallocate(pbl_indices) +end if + call check('nf90_close', nf90_close(ncid)) deallocate(ncid_var) end subroutine ufo_geovals_write_netcdf +!subroutine ufo_geovals_write_netcdf(self, filename) +!use netcdf +!implicit none +!type(ufo_geovals), intent(inout) :: self +!character(max_string), intent(in) :: filename +! +!integer :: i +!integer :: ncid, dimid_nlocs, dimid_nval, dims(2) +!integer, allocatable :: ncid_var(:) +! +!allocate(ncid_var(self%nvar)) +! +!call check('nf90_create', nf90_create(trim(filename),nf90_hdf5,ncid)) +!! TODO(wsmigaj): define a new format with nlocs replaced by nprofiles and stored +!! separately for each variable +!call check('nf90_def_dim', nf90_def_dim(ncid,'nlocs',self%nlocs, dimid_nlocs)) +!dims(2) = dimid_nlocs +! +!do i = 1, self%nvar +! call check('nf90_def_dim', & +! nf90_def_dim(ncid,trim(self%variables(i))//"_nval",self%geovals(i)%nval, dimid_nval)) +! dims(1) = dimid_nval +! call check('nf90_def_var', & +! nf90_def_var(ncid,trim(self%variables(i)),nf90_float,dims,ncid_var(i))) +!enddo +! +!call check('nf90_enddef', nf90_enddef(ncid)) +! +!do i = 1, self%nvar +! call check('nf90_put_var', nf90_put_var(ncid,ncid_var(i),self%geovals(i)%vals(:,:))) +!enddo +! +!call check('nf90_close', nf90_close(ncid)) +!deallocate(ncid_var) +! +!end subroutine ufo_geovals_write_netcdf + ! ------------------------------------------------------------------------------ subroutine ufo_geovals_fill(self, varname, c_nprofiles, c_indx, c_nlev, c_vals, levels_top_down) diff --git a/src/ufo/variabletransforms/CMakeLists.txt b/src/ufo/variabletransforms/CMakeLists.txt index ef9ccd424..e59c36fcf 100644 --- a/src/ufo/variabletransforms/CMakeLists.txt +++ b/src/ufo/variabletransforms/CMakeLists.txt @@ -8,6 +8,7 @@ set ( variabletransforms_files Cal_AdjustedHeightCoordinate.h Cal_AdjustedHeightCoordinate.cc + Cal_GnssroRefractivityGradient.cc Cal_HeightFromPressure.h Cal_HeightFromPressure.cc Cal_Humidity.h @@ -16,6 +17,7 @@ set ( variabletransforms_files Cal_IceThicknessFromFreeboard.cc Cal_Logarithm.h Cal_Logarithm.cc + Cal_NonLocalPseudoExcessPhase.cc Cal_PotentialTFromT.h Cal_PotentialTFromT.cc Cal_PressureFromHeight.h diff --git a/src/ufo/variabletransforms/Cal_GnssroRefractivityGradient.cc b/src/ufo/variabletransforms/Cal_GnssroRefractivityGradient.cc new file mode 100644 index 000000000..15e0585ea --- /dev/null +++ b/src/ufo/variabletransforms/Cal_GnssroRefractivityGradient.cc @@ -0,0 +1,184 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#include +#include "ufo/filters/VariableTransformParametersBase.h" +#include "ufo/operators/gnssro/utils/RefractivityCalculator.h" +#include "ufo/utils/Constants.h" +#include "ufo/utils/VertInterp.interface.h" +#include "ufo/variabletransforms/Cal_GnssroRefractivityGradient.h" + +namespace ufo { + +/************************************************************************************/ +// Cal_GnssroRefractivityGradient +/************************************************************************************/ + +// Register the variable transform with ufo::TransformFactory +static TransformMaker + makerCal_GnssroRefractivityGradient("GnssroRefractivityGradient"); + +Cal_GnssroRefractivityGradient::Cal_GnssroRefractivityGradient( + const Parameters_ &options, + const ObsFilterData &data, + const std::shared_ptr> &flags, + const std::shared_ptr> &obserr) + : TransformBase(options, data, flags, obserr) + , group_(options.group) + , refrvariable_(options.refrVariable) + , heightvariable_(options.heightVariable) + , gradvariable_(options.gradVariable) + , minSuperRefraction_(options.minSuperRefraction) + , minSRThreshold_(minSuperRefraction_ * Constants::superRefractionCritVal) // N/m + , nearThreshold_(0.6 * Constants::superRefractionCritVal) // N/m + , quiteNearThreshold_(0.7 * Constants::superRefractionCritVal) // N/m + , veryNearThreshold_(0.8 * Constants::superRefractionCritVal) // N/m + , extremelyNearThreshold_(0.9 * Constants::superRefractionCritVal) // N/m + , calcDuctingFlag_(options.calculateDuctingFlag) + , calcProfileDuctingFlag_(options.calculateProfileDuctingFlag) + , extractor_(obsdb_) +{ + // Validate configured setting of min super refraction + if (minSuperRefraction_ < 0.0 || minSuperRefraction_ >= 0.6) { + throw eckit::BadValue("GnssroRefractivityGradient: bad min super refraction value " + + std::to_string(minSuperRefraction_) + + ": must be >= 0 and < near ducting threshold of 0.6"); + } +} + +/************************************************************************************/ + +void Cal_GnssroRefractivityGradient::runTransform(const std::vector &apply) { + const char funcName[] = "Cal_GnssroRefractivityGradient::runTransform"; + oops::Log::trace() << funcName << ": START: obsName: " << obsName() << std::endl; + + // Allocate space for refractivity gradient and optional flag variables. + const size_t nlocs = obsdb_.nlocs(); + std::vector refrGradients(nlocs, missingValueFloat); + + std::vector ductingFlags; + if (calcDuctingFlag_) { + ductingFlags.assign(nlocs, NONDUCTING); + } + std::vector profileDuctingFlags; + if (calcProfileDuctingFlag_) { + profileDuctingFlags.assign(nlocs, NONDUCTING); + } + + // Get observed refractivity and heights. + std::vector refractivities; + std::vector geomHeights; + getObservation(group_, refrvariable_, refractivities, true); + getObservation("MetaData", heightvariable_, geomHeights, true); + + double dmiss = util::missingValue(); + + // Iterate over GNSSRO profiles found by the GnssroProfileExtractor. + for (GnssroProfileExtractor::const_iterator itr = extractor_.cbegin(); + itr != extractor_.cend(); ++itr) + { + oops::Log::debug() << funcName << ": Profile slice " + << itr->seqNum() << " from " << itr->start() + << " to " << itr->end() << std::endl; + + // Iterate over all the levels in the profile, excluding the topmost. + // dN/dz for a given level will be determined by the differences + // between that level and the level above it. + double maxDuctingFlag = NONDUCTING; + std::size_t topObIdx = itr->end() - 1; + for (std::size_t obIdx = itr->start(); obIdx < topObIdx; ++obIdx) + { + // Do not compute derived values if ob has been excluded by the where statement + if (!apply[obIdx]) continue; + + // If the refractivity or height is missing, skip this ob. + if (refractivities[obIdx] == dmiss) continue; + if (geomHeights[obIdx] == dmiss) continue; + + // Compute dN/dZ in units N-units / meter + refrGradients[obIdx] = (refractivities[obIdx+1] - refractivities[obIdx]) + / (geomHeights[obIdx+1] - geomHeights[obIdx]); + + oops::Log::debug() << funcName << ", seqNum=" << itr->seqNum() + << ", obIdx=" << obIdx << ", Ntop=" << refractivities[obIdx+1] + << ", Nbot=" << refractivities[obIdx] + << ", Htop=" << geomHeights[obIdx+1] + << ", Hbot=" << geomHeights[obIdx] + << ", dN/dz=" << refrGradients[obIdx] + << ", minSRThreshold=" << minSRThreshold_ << std::endl; + + // Set optional flags indicating non-ducting, near-ducting, or ducting + if (calcDuctingFlag_ || calcProfileDuctingFlag_) { + double absRefGradient = std::abs(refrGradients[obIdx]); + int ductingFlag = NONDUCTING; + if (absRefGradient >= Constants::superRefractionCritVal) { + ductingFlag = DUCTING; + } else if (absRefGradient >= extremelyNearThreshold_) { + ductingFlag = EXTREMELY_NEAR_DUCTING; + } else if (absRefGradient >= veryNearThreshold_) { + ductingFlag = VERY_NEAR_DUCTING; + } else if (absRefGradient >= quiteNearThreshold_) { + ductingFlag = QUITE_NEAR_DUCTING; + } else if (absRefGradient >= nearThreshold_) { + ductingFlag = NEAR_DUCTING; + } else if (absRefGradient >= minSRThreshold_) { + ductingFlag = MIN_SUPER_REFRACTION; + } + + if (ductingFlag != 0) { + oops::Log::debug() << funcName << ", seqNum=" << itr->seqNum() + << ", obIdx=" << obIdx << ", detected ductingFlag=" + << ductingFlag << std::endl; + } + + if (calcDuctingFlag_) { + ductingFlags[obIdx] = ductingFlag; + } + if (calcProfileDuctingFlag_) { + if (ductingFlag > maxDuctingFlag) { + maxDuctingFlag = ductingFlag; + } + } + } // end logic for computing ducting flags. + } // end iteration over observations in a single RO profile + + // Set ducting flag for the entire profile based on the level with the + // more extreme ducting properties. + if (calcProfileDuctingFlag_ && maxDuctingFlag != NONDUCTING) { + std::fill(profileDuctingFlags.begin() + itr->start(), + profileDuctingFlags.begin() + itr->end(), + maxDuctingFlag); + } + oops::Log::debug() << funcName << ", seqNum=" << itr->seqNum() + << ", profileDuctingFlag=" << maxDuctingFlag << std::endl; + } // end iteration over RO profiles + + // put new variable at existing locations + putObservation(gradvariable_, refrGradients, getDerivedGroup(group_)); + if (calcDuctingFlag_) { + putObservation("sampleDuctingFlag", ductingFlags, getDerivedGroup(group_)); + } + if (calcProfileDuctingFlag_) { + putObservation("profileDuctingFlag", profileDuctingFlags, getDerivedGroup(group_)); + } + oops::Log::trace() << funcName << ": DONE: obsName: " << obsName() << ", added " + << refrGradients.size() << " obs to " << gradvariable_ << " in group " + << getDerivedGroup(group_) << std::endl; + return; +} + +} // namespace ufo diff --git a/src/ufo/variabletransforms/Cal_GnssroRefractivityGradient.h b/src/ufo/variabletransforms/Cal_GnssroRefractivityGradient.h new file mode 100644 index 000000000..13e6d4b57 --- /dev/null +++ b/src/ufo/variabletransforms/Cal_GnssroRefractivityGradient.h @@ -0,0 +1,100 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#ifndef UFO_VARIABLETRANSFORMS_CAL_GNSSROREFRACTIVITYGRADIENT_H_ +#define UFO_VARIABLETRANSFORMS_CAL_GNSSROREFRACTIVITYGRADIENT_H_ + +#include +#include +#include +#include +#include + +#include "ufo/operators/gnssro/utils/GnssroProfileExtractor.h" +#include "ufo/variabletransforms/TransformBase.h" + +namespace ufo { + +/// Configuration parameters for the wind components transform. +class Cal_GnssroRefractivityGradientParameters: public VariableTransformParametersBase { + OOPS_CONCRETE_PARAMETERS(Cal_GnssroRefractivityGradientParameters, + VariableTransformParametersBase); + + public: + /// Observation group name. Default is ObsValue. + oops::Parameter group{"group", "ObsValue", this}; + oops::Parameter refrVariable{"refractivity variable", + "atmosphericRefractivity", this}; + oops::Parameter heightVariable{"height variable", + "height", this}; + oops::Parameter gradVariable{"gradient variable", + "atmosphericRefractivityGradient", this}; + oops::Parameter minSuperRefraction{"min super refraction", 0.5, this}; + oops::Parameter calculateDuctingFlag{"calculate ducting flag", false, this}; + oops::Parameter calculateProfileDuctingFlag{"calculate profile ducting flag", + false, this}; +}; + +/*! +* \brief GNSSRO refractivity gradient filter +* +* \details Performs the vertical derivative of refractivity wrt geometric height dN/dz. +* The newly calculated variables are included in the same obs space in a group with +* "Derived" as the prefix; default is DerivedObsValue. +/// +* +* See VariableTransformParametersBase for filter setup. +*/ +class Cal_GnssroRefractivityGradient: public TransformBase { + public: + static constexpr int NONDUCTING = 0; + static constexpr int MIN_SUPER_REFRACTION = 1; + static constexpr int NEAR_DUCTING = 2; + static constexpr int QUITE_NEAR_DUCTING = 3; + static constexpr int VERY_NEAR_DUCTING = 4; + static constexpr int EXTREMELY_NEAR_DUCTING = 5; + static constexpr int DUCTING = 6; + + typedef Cal_GnssroRefractivityGradientParameters Parameters_; + + Cal_GnssroRefractivityGradient(const Parameters_ &options, + const ObsFilterData &data, + const std::shared_ptr> &flags, + const std::shared_ptr> &obserr); + // Run variable conversion + void runTransform(const std::vector &apply) override; + + private: + // Member data + std::string refrvariable_; + std::string heightvariable_; + std::string gradvariable_; + std::string group_; + float minSuperRefraction_; // Smallest fraction of critical value tracked. + float minSRThreshold_; // units: N-units / meter + float nearThreshold_; + float quiteNearThreshold_; + float veryNearThreshold_; + float extremelyNearThreshold_; + bool calcDuctingFlag_; + bool calcProfileDuctingFlag_; + GnssroProfileExtractor extractor_; +}; + +} // namespace ufo + +#endif // UFO_VARIABLETRANSFORMS_CAL_GNSSROREFRACTIVITYGRADIENT_H_ diff --git a/src/ufo/variabletransforms/Cal_NonLocalPseudoExcessPhase.cc b/src/ufo/variabletransforms/Cal_NonLocalPseudoExcessPhase.cc new file mode 100644 index 000000000..cf2f086b7 --- /dev/null +++ b/src/ufo/variabletransforms/Cal_NonLocalPseudoExcessPhase.cc @@ -0,0 +1,165 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#include "ufo/filters/VariableTransformParametersBase.h" +#include "ufo/operators/gnssro/utils/RefractivityCalculator.h" +#include "ufo/utils/Constants.h" +#include "ufo/utils/VertInterp.interface.h" +#include "ufo/variabletransforms/Cal_NonLocalPseudoExcessPhase.h" + +namespace ufo { + +/************************************************************************************/ +// Cal_NonLocalPseudoExcessPhase +/************************************************************************************/ + +// Register the variable transform with ufo::TransformFactory +static TransformMaker + makerCal_NonLocalPseudoExcessPhase("NonLocalPseudoExcessPhase"); + +Cal_NonLocalPseudoExcessPhase::Cal_NonLocalPseudoExcessPhase( + const Parameters_ &options, + const ObsFilterData &data, + const std::shared_ptr> &flags, + const std::shared_ptr> &obserr) + : TransformBase(options, data, flags, obserr) + , group_(options.group) + , refrvariable_(options.refrVariable) + , nlpepvariable_(options.nlpepVariable) + , rpParams_(options.rayPathGenType, options.approxRayLength, + options.top2D, options.res, options.nHoriz) + , orchestrator_(obsdb_, rpParams_) +{} + +/************************************************************************************/ + +void Cal_NonLocalPseudoExcessPhase::runTransform(const std::vector &apply) { + const char funcName[] = "Cal_NonLocalPseudoExcessPhase::runTransform"; + oops::Log::trace() << funcName << ": START: obsName: " << obsName() + << ", numProfiles=" << orchestrator_.profiles().size() + << ", nlocs=" << obsdb_.nlocs() << std::endl; + + // Confirm our sizing for number of obs is what we expect. + const size_t nlocs = obsdb_.nlocs(); + std::size_t totalNumRays = orchestrator_.totalNumRays(); + if (totalNumRays != nlocs) + { + oops::Log::error() << funcName << ": size mismatch between obsdb nlocs (" + << nlocs << ") and Orchestrator's total number of rays (" + << totalNumRays << ")" << std::endl; + throw eckit::BadValue("size mismatch between obsdb nlocs and Orchestrator in " + "Cal_NonLocalPseudoExcessPhase", Here()); + } + + // Allocate space for nonLocalPseudoExcessPhase + std::vector nlpep(nlocs, missingValueFloat); + + // Get observed refractivity + std::vector refractivities; + std::vector geomHeights; + getObservation(group_, refrvariable_, refractivities, true); + getObservation("MetaData", "height", geomHeights, true); + + double dmiss = util::missingValue(); + + std::size_t ovecIdx = 0; // Index into locations in ObsVector (nlocs) + // Iterate over RO profiles. + const GnssroRayPathOrchestrator::Profiles_& profiles = orchestrator_.profiles(); + for (const auto& roProfile : profiles) + { + int seqNum = roProfile->seqNum(); + oops::Log::debug() << funcName << ": processing profile seqNum=" << seqNum + << " with " << roProfile->numRays() << " rays" << std::endl; + + // Get subset of the observation arrays associated with this RO profile. + const double * profileGeomHeights = geomHeights.data() + ovecIdx; + const double * profileRefractivities = refractivities.data() + ovecIdx; + + // Iterate over rays within this profile. + for (std::size_t rayIdx = 0; rayIdx < roProfile->numRays(); ++rayIdx, ++ovecIdx) + { + const GnssroRayPathOrchestrator::Ray_& ray = roProfile->ray(rayIdx); + + // Do not compute NLPEP if this obs has been excluded by the where statement + if (!apply[ovecIdx]) continue; + + // If the tangent point refractivity is missing, skip this ob. + if (refractivities[ovecIdx] == dmiss) continue; + + std::size_t tpNodeIdx = ray->tpNodeIdx(); + double accumNLPEP = 0.0; + for (std::size_t nodeIdx = 0; nodeIdx < ray->numNodes(); ++nodeIdx) + { + // Get refractivity at the average node height, assuming a + // spherically symmetric atmosphere. This allows us to treat + // the observations in the profile as though they formed a + // vertical column. + double refr = dmiss; + double geomHgt = dmiss; + if (nodeIdx == tpNodeIdx) { + // Tangent point of ray: use the observed values. + geomHgt = geomHeights[ovecIdx]; + refr = refractivities[ovecIdx]; + } else { // Average of heights of edges of this ray segment. + float thisGeomHgt = ray->geomHgt(nodeIdx); + float nextGeomHgt = ray->geomHgt(nodeIdx + 1); + geomHgt = static_cast(0.5 * (thisGeomHgt + nextGeomHgt)); + + // Get vertical interpolation weights using geometric height. + // Implementation note: If the height is out-of-bounds, this will give us + // weights for 100% of the topmost or bottommost value. + int wi; // Base vertical index for vertical interpolation (1-based). + double wf; // Weight for base index (rest of weight goes to idx wi + 1). + vert_interp_weights_f90(roProfile->numRays(), geomHgt, profileGeomHeights, wi, wf); + + // Vertically interpolate refractivity to target geometric height + vert_interp_apply_f90(roProfile->numRays(), profileRefractivities, refr, wi, wf); + } + + double segmentLength = static_cast(ray->segLen(nodeIdx)); + if (refr != dmiss) { + accumNLPEP += (refr * segmentLength); + } else { + oops::Log::warning() << funcName << ": seqNum=" << seqNum << ", rayIdx=" << rayIdx + << ", ovecIdx=" << ovecIdx << ", nodeIdx=" << nodeIdx << ", tpNodeIdx=" + << tpNodeIdx << ", numNodes=" << ray->numNodes() << ", geomHgt=" << geomHgt + << ", tpRefr=" << refractivities[ovecIdx] + << ", refr is missing for this node, so NLPEP set to missing" + << std::endl; + accumNLPEP = dmiss; + break; + } + } // end of iteration over nodes in a single ray + + if (accumNLPEP != dmiss) { + // All the nodes had non-missing data, so we could get a complete accumulation. + // The factor converts units from Refractivity (N) to excess index of refraction. + // The resulting units are in meters. + nlpep[ovecIdx] = accumNLPEP * RefractivityCalculator::N_TO_EXCESS_IOR; + } + } // end iteration over rays in a single RO profile + } // end iteration over RO profiles + + // put new variable at existing locations + putObservation(nlpepvariable_, nlpep, getDerivedGroup(group_)); + oops::Log::trace() << funcName << ": DONE: obsName: " << obsName() << ", added " + << nlpep.size() << " obs to " << nlpepvariable_ << " in group " + << getDerivedGroup(group_) << std::endl; + return; +} + +} // namespace ufo diff --git a/src/ufo/variabletransforms/Cal_NonLocalPseudoExcessPhase.h b/src/ufo/variabletransforms/Cal_NonLocalPseudoExcessPhase.h new file mode 100644 index 000000000..0dc021ca5 --- /dev/null +++ b/src/ufo/variabletransforms/Cal_NonLocalPseudoExcessPhase.h @@ -0,0 +1,92 @@ +/* + * (C) Copyright 2025 Space Sciences and Engineering, LLC (dba PlanetiQ). + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * author Steve Marshall (smarshall@planetiq.com) + */ +#ifndef UFO_VARIABLETRANSFORMS_CAL_NONLOCALPSEUDOEXCESSPHASE_H_ +#define UFO_VARIABLETRANSFORMS_CAL_NONLOCALPSEUDOEXCESSPHASE_H_ + +#include +#include +#include +#include +#include + +#include "ufo/operators/gnssro/utils/GnssroProfileRayPaths.h" +#include "ufo/operators/gnssro/utils/GnssroRayPathOrchestrator.h" +#include "ufo/operators/gnssro/utils/GnssroRayPathParameters.h" +#include "ufo/variabletransforms/TransformBase.h" + +namespace ufo { + +/// Configuration parameters for the wind components transform. +class Cal_NonLocalPseudoExcessPhaseParameters: public VariableTransformParametersBase { + OOPS_CONCRETE_PARAMETERS(Cal_NonLocalPseudoExcessPhaseParameters, + VariableTransformParametersBase); + + public: + /// Observation group name. Default is ObsValue. + oops::Parameter group{"group", "ObsValue", this}; + oops::Parameter refrVariable{"refractivity variable", + "atmosphericRefractivity", this}; + oops::Parameter nlpepVariable{"nlpep variable", + "nonLocalPseudoExcessPhase", this}; + oops::Parameter rayPathGenType{"ray_path_gen_type", + GnssroRayPathParameters::DEFAULT_RAY_PATH_GEN_TYPE, this}; + oops::Parameter approxRayLength{"ray_length", + GnssroRayPathParameters::DEFAULT_APPROX_RAY_LENGTH_KM, this}; + oops::Parameter res{"res", + GnssroRayPathParameters::DEFAULT_HORIZONTAL_RES_KM, this}; + oops::Parameter top2D{"top_2d", + GnssroRayPathParameters::DEFAULT_TOP2D_KM, this}; + oops::Parameter nHoriz{"n_horiz", + GnssroRayPathParameters::DEFAULT_NHORIZ, this}; +}; + +/*! +* \brief nonLocalPseudoExcessPhase filter +* +* \details Performs a variable conversion from atmosphericRefractivity to nonLocalPseudoExcessPhase. +* nonLocalPseudoExcessPhase is a path-integrated form of refractivity, dependent on +* a pre-defined ray path. This ray path must match the ray path used by an associated +* GNSSRO forward operator, e.g. RefNLPEP2D. +* The newly calculated variables are included in the same obs space. +/// +* +* See VariableTransformParametersBase for filter setup. +*/ +class Cal_NonLocalPseudoExcessPhase: public TransformBase { + public: + typedef Cal_NonLocalPseudoExcessPhaseParameters Parameters_; + + Cal_NonLocalPseudoExcessPhase(const Parameters_ &options, + const ObsFilterData &data, + const std::shared_ptr> &flags, + const std::shared_ptr> &obserr); + // Run variable conversion + void runTransform(const std::vector &apply) override; + + private: + // Member data + std::string refrvariable_; + std::string nlpepvariable_; + std::string group_; + GnssroRayPathParameters rpParams_; + GnssroRayPathOrchestrator orchestrator_; +}; + +} // namespace ufo + +#endif // UFO_VARIABLETRANSFORMS_CAL_NONLOCALPSEUDOEXCESSPHASE_H_ diff --git a/test/testinput/unit_tests/operators/CMakeLists.txt b/test/testinput/unit_tests/operators/CMakeLists.txt index ef4c56b79..7050848a2 100644 --- a/test/testinput/unit_tests/operators/CMakeLists.txt +++ b/test/testinput/unit_tests/operators/CMakeLists.txt @@ -395,7 +395,26 @@ ufo_add_test( NAME test_ufo_linopr_groundgnssmetoffice LIBS ufo LABELS operators WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/../../../ ) - +# + ufo_add_test( NAME test_ufo_opr_gnssrorefnlpep2d + TIER 1 + ECBUILD + COMMAND ${CMAKE_BINARY_DIR}/bin/test_ObsOperator.x + ARGS "${CMAKE_CURRENT_SOURCE_DIR}/gnssrorefnlpep2d.yaml" + MPI 1 + LIBS ufo + LABELS operators + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/../../../ ) +# + ufo_add_test( NAME test_ufo_linopr_gnssrorefnlpep2d + TIER 1 + ECBUILD + COMMAND ${CMAKE_BINARY_DIR}/bin/test_ObsOperatorTLAD.x + ARGS "${CMAKE_CURRENT_SOURCE_DIR}/gnssrorefnlpep2d.yaml" + MPI 1 + LIBS ufo + LABELS operators + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/../../../ ) # Tests specific to ROPP if( ${ropp-ufo_FOUND} ) diff --git a/test/testinput/unit_tests/operators/gnssrorefnlpep2d.yaml b/test/testinput/unit_tests/operators/gnssrorefnlpep2d.yaml new file mode 100644 index 000000000..8e66ac06d --- /dev/null +++ b/test/testinput/unit_tests/operators/gnssrorefnlpep2d.yaml @@ -0,0 +1,35 @@ +time window: + begin: 2025-07-03T21:00:00Z + end: 2025-07-04T03:00:00Z +observations: +- obs operator: + name: GnssroRefNLPEP2D + obs options: + ray_path_gen_type: "StraightLine" + res: 11.0 + top_2d: 22.0 + n_horiz: 19 + refr_algo: "RuegerBevis" + use_compress: true + obs space: + name: GnssroRef + obsdatain: + engine: + type: H5File + obsfile: Data/ufo/testinput_tier_1/gnssro_obs_2025070400_refnlpep2d.nc4 + observed variables: [] + derived variables: [nonLocalPseudoExcessPhase] + simulated variables: [nonLocalPseudoExcessPhase] +# obsdataout: +# engine: +# type: H5File +# obsfile: Data/gnssro_obs_2025070400_refnlpep2d_nlpep_opr_out.nc4 + geovals: + filename: Data/ufo/testinput_tier_1/gnssro_geoval_2025070400_refnlpep2d.nc4 + rms ref: 35.181904131227647 + tolerance: 1.0e-11 + linear obs operator test: + coef TL: 0.1 + tolerance TL: 6.0e-4 + tolerance AD: 1.0e-13 + iterations TL: 10 diff --git a/test/testinput/unit_tests/variabletransforms/CMakeLists.txt b/test/testinput/unit_tests/variabletransforms/CMakeLists.txt index 8565a033d..dfde13d0a 100644 --- a/test/testinput/unit_tests/variabletransforms/CMakeLists.txt +++ b/test/testinput/unit_tests/variabletransforms/CMakeLists.txt @@ -337,4 +337,26 @@ ufo_add_test( NAME test_ufo_variabletransforms_radarbeamgeometry LABELS unit_tests variabletransforms WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/../../../ TEST_DEPENDS ufo_get_ufo_test_data ) +# +ufo_add_test( NAME test_ufo_variabletransforms_nonlocalpseudoexcessphase + TIER 1 + ECBUILD + COMMAND ${CMAKE_BINARY_DIR}/bin/test_ObsFilters.x + ARGS "${CMAKE_CURRENT_SOURCE_DIR}/variabletransforms_nonlocalpseudoexcessphase.yaml" + MPI 1 + LIBS ufo + LABELS unit_tests variabletransforms + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/../../../ + TEST_DEPENDS ufo_get_ufo_test_data ) +# +ufo_add_test( NAME test_ufo_variabletransforms_gnssrorefractivitygradient + TIER 1 + ECBUILD + COMMAND ${CMAKE_BINARY_DIR}/bin/test_ObsFilters.x + ARGS "${CMAKE_CURRENT_SOURCE_DIR}/variabletransforms_gnssrorefractivitygradient.yaml" + MPI 1 + LIBS ufo + LABELS unit_tests variabletransforms + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/../../../ + TEST_DEPENDS ufo_get_ufo_test_data ) diff --git a/test/testinput/unit_tests/variabletransforms/variabletransforms_gnssrorefractivitygradient.yaml b/test/testinput/unit_tests/variabletransforms/variabletransforms_gnssrorefractivitygradient.yaml new file mode 100644 index 000000000..de32ab2c7 --- /dev/null +++ b/test/testinput/unit_tests/variabletransforms/variabletransforms_gnssrorefractivitygradient.yaml @@ -0,0 +1,81 @@ +time window: + begin: 2025-07-03T21:00:00Z + end: 2025-07-04T03:00:00Z +observations: +- obs space: + name: GnssroRefGradient with ducting flags + obsdatain: + engine: + type: H5File + obsfile: Data/ufo/testinput_tier_1/gnssrorefractivitygradient_transforms.nc4 + observed variables: [atmosphericRefractivity] + derived variables: [atmosphericRefractivityGradient] + simulated variables: [atmosphericRefractivityGradient] + obsdataout: + engine: + type: H5File + obsfile: Data/gnssrorefractivitygradient_transforms_output_with_ducting_flags.nc4 + obs filters: + - filter: Variable Transforms + Transform: GnssroRefractivityGradient + UseValidDataOnly: false + SkipWhenNoObs: false + group: "ObsValue" + min super refraction: 0.45 + calculate ducting flag: true + calculate profile ducting flag: true + - filter: Perform Action + action: + name: assign error + error parameter: 1 + passedBenchmark: 48 # Topmost levels of 2 RO profiles have missing atmosphericRefractivityGradient + compareVariables: # test output matches precalculated values + - reference: + name: TestReference/atmosphericRefractivityGradient + test: + name: DerivedObsValue/atmosphericRefractivityGradient + relTol: 5.0e-6 + - reference: + name: TestReference/sampleDuctingFlag + test: + name: DerivedObsValue/sampleDuctingFlag + - reference: + name: TestReference/profileDuctingFlag + test: + name: DerivedObsValue/profileDuctingFlag + +- obs space: + name: GnssroRefGradient without ducting flags + obsdatain: + engine: + type: H5File + obsfile: Data/ufo/testinput_tier_1/gnssrorefractivitygradient_transforms.nc4 + observed variables: [atmosphericRefractivity] + derived variables: [atmosphericRefractivityGradient] + simulated variables: [atmosphericRefractivityGradient] + obsdataout: + engine: + type: H5File + obsfile: Data/gnssrorefractivitygradient_transforms_output_wo_ducting_flags.nc4 + obs filters: + - filter: Variable Transforms + Transform: GnssroRefractivityGradient + UseValidDataOnly: false + SkipWhenNoObs: false + group: "ObsValue" + calculate ducting flag: false + calculate profile ducting flag: false + - filter: Perform Action + action: + name: assign error + error parameter: 1 + passedBenchmark: 48 # Topmost levels of 2 RO profiles have missing atmosphericRefractivityGradient + compareVariables: # test output matches precalculated values + - reference: + name: TestReference/atmosphericRefractivityGradient + test: + name: DerivedObsValue/atmosphericRefractivityGradient + relTol: 5.0e-6 + expectVariablesNotToExist: + - name: DerivedObsValue/sampleDuctingFlag + - name: DerivedObsValue/profileDuctingFlag diff --git a/test/testinput/unit_tests/variabletransforms/variabletransforms_nonlocalpseudoexcessphase.yaml b/test/testinput/unit_tests/variabletransforms/variabletransforms_nonlocalpseudoexcessphase.yaml new file mode 100644 index 000000000..f9b0a474e --- /dev/null +++ b/test/testinput/unit_tests/variabletransforms/variabletransforms_nonlocalpseudoexcessphase.yaml @@ -0,0 +1,38 @@ +time window: + begin: 2025-07-03T21:00:00Z + end: 2025-07-04T03:00:00Z +observations: +- obs space: + name: GnssroRef + obsdatain: + engine: + type: H5File + obsfile: Data/ufo/testinput_tier_1/nonlocalpseudoexcessphase_transforms.nc4 + observed variables: [atmosphericRefractivity] + derived variables: [nonLocalPseudoExcessPhase] + simulated variables: [nonLocalPseudoExcessPhase] +# obsdataout: +# engine: +# type: H5File +# obsfile: Data/nonlocalpseudoexcessphase_transforms_output.nc4 + obs filters: + - filter: Variable Transforms + Transform: NonLocalPseudoExcessPhase + UseValidDataOnly: true + res: 11.0 + top_2d: 60.0 + n_horiz: 19 + group: "ObsValue" + refractivity variable: "atmosphericRefractivity" + nlpep variable: "nonLocalPseudoExcessPhase" + - filter: Perform Action + action: + name: assign error + error parameter: 1 + passedBenchmark: 50 # All 50 obs pass for nonLocalPseudoExcessPhase + compareVariables: # test output matches precalculated values + - reference: + name: TestReference/nonLocalPseudoExcessPhase + test: + name: DerivedObsValue/nonLocalPseudoExcessPhase + relTol: 5.0e-6