diff --git a/include/gridpp.h b/include/gridpp.h index 331446f..bb01d09 100644 --- a/include/gridpp.h +++ b/include/gridpp.h @@ -325,161 +325,159 @@ namespace gridpp { /** Optimal interpolation for an ensemble gridded field (alternative version). This is an experimental method. * @param bgrid Grid of background field - * @param background_l background 3D vector of (left) background values to update (Y, X, E) - * @param background_L background 3D vector of (LEFT) background values (Y, X, E) used to compute correlations + * @param bratios 2D vector (Y, X) representing the ratio of background error standard deviation at grid points to that at observation points. The background at grid points is the value being updated, while the background at station points shares the units and time of the observations, which may differ from the background at grid points. This vector contains coefficients (0-1) that adjust the analysis at grid points, accounting for differences in units and variability between the innovations (observation minus background) and the background at grid points. For example, if trusting the ensemble spread, bratios can be set as the ratio between the ensemble spread of the background to be updated and that used to compute the innovations. If the ensemble spread is not trusted at specific times or grid points, the bratios can be based on a typical expected ratio of spreads from multiple ensemble background realizations. + * @param background 3D vector (Y, X, E) representing the background values at grid points to be updated. + * @param background_corr 3D vector (Y, X, E) representing the background values used to compute dynamic correlations. * @param obs_points observation points. S = num. observations - * @param obs 2D vector of perturbed observations (S, E) - * @param pbackground_r background 2D vector of (right) background values used to compute innovations (S, E) - * @param pbackground_R background 2D vector of (RIGHT) background values (S, E) used to compute correlations + * @param pobs 2D vector of perturbed observations (S, E) + * @param pratios 1D vector (S) representing the ratio of observation to background error variance. These coefficients (0-1) indicate the relative trust in observations versus the background. A value of 1 means equal trust in both, while values close to 0 indicate greater trust in the observations. For example, a value of 0.1 means the observations are trusted 10 times more than the background. + * @param pbackground 2D vector (S, E) representing the background values at observation points used to compute innovations. + * @param pbackground_corr 2D vector (S, E) representing the background values at observation points used to compute dynamic correlations. * @param structure Structure function - * @param var_ratios_or variance_ratio (ratio of observation to right background error variance) - * @param std_ratios_lr standard deviation ratio (ratio of left to right background error standard deviation) - * @param weight given to the analysis increment + * @param bweights 2D vector (Y, X) representing user-defined weights at grid points. The analysis at each grid point is calculated as the background value plus the analysis increment, multiplied by the weight at that point. These weights are useful when iterating over multiple observation times, for instance, allowing the user to prioritize observations that match the time of the original background. For example, if running the function over three observation times (the same as the original background, one hour prior, and two hours prior), setting the weight for the original background time to 0.8 and 0.1 for the other two times ensures that most of the modification comes from the analysis at the original background time. * @param max_points Maximum number of observations to use inside localization zone; Use 0 to disable - * @param dynamic_correlations Use ensemble-based correlations. If true, the structure function is used to define localization functions. If false, the structure function is used to define the correlations. + * @param dynamic_correlations Determines whether to use flow-dependent correlations derived from the ensembles. If true, the structure function defines localization functions. If false, the structure function defines static (non-flow-dependent) correlations. * @param allow_extrapolation Allow EnSI to extrapolate increments outside increments at observations * @returns 3D vector of analised values (Y, X, E) */ - vec3 optimal_interpolation_ensi_lr(const Grid& bgrid, - const vec3& background_l, - const vec3& background_L, + vec3 optimal_interpolation_ensi_multi(const Grid& bgrid, + const vec2& bratios, + const vec3& background, + const vec3& background_corr, const Points& obs_points, - const vec2& obs, - const vec2& pbackground_r, - const vec2& pbackground_R, + const vec2& pobs, + const vec& pratios, + const vec2& pbackground, + const vec2& pbackground_corr, const StructureFunction& structure, - float var_ratios_or, - float std_ratios_lr, - float weight, + const vec2& bweights, int max_points, bool dynamic_correlations=true, bool allow_extrapolation=true); /** Optimal interpolation for an ensemble gridded field (alternative version) with ensemble-based correlations. This is an experimental method. - * @param bpoints Points of background field - * @param background_l background 2D vector of (left) background values to update (M, E) M=num. grid points - * @param background_L 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 pbackground_r background 2D vector of (right) background values used to compute innovations (S, E) - * @param pbackground_R background 2D vector of (RIGHT) background values (S, E) used to compute correlations + * @param bpoints Points of background field (L=num. grid points) + * @param bratios 1D vector (L) representing the ratio of background error standard deviation at grid points to that at observation points. The background at grid points is the value being updated, while the background at station points shares the units and time of the observations, which may differ from the background at grid points. This vector contains coefficients (0-1) that adjust the analysis at grid points, accounting for differences in units and variability between the innovations (observation minus background) and the background at grid points. For example, if trusting the ensemble spread, bratios can be set as the ratio between the ensemble spread of the background to be updated and that used to compute the innovations. If the ensemble spread is not trusted at specific times or grid points, the bratios can be based on a typical expected ratio of spreads from multiple ensemble background realizations. + * @param background 2D vector (L, E) representing the background values at grid points to be updated. + * @param background_corr 2D vector (L, E) representing the background values used to compute dynamic correlations. + * @param obs_points Observation points (S = num. observations) + * @param pobs 2D vector of perturbed observations (S, E) + * @param pratios 1D vector (S) representing the ratio of observation to background error variance. These coefficients (0-1) indicate the relative trust in observations versus the background. A value of 1 means equal trust in both, while values close to 0 indicate greater trust in the observations. For example, a value of 0.1 means the observations are trusted 10 times more than the background. + * @param pbackground 2D vector (S, E) representing the background values at observation points used to compute innovations. + * @param pbackground_corr 2D vector (S, E) representing the background values at observation points used to compute dynamic correlations. * @param structure Structure function for the localization function - * @param var_ratios_or variance_ratio (ratio of observation to right background error variance) - * @param std_ratios_lr standard deviation ratio (ratio of left to right background error standard deviation) - * @param weight given to the analysis increment + * @param bweights 1D vector (L) representing user-defined weights at grid points. The analysis at each grid point is calculated as the background value plus the analysis increment, multiplied by the weight at that point. These weights are useful when iterating over multiple observation times, for instance, allowing the user to prioritize observations that match the time of the original background. For example, if running the function over three observation times (the same as the original background, one hour prior, and two hours prior), setting the weight for the original background time to 0.8 and 0.1 for the other two times ensures that most of the modification comes from the analysis at the original background time. * @param max_points Maximum number of observations to use inside localization zone; Use 0 to disable * @param allow_extrapolation Allow EnSI to extrapolate increments outside increments at observations - * @returns 2D vector of analised values (M, E) + * @returns 2D vector of analised values (L, E) */ - vec2 optimal_interpolation_ensi_lr(const Points& bpoints, - const vec2& background_l, - const vec2& background_L, + vec2 optimal_interpolation_ensi_multi(const Points& bpoints, + const vec& bratios, + const vec2& background, + const vec2& background_corr, const Points& obs_points, - const vec2& obs, - const vec2& pbackground_r, - const vec2& pbackground_R, + const vec2& pobs, + const vec& pratios, + const vec2& pbackground, + const vec2& pbackground_corr, const StructureFunction& structure, - float var_ratios_or, - float std_ratios_lr, - float weight, + const vec& bweights, int max_points, bool allow_extrapolation=true); /** Optimal interpolation for an ensemble gridded field (alternative version) with static correlations. This is an experimental method. - * @param bpoints Points of background field - * @param background_l background 2D vector of (left) background values to update (M, E) M=num. grid points - * @param obs_points Observation points - * @param obs 2D vector of perturbed observations (S, E) S = num. obs points - * @param pbackground_r background 2D vector of (right) background values used to compute innovations (S, E) + * @param bpoints Points of background field (L=num. grid points) + * @param bratios 1D vector (L) representing the ratio of background error standard deviation at grid points to that at observation points. The background at grid points is the value being updated, while the background at station points shares the units and time of the observations, which may differ from the background at grid points. This vector contains coefficients (0-1) that adjust the analysis at grid points, accounting for differences in units and variability between the innovations (observation minus background) and the background at grid points. For example, if trusting the ensemble spread, bratios can be set as the ratio between the ensemble spread of the background to be updated and that used to compute the innovations. If the ensemble spread is not trusted at specific times or grid points, the bratios can be based on a typical expected ratio of spreads from multiple ensemble background realizations. + * @param background 2D vector (L, E) representing the background values at grid points to be updated. + * @param obs_points Observation points (S = num. observations) + * @param pobs 2D vector of perturbed observations (S, E) + * @param pratios 1D vector (S) representing the ratio of observation to background error variance. These coefficients (0-1) indicate the relative trust in observations versus the background. A value of 1 means equal trust in both, while values close to 0 indicate greater trust in the observations. For example, a value of 0.1 means the observations are trusted 10 times more than the background. + * @param pbackground 2D vector (S, E) representing the background values at observation points used to compute innovations. * @param structure Structure function for the static correlations - * @param var_ratios_or variance_ratio (ratio of observation to right background error variance) - * @param std_ratios_lr 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 bweights 1D vector (L) representing user-defined weights at grid points. The analysis at each grid point is calculated as the background value plus the analysis increment, multiplied by the weight at that point. These weights are useful when iterating over multiple observation times, for instance, allowing the user to prioritize observations that match the time of the original background. For example, if running the function over three observation times (the same as the original background, one hour prior, and two hours prior), setting the weight for the original background time to 0.8 and 0.1 for the other two times ensures that most of the modification comes from the analysis at the original background time. + * @param max_points Maximum number of observations to use inside the zone defined by the static correlation; Use 0 to disable * @param allow_extrapolation Allow EnSI to extrapolate increments outside increments at observations - * @returns 2D vector of analised values (M, E) + * @returns 2D vector of analised values (L, E) */ - vec2 optimal_interpolation_ensi_staticcorr_lr(const Points& bpoints, - const vec2& background_l, + vec2 optimal_interpolation_ensi_staticcorr_multi(const Points& bpoints, + const vec& bratios, + const vec2& background, const Points& obs_points, - const vec2& obs, - const vec2& pbackground_r, + const vec2& pobs, + const vec& pratios, + const vec2& pbackground, const StructureFunction& structure, - float var_ratios_or, - float std_ratios_lr, - float weight, + const vec& bweights, int max_points, bool allow_extrapolation=true); /** Optimal interpolation for an ensemble gridded field (alternative version) with ensemble-based correlations. This is an experimental method. (version that works with R bindings) - * @param bpoints Points of background field - * @param background_l background 2D vector of (left) background values to update (M, E) M=num. grid points - * @param background_L 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 pbackground_r background 2D vector of (right) background values used to compute innovations (S, E) - * @param pbackground_R background 2D vector of (RIGHT) background values (S, E) used to compute correlations + * @param bpoints Points of background field (L=num. grid points) + * @param bratios 1D vector (L) representing the ratio of background error standard deviation at grid points to that at observation points. The background at grid points is the value being updated, while the background at station points shares the units and time of the observations, which may differ from the background at grid points. This vector contains coefficients (0-1) that adjust the analysis at grid points, accounting for differences in units and variability between the innovations (observation minus background) and the background at grid points. For example, if trusting the ensemble spread, bratios can be set as the ratio between the ensemble spread of the background to be updated and that used to compute the innovations. If the ensemble spread is not trusted at specific times or grid points, the bratios can be based on a typical expected ratio of spreads from multiple ensemble background realizations. + * @param background 2D vector (L, E) representing the background values at grid points to be updated. + * @param background_corr 2D vector (L, E) representing the background values used to compute dynamic correlations. + * @param obs_points Observation points (S = num. observations) + * @param pobs 2D vector of perturbed observations (S, E) + * @param pratios 1D vector (S) representing the ratio of observation to background error variance. These coefficients (0-1) indicate the relative trust in observations versus the background. A value of 1 means equal trust in both, while values close to 0 indicate greater trust in the observations. For example, a value of 0.1 means the observations are trusted 10 times more than the background. + * @param pbackground 2D vector (S, E) representing the background values at observation points used to compute innovations. + * @param pbackground_corr 2D vector (S, E) representing the background values at observation points used to compute dynamic 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 var_ratios_or variance_ratio (ratio of observation to right background error variance) - * @param std_ratios_lr standard deviation ratio (ratio of left to right background error standard deviation) - * @param weight given to the analysis increment + * @param bweights 1D vector (L) representing user-defined weights at grid points. The analysis at each grid point is calculated as the background value plus the analysis increment, multiplied by the weight at that point. These weights are useful when iterating over multiple observation times, for instance, allowing the user to prioritize observations that match the time of the original background. For example, if running the function over three observation times (the same as the original background, one hour prior, and two hours prior), setting the weight for the original background time to 0.8 and 0.1 for the other two times ensures that most of the modification comes from the analysis at the original background time. * @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 (M, E) + * @param allow_extrapolation Allow EnSI to extrapolate increments outside increments at observations + * @returns 2D vector of analised values (L, E) */ - vec2 R_optimal_interpolation_ensi_lr(const Points& bpoints, - const vec2& background_l, - const vec2& background_L, + vec2 R_optimal_interpolation_ensi_multi(const Points& bpoints, + const vec& bratios, + const vec2& background, + const vec2& background_corr, const Points& obs_points, - const vec2& obs, - const vec2& pbackground_r, - const vec2& pbackground_R, + const vec2& pobs, + const vec& pratios, + const vec2& pbackground, + const vec2& pbackground_corr, /* const StructureFunction& structure, */ int which_structfun, float dh, float dz, float dw, - float var_ratios_or, - float std_ratios_lr, - float weight, + const vec& bweights, int max_points, bool allow_extrapolation=true); /** Optimal interpolation for an ensemble gridded field (alternative version) with static correlations. This is an experimental method. (version that works with R bindings) - * @param bpoints Points of background field - * @param background_l background 2D vector of (left) background values to update (M, E) M=num. grid points - * @param background_L 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 pbackground_r background 2D vector of (right) background values used to compute innovations (S, E) - * @param pbackground_R background 2D vector of (RIGHT) background values (S, E) used to compute correlations + * @param bpoints Points of background field (L=num. grid points) + * @param bratios 1D vector (L) representing the ratio of background error standard deviation at grid points to that at observation points. The background at grid points is the value being updated, while the background at station points shares the units and time of the observations, which may differ from the background at grid points. This vector contains coefficients (0-1) that adjust the analysis at grid points, accounting for differences in units and variability between the innovations (observation minus background) and the background at grid points. For example, if trusting the ensemble spread, bratios can be set as the ratio between the ensemble spread of the background to be updated and that used to compute the innovations. If the ensemble spread is not trusted at specific times or grid points, the bratios can be based on a typical expected ratio of spreads from multiple ensemble background realizations. + * @param background 2D vector (L, E) representing the background values at grid points to be updated. + * @param obs_points Observation points (S = num. observations) + * @param pobs 2D vector of perturbed observations (S, E) + * @param pratios 1D vector (S) representing the ratio of observation to background error variance. These coefficients (0-1) indicate the relative trust in observations versus the background. A value of 1 means equal trust in both, while values close to 0 indicate greater trust in the observations. For example, a value of 0.1 means the observations are trusted 10 times more than the background. + * @param pbackground 2D vector (S, E) representing the background values at observation points used to compute innovations. * @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 var_ratios_or variance_ratio (ratio of observation to right background error variance) - * @param std_ratios_lr 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 2D vector of analised values (M, E) + * @param bweights 1D vector (L) representing user-defined weights at grid points. The analysis at each grid point is calculated as the background value plus the analysis increment, multiplied by the weight at that point. These weights are useful when iterating over multiple observation times, for instance, allowing the user to prioritize observations that match the time of the original background. For example, if running the function over three observation times (the same as the original background, one hour prior, and two hours prior), setting the weight for the original background time to 0.8 and 0.1 for the other two times ensures that most of the modification comes from the analysis at the original background time. + * @param max_points Maximum number of observations to use inside the zone defined by the static correlation; Use 0 to disable + * @param allow_extrapolation Allow EnSI to extrapolate increments outside increments at observations + * @returns 2D vector of analised values (L, E) */ - vec2 R_optimal_interpolation_ensi_staticcorr_lr(const Points& bpoints, - const vec2& background_l, + vec2 R_optimal_interpolation_ensi_staticcorr_multi(const Points& bpoints, + const vec& bratios, + const vec2& background, const Points& obs_points, - const vec2& obs, - const vec2& pbackground_r, + const vec2& pobs, + const vec& pratios, + const vec2& pbackground, /* const StructureFunction& structure, */ int which_structfun, float dh, float dz, float dw, - float var_ratios_or, - float std_ratios_lr, - float weight, + const vec& bweights, int max_points, bool allow_extrapolation=true); @@ -2150,7 +2148,7 @@ namespace gridpp { * @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 linear_rho(float dist, float length, float min) const; float m_localization_distance; }; class MultipleStructure: public StructureFunction { @@ -2202,6 +2200,143 @@ namespace gridpp { float m_min_rho; bool m_is_spatial; }; + /** SOAR structure function based on distance, elevation, and land area fraction */ + class SoarStructure: public StructureFunction { + public: + /** Exponential structure function + * @param h: Horizontal decorrelation length >=0 [m] + * @param v: Vertical decorrelation length >=0 [m]. If 0, disable decorrelation. + * @param w: Land/sea decorrelation length >=0 [1]. If 0, disable decorrelation. + * @param hmax: Truncate horizontal correlation beyond this length [m]. If undefined, 3.64 * h. + */ + SoarStructure(float h, float v=0, float w=0, float hmax=MV); + + /** Soar structure function where decorrelation varyies spatially + * @param grid: Grid of decorrelation field + * @param h: 2D vector of horizontal decorrelation lengths >=0, same size as grid [m] + * @param v: 2D vector of Vertical decorrelation lengths >=0 [m]. Set all to 0 to disable decorrelation. + * @param w: 2D vector of land/sea decorrelation lengths >=0 [1]. Set all to 0 to disable decorrelation. + * @param min_rho: Truncate horizontal correlation when rho is less than this value [m]. + */ + SoarStructure(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; + }; + /** TOAR structure function based on distance, elevation, and land area fraction */ + class ToarStructure: public StructureFunction { + public: + /** Exponential structure function + * @param h: Horizontal decorrelation length >=0 [m] + * @param v: Vertical decorrelation length >=0 [m]. If 0, disable decorrelation. + * @param w: Land/sea decorrelation length >=0 [1]. If 0, disable decorrelation. + * @param hmax: Truncate horizontal correlation beyond this length [m]. If undefined, 3.64 * h. + */ + ToarStructure(float h, float v=0, float w=0, float hmax=MV); + + /** Toar structure function where decorrelation varyies spatially + * @param grid: Grid of decorrelation field + * @param h: 2D vector of horizontal decorrelation lengths >=0, same size as grid [m] + * @param v: 2D vector of Vertical decorrelation lengths >=0 [m]. Set all to 0 to disable decorrelation. + * @param w: 2D vector of land/sea decorrelation lengths >=0 [1]. Set all to 0 to disable decorrelation. + * @param min_rho: Truncate horizontal correlation when rho is less than this value [m]. + */ + ToarStructure(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; + }; + /** Powerlaw structure function based on distance, elevation, and land area fraction */ + class PowerlawStructure: public StructureFunction { + public: + /** Exponential structure function + * @param h: Horizontal decorrelation length >=0 [m] + * @param v: Vertical decorrelation length >=0 [m]. If 0, disable decorrelation. + * @param w: Land/sea decorrelation length >=0 [1]. If 0, disable decorrelation. + * @param hmax: Truncate horizontal correlation beyond this length [m]. If undefined, 3.64 * h. + */ + PowerlawStructure(float h, float v=0, float w=0, float hmax=MV); + + /** Powerlaw structure function where decorrelation varyies spatially + * @param grid: Grid of decorrelation field + * @param h: 2D vector of horizontal decorrelation lengths >=0, same size as grid [m] + * @param v: 2D vector of Vertical decorrelation lengths >=0 [m]. Set all to 0 to disable decorrelation. + * @param w: 2D vector of land/sea decorrelation lengths >=0 [1]. Set all to 0 to disable decorrelation. + * @param min_rho: Truncate horizontal correlation when rho is less than this value [m]. + */ + PowerlawStructure(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; + }; + /** Linear structure function based on distance, elevation, and land area fraction */ + class LinearStructure: public StructureFunction { + public: + /** Exponential structure function + * @param h: Horizontal decorrelation length >=0 [m] + * @param v: Vertical decorrelation length >=0 [m]. If 0, disable decorrelation. + * @param w: Land/sea decorrelation length >=0 [1]. If 0, disable decorrelation. + * @param hmin: Horizontal decorrelation length >=0 [m] + * @param vmin: Vertical decorrelation length >=0 [m]. If 0, disable decorrelation. + * @param wmin: Land/sea decorrelation length >=0 [1]. If 0, disable decorrelation. + * @param hmax: Truncate horizontal correlation beyond this length [m]. If undefined, 3.64 * h. + */ + LinearStructure(float h, float v=0, float w=0, float hmin=0, float vmin=0, float wmin=0, float hmax=MV); + + /** Linear structure function where decorrelation varyies spatially + * @param grid: Grid of decorrelation field + * @param h: 2D vector of horizontal decorrelation lengths >=0, same size as grid [m] + * @param v: 2D vector of Vertical decorrelation lengths >=0 [m]. Set all to 0 to disable decorrelation. + * @param w: 2D vector of land/sea decorrelation lengths >=0 [1]. Set all to 0 to disable decorrelation. + * @param hmin: 2D vector of horizontal decorrelation lengths >=0, same size as grid [m] + * @param vmin: 2D vector of Vertical decorrelation lengths >=0 [m]. Set all to 0 to disable decorrelation. + * @param wmin: 2D vector of land/sea decorrelation lengths >=0 [1]. Set all to 0 to disable decorrelation. + * @param min_rho: Truncate horizontal correlation when rho is less than this value [m]. + */ + LinearStructure(Grid grid, vec2 h, vec2 v, vec2 w, vec2 hmin, vec2 vmin, vec2 wmin, 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, float hmin) const; + Grid m_grid; + vec2 mH; + vec2 mV; + vec2 mW; + vec2 mHmin; + vec2 mVmin; + vec2 mWmin; + float m_min_rho; + bool m_is_spatial; + }; /** Mix of structure functions based on distance(SOAR), elevation(Gauss), and land area fraction(linear) */ class MixAStructure: public StructureFunction { diff --git a/src/api/oi_ensi_lr.cpp b/src/api/oi_ensi_multi.cpp similarity index 85% rename from src/api/oi_ensi_lr.cpp rename to src/api/oi_ensi_multi.cpp index 047fa7a..85fe107 100644 --- a/src/api/oi_ensi_lr.cpp +++ b/src/api/oi_ensi_multi.cpp @@ -30,17 +30,17 @@ void print_matrix(Matrix matrix) { template void print_matrix< ::mattype>(::mattype matrix); template void print_matrix< ::cxtype>(::cxtype matrix); -vec3 gridpp::optimal_interpolation_ensi_lr(const gridpp::Grid& bgrid, - const vec3& background_l, - const vec3& background_L, +vec3 gridpp::optimal_interpolation_ensi_multi(const gridpp::Grid& bgrid, + const vec2& bratios, + const vec3& background, + const vec3& background_corr, const gridpp::Points& points, const vec2& pobs, - const vec2& pbackground_r, - const vec2& pbackground_R, + const vec& pratios, + const vec2& pbackground, + const vec2& pbackground_corr, const gridpp::StructureFunction& structure, - float var_ratios_or, - float std_ratios_lr, - float weight, + const vec2& bweights, int max_points, bool dynamic_correlations, bool allow_extrapolation) { @@ -52,7 +52,7 @@ vec3 gridpp::optimal_interpolation_ensi_lr(const gridpp::Grid& bgrid, int nS = points.size(); if(nS == 0) - return background_l; + return background; int nY = bgrid.size()[0]; int nX = bgrid.size()[1]; @@ -67,54 +67,71 @@ vec3 gridpp::optimal_interpolation_ensi_lr(const gridpp::Grid& bgrid, throw std::invalid_argument("Both background grid and observations points must be of same coordinate type (lat/lon or x/y)"); } // Check ensembles have consistent sizes - int nE = background_l[0][0].size(); - if(background_l.size() != nY || background_l[0].size() != nX) { + int nE = background[0][0].size(); + if(background.size() != nY || background[0].size() != nX) { std::stringstream ss; - ss << "Input left field (" << background_l.size() << "," << background_l[0].size() << "," << background_l[0][0].size() << ") is not the same size as the grid (" << nY << "," << nX << "," << nE << ")"; + ss << "Input left field (" << background.size() << "," << background[0].size() << "," << background[0][0].size() << ") is not the same size as the grid (" << nY << "," << nX << "," << nE << ")"; throw std::invalid_argument(ss.str()); } - if(background_L.size() != nY || background_L[0].size() != nX || background_L[0][0].size() != nE) { + if(background_corr.size() != nY || background_corr[0].size() != nX || background_corr[0][0].size() != nE) { std::stringstream ss; - ss << "Input LEFT field (" << background_L.size() << "," << background_L[0].size() << "," << background_L[0][0].size() << ") is not the same size as the grid (" << nY << "," << nX << "," << nE << ")"; + ss << "Input LEFT field (" << background_corr.size() << "," << background_corr[0].size() << "," << background_corr[0][0].size() << ") is not the same size as the grid (" << nY << "," << nX << "," << nE << ")"; throw std::invalid_argument(ss.str()); } - if(pbackground_r.size() != nS || pbackground_r[0].size() != nE) { + if(bratios.size() != nY || bratios[0].size() != nX) { std::stringstream ss; - ss << "Input right field at observation location (" << pbackground_r.size() << "," << pbackground_r[0].size() << ") and points (" << nS << "," << nE << ") size mismatch"; + ss << "Input bratios field (" << bratios.size() << "," << bratios[0].size() << ") is not the same size as the grid (" << nY << "," << nX << ")"; throw std::invalid_argument(ss.str()); } - - if(pbackground_R.size() != nS || pbackground_R[0].size() != nE) { + if(bweights.size() != nY || bweights[0].size() != nX) { + std::stringstream ss; + ss << "Input bweights field (" << bweights.size() << "," << bweights[0].size() << ") is not the same size as the grid (" << nY << "," << nX << ")"; + throw std::invalid_argument(ss.str()); + } + if(pbackground.size() != nS || pbackground[0].size() != nE) { std::stringstream ss; - ss << "Input RIGHT field at observation location (" << pbackground_R.size() << "," << pbackground_R[0].size() << ") and points (" << nS << "," << nE << ") size mismatch"; + ss << "Input right field at observation location (" << pbackground.size() << "," << pbackground[0].size() << ") and points (" << nS << "," << nE << ") size mismatch"; + throw std::invalid_argument(ss.str()); + } + if(pbackground_corr.size() != nS || pbackground_corr[0].size() != nE) { + std::stringstream ss; + ss << "Input RIGHT field at observation location (" << pbackground_corr.size() << "," << pbackground_corr[0].size() << ") and points (" << nS << "," << nE << ") size mismatch"; throw std::invalid_argument(ss.str()); } - // Check observations have consistent size if(pobs.size() != nS || pobs[0].size() != nE) { std::stringstream ss; ss << "Observations (" << pobs.size() << "," << pobs[0].size() << ") and points (" << nS << "," << nE << ") size mismatch"; throw std::invalid_argument(ss.str()); } - + if(pratios.size() != nS) { + std::stringstream ss; + ss << "Ratios (" << pratios.size() << ") and points (" << nS << ") size mismatch"; + throw std::invalid_argument(ss.str()); + } + gridpp::Points bpoints = bgrid.to_points(); - vec2 background_l1 = gridpp::init_vec2(nY * nX, nE); - vec2 background_L1 = gridpp::init_vec2(nY * nX, nE); + vec2 background1 = gridpp::init_vec2(nY * nX, nE); + vec2 background_corr1 = gridpp::init_vec2(nY * nX, nE); + vec bratios1(nY * nX); + vec bweights1(nY * nX); int count = 0; for(int y = 0; y < nY; y++) { for(int x = 0; x < nX; x++) { + bratios1[count] = bratios[y][x]; + bweights1[count] = bweights[y][x]; for(int e = 0; e < nE; e++) { - background_l1[count][e] = background_l[y][x][e]; - background_L1[count][e] = background_L[y][x][e]; + background1[count][e] = background[y][x][e]; + background_corr1[count][e] = background_corr[y][x][e]; } count++; } } vec2 output1 = gridpp::init_vec2(nY * nX, nE); if(dynamic_correlations) { - output1 = optimal_interpolation_ensi_lr(bpoints, background_l1, background_L1, points, pobs, pbackground_r, pbackground_R, structure, var_ratios_or, std_ratios_lr, weight, max_points, allow_extrapolation); + output1 = optimal_interpolation_ensi_multi(bpoints, bratios1, background1, background_corr1, points, pobs, pratios, pbackground, pbackground_corr, structure, bweights1, max_points, allow_extrapolation); } else { - output1 = optimal_interpolation_ensi_staticcorr_lr(bpoints, background_l1, points, pobs, pbackground_r, structure, var_ratios_or, std_ratios_lr, weight, max_points, allow_extrapolation); + output1 = optimal_interpolation_ensi_staticcorr_multi(bpoints, bratios1, background1, points, pobs, pratios, pbackground, structure, bweights1, max_points, allow_extrapolation); } vec3 output = gridpp::init_vec3(nY, nX, nE); count = 0; @@ -129,17 +146,17 @@ vec3 gridpp::optimal_interpolation_ensi_lr(const gridpp::Grid& bgrid, return output; } -vec2 gridpp::optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, - const vec2& background_l, - const vec2& background_L, +vec2 gridpp::optimal_interpolation_ensi_multi(const gridpp::Points& bpoints, + const vec& bratios, + const vec2& background, + const vec2& background_corr, const gridpp::Points& points, - const vec2& pobs, - const vec2& pbackground_r, - const vec2& pbackground_R, + const vec2& pobs, + const vec& pratios, + const vec2& pbackground, + const vec2& pbackground_corr, const gridpp::StructureFunction& structure, - float var_ratios_or, - float std_ratios_lr, - float weight, + const vec& bweights, int max_points, bool allow_extrapolation) { if(max_points < 0) @@ -147,32 +164,32 @@ vec2 gridpp::optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, 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_l.size() != bpoints.size()) + if(background.size() != bpoints.size()) throw std::invalid_argument("Input left field is not the same size as the grid"); - if(background_L.size() != bpoints.size()) + if(background_corr.size() != bpoints.size()) throw std::invalid_argument("Input LEFT field is not the same size as the grid"); if(pobs.size() != points.size()) throw std::invalid_argument("Observations and points exception mismatch"); - if(pbackground_r.size() != points.size()) + if(pbackground.size() != points.size()) throw std::invalid_argument("Background rigth and points size mismatch"); - if(pbackground_R.size() != points.size()) + if(pbackground_corr.size() != points.size()) throw std::invalid_argument("Background RIGTH and points size mismatch"); float default_min_std = 0.0013; int nS = points.size(); if(nS == 0) - return background_l; + return background; int mY = -1; // Write debug information for this station index bool diagnose = false; - int nY = background_l.size(); - int nEns = background_l[0].size(); + int nY = background.size(); + int nEns = background[0].size(); // Prepare output matrix float missing_value = -99999.999; /* vec2 output = gridpp::init_vec2(nY, nEns, missing_value); */ - vec2 output = background_l; + vec2 output = background; vec blats = bpoints.get_lats(); vec blons = bpoints.get_lons(); @@ -199,14 +216,14 @@ vec2 gridpp::optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, for(int e = 0; e < nEns; e++) { int numInvalid = 0; for(int y = 0; y < nY; y++) { - float value_l = background_l[y][e]; - float value_L = background_L[y][e]; + float value_l = background[y][e]; + float value_L = background_corr[y][e]; if(!gridpp::is_valid(value_l) || !gridpp::is_valid(value_L)) numInvalid++; } for(int i = 0; i < nS; i++) { - float value_r = pbackground_r[i][e]; - float value_R = pbackground_R[i][e]; + float value_r = pbackground[i][e]; + float value_R = pbackground_corr[i][e]; if(!gridpp::is_valid(value_r) || !gridpp::is_valid(value_R)) numInvalid++; } @@ -216,7 +233,7 @@ vec2 gridpp::optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, } } if(nValidEns == 0) - return background_l; + return background; // gZ_R(nY, nValidEns): used to compute ensemble-based background correlations i) between yth gridpoint and observations ii) among observations vec2 gZ_R = gridpp::init_vec2(nY, nValidEns); // useful to compute dynamical correlations @@ -224,7 +241,7 @@ vec2 gridpp::optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, vec pbackgroundValid_R(nValidEns); for(int e = 0; e < nValidEns; e++) { int ei = validEns[e]; - pbackgroundValid_R[e] = pbackground_R[i][ei]; + pbackgroundValid_R[e] = pbackground_corr[i][ei]; } float mean = gridpp::calc_statistic(pbackgroundValid_R, gridpp::Mean); float std = gridpp::calc_statistic(pbackgroundValid_R, gridpp::Std); @@ -249,6 +266,8 @@ vec2 gridpp::optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, // but it doesnt seem to help // #pragma omp parallel for for(int y = 0; y < nY; y++) { + float std_ratios_lr = bratios[y]; + float weight = bweights[y]; float lat = blats[y]; float lon = blons[y]; float elev = belevs[y]; @@ -325,7 +344,7 @@ vec2 gridpp::optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, vec backgroundValid_L(nValidEns); for(int e = 0; e < nValidEns; e++) { int ei = validEns[e]; - backgroundValid_L[e] = background_L[y][ei]; + backgroundValid_L[e] = background_corr[y][ei]; } float mean = gridpp::calc_statistic(backgroundValid_L, gridpp::Mean); float std = gridpp::calc_statistic(backgroundValid_L, gridpp::Std); @@ -335,7 +354,7 @@ vec2 gridpp::optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, } // lZ_R: used to compute ensemble-based background correlations i) between yth gridpoint and observations ii) among observations mattype lZ_R(lS, nValidEns); - // Innovation: Observation - Background_r + // Innovation: Observation - background mattype lInnov(lS, nValidEns, arma::fill::zeros); // lR_dd: Observation error correlation matrix mattype lR_dd(lS, lS, arma::fill::zeros); @@ -344,16 +363,16 @@ vec2 gridpp::optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, // lLoc2D: localization for ensemble-based background correlations among observations mattype lLoc2D(lS, lS, arma::fill::zeros); for(int i = 0; i < lS; i++) { + int index = lLocIndices[i]; // lR_dd: diagonal observation error correlation matrix - lR_dd(i, i) = var_ratios_or; + lR_dd(i, i) = pratios[index]; // lLoc1D between yth gridpoint and ith observation computed before lLoc1D(0, i) = lRhos(i); // compute lZ_R and lInnov - int index = lLocIndices[i]; for(int e = 0; e < nValidEns; e++) { lZ_R(i, e) = gZ_R[index][e]; int ei = validEns[e]; - lInnov(i, ei) = pobs[index][ei] - pbackground_r[index][ei]; + lInnov(i, ei) = pobs[index][ei] - pbackground[index][ei]; } // compute lLoc2D Point p1 = point_vec[index]; @@ -405,35 +424,35 @@ vec2 gridpp::optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, } } - // Compute the analysis as the updated background_l + // Compute the analysis as the updated background for(int e = 0; e < nValidEns; e++) { int ei = validEns[e]; - output[y][ei] = background_l[y][ei] + dx[e]; + output[y][ei] = background[y][ei] + dx[e]; } // debug /* for(int i = 0; i < lS; i++) { // compute lZ_R and lInnov int index = lLocIndices[i]; - std::cout << i << " backg_r obs " << pbackground_r[index][0] << " " << pobs[index][0] << std::endl; + std::cout << i << " backg_r obs " << pbackground[index][0] << " " << pobs[index][0] << std::endl; } // end loop over closer observations for(int e = 0; e < nValidEns; e++) { int ei = validEns[e]; - std::cout << ei << " backg_l analysis " << background_l[y][ei] << " " << output[y][ei] << std::endl; + std::cout << ei << " backg_l analysis " << background[y][ei] << " " << output[y][ei] << std::endl; } */ } // end loop over gridpoint return output; -} // end of optimal_interpolation_ensi_lr +} // end of optimal_interpolation_ensi_multi -vec2 gridpp::optimal_interpolation_ensi_staticcorr_lr(const gridpp::Points& bpoints, - const vec2& background_l, +vec2 gridpp::optimal_interpolation_ensi_staticcorr_multi(const gridpp::Points& bpoints, + const vec& bratios, + const vec2& background, const gridpp::Points& points, const vec2& pobs, - const vec2& pbackground_r, + const vec& pratios, + const vec2& pbackground, const gridpp::StructureFunction& structure, - float var_ratios_or, - float std_ratios_lr, - float weight, + const vec& bweights, int max_points, bool allow_extrapolation) { if(max_points < 0) @@ -441,26 +460,26 @@ vec2 gridpp::optimal_interpolation_ensi_staticcorr_lr(const gridpp::Points& bpoi 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_l.size() != bpoints.size()) + if(background.size() != bpoints.size()) throw std::invalid_argument("Input left field is not the same size as the grid"); if(pobs.size() != points.size()) throw std::invalid_argument("Observations and points exception mismatch"); - if(pbackground_r.size() != points.size()) + if(pbackground.size() != points.size()) throw std::invalid_argument("Background rigth and points size mismatch"); int nS = points.size(); if(nS == 0) - return background_l; + return background; int mY = -1; // Write debug information for this station index bool diagnose = false; - int nY = background_l.size(); - int nEns = background_l[0].size(); + int nY = background.size(); + int nEns = background[0].size(); // Prepare output matrix float missing_value = -99999.999; - vec2 output = background_l; + vec2 output = background; vec blats = bpoints.get_lats(); vec blons = bpoints.get_lons(); @@ -486,12 +505,12 @@ vec2 gridpp::optimal_interpolation_ensi_staticcorr_lr(const gridpp::Points& bpoi for(int e = 0; e < nEns; e++) { int numInvalid = 0; for(int y = 0; y < nY; y++) { - float value_l = background_l[y][e]; + float value_l = background[y][e]; if(!gridpp::is_valid(value_l)) numInvalid++; } for(int i = 0; i < nS; i++) { - float value_r = pbackground_r[i][e]; + float value_r = pbackground[i][e]; if(!gridpp::is_valid(value_r)) numInvalid++; } @@ -501,13 +520,15 @@ vec2 gridpp::optimal_interpolation_ensi_staticcorr_lr(const gridpp::Points& bpoi } } if(nValidEns == 0) - return background_l; + return background; // 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 std_ratios_lr = bratios[y]; + float weight = bweights[y]; float lat = blats[y]; float lon = blons[y]; float elev = belevs[y]; @@ -580,7 +601,7 @@ vec2 gridpp::optimal_interpolation_ensi_staticcorr_lr(const gridpp::Points& bpoi lLafs[i] = plafs[index]; } - // Innovation: Observation - Background_r + // Innovation: Observation - background mattype lInnov(lS, nValidEns, arma::fill::zeros); // lR_dd: Observation error correlation matrix mattype lR_dd(lS, lS, arma::fill::zeros); @@ -589,15 +610,15 @@ vec2 gridpp::optimal_interpolation_ensi_staticcorr_lr(const gridpp::Points& bpoi // lCorr2D: background correlations among observations mattype lCorr2D(lS, lS, arma::fill::zeros); for(int i = 0; i < lS; i++) { + int index = lLocIndices[i]; // lR_dd: diagonal observation error correlation matrix - lR_dd(i, i) = var_ratios_or; + lR_dd(i, i) = pratios[index]; // lCorr1D between yth gridpoint and ith observation computed before lCorr1D(0, i) = lRhos(i); // compute lInnov - int index = lLocIndices[i]; for(int e = 0; e < nValidEns; e++) { int ei = validEns[e]; - lInnov(i, ei) = pobs[index][ei] - pbackground_r[index][ei]; + lInnov(i, ei) = pobs[index][ei] - pbackground[index][ei]; } // compute lCorr2D Point p1 = point_vec[index]; @@ -643,35 +664,35 @@ vec2 gridpp::optimal_interpolation_ensi_staticcorr_lr(const gridpp::Points& bpoi } } - // Compute the analysis as the updated background_l + // Compute the analysis as the updated background for(int e = 0; e < nValidEns; e++) { int ei = validEns[e]; - output[y][ei] = background_l[y][ei] + dx[e]; + output[y][ei] = background[y][ei] + dx[e]; } } // end loop over gridpoint return output; -} // end of R_optimal_interpolation_ensi_staticcorr_lr +} // end of R_optimal_interpolation_ensi_staticcorr_multi //---------------------------------------------------------------------------- // R-specifics //---------------------------------------------------------------------------- /* command to fix the gridpp.R file: for i in {1..50}; do sed -i "s/all(sapply(argv\[\[$i\]\] , is.integer) || sapply(argv\[\[$i\]\], is.numeric))/all(sapply(argv\[\[$i\]\] , is.integer) | sapply(argv\[\[$i\]\], is.numeric))/g" gridpp.R; done */ -vec2 gridpp::R_optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, - const vec2& background_l, - const vec2& background_L, +vec2 gridpp::R_optimal_interpolation_ensi_multi(const gridpp::Points& bpoints, + const vec& bratios, + const vec2& background, + const vec2& background_corr, const gridpp::Points& points, const vec2& pobs, - const vec2& pbackground_r, - const vec2& pbackground_R, + const vec& pratios, + const vec2& pbackground, + const vec2& pbackground_corr, /* const gridpp::StructureFunction& structure, it creates problem with R bindings */ int which_structfun, float dh, float dz, float dw, - float var_ratios_or, - float std_ratios_lr, - float weight, + const vec& bweights, int max_points, bool allow_extrapolation) { if(max_points < 0) @@ -679,15 +700,15 @@ vec2 gridpp::R_optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, 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_l.size() != bpoints.size()) + if(background.size() != bpoints.size()) throw std::invalid_argument("Input left field is not the same size as the grid"); - if(background_L.size() != bpoints.size()) + if(background_corr.size() != bpoints.size()) throw std::invalid_argument("Input LEFT field is not the same size as the grid"); if(pobs.size() != points.size()) throw std::invalid_argument("Observations and points exception mismatch"); - if(pbackground_r.size() != points.size()) + if(pbackground.size() != points.size()) throw std::invalid_argument("Background rigth and points size mismatch"); - if(pbackground_R.size() != points.size()) + if(pbackground_corr.size() != points.size()) throw std::invalid_argument("Background RIGTH and points size mismatch"); float hmax = 7 * dh; @@ -706,18 +727,18 @@ vec2 gridpp::R_optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, int nS = points.size(); if(nS == 0) - return background_l; + return background; int mY = -1; // Write debug information for this station index bool diagnose = false; - int nY = background_l.size(); - int nEns = background_l[0].size(); + int nY = background.size(); + int nEns = background[0].size(); // Prepare output matrix float missing_value = -99999.999; /* vec2 output = gridpp::init_vec2(nY, nEns, missing_value); */ - vec2 output = background_l; + vec2 output = background; vec blats = bpoints.get_lats(); vec blons = bpoints.get_lons(); @@ -744,14 +765,14 @@ vec2 gridpp::R_optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, for(int e = 0; e < nEns; e++) { int numInvalid = 0; for(int y = 0; y < nY; y++) { - float value_l = background_l[y][e]; - float value_L = background_L[y][e]; + float value_l = background[y][e]; + float value_L = background_corr[y][e]; if(!gridpp::is_valid(value_l) || !gridpp::is_valid(value_L)) numInvalid++; } for(int i = 0; i < nS; i++) { - float value_r = pbackground_r[i][e]; - float value_R = pbackground_R[i][e]; + float value_r = pbackground[i][e]; + float value_R = pbackground_corr[i][e]; if(!gridpp::is_valid(value_r) || !gridpp::is_valid(value_R)) numInvalid++; } @@ -761,7 +782,7 @@ vec2 gridpp::R_optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, } } if(nValidEns == 0) - return background_l; + return background; // gZ_R(nY, nValidEns): used to compute ensemble-based background correlations i) between yth gridpoint and observations ii) among observations vec2 gZ_R = gridpp::init_vec2(nY, nValidEns); // useful to compute dynamical correlations @@ -769,7 +790,7 @@ vec2 gridpp::R_optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, vec pbackgroundValid_R(nValidEns); for(int e = 0; e < nValidEns; e++) { int ei = validEns[e]; - pbackgroundValid_R[e] = pbackground_R[i][ei]; + pbackgroundValid_R[e] = pbackground_corr[i][ei]; } float mean = gridpp::calc_statistic(pbackgroundValid_R, gridpp::Mean); float std = gridpp::calc_statistic(pbackgroundValid_R, gridpp::Std); @@ -794,6 +815,8 @@ vec2 gridpp::R_optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, // but it doesnt seem to help // #pragma omp parallel for for(int y = 0; y < nY; y++) { + float std_ratios_lr = bratios[y]; + float weight = bweights[y]; float lat = blats[y]; float lon = blons[y]; float elev = belevs[y]; @@ -893,7 +916,7 @@ vec2 gridpp::R_optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, vec backgroundValid_L(nValidEns); for(int e = 0; e < nValidEns; e++) { int ei = validEns[e]; - backgroundValid_L[e] = background_L[y][ei]; + backgroundValid_L[e] = background_corr[y][ei]; } float mean = gridpp::calc_statistic(backgroundValid_L, gridpp::Mean); float std = gridpp::calc_statistic(backgroundValid_L, gridpp::Std); @@ -903,7 +926,7 @@ vec2 gridpp::R_optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, } // lZ_R: used to compute ensemble-based background correlations i) between yth gridpoint and observations ii) among observations mattype lZ_R(lS, nValidEns); - // Innovation: Observation - Background_r + // Innovation: Observation - background mattype lInnov(lS, nValidEns, arma::fill::zeros); // lR_dd: Observation error correlation matrix mattype lR_dd(lS, lS, arma::fill::zeros); @@ -912,16 +935,16 @@ vec2 gridpp::R_optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, // lLoc2D: localization for ensemble-based background correlations among observations mattype lLoc2D(lS, lS, arma::fill::zeros); for(int i = 0; i < lS; i++) { + int index = lLocIndices[i]; // lR_dd: diagonal observation error correlation matrix - lR_dd(i, i) = var_ratios_or; + lR_dd(i, i) = pratios[index]; // lLoc1D between yth gridpoint and ith observation computed before lLoc1D(0, i) = lRhos(i); // compute lZ_R and lInnov - int index = lLocIndices[i]; for(int e = 0; e < nValidEns; e++) { lZ_R(i, e) = gZ_R[index][e]; int ei = validEns[e]; - lInnov(i, ei) = pobs[index][ei] - pbackground_r[index][ei]; + lInnov(i, ei) = pobs[index][ei] - pbackground[index][ei]; } // compute lLoc2D Point p1 = point_vec[index]; @@ -984,41 +1007,41 @@ vec2 gridpp::R_optimal_interpolation_ensi_lr(const gridpp::Points& bpoints, } } - // Compute the analysis as the updated background_l + // Compute the analysis as the updated background for(int e = 0; e < nValidEns; e++) { int ei = validEns[e]; - output[y][ei] = background_l[y][ei] + dx[e]; + output[y][ei] = background[y][ei] + dx[e]; } // debug /* for(int i = 0; i < lS; i++) { // compute lZ_R and lInnov int index = lLocIndices[i]; - std::cout << i << " backg_r obs " << pbackground_r[index][0] << " " << pobs[index][0] << std::endl; + std::cout << i << " backg_r obs " << pbackground[index][0] << " " << pobs[index][0] << std::endl; } // end loop over closer observations for(int e = 0; e < nValidEns; e++) { int ei = validEns[e]; - std::cout << ei << " backg_l analysis " << background_l[y][ei] << " " << output[y][ei] << std::endl; + std::cout << ei << " backg_l analysis " << background[y][ei] << " " << output[y][ei] << std::endl; } */ } // end loop over gridpoint return output; -} // end of R_optimal_interpolation_ensi_lr +} // end of R_optimal_interpolation_ensi_multi /* command to fix the gridpp.R file: for i in {1..50}; do sed -i "s/all(sapply(argv\[\[$i\]\] , is.integer) || sapply(argv\[\[$i\]\], is.numeric))/all(sapply(argv\[\[$i\]\] , is.integer) | sapply(argv\[\[$i\]\], is.numeric))/g" gridpp.R; done */ -vec2 gridpp::R_optimal_interpolation_ensi_staticcorr_lr(const gridpp::Points& bpoints, - const vec2& background_l, +vec2 gridpp::R_optimal_interpolation_ensi_staticcorr_multi(const gridpp::Points& bpoints, + const vec& bratios, + const vec2& background, const gridpp::Points& points, const vec2& pobs, - const vec2& pbackground_r, + const vec& pratios, + const vec2& pbackground, /* const gridpp::StructureFunction& structure, it creates problem with R bindings */ int which_structfun, float dh, float dz, float dw, - float var_ratios_or, - float std_ratios_lr, - float weight, + const vec& bweights, int max_points, bool allow_extrapolation) { if(max_points < 0) @@ -1026,11 +1049,11 @@ vec2 gridpp::R_optimal_interpolation_ensi_staticcorr_lr(const gridpp::Points& bp 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_l.size() != bpoints.size()) + if(background.size() != bpoints.size()) throw std::invalid_argument("Input left field is not the same size as the grid"); if(pobs.size() != points.size()) throw std::invalid_argument("Observations and points exception mismatch"); - if(pbackground_r.size() != points.size()) + if(pbackground.size() != points.size()) throw std::invalid_argument("Background rigth and points size mismatch"); float hmax = 7 * dh; @@ -1038,17 +1061,17 @@ vec2 gridpp::R_optimal_interpolation_ensi_staticcorr_lr(const gridpp::Points& bp int nS = points.size(); if(nS == 0) - return background_l; + return background; int mY = -1; // Write debug information for this station index bool diagnose = false; - int nY = background_l.size(); - int nEns = background_l[0].size(); + int nY = background.size(); + int nEns = background[0].size(); // Prepare output matrix float missing_value = -99999.999; - vec2 output = background_l; + vec2 output = background; vec blats = bpoints.get_lats(); vec blons = bpoints.get_lons(); @@ -1074,12 +1097,12 @@ vec2 gridpp::R_optimal_interpolation_ensi_staticcorr_lr(const gridpp::Points& bp for(int e = 0; e < nEns; e++) { int numInvalid = 0; for(int y = 0; y < nY; y++) { - float value_l = background_l[y][e]; + float value_l = background[y][e]; if(!gridpp::is_valid(value_l)) numInvalid++; } for(int i = 0; i < nS; i++) { - float value_r = pbackground_r[i][e]; + float value_r = pbackground[i][e]; if(!gridpp::is_valid(value_r)) numInvalid++; } @@ -1089,13 +1112,15 @@ vec2 gridpp::R_optimal_interpolation_ensi_staticcorr_lr(const gridpp::Points& bp } } if(nValidEns == 0) - return background_l; + return background; // 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 std_ratios_lr = bratios[y]; + float weight = bweights[y]; float lat = blats[y]; float lon = blons[y]; float elev = belevs[y]; @@ -1189,7 +1214,7 @@ vec2 gridpp::R_optimal_interpolation_ensi_staticcorr_lr(const gridpp::Points& bp lLafs[i] = plafs[index]; } - // Innovation: Observation - Background_r + // Innovation: Observation - background mattype lInnov(lS, nValidEns, arma::fill::zeros); // lR_dd: Observation error correlation matrix mattype lR_dd(lS, lS, arma::fill::zeros); @@ -1198,15 +1223,15 @@ vec2 gridpp::R_optimal_interpolation_ensi_staticcorr_lr(const gridpp::Points& bp // lCorr2D: background correlations among observations mattype lCorr2D(lS, lS, arma::fill::zeros); for(int i = 0; i < lS; i++) { + int index = lLocIndices[i]; // lR_dd: diagonal observation error correlation matrix - lR_dd(i, i) = var_ratios_or; + lR_dd(i, i) = pratios[index]; // lCorr1D between yth gridpoint and ith observation computed before lCorr1D(0, i) = lRhos(i); // compute lInnov - int index = lLocIndices[i]; for(int e = 0; e < nValidEns; e++) { int ei = validEns[e]; - lInnov(i, ei) = pobs[index][ei] - pbackground_r[index][ei]; + lInnov(i, ei) = pobs[index][ei] - pbackground[index][ei]; } // compute lCorr2D Point p1 = point_vec[index]; @@ -1263,11 +1288,11 @@ vec2 gridpp::R_optimal_interpolation_ensi_staticcorr_lr(const gridpp::Points& bp } } - // Compute the analysis as the updated background_l + // Compute the analysis as the updated background for(int e = 0; e < nValidEns; e++) { int ei = validEns[e]; - output[y][ei] = background_l[y][ei] + dx[e]; + output[y][ei] = background[y][ei] + dx[e]; } } // end loop over gridpoint return output; -} // end of R_optimal_interpolation_ensi_staticcorr_lr +} // end of R_optimal_interpolation_ensi_staticcorr_multi diff --git a/src/api/structure.cpp b/src/api/structure.cpp index 4a7db60..cad5c18 100644 --- a/src/api/structure.cpp +++ b/src/api/structure.cpp @@ -73,13 +73,16 @@ float gridpp::StructureFunction::powerlaw_rho(float dist, float length) const { 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) +float gridpp::StructureFunction::linear_rho(float dist, float length, float min) const { + if(!gridpp::is_valid(min) || min < 0 || !gridpp::is_valid(length) || length < 0) // Disabled return 1; - if(!gridpp::is_valid(normdist)) + if(!gridpp::is_valid(dist)) return 0; - return (1 - (1-min) * abs(normdist)); + float normdist = abs(dist) / length; + if(normdist > 1) + normdist = 1; + return (1 - (1-min) * normdist); } float gridpp::StructureFunction::localization_distance(const Point& p) const { return m_localization_distance; @@ -134,7 +137,9 @@ gridpp::StructureFunctionPtr gridpp::MultipleStructure::clone() const { return std::make_shared(*m_structure_h, *m_structure_v, *m_structure_w); } +//============================================================================= /** Barnes */ +//============================================================================= gridpp::BarnesStructure::BarnesStructure(float h, float v, float w, float hmax) : m_is_spatial(false) { if(gridpp::is_valid(hmax) && hmax < 0) @@ -276,7 +281,9 @@ float gridpp::BarnesStructure::localization_distance(float h) const { return sqrt(-2*log(m_min_rho)) * h; } +//============================================================================= /** Cressman */ +//============================================================================= gridpp::CressmanStructure::CressmanStructure(float h, float v, float w) : gridpp::StructureFunction(h) { if(!gridpp::is_valid(v) || v < 0) @@ -303,7 +310,9 @@ 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) @@ -373,7 +382,7 @@ float gridpp::MixAStructure::corr(const Point& p1, const Point& p2) const { } if(gridpp::is_valid(p1.laf) && gridpp::is_valid(p2.laf)) { float lafdist = p1.laf - p2.laf; - rho *= gridpp::StructureFunction::linear_rho(lafdist, w); + rho *= gridpp::StructureFunction::linear_rho(lafdist, 1., w); } } else { @@ -387,7 +396,7 @@ float gridpp::MixAStructure::corr(const Point& p1, const Point& p2) const { } 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]); + rho *= gridpp::StructureFunction::linear_rho(lafdist, 1., mW[0][0]); } } return rho; @@ -416,7 +425,7 @@ vec gridpp::MixAStructure::corr(const Point& p1, const std::vector& p2) c } 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); + rho *= gridpp::StructureFunction::linear_rho(lafdist, 1., w); } } output[i] = rho; @@ -445,6 +454,631 @@ float gridpp::MixAStructure::localization_distance(float h) const { return sqrt(-2*log(m_min_rho)) * h; } +//============================================================================= +/** SOAR */ +//============================================================================= +gridpp::SoarStructure::SoarStructure(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::SoarStructure::SoarStructure(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::SoarStructure::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::soar_rho(vdist, v); + } + if(gridpp::is_valid(p1.laf) && gridpp::is_valid(p2.laf)) { + float lafdist = p1.laf - p2.laf; + rho *= gridpp::StructureFunction::soar_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::soar_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::soar_rho(lafdist, mW[0][0]); + } + } + return rho; +} +vec gridpp::SoarStructure::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::soar_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::soar_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::SoarStructure::clone() const { + return std::make_shared(m_grid, mH, mV, mW, m_min_rho); +} +float gridpp::SoarStructure::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::SoarStructure::localization_distance(float h) const { +// To find an approximate analytical solution for the inverse of the function y=(1+x)exp(-x) when y is between 0 and 1, we can follow a similar approach used in approximations for transcendental functions. + if(m_min_rho > 0.5) { + return sqrt( 2 * (1-m_min_rho)) * h; + } + else { + return ( -log(m_min_rho) - log(log(m_min_rho))) * h; + } +} +//----------------------------------------------------------------------------- +/** End SOAR */ +//----------------------------------------------------------------------------- + +//============================================================================= +/** TOAR */ +//============================================================================= +gridpp::ToarStructure::ToarStructure(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 + pow(hmax / h, 2)/3) * 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::ToarStructure::ToarStructure(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::ToarStructure::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::toar_rho(hdist, h); + if(gridpp::is_valid(p1.elev) && gridpp::is_valid(p2.elev)) { + float vdist = p1.elev - p2.elev; + rho *= gridpp::StructureFunction::toar_rho(vdist, v); + } + if(gridpp::is_valid(p1.laf) && gridpp::is_valid(p2.laf)) { + float lafdist = p1.laf - p2.laf; + rho *= gridpp::StructureFunction::toar_rho(lafdist, w); + } + } + else { + if(hdist > localization_distance(p1)) + return 0; + + rho = gridpp::StructureFunction::toar_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::toar_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::toar_rho(lafdist, mW[0][0]); + } + } + return rho; +} +vec gridpp::ToarStructure::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::toar_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::toar_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::toar_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::ToarStructure::clone() const { + return std::make_shared(m_grid, mH, mV, mW, m_min_rho); +} +float gridpp::ToarStructure::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::ToarStructure::localization_distance(float h) const { +// To find an approximate analytical solution for the inverse of the function y=(1+x+1/3*x^2)exp(-x) when y is between 0 and 1, we can follow a similar approach used in approximations for transcendental functions. + if(m_min_rho > 0.5) { + return sqrt( 6 * (1-m_min_rho)) * h; + } + else { + return ( -log(m_min_rho)) * h; + } +} +//----------------------------------------------------------------------------- +/** End TOAR */ +//----------------------------------------------------------------------------- + +//============================================================================= +/** Powerlaw */ +//============================================================================= +gridpp::PowerlawStructure::PowerlawStructure(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 / (1 + 0.5 * pow(hmax / h, 2)); + 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::PowerlawStructure::PowerlawStructure(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::PowerlawStructure::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::powerlaw_rho(hdist, h); + if(gridpp::is_valid(p1.elev) && gridpp::is_valid(p2.elev)) { + float vdist = p1.elev - p2.elev; + rho *= gridpp::StructureFunction::powerlaw_rho(vdist, v); + } + if(gridpp::is_valid(p1.laf) && gridpp::is_valid(p2.laf)) { + float lafdist = p1.laf - p2.laf; + rho *= gridpp::StructureFunction::powerlaw_rho(lafdist, w); + } + } + else { + if(hdist > localization_distance(p1)) + return 0; + + rho = gridpp::StructureFunction::powerlaw_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::powerlaw_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::powerlaw_rho(lafdist, mW[0][0]); + } + } + return rho; +} +vec gridpp::PowerlawStructure::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::powerlaw_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::powerlaw_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::powerlaw_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::PowerlawStructure::clone() const { + return std::make_shared(m_grid, mH, mV, mW, m_min_rho); +} +float gridpp::PowerlawStructure::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::PowerlawStructure::localization_distance(float h) const { + return sqrt( 2 * (1-m_min_rho) / m_min_rho) * h; +} +//----------------------------------------------------------------------------- +/** End Powerlaw */ +//----------------------------------------------------------------------------- + +//============================================================================= +/** Linear */ +//============================================================================= +gridpp::LinearStructure::LinearStructure(float h, float v, float w, float hmin, float vmin, float wmin, float hmax) : + m_is_spatial(false) { + if(gridpp::is_valid(hmax) && hmax < 0) + throw std::invalid_argument("min 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(hmin) || hmin < 0) + throw std::invalid_argument("hmin must be >= 0"); + if(!gridpp::is_valid(vmin) || vmin < 0) + throw std::invalid_argument("vmin must be >= 0"); + if(!gridpp::is_valid(wmin) || wmin < 0) + throw std::invalid_argument("wmin must be >= 0"); + + if(gridpp::is_valid(hmax)) + if(hmax > h) { + m_min_rho = hmin; + } + else { + m_min_rho = (1 - (1-hmin) * 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::LinearStructure::LinearStructure(Grid grid, vec2 h, vec2 v, vec2 w, vec2 hmin, vec2 vmin, vec2 wmin, float min_rho) : + m_grid(grid), + m_min_rho(min_rho), + mH(h), + mV(v), + mW(w), + mHmin(h), + mVmin(v), + mWmin(w) { + if(mH.size() == 1 && mH[0].size() == 1 && mV.size() == 1 && mV[0].size() == 1 && mW.size() == 1 && mW[0].size() == 1 && mHmin.size() == 1 && mHmin[0].size() == 1 && mVmin.size() == 1 && mVmin[0].size() == 1 && mWmin.size() == 1 && mWmin[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() || grid.size()[0] != hmin.size() || grid.size()[0] != vmin.size() || grid.size()[0] != wmin.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() || grid.size()[1] != hmin[0].size() || grid.size()[1] != vmin[0].size() || grid.size()[1] != wmin[0].size()) + throw std::invalid_argument("Grid size not the same as scale size"); + } +} +float gridpp::LinearStructure::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]]; + float hmin = mHmin[I[0]][I[1]]; + float vmin = mVmin[I[0]][I[1]]; + float wmin = mWmin[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, hmin); + if(hdist > loc_dist) + return 0; + + rho = gridpp::StructureFunction::linear_rho(hdist, h, hmin); + if(gridpp::is_valid(p1.elev) && gridpp::is_valid(p2.elev)) { + float vdist = p1.elev - p2.elev; + rho *= gridpp::StructureFunction::linear_rho(vdist, v, vmin); + } + if(gridpp::is_valid(p1.laf) && gridpp::is_valid(p2.laf)) { + float lafdist = p1.laf - p2.laf; + rho *= gridpp::StructureFunction::linear_rho(lafdist, w, wmin); + } + } + else { + if(hdist > localization_distance(p1)) + return 0; + + rho = gridpp::StructureFunction::linear_rho(hdist, mH[0][0], mHmin[0][0]); + if(gridpp::is_valid(p1.elev) && gridpp::is_valid(p2.elev)) { + float vdist = p1.elev - p2.elev; + rho *= gridpp::StructureFunction::linear_rho(vdist, mV[0][0], mVmin[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], mWmin[0][0]); + } + } + return rho; +} +vec gridpp::LinearStructure::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 hmin = mHmin[I[0]][I[1]]; + float vmin = mVmin[I[0]][I[1]]; + float wmin = mWmin[I[0]][I[1]]; + float loc_dist = localization_distance(h, hmin); + 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::linear_rho(hdist, h, hmin); + if(gridpp::is_valid(p1.elev) && gridpp::is_valid(p2[i].elev)) { + float vdist = p1.elev - p2[i].elev; + rho *= gridpp::StructureFunction::linear_rho(vdist, v, vmin); + } + 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, wmin); + } + } + output[i] = rho; + } + } + else { + for(int i = 0; i < p2.size(); i++) { + output[i] = corr(p1, p2[i]); + } + } + return output; +} +gridpp::StructureFunctionPtr gridpp::LinearStructure::clone() const { + return std::make_shared(m_grid, mH, mV, mW, mHmin, mVmin, mWmin, m_min_rho); +} +float gridpp::LinearStructure::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]], mHmin[I[0]][I[1]]); + } + else { + return localization_distance(mH[0][0], mHmin[0][0]); + } +} +float gridpp::LinearStructure::localization_distance(float h, float hmin) const { + if( m_min_rho < hmin) { + return h; + } + else { + return (1 - m_min_rho) / (1 - hmin) * h; + } +} +//----------------------------------------------------------------------------- +/** End Linear */ +//----------------------------------------------------------------------------- + /** CrossValidation */ gridpp::CrossValidation::CrossValidation(StructureFunction& structure, float dist) : StructureFunction(0){ diff --git a/swig/gridpp.i b/swig/gridpp.i index 192f233..3a6b7e4 100644 --- a/swig/gridpp.i +++ b/swig/gridpp.i @@ -46,6 +46,10 @@ %shared_ptr(gridpp::StructureFunction) %shared_ptr(gridpp::MultipleStructure) %shared_ptr(gridpp::BarnesStructure) +%shared_ptr(gridpp::SoarStructure) +%shared_ptr(gridpp::ToarStructure) +%shared_ptr(gridpp::PowerlawStructure) +%shared_ptr(gridpp::LinearStructure) %shared_ptr(gridpp::MixAStructure) %shared_ptr(gridpp::CressmanStructure) %shared_ptr(gridpp::CrossValidation)