From 1a3bb800b8f17812de6ae26917cd6dc198478a34 Mon Sep 17 00:00:00 2001 From: Cristian Lussana Date: Thu, 12 Sep 2024 19:22:56 +0200 Subject: [PATCH] added structure functions and R_optimal_... functions --- include/gridpp.h | 119 ++++++++-- src/api/oi_ensi.cpp | 489 +++++++++++++++++++++++++++++++++++++++++ src/api/oi_ensi_lr.cpp | 50 ++++- src/api/structure.cpp | 180 +++++++++++++++ swig/gridpp.i | 1 + 5 files changed, 820 insertions(+), 19 deletions(-) diff --git a/include/gridpp.h b/include/gridpp.h index 5c7542c8..05f7c126 100644 --- a/include/gridpp.h +++ b/include/gridpp.h @@ -293,6 +293,36 @@ namespace gridpp { int max_points, bool allow_extrapolation=true); + /** Optimal interpolation for an ensemble point field (R bindings) + * See Lussana et al 2019 (DOI: 10.1002/qj.3646) + * @param bpoints Points of background field + * @param background 2D vector of background values (L, E) + * @param obs_points Observation points + * @param obs 1D vector of observations + * @param obs_standard_deviations 1D vector of observation standard deviations + * @param background_at_points 2D vector of background at observation points (L, E) + * @param which_structfun structure function to use (0=Barnes;1=MixA) + * @param dh length scale for the horizontal structure function + * @param dz length scale for the vertical structure function + * @param dw minimum value of the correlation coefficient for laf structure function + * @param max_points Maximum number of observations to use inside localization zone; Use 0 to disable + * @param allow_extrapolation Allow OI to extrapolate increments outside increments at observations + * @returns 2D vector of analised values (L, E) + */ + vec2 R_optimal_interpolation_ensi(const Points& bpoints, + const vec2& background, + const Points& obs_points, + const vec& obs, + const vec& obs_standard_deviations, + const vec2& background_at_points, +/* const StructureFunction& structure, */ + int which_structfun, + float dh, + float dz, + float dw, + int max_points, + bool allow_extrapolation=true); + /** Optimal interpolation for an ensemble gridded field (alternative version) * Work in progress * @param bgrid Grid of background field @@ -324,26 +354,27 @@ namespace gridpp { int max_points, bool allow_extrapolation=true); */ - /** Optimal interpolation for an ensemble gridded field (alternative version) + /** Optimal interpolation for an ensemble gridded field (alternative version that works with R bindings) * Work in progress - * @param bgrid Grid of background field - * @param background 3D vector of (left) background values to update (Y, X, E) - * @param background 3D vector of (LEFT) background values (Y, X, E) used to compute correlations - * @param obs_points observation points - * @param obs 2D vector of perturbed observations (S, E) - * @param background 3D vector of (right) background values used to compute innovations (Y, X, E) - * @param background 3D vector of (RIGHT) background values (Y, X, E) used to compute correlations - * @param dh structure Structure function par 1 - * @param dz structure Structure function par 2 - * @param dw structure Structure function par 3 + * @param bpoints Points of background field + * @param background 2D vector of (left) background values to update (M, E) M=num. grid points + * @param background 2D vector of (LEFT) background values (M, E) used to compute correlations + * @param obs_points Observation points + * @param obs 2D vector of perturbed observations (S, E) S=num. obs points + * @param background 2D vector of (right) background values used to compute innovations (S, E) + * @param background 2D vector of (RIGHT) background values (S, E) used to compute correlations + * @param which_structfun structure function to use (0=Barnes;1=MixA) + * @param dh length scale for the horizontal structure function + * @param dz length scale for the vertical structure function + * @param dw minimum value of the correlation coefficient for laf structure function * @param variance_ratio (ratio of observation to right background error variance) * @param standard deviation ratio (ratio of left to right background error standard deviation) * @param weight given to the analysis increment * @param max_points Maximum number of observations to use inside localization zone; Use 0 to disable * @param allow_extrapolation Allow OI to extrapolate increments outside increments at observations - * @returns 3D vector of analised values (Y, X, E) + * @returns 2D vector of analised values (M, E) */ - vec2 optimal_interpolation_ensi_lr(const Points& bpoints, + vec2 R_optimal_interpolation_ensi_lr(const Points& bpoints, const vec2& background_l, const vec2& background_L, const Points& obs_points, @@ -351,6 +382,7 @@ namespace gridpp { const vec2& pbackground_r, const vec2& pbackground_R, /* const StructureFunction& structure, */ + int which_structfun, float dh, float dz, float dw, @@ -2000,6 +2032,34 @@ namespace gridpp { * @returns Cressman rho */ float cressman_rho(float dist, float length) const; + + /** Compactly supported second-order autoregressive correlation function + * @param dist Distance between points. Same units as 'length' + * @param length Length scale + * @returns SOAR rho + */ + float soar_rho(float dist, float length) const; + + /** Compactly supported third-order autoregressive correlation function + * @param dist Distance between points. Same units as 'length' + * @param length Length scale + * @returns TOAR rho + */ + float toar_rho(float dist, float length) const; + + /** Powerlaw correlation function + * @param dist Distance between points. Same units as 'length' + * @param length Length scale + * @returns powerlaw rho + */ + float powerlaw_rho(float dist, float length) const; + + /** Linear correlation function + * @param normdist Normalized distance between points. Must be in the range -1 to 1. + * @param min Minimum allowed value for the correlation (if less than 0, the return 1) + * @returns linear rho + */ + float linear_rho(float normdist, float min) const; float m_localization_distance; }; class MultipleStructure: public StructureFunction { @@ -2052,6 +2112,39 @@ namespace gridpp { bool m_is_spatial; }; + /** Mix of structure functions based on distance(SOAR), elevation(Gauss), and land area fraction(linear) */ + class MixAStructure: public StructureFunction { + public: + /** Structure function SOAR - Gaussian - Linear + * @param h: Horizontal decorrelation length >=0 [m] (SOAR) + * @param v: Vertical decorrelation length >=0 [m] (Gaussian). If 0, disable decorrelation. + * @param w: Land/sea decorrelation length >=0 [1] (Linear). If 0, disable decorrelation. + * @param hmax: Truncate horizontal correlation beyond this length [m]. If undefined, 3.64 * h. + */ + MixAStructure(float h, float v=0, float w=0, float hmax=MV); + + /** MixA structure function where decorrelation varyies spatially + * @param grid: Grid of decorrelation field + * @param h: 2D vector of horizontal decorrelation lengths >=0 (SOAR) , same size as grid [m] + * @param v: 2D vector of Vertical decorrelation lengths >=0 [m] (Gaussian). Set all to 0 to disable decorrelation. + * @param w: 2D vector of land/sea decorrelation lengths >=0 [1] (Linear). Set all to 0 to disable decorrelation. + * @param min_rho: Truncate horizontal correlation when rho is less than this value [m]. + */ + MixAStructure(Grid grid, vec2 h, vec2 v, vec2 w, float min_rho=StructureFunction::default_min_rho); + float corr(const Point& p1, const Point& p2) const; + vec corr(const Point& p1, const std::vector& p2) const; + StructureFunctionPtr clone() const; + float localization_distance(const Point& p) const; + private: + float localization_distance(float h) const; + Grid m_grid; + vec2 mH; + vec2 mV; + vec2 mW; + float m_min_rho; + bool m_is_spatial; + }; + /** Simple structure function based on distance, elevation, and land area fraction */ class CressmanStructure: public StructureFunction { public: diff --git a/src/api/oi_ensi.cpp b/src/api/oi_ensi.cpp index 4f850088..e75afaa9 100644 --- a/src/api/oi_ensi.cpp +++ b/src/api/oi_ensi.cpp @@ -110,6 +110,7 @@ vec3 gridpp::optimal_interpolation_ensi(const gridpp::Grid& bgrid, } return output; } + vec2 gridpp::optimal_interpolation_ensi(const gridpp::Points& bpoints, const vec2& background, const gridpp::Points& points, @@ -565,3 +566,491 @@ vec2 gridpp::optimal_interpolation_ensi(const gridpp::Points& bpoints, } return output; } + +vec2 gridpp::R_optimal_interpolation_ensi(const gridpp::Points& bpoints, + const vec2& background, + const gridpp::Points& points, + const vec& pobs, // gObs + const vec& psigmas, // pci + const vec2& pbackground, +/* const gridpp::StructureFunction& structure, */ + int which_structfun, + float dh, + float dz, + float dw, + int max_points, + bool allow_extrapolation) { + if(max_points < 0) + throw std::invalid_argument("max_points must be >= 0"); + if(bpoints.get_coordinate_type() != points.get_coordinate_type()) { + throw std::invalid_argument("Both background and observations points must be of same coorindate type (lat/lon or x/y)"); + } + if(background.size() != bpoints.size()) + throw std::invalid_argument("Input field is not the same size as the grid"); + if(pobs.size() != points.size()) + throw std::invalid_argument("Observations and points exception mismatch"); + if(psigmas.size() != points.size()) + throw std::invalid_argument("Sigmas and points size mismatch"); + if(pbackground.size() != points.size()) + throw std::invalid_argument("Background and points size mismatch"); + + float hmax = 7 * dh; +/* + if(which_structfun == 0) { + BarnesStructure structure = BarnesStructure( dh, dz, dw, hmax); + } + else if(which_structfun == 1) { + MixAStructure structure = MixAStructure( dh, dz, dw, hmax); + } + else { + BarnesStructure structure = BarnesStructure( dh, dz, dw, hmax); + } */ + + int nS = points.size(); + if(nS == 0) + return background; + + int mY = -1; // Write debug information for this station index + int num_parameters = 2; + float sigmac = 0.5; + float delta = 1; + bool diagnose = false; + + int nY = background.size(); + int nEns = background[0].size(); + + // Prepare output matrix + vec2 output = background; + + int num_condition_warning = 0; + int num_real_part_warning = 0; + + vec blats = bpoints.get_lats(); + vec blons = bpoints.get_lons(); + vec belevs = bpoints.get_elevs(); + vec blafs = bpoints.get_lafs(); + + vec plats = points.get_lats(); + vec plons = points.get_lons(); + vec pelevs = points.get_elevs(); + vec plafs = points.get_lafs(); + + // Compute Y + vec2 gY = pbackground; + vec gYhat(nS); + for(int i = 0; i < nS; i++) { + float mean = gridpp::calc_statistic(gY[i], gridpp::Mean); + for(int e = 0; e < nEns; e++) { + float value = gY[i][e]; + if(gridpp::is_valid(value) && gridpp::is_valid(mean)) { + gY[i][e] -= mean; + } + } + gYhat[i] = mean; + } + + // Create all objects of type Point (to save time on creation later) + std::vector point_vec; + point_vec.reserve(nS); + for(int i = 0; i < nS; i++) { + point_vec.push_back(points.get_point(i)); + } + + // Calculate number of valid members + int nValidEns = 0; + ivec validEns; + for(int e = 0; e < nEns; e++) { + int numInvalid = 0; + for(int y = 0; y < nY; y++) { + float value = background[y][e]; + if(!gridpp::is_valid(value)) + numInvalid++; + } + if(numInvalid == 0) { + validEns.push_back(e); + nValidEns++; + } + } + + // This causes segmentation fault when building the gridpp pypi package + // 1) Tested removing num_condition_warning and num_real_part_warning from the parallel loop + // but it doesnt seem to help + // #pragma omp parallel for + for(int y = 0; y < nY; y++) { + float lat = blats[y]; + float lon = blons[y]; + float elev = belevs[y]; + float laf = blafs[y]; + Point p1 = bpoints.get_point(y); + float localizationRadius = 0; + if(which_structfun == 0) { + BarnesStructure structure = BarnesStructure( dh, dz, dw, hmax); + localizationRadius = structure.localization_distance(p1); + } + else if(which_structfun == 1) { + MixAStructure structure = MixAStructure( dh, dz, dw, hmax); + localizationRadius = structure.localization_distance(p1); + } + + // Create list of locations for this gridpoint + ivec lLocIndices0 = points.get_neighbours(lat, lon, localizationRadius); + if(lLocIndices0.size() == 0) { + // If we have too few observations though, then use the background + continue; + } + ivec lLocIndices; + lLocIndices.reserve(lLocIndices0.size()); + std::vector > lRhos0; + // Calculate gridpoint to observation rhos + lRhos0.reserve(lLocIndices0.size()); + std::vector p2; + p2.reserve(lLocIndices0.size()); + for(int i = 0; i < lLocIndices0.size(); i++) { + int index = lLocIndices0[i]; + p2.push_back(point_vec[index]); + } + vec rhos(lLocIndices0.size()); + if(which_structfun == 0) { + BarnesStructure structure = BarnesStructure( dh, dz, dw, hmax); + rhos = structure.corr_background(p1, p2); + } + else if(which_structfun == 1) { + MixAStructure structure = MixAStructure( dh, dz, dw, hmax); + rhos = structure.corr_background(p1, p2); + } + for(int i = 0; i < lLocIndices0.size(); i++) { + int index = lLocIndices0[i]; + if(gridpp::is_valid(pobs[index])) { + if(rhos[i] > 0) { + lRhos0.push_back(std::pair(rhos[i], i)); + } + } + } + + // Make sure we don't use too many observations + arma::vec lRhos; + if(max_points > 0 && lRhos0.size() > max_points) { + // If sorting is enabled and we have too many locations, then only keep the best ones based on rho. + // Otherwise, just use the last locations added + lRhos = arma::vec(max_points); + std::sort(lRhos0.begin(), lRhos0.end(), ::sort_pair_first()); + for(int i = 0; i < max_points; i++) { + // The best values start at the end of the array + int index = lRhos0[lRhos0.size() - 1 - i].second; + lLocIndices.push_back(lLocIndices0[index]); + lRhos(i) = lRhos0[lRhos0.size() - 1 - i].first; + } + } + else { + lRhos = arma::vec(lRhos0.size()); + for(int i = 0; i < lRhos0.size(); i++) { + int index = lRhos0[i].second; + lLocIndices.push_back(lLocIndices0[index]); + lRhos(i) = lRhos0[i].first; + } + } + + int lS = lLocIndices.size(); + if(lS == 0) { + // If we have too few observations though, then use the background + continue; + } + + vectype lObs(lS); + vectype lElevs(lS); + vectype lLafs(lS); + for(int i = 0; i < lLocIndices.size(); i++) { + int index = lLocIndices[i]; + lObs[i] = pobs[index]; + lElevs[i] = pelevs[index]; + lLafs[i] = plafs[index]; + } + + // Compute Y (model at obs-locations) + mattype lY(lS, nValidEns); + vectype lYhat(lS); + + for(int i = 0; i < lS; i++) { + // Use the nearest neighbour for this location + int index = lLocIndices[i]; + for(int e = 0; e < nValidEns; e++) { + int ei = validEns[e]; + lY(i, e) = gY[index][ei]; + } + lYhat[i] = gYhat[index]; + } + + // Compute Rinv + mattype Rinv(lS, lS, arma::fill::zeros); + if(num_parameters == 2) { + for(int i = 0; i < lS; i++) { + int index = lLocIndices[i]; + Rinv(i, i) = lRhos[i] / (psigmas[index] * psigmas[index]); + } + } + else if(num_parameters == 3) { + /* + // Inverting the matrix is more complicated, since the radar observations + // have covariances. Therefore invert the covariance matrix for the radar part and + // insert the values into the big inverse matrix. + // std::cout << "Computing R matrix" << std::endl; + // R = get_precipitation_r(gRadarL, pci, lLocIndices, lRhos); + // Compute little R + ivec gRadarIndices; + gRadarIndices.reserve(lS); + ivec lRadarIndices; + lRadarIndices.reserve(lS); + for(int i = 0; i < lS; i++) { + int index = lLocIndices[i]; + if(gRadarL[index] > 0) { + gRadarIndices.push_back(index); + lRadarIndices.push_back(i); + } + } + int lNumRadar = gRadarIndices.size(); + + // Compute R tilde r + mattype radarR(lNumRadar, lNumRadar, arma::fill::zeros); + for(int i = 0; i < lNumRadar; i++) { + for(int j = 0; j < lNumRadar; j++) { + int gIndex_i = gRadarIndices[i]; + int gIndex_j = gRadarIndices[j]; + int lIndex_i = lRadarIndices[i]; + int lIndex_j = lRadarIndices[j]; + if(i == j) { + radarR(i, i) = 1; + } + else { + // Equation 5 + float dist = Util::getDistance(gLocations[gIndex_i].lat(), gLocations[gIndex_i].lon(), gLocations[gIndex_j].lat(), gLocations[gIndex_j].lon(), true); + float h = dist / mHLengthC; + float rho = (1 + h) * exp(-h); + radarR(i, j) = rho; + } + } + } + + float cond = arma::rcond(radarR); + if(cond <= 0) { + std::stringstream ss; + ss << "Condition number of " << cond << " for radar values. Using raw values"; + gridpp::warning(ss.str()); + for(int e = 0; e < nEns; e++) { + (*output)(y, x, e) = (*field)(y, x, e); // Util::MV; + } + continue; + } + + mattype radarRinv(lNumRadar, lNumRadar, arma::fill::zeros); + radarRinv = arma::inv(radarR); + + for(int i = 0; i < lS; i++) { + int index = lLocIndices[i]; + Rinv(i, i) = lRhos[i] / (sigma * sigma * pci[index]); + } + // Overwrite where we have radar pixels + for(int i = 0; i < lNumRadar; i++) { + int ii = lRadarIndices[i]; + for(int j = 0; j < lNumRadar; j++) { + int jj = lRadarIndices[j]; + Rinv(ii, jj) = sqrt(lRhos[ii] * lRhos[jj]) / (sigmac * sigmac) * radarRinv(i, j); + } + } + */ + } + else { + abort(); + } + + // Compute C matrix + // k x nS * nS x nS + mattype C(nValidEns, lS); + C = lY.t() * Rinv; + + mattype Pinv(nValidEns, nValidEns); + float diag = 1 / delta * (nValidEns - 1); + + Pinv = C * lY + diag * arma::eye(nValidEns, nValidEns); + float cond = arma::rcond(Pinv); + if(cond <= 0) { + num_condition_warning++; + continue; + } + + // Compute sqrt of matrix. Armadillo 6.6 has this function, but on many systems this + // is not available. Therefore, compute sqrt using the method found in 6.6 + // cxtype Wcx(nValidEns, nValidEns); + // status = arma::sqrtmat(Wcx, (nValidEns - 1) * P); + // mattype W = arma::real(Wcx); + + mattype P = arma::inv(Pinv); + vectype eigval; + mattype eigvec; + bool status = arma::eig_sym(eigval, eigvec, (nValidEns - 1) * P); + if(!status) { + std::cout << "Cannot find eigenvector:" << std::endl; + std::cout << "Lat: " << lat << std::endl; + std::cout << "Lon: " << lon << std::endl; + std::cout << "Elev: " << elev << std::endl; + std::cout << "Laf: " << laf << std::endl; + std::cout << "Pinv" << std::endl; + print_matrix(Pinv); + std::cout << "P" << std::endl; + print_matrix(P); + std::cout << "Y:" << std::endl; + print_matrix(lY); + std::cout << "lObs:" << std::endl; + print_matrix(lObs); + std::cout << "Yhat" << std::endl; + print_matrix(lYhat); + } + eigval = sqrt(eigval); + mattype Wcx = eigvec * arma::diagmat(eigval) * eigvec.t(); + mattype W = arma::real(Wcx); + + if(W.n_rows == 0) { + num_real_part_warning++; + continue; + } + + // Compute PC + mattype PC(nValidEns, lS); + PC = P * C; + + // Compute w + vectype w(nValidEns); + if(diagnose) + w = PC * (arma::ones(lS)); + else + w = PC * (lObs - lYhat); + + // Add w to W + for(int e = 0; e < nValidEns; e++) { + for(int e2 = 0; e2 < nValidEns; e2 ++) { + W(e, e2) = W(e, e2) + w(e) ; + } + } + + // Compute X (perturbations about model mean) + vectype X(nValidEns); + float total = 0; + int count = 0; + for(int e = 0; e < nValidEns; e++) { + int ei = validEns[e]; + float value = background[y][ei]; + if(gridpp::is_valid(value)) { + X(e) = value; + total += value; + count++; + } + } + float ensMean = total / count; + for(int e = 0; e < nValidEns; e++) { + X(e) -= ensMean; + } + + // Write debugging information + if(y == mY) { + std::cout << "Lat: " << lat << std::endl; + std::cout << "Lon: " << lon << " " << lat << " " << std::endl; + std::cout << "Elev: " << elev << std::endl; + std::cout << "Laf: " << laf << std::endl; + std::cout << "Num obs: " << lS << std::endl; + std::cout << "Num ens: " << nValidEns << std::endl; + std::cout << "rhos" << std::endl; + print_matrix(lRhos); + std::cout << "P" << std::endl; + print_matrix(P); + std::cout << "C" << std::endl; + print_matrix(C); + std::cout << "C * lY" << std::endl; + print_matrix(C * lY); + std::cout << "PC" << std::endl; + print_matrix(PC); + std::cout << "W" << std::endl; + print_matrix(W); + std::cout << "w" << std::endl; + print_matrix(w); + std::cout << "Y:" << std::endl; + print_matrix(lY); + std::cout << "Yhat" << std::endl; + print_matrix(lYhat); + std::cout << "lObs" << std::endl; + print_matrix(lObs); + std::cout << "lObs - Yhat" << std::endl; + print_matrix(lObs - lYhat); + std::cout << "X" << std::endl; + print_matrix(X); + std::cout << "elevs" << std::endl; + print_matrix(lElevs); + std::cout << "lafs" << std::endl; + print_matrix(lLafs); + std::cout << "Analysis increment:" << std::endl; + print_matrix(X.t() * W); + std::cout << "My: " << arma::mean(arma::dot(lObs - lYhat, lRhos) / lS) << std::endl; + } + + // Compute analysis + for(int e = 0; e < nValidEns; e++) { + int ei = validEns[e]; + float total = 0; + for(int k = 0; k < nValidEns; k++) { + total += X(k) * W(k, e); + } + + float currIncrement = total; + + float raw = ensMean; + + /////////////////////////////// + // Anti-extrapolation filter // + /////////////////////////////// + if(!allow_extrapolation) { + // Don't allow a final increment that is larger than any increment + // at station points + float maxInc = arma::max(lObs - (lY[e] + lYhat)); + float minInc = arma::min(lObs - (lY[e] + lYhat)); + if(y == mY) { + std::cout << "Increments: " << maxInc << " " << minInc << " " << currIncrement << std::endl; + } + + // The increment for this member. currIncrement is the increment relative to + // ensemble mean + float memberIncrement = currIncrement - X(e); + // Adjust increment if it gives a member increment that is outside the range + // of the observation increments + if(y == mY) { + std::cout << "Analysis increment: " << memberIncrement << " " << ensMean << " " << currIncrement << " " << X(e) << std::endl; + } + if(maxInc > 0 && memberIncrement > maxInc) { + currIncrement = maxInc + X(e); + } + else if(maxInc < 0 && memberIncrement > 0) { + currIncrement = 0 + X(e); + } + else if(minInc < 0 && memberIncrement < minInc) { + currIncrement = minInc + X(e); + } + else if(minInc > 0 && memberIncrement < 0) { + currIncrement = 0 + X(e); + } + if(y == mY) { + std::cout << "Final increment: " << currIncrement << " " << currIncrement - X(e) << std::endl; + } + } + output[y][ei] = ensMean + currIncrement; + } + } + + if(num_condition_warning > 0) { + std::stringstream ss; + ss << "Condition number error in " << num_condition_warning << " points. Using raw values in those points."; + gridpp::warning(ss.str()); + } + if(num_real_part_warning > 0) { + std::stringstream ss; + ss << "Could not find the real part of W in " << num_real_part_warning << " points. Using raw values in those points."; + gridpp::warning(ss.str()); + } + return output; +} diff --git a/src/api/oi_ensi_lr.cpp b/src/api/oi_ensi_lr.cpp index c667b925..a9d1c193 100644 --- a/src/api/oi_ensi_lr.cpp +++ b/src/api/oi_ensi_lr.cpp @@ -121,7 +121,7 @@ template void print_matrix< ::cxtype>(::cxtype matrix); } return output; } */ -vec2 gridpp::optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, +vec2 gridpp::R_optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, const vec2& background_l, const vec2& background_L, const gridpp::Points& points, @@ -129,6 +129,7 @@ vec2 gridpp::optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, const vec2& pbackground_r, const vec2& pbackground_R, /* const gridpp::StructureFunction& structure, it creates problem with R bindings */ + int which_structfun, float dh, float dz, float dw, @@ -154,7 +155,17 @@ vec2 gridpp::optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, throw std::invalid_argument("Background RIGTH and points size mismatch"); float hmax = 7 * dh; - BarnesStructure structure = BarnesStructure( dh, dz, dw, hmax); +/* + BarnesStructure structure = BarnesStructure( dh, dz, dw, hmax)t; + if(which_structfun == 0) { + BarnesStructure structure = BarnesStructure( dh, dz, dw, hmax); + } + else if(which_structfun == 1) { + MixAStructure structure = MixAStructure( dh, dz, dw, hmax); + } + else { + BarnesStructure structure = BarnesStructure( dh, dz, dw, hmax); + } */ int nS = points.size(); if(nS == 0) @@ -251,8 +262,17 @@ vec2 gridpp::optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, float elev = belevs[y]; float laf = blafs[y]; Point p1 = bpoints.get_point(y); - float localizationRadius = structure.localization_distance(p1); - + /* float localizationRadius = structure.localization_distance(p1); */ + float localizationRadius = 0; + if(which_structfun == 0) { + BarnesStructure structure = BarnesStructure( dh, dz, dw, hmax); + localizationRadius = structure.localization_distance(p1); + } + else if(which_structfun == 1) { + MixAStructure structure = MixAStructure( dh, dz, dw, hmax); + localizationRadius = structure.localization_distance(p1); + } + // Create list of locations for this gridpoint ivec lLocIndices0 = points.get_neighbours(lat, lon, localizationRadius); if(lLocIndices0.size() == 0) { @@ -271,7 +291,16 @@ vec2 gridpp::optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, int index = lLocIndices0[i]; p2.push_back(point_vec[index]); } - vec rhos = structure.corr_background(p1, p2); + vec rhos(lLocIndices0.size()); + if(which_structfun == 0) { + BarnesStructure structure = BarnesStructure( dh, dz, dw, hmax); + rhos = structure.corr_background(p1, p2); + } + else if(which_structfun == 1) { + MixAStructure structure = MixAStructure( dh, dz, dw, hmax); + rhos = structure.corr_background(p1, p2); + } +/* vec rhos = structure.corr_background(p1, p2); */ for(int i = 0; i < lLocIndices0.size(); i++) { int index = lLocIndices0[i]; if(gridpp::is_valid(pobs[index][0])) { @@ -360,7 +389,16 @@ vec2 gridpp::optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, int index_j = lLocIndices[j]; p2[j] = point_vec[index_j]; } - vec corr = structure.corr(p1, p2); + vec corr(lS); + if(which_structfun == 0) { + BarnesStructure structure = BarnesStructure( dh, dz, dw, hmax); + corr = structure.corr_background(p1, p2); + } + else if(which_structfun == 1) { + MixAStructure structure = MixAStructure( dh, dz, dw, hmax); + corr = structure.corr_background(p1, p2); + } + /* vec corr = structure.corr(p1, p2); */ for(int j = 0; j < lS; j++) lLoc2D(i, j) = corr[j]; } // end loop over closer observations diff --git a/src/api/structure.cpp b/src/api/structure.cpp index f198e044..4a7db600 100644 --- a/src/api/structure.cpp +++ b/src/api/structure.cpp @@ -42,6 +42,45 @@ float gridpp::StructureFunction::cressman_rho(float dist, float length) const { return 0; return (length * length - dist * dist) / (length * length + dist * dist); } +/* second order autoregressive function */ +float gridpp::StructureFunction::soar_rho(float dist, float length) const { + if(!gridpp::is_valid(length) || length == 0) + // Disabled + return 1; + if(!gridpp::is_valid(dist)) + return 0; + float v = dist / length; + return (1 + v) * exp(-v); +} +/* third order autoregressive function */ +float gridpp::StructureFunction::toar_rho(float dist, float length) const { + if(!gridpp::is_valid(length) || length == 0) + // Disabled + return 1; + if(!gridpp::is_valid(dist)) + return 0; + float v = dist / length; + return (1 + v + (v * v) / 3) * exp(-v); +} +/* powerlaw function */ +float gridpp::StructureFunction::powerlaw_rho(float dist, float length) const { + if(!gridpp::is_valid(length) || length == 0) + // Disabled + return 1; + if(!gridpp::is_valid(dist)) + return 0; + float v = dist / length; + return 1 / (1 + 0.5 * v * v); +} +/* linear function returns correlation between min and 1. normdist between -1 and 1 */ +float gridpp::StructureFunction::linear_rho(float normdist, float min) const { + if(!gridpp::is_valid(min) || min < 0) + // Disabled + return 1; + if(!gridpp::is_valid(normdist)) + return 0; + return (1 - (1-min) * abs(normdist)); +} float gridpp::StructureFunction::localization_distance(const Point& p) const { return m_localization_distance; } @@ -264,6 +303,147 @@ float gridpp::CressmanStructure::corr(const Point& p1, const Point& p2) const { gridpp::StructureFunctionPtr gridpp::CressmanStructure::clone() const { return std::make_shared(mH, mV, mW); } +/** MixA */ +gridpp::MixAStructure::MixAStructure(float h, float v, float w, float hmax) : + m_is_spatial(false) { + if(gridpp::is_valid(hmax) && hmax < 0) + throw std::invalid_argument("hmax must be >= 0"); + if(!gridpp::is_valid(h) || h < 0) + throw std::invalid_argument("h must be >= 0"); + if(!gridpp::is_valid(v) || v < 0) + throw std::invalid_argument("v must be >= 0"); + if(!gridpp::is_valid(w) || w < 0) + throw std::invalid_argument("w must be >= 0"); + + if(gridpp::is_valid(hmax)) + m_min_rho = (1 + hmax / h) * exp(-hmax / h); + else + m_min_rho = default_min_rho; + vec2 h2(1); + h2[0].push_back(h); + vec2 v2(1); + v2[0].push_back(v); + vec2 w2(1); + w2[0].push_back(w); + mH = h2; + mV = v2; + mW = w2; +} +gridpp::MixAStructure::MixAStructure(Grid grid, vec2 h, vec2 v, vec2 w, float min_rho) : + m_grid(grid), + m_min_rho(min_rho), + mH(h), + mV(v), + mW(w) { + if(mH.size() == 1 && mH[0].size() == 1 && mV.size() == 1 && mV[0].size() == 1 && mW.size() == 1 && mW[0].size() == 1) { + m_is_spatial = false; + } + else { + m_is_spatial = true; + if(grid.size()[0] != h.size() || grid.size()[0] != v.size() || grid.size()[0] != w.size()) + throw std::invalid_argument("Grid size not the same as scale size"); + if(grid.size()[1] != h[0].size() || grid.size()[1] != v[0].size() || grid.size()[1] != w[0].size()) + throw std::invalid_argument("Grid size not the same as scale size"); + } +} +float gridpp::MixAStructure::corr(const Point& p1, const Point& p2) const { + float hdist = gridpp::KDTree::calc_straight_distance(p1, p2); + float rho = 1; + if(m_is_spatial) { + // This part is slower because of the nearest neighbour lookup + ivec I = m_grid.get_nearest_neighbour(p1.lat, p1.lon); + if(I[0] > mH.size()) + throw std::runtime_error("Invalid I[0]"); + if(I[1] > mH[I[0]].size()) + throw std::runtime_error("Invalid I[1]"); + + float h = mH[I[0]][I[1]]; + float v = mV[I[0]][I[1]]; + float w = mW[I[0]][I[1]]; + + // Use h to compute this, so we don't have to get nearest neighbour twice + float loc_dist = localization_distance(h); + if(hdist > loc_dist) + return 0; + + rho = gridpp::StructureFunction::soar_rho(hdist, h); + if(gridpp::is_valid(p1.elev) && gridpp::is_valid(p2.elev)) { + float vdist = p1.elev - p2.elev; + rho *= gridpp::StructureFunction::barnes_rho(vdist, v); + } + if(gridpp::is_valid(p1.laf) && gridpp::is_valid(p2.laf)) { + float lafdist = p1.laf - p2.laf; + rho *= gridpp::StructureFunction::linear_rho(lafdist, w); + } + } + else { + if(hdist > localization_distance(p1)) + return 0; + + rho = gridpp::StructureFunction::soar_rho(hdist, mH[0][0]); + if(gridpp::is_valid(p1.elev) && gridpp::is_valid(p2.elev)) { + float vdist = p1.elev - p2.elev; + rho *= gridpp::StructureFunction::barnes_rho(vdist, mV[0][0]); + } + if(gridpp::is_valid(p1.laf) && gridpp::is_valid(p2.laf)) { + float lafdist = p1.laf - p2.laf; + rho *= gridpp::StructureFunction::linear_rho(lafdist, mW[0][0]); + } + } + return rho; +} +vec gridpp::MixAStructure::corr(const Point& p1, const std::vector& p2) const { + vec output(p2.size()); + if(m_is_spatial) { + ivec I = m_grid.get_nearest_neighbour(p1.lat, p1.lon); + if(I[0] > mH.size()) + throw std::runtime_error("Invalid I[0]"); + if(I[1] > mH[I[0]].size()) + throw std::runtime_error("Invalid I[1]"); + + float h = mH[I[0]][I[1]]; + float v = mV[I[0]][I[1]]; + float w = mW[I[0]][I[1]]; + float loc_dist = localization_distance(h); + for(int i = 0; i < p2.size(); i++) { + float hdist = gridpp::KDTree::calc_straight_distance(p1, p2[i]); + float rho = 0; + if(hdist <= loc_dist) { + rho = gridpp::StructureFunction::soar_rho(hdist, h); + if(gridpp::is_valid(p1.elev) && gridpp::is_valid(p2[i].elev)) { + float vdist = p1.elev - p2[i].elev; + rho *= gridpp::StructureFunction::barnes_rho(vdist, v); + } + if(gridpp::is_valid(p1.laf) && gridpp::is_valid(p2[i].laf)) { + float lafdist = p1.laf - p2[i].laf; + rho *= gridpp::StructureFunction::linear_rho(lafdist, w); + } + } + output[i] = rho; + } + } + else { + for(int i = 0; i < p2.size(); i++) { + output[i] = corr(p1, p2[i]); + } + } + return output; +} +gridpp::StructureFunctionPtr gridpp::MixAStructure::clone() const { + return std::make_shared(m_grid, mH, mV, mW, m_min_rho); +} +float gridpp::MixAStructure::localization_distance(const Point& p) const { + if(m_is_spatial) { + ivec I = m_grid.get_nearest_neighbour(p.lat, p.lon); + return localization_distance(mH[I[0]][I[1]]); + } + else { + return localization_distance(mH[0][0]); + } +} +float gridpp::MixAStructure::localization_distance(float h) const { + return sqrt(-2*log(m_min_rho)) * h; +} /** CrossValidation */ gridpp::CrossValidation::CrossValidation(StructureFunction& structure, float dist) : diff --git a/swig/gridpp.i b/swig/gridpp.i index 4c6ee692..192f233e 100644 --- a/swig/gridpp.i +++ b/swig/gridpp.i @@ -46,6 +46,7 @@ %shared_ptr(gridpp::StructureFunction) %shared_ptr(gridpp::MultipleStructure) %shared_ptr(gridpp::BarnesStructure) +%shared_ptr(gridpp::MixAStructure) %shared_ptr(gridpp::CressmanStructure) %shared_ptr(gridpp::CrossValidation)