KJB
|
a class to interpolate elevation values using Gaussian processes More...
#include <kriging.h>
Public Member Functions | |
Kriging_interpolator (double, const Ned13_one_degree_grid::IntegralLL &, Ned13_grid_cache *, Pthread_mutex *=0) | |
Create object that can interpolate points in a small neighborhood. More... | |
double | operator() (const TopoFusion::pt &utm) const |
interpolate (krige) the elevation at the given query point More... | |
double | diff_e (const TopoFusion::pt &) const |
compute the east-west slope of the terrain at query pt More... | |
double | diff_n (const TopoFusion::pt &) const |
compute the north-south slope of the terrain at query pt More... | |
a class to interpolate elevation values using Gaussian processes
This class performs Gaussian process interpolation of elevation values using a technique commonly known as "kriging." It is a straighforward application of Gaussian processes but the data structures and interface is customized for GIS data, specifically elevation.
Objects of this class have a substantial memory footprint. They are copyable but you should avoid copying them unnecessarily.
Basic usage: the idea is that you have access to a sparse grid of known elevation points, in the form of NED13 data, and you wish to densely interpolate points in the center grid cell. At least, this is how I use it. To do so, you instantiate this class with a pointer to a NED 13 grid cache, the characteristic length of the Gaussian Process (which determines the smoothness of the interpolation) and the center point of the grid, i.e., a representative point nearby where all other interpolations will be performed. Then you can use the function operator to fetch the estimated elevation. Also the diff_e and diff_n functions let you fetch the east-west and north-south gradient at a point.
If you have multiple threads all sharing a single NED grid cache and trying to build these objects, you must help serialize cache access, because the NED grid cache fills itself from the file system, and you don't want two or more parallel attempts to enlarge the cache that way, especially for the same grid file. In that case, just create a kjb::Pthread_mutex object elsewhere, and all threads should provide to the ctor that mutex pointer. If the above conditions do not apply (e.g., single-threaded or the cache is not shared) then omit the pointer or pass in a NULL value.
If you wish to have all three output values, elevation and the two gradient components, then be sure to call all these functions one after another at the same point, without any intervening queries at a different point, because the object caches its internal state. This makes the gradient computation incrementally cheap when it immediately precedes or follows the elevation results. If you fail to maintain this "sequential-spatial coherency," performance will degrade.
kjb::Kriging_interpolator::Kriging_interpolator | ( | double | , |
const Ned13_one_degree_grid::IntegralLL & | , | ||
Ned13_grid_cache * | , | ||
Pthread_mutex * | = 0 |
||
) |
Create object that can interpolate points in a small neighborhood.
The "small" neighborhood size is roughly a 10x10 meters. You create the interpolator object using this ctor, then you can interpolate any point about 7 meters away or less, from 'training_center_loc' the center point, and this will give you a well-behaved elevation estimate, and incline (gradient) estimate too if you want it. The interpolation method is called "kriging" in the geostatistics literature, but in the machine learning community we would call that "using a Gaussian process for regression." The GP uses the popular squared-exponential kernel (a/k/a covariance).
characteristic_spatial_frequency_squared | This value represents the spatial variability of the correlation of the GP model. Smaller means slower-varying, and a smoother output. Don't make it too smooth. For NED13 data, 1/400 is a nice value, 1/100 is probably too quick, and 1/900 is probably too smooth. The inverse square root is a value known as the "characteristic length." |
training_center_loc | The centerpoint, as a NED13 grid location (i.e., an earth coordinate using latitude and longitude on a one-third arcsecond grid.) It does not have to be the exact center, but for best results it should be closer to your query points (which are sent to the functor operator). The value in 'training_center_loc' represents the center of the square of training data used to train the Gaussian process model for this object. |
cache | A NED13 grid caching object, able to read NED13 elevation data. |
pmt_grid_cache_serializer | Optional mutex pointer – set equal to NULL in single-threaded applications. If multiple threads will be building these objects, and if the threads are sharing a common NED grid cache, you should pass in a pointer to a common mutex – |
double kjb::Kriging_interpolator::diff_e | ( | const TopoFusion::pt & | ) | const |
compute the east-west slope of the terrain at query pt
If you plan to query both gradients, or elevation, please
double kjb::Kriging_interpolator::diff_n | ( | const TopoFusion::pt & | ) | const |
compute the north-south slope of the terrain at query pt
If you plan to query both gradients, or elevation, please
|
inline |
interpolate (krige) the elevation at the given query point
If you plan to query gradients too, please