/*
* This file is part of ELKI:
* Environment for Developing KDD-Applications Supported by Index-Structures
*
* Copyright (C) 2017
* ELKI Development Team
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
package de.lmu.ifi.dbs.elki.math.geodesy;
import static org.junit.Assert.assertEquals;
import static org.junit.Assert.assertFalse;
import org.junit.Test;
/**
* Unit test that cross-validates some distance computations with different
* implementations. Note that it is common to see some difference in these
* computations, unfortunately - even within R: "sp", "SDMTools" and "geosphere"
* give different values. But in the end, it's always just an approximation.
*
* @author Erich Schubert
* @since 0.6.0
*/
public class EarthModelsTest {
// New York City
final double[] NEW_YORK = new double[] { 40.67, -73.94 };
// Munich City
final double[] MUNICH = new double[] { 48.133333, 11.566667 };
// Close to Munich city.
final double[] MUNICH_AIRPORT = new double[] { 48.353889, 11.786111 };
// Beijing
final double[] BEIJING = new double[] { 39.913889, 116.391667 };
// Sydney
final double[] SYDNEY = new double[] { -33.859972, 151.211111 };
// South pole
final double[] SOUTH = new double[] { -90., 123. };
// North pole
final double[] NORTH = new double[] { 90., -11. };
// Null
final double[] NULL = new double[] { 0., 0. };
// Antipodal of New York
final double[] ANEW_YORK = new double[] { -40.67, 180. + -73.94 };
// The full data set
final double[][] DATA = new double[][] { NEW_YORK, MUNICH, MUNICH_AIRPORT, BEIJING, SYDNEY, SOUTH, NORTH, ANEW_YORK, NULL };
// Names, for reporting
final String[] NAMES = new String[] { "New York", "Munich", "Munich Airport", "Beijing", "Sydney", "South Pole", "North Pole", "AnitpodalNY", "Null" };
// Distance matrix, as computed by R "sp" package.
final double[][] R_SP_WGS84 = new double[][] { //
{ 0.00000000e+00, 6.50366148e+06, 6.50686263e+06, //
1.10191059e+07, 1.59817953e+07, 1.44857385e+07, //
5.49643716e+06, 2.00089752e+07, 8.65713703e+06, //
}, //
{ 6.50366148e+06, 0.00000000e+00, 2.94317181e+04, //
7.77130051e+06, 1.63129394e+07, 1.53127242e+07, //
4.66770838e+06, 1.34974501e+07, 5.44172321e+06, //
}, //
{ 6.50686263e+06, 2.94317181e+04, 0.00000000e+00, //
7.74340170e+06, 1.62922795e+07, 1.53371802e+07, //
4.64319863e+06, 1.34940877e+07, 5.46970119e+06, //
}, //
{ 1.10191059e+07, 7.77130051e+06, 7.74340170e+06, //
0.00000000e+00, 8.90595105e+06, 1.44020175e+07, //
5.58032509e+06, 8.97021429e+06, 1.22298940e+07, //
}, //
{ 1.59817953e+07, 1.63129394e+07, 1.62922795e+07, //
8.90595105e+06, 0.00000000e+00, 6.25152406e+06, //
1.37320772e+07, 4.03034506e+06, 1.52087327e+07, //
}, //
{ 1.44857385e+07, 1.53127242e+07, 1.53371802e+07, //
1.44020175e+07, 6.25152406e+06, 0.00000000e+00, //
1.99703264e+07, 5.49643716e+06, 9.99356093e+06, //
}, //
{ 5.49643716e+06, 4.66770838e+06, 4.64319863e+06, //
5.58032509e+06, 1.37320772e+07, 1.99703264e+07, //
0.00000000e+00, 1.44857385e+07, 9.99356093e+06, //
}, //
{ 2.00089752e+07, 1.34974501e+07, 1.34940877e+07, //
8.97021429e+06, 4.03034506e+06, 5.49643716e+06, //
1.44857385e+07, 0.00000000e+00, 1.13581966e+07, //
}, //
{ 8.65713703e+06, 5.44172321e+06, 5.46970119e+06, //
1.22298940e+07, 1.52087327e+07, 9.99356093e+06, //
9.99356093e+06, 1.13581966e+07, 0.00000000e+00, //
}, //
};
// Distance matrix, as computed by R "SDMTools" with avoiding segfaults
final double[][] SDM_WGS84 = new double[][] { //
{ 0.00000000e+00, 6.49607649e+06, 6.49934560e+06, //
1.10186899e+07, 1.59618991e+07, 1.45058923e+07, //
5.49803918e+06, 1.99845542e+07, 8.64487707e+06, //
}, //
{ 6.49607649e+06, 0.00000000e+00, 2.94333400e+04, //
7.76413413e+06, 1.63019586e+07, 1.53352188e+07, //
4.66871262e+06, 1.34959988e+07, 5.44849943e+06, //
}, //
{ 6.49934560e+06, 2.94333400e+04, 0.00000000e+00, //
7.73629144e+06, 1.62815962e+07, 1.53597436e+07, //
4.64418788e+06, 1.34928316e+07, 5.47647977e+06, //
}, //
{ 1.10186899e+07, 7.76413413e+06, 7.73629144e+06, //
0.00000000e+00, 8.91481946e+06, 1.44219335e+07, //
5.58199793e+06, 8.98486478e+06, 1.22110932e+07, //
}, //
{ 1.59618991e+07, 1.63019586e+07, 1.62815962e+07, //
8.91481946e+06, 0.00000000e+00, 6.25383635e+06, //
1.37500951e+07, 4.02274801e+06, 1.51936337e+07, //
}, //
{ 1.45058923e+07, 1.53352188e+07, 1.53597436e+07, //
1.44219335e+07, 6.25383635e+06, 0.00000000e+00, //
2.00039315e+07, 5.49803918e+06, 1.00019657e+07, //
}, //
{ 5.49803918e+06, 4.66871262e+06, 4.64418788e+06, //
5.58199793e+06, 1.37500951e+07, 2.00039315e+07, //
0.00000000e+00, 1.45058923e+07, 1.00019657e+07, //
}, //
{ 1.99845542e+07, 1.34959988e+07, 1.34928316e+07, //
8.98486478e+06, 4.02274801e+06, 5.49803918e+06, //
1.45058923e+07, 0.00000000e+00, 1.13403339e+07, //
}, //
{ 8.64487707e+06, 5.44849943e+06, 5.47647977e+06, //
1.22110932e+07, 1.51936337e+07, 1.00019657e+07, //
1.00019657e+07, 1.13403339e+07, 0.00000000e+00, //
}, //
};
// Distance matrix, cosine distance in "geosphere" package:
final double[][] GEOSPHERE_COSINE = new double[][] { //
{ 0.000000000000000e+00, 6.486356338365966e+06, 6.489535193201737e+06, //
1.099371420315689e+07, 1.599344217645813e+07, 1.452986159024009e+07, //
5.485253480114362e+06, 2.001511507035445e+07, 8.660734833379321e+06, //
}, //
{ 6.486356338365966e+06, 0.000000000000000e+00, 2.941993519886598e+04, //
7.750900176348899e+06, 1.632910711431917e+07, 1.535974752803661e+07, //
4.655367542317842e+06, 1.352875873198849e+07, 5.467217714777914e+06, //
}, //
{ 6.489535193201737e+06, 2.941993519886598e+04, 0.000000000000000e+00, //
7.723050836613588e+06, 1.630864631433774e+07, 1.538427227092249e+07, //
4.630842799431969e+06, 1.352557987715272e+07, 5.495190660363145e+06, //
}, //
{ 1.099371420315689e+07, 7.750900176348899e+06, 7.723050836613588e+06, //
0.000000000000000e+00, 8.948985711203361e+06, 1.444578576429031e+07, //
5.569329306064145e+06, 9.021400867197569e+06, 1.222413440358347e+07, //
}, //
{ 1.599344217645813e+07, 1.632910711431917e+07, 1.630864631433774e+07, //
8.948985711203361e+06, 0.000000000000000e+00, 6.242495113738451e+06, //
1.377261995661600e+07, 4.021672893896329e+06, 1.520028800598781e+07, //
}, //
{ 1.452986159024009e+07, 1.535974752803661e+07, 1.538427227092249e+07, //
1.444578576429031e+07, 6.242495113738451e+06, 0.000000000000000e+00, //
2.001511507035445e+07, 5.485253480114362e+06, 1.000755753517723e+07, //
}, //
{ 5.485253480114362e+06, 4.655367542317842e+06, 4.630842799431969e+06, //
5.569329306064145e+06, 1.377261995661600e+07, 2.001511507035445e+07, //
0.000000000000000e+00, 1.452986159024009e+07, 1.000755753517723e+07, //
}, //
{ 2.001511507035445e+07, 1.352875873198849e+07, 1.352557987715272e+07, //
9.021400867197569e+06, 4.021672893896329e+06, 5.485253480114362e+06, //
1.452986159024009e+07, 0.000000000000000e+00, 1.135438023697513e+07, //
}, //
{ 8.660734833379321e+06, 5.467217714777914e+06, 5.495190660363145e+06, //
1.222413440358347e+07, 1.520028800598781e+07, 1.000755753517723e+07, //
1.000755753517723e+07, 1.135438023697513e+07, 0.000000000000000e+00, //
}, //
};
// Distance matrix, cosine distance in "geosphere" package:
final double[][] GEOSPHERE_HAVERSINE = new double[][] { //
{ 0.000000000000000e+00, 6.486356338365965e+06, 6.489535193201736e+06, //
1.099371420315688e+07, 1.599344217645813e+07, 1.452986159024009e+07, //
5.485253480114363e+06, 2.001511507035445e+07, 8.660734833379321e+06, //
}, //
{ 6.486356338365965e+06, 0.000000000000000e+00, 2.941993519893691e+04, //
7.750900176348901e+06, 1.632910711431917e+07, 1.535974752803662e+07, //
4.655367542317841e+06, 1.352875873198849e+07, 5.467217714777914e+06, //
}, //
{ 6.489535193201736e+06, 2.941993519893691e+04, 0.000000000000000e+00, //
7.723050836613588e+06, 1.630864631433774e+07, 1.538427227092249e+07, //
4.630842799431968e+06, 1.352557987715272e+07, 5.495190660363145e+06, //
}, //
{ 1.099371420315688e+07, 7.750900176348901e+06, 7.723050836613588e+06, //
0.000000000000000e+00, 8.948985711203359e+06, 1.444578576429031e+07, //
5.569329306064145e+06, 9.021400867197569e+06, 1.222413440358347e+07, //
}, //
{ 1.599344217645813e+07, 1.632910711431917e+07, 1.630864631433774e+07, //
8.948985711203359e+06, 0.000000000000000e+00, 6.242495113738450e+06, //
1.377261995661600e+07, 4.021672893896329e+06, 1.520028800598781e+07, //
}, //
{ 1.452986159024009e+07, 1.535974752803662e+07, 1.538427227092249e+07, //
1.444578576429031e+07, 6.242495113738450e+06, 0.000000000000000e+00, //
2.001511507035445e+07, 5.485253480114363e+06, 1.000755753517723e+07, //
}, //
{ 5.485253480114363e+06, 4.655367542317841e+06, 4.630842799431968e+06, //
5.569329306064145e+06, 1.377261995661600e+07, 2.001511507035445e+07, //
0.000000000000000e+00, 1.452986159024009e+07, 1.000755753517723e+07, //
}, //
{ 2.001511507035445e+07, 1.352875873198849e+07, 1.352557987715272e+07, //
9.021400867197569e+06, 4.021672893896329e+06, 5.485253480114363e+06, //
1.452986159024009e+07, 0.000000000000000e+00, 1.135438023697513e+07, //
}, //
{ 8.660734833379321e+06, 5.467217714777914e+06, 5.495190660363145e+06, //
1.222413440358347e+07, 1.520028800598781e+07, 1.000755753517723e+07, //
1.000755753517723e+07, 1.135438023697513e+07, 0.000000000000000e+00, //
}, //
};
// Distance matrix, cosine distance in "geosphere" package:
final double[][] GEOSPHERE_VINCENTY_SPHERE = new double[][] { //
{ 0.000000000000000e+00, 6.486356338365965e+06, 6.489535193201737e+06, //
1.099371420315689e+07, 1.599344217645813e+07, 1.452986159024009e+07, //
5.485253480114362e+06, 2.001511507035445e+07, 8.660734833379321e+06, //
}, //
{ 6.486356338365966e+06, 0.000000000000000e+00, 2.941993519893652e+04, //
7.750900176348899e+06, 1.632910711431917e+07, 1.535974752803661e+07, //
4.655367542317842e+06, 1.352875873198849e+07, 5.467217714777914e+06, //
}, //
{ 6.489535193201737e+06, 2.941993519893681e+04, 0.000000000000000e+00, //
7.723050836613588e+06, 1.630864631433774e+07, 1.538427227092249e+07, //
4.630842799431969e+06, 1.352557987715272e+07, 5.495190660363145e+06, //
}, //
{ 1.099371420315689e+07, 7.750900176348899e+06, 7.723050836613588e+06, //
0.000000000000000e+00, 8.948985711203361e+06, 1.444578576429031e+07, //
5.569329306064145e+06, 9.021400867197569e+06, 1.222413440358347e+07, //
}, //
{ 1.599344217645813e+07, 1.632910711431917e+07, 1.630864631433774e+07, //
8.948985711203361e+06, 0.000000000000000e+00, 6.242495113738451e+06, //
1.377261995661600e+07, 4.021672893896329e+06, 1.520028800598781e+07, //
}, //
{ 1.452986159024009e+07, 1.535974752803661e+07, 1.538427227092249e+07, //
1.444578576429031e+07, 6.242495113738451e+06, 0.000000000000000e+00, //
2.001511507035445e+07, 5.485253480114362e+06, 1.000755753517723e+07, //
}, //
{ 5.485253480114362e+06, 4.655367542317842e+06, 4.630842799431969e+06, //
5.569329306064145e+06, 1.377261995661600e+07, 2.001511507035445e+07, //
0.000000000000000e+00, 1.452986159024009e+07, 1.000755753517723e+07, //
}, //
{ 2.001511507035445e+07, 1.352875873198849e+07, 1.352557987715272e+07, //
9.021400867197569e+06, 4.021672893896328e+06, 5.485253480114362e+06, //
1.452986159024009e+07, 0.000000000000000e+00, 1.135438023697513e+07, //
}, //
{ 8.660734833379321e+06, 5.467217714777914e+06, 5.495190660363145e+06, //
1.222413440358347e+07, 1.520028800598781e+07, 1.000755753517723e+07, //
1.000755753517723e+07, 1.135438023697513e+07, 0.000000000000000e+00, //
}, //
};
// Distance matrix, cosine distance in "geosphere" package:
final double[][] GEOSPHERE_VINCENTY_WGS84 = new double[][] { //
{ 0.000000000000000e+00, 6.503767848090770e+06, 6.506974602807214e+06, //
1.101910369872154e+07, 1.599270764104164e+07, 1.450589228275554e+07, //
5.498039175797141e+06, Double.NaN, 8.661075093360968e+06, //
}, //
{ 6.503767848090770e+06, 0.000000000000000e+00, 2.944678565227988e+04, //
7.771413604903577e+06, 1.632562828361734e+07, 1.533521883690098e+07, //
4.668712621651703e+06, 1.351202375815323e+07, 5.449076622477260e+06, //
}, //
{ 6.506974602807214e+06, 2.944678565227961e+04, 0.000000000000000e+00, //
7.743520548885130e+06, 1.630501450760304e+07, 1.535974357438163e+07, //
4.644187884171052e+06, 1.350871506702908e+07, 5.477072159975932e+06, //
}, //
{ 1.101910369872154e+07, 7.771413604903577e+06, 7.743520548885130e+06, //
0.000000000000000e+00, 8.918923554119145e+06, 1.442193352759095e+07, //
5.581997930961725e+06, 8.985204289439090e+06, 1.223306246456788e+07, //
}, //
{ 1.599270764104165e+07, 1.632562828361734e+07, 1.630501450760304e+07, //
8.918923554119145e+06, 0.000000000000000e+00, 6.253836350029638e+06, //
1.375009510852304e+07, 4.030502513001381e+06, 1.521100587039189e+07, //
}, //
{ 1.450589228275554e+07, 1.533521883690098e+07, 1.535974357438163e+07, //
1.442193352759095e+07, 6.253836350029638e+06, 0.000000000000000e+00, //
2.000393145855268e+07, 5.498039175797141e+06, 1.000196572927634e+07, //
}, //
{ 5.498039175797141e+06, 4.668712621651703e+06, 4.644187884171051e+06, //
5.581997930961725e+06, 1.375009510852304e+07, 2.000393145855268e+07, //
0.000000000000000e+00, 1.450589228275554e+07, 1.000196572927634e+07, //
}, //
{ Double.NaN, 1.351202375815323e+07, 1.350871506702908e+07, //
8.985204289439090e+06, 4.030502513001381e+06, 5.498039175797141e+06, //
1.450589228275554e+07, 0.000000000000000e+00, 1.136156548644650e+07, //
}, //
{ 8.661075093360968e+06, 5.449076622477260e+06, 5.477072159975932e+06, //
1.223306246456788e+07, 1.521100587039189e+07, 1.000196572927634e+07, //
1.000196572927634e+07, 1.136156548644650e+07, 0.000000000000000e+00, //
}, //
};
@Test
public void testWGS84SpheroidEarth() {
// WGS84 Vincenty to WGS84 Haversine: .3% error on test set.
testEarthModel(WGS84SpheroidEarthModel.STATIC, R_SP_WGS84, .00316, 1e-12);
// WGS84 Vincenty to WGS84 Vincenty: also .2% error on test set.
testEarthModel(WGS84SpheroidEarthModel.STATIC, SDM_WGS84, .001936, 1e-12);
// WGS84 Vincenty to WGS84 Vincenty: with "geosphere" we have a high
// agreement.
testEarthModel(WGS84SpheroidEarthModel.STATIC, GEOSPHERE_VINCENTY_WGS84, 6.1763e-12, 1e-12);
}
@Test
public void testHaversineEarth() {
// Spherical Haversine to WGS84 Haversine: .6% error on test set.
testEarthModel(SphericalHaversineEarthModel.STATIC, R_SP_WGS84, .005674, 1e-12);
// Spherical Haversine to WGS84 Vincenty: .4% error on test set.
testEarthModel(SphericalHaversineEarthModel.STATIC, SDM_WGS84, .00405, 1e-12);
// WGS84 Vincenty to WGS84 Vincenty: with "geosphere" we have a high
// agreement.
testEarthModel(SphericalHaversineEarthModel.STATIC, GEOSPHERE_HAVERSINE, 6.662e-14, 1e-12);
}
@Test
public void testCosineEarth() {
// Spherical Cosine to WGS84 Haversine: 0% error on test set.
testEarthModel(SphericalCosineEarthModel.STATIC, R_SP_WGS84, .005674, 0);
// Spherical Cosine to WGS84 Vincenty: 0% error on test set.
testEarthModel(SphericalCosineEarthModel.STATIC, SDM_WGS84, .00405, 0.);
// Spherical Cosine to Cosine: with "geosphere" we have a high agreement.
testEarthModel(SphericalCosineEarthModel.STATIC, GEOSPHERE_COSINE, 4.4409e-16, 3.060705e-7);
}
@Test
public void testSphericalEarth() {
// Spherical Vincenty to WGS84 Haversine: .6% error on test set.
testEarthModel(SphericalVincentyEarthModel.STATIC, R_SP_WGS84, .005674, 1e-12);
// Spherical Vincenty to WGS84 Vincenty: .4% error on test set.
testEarthModel(SphericalVincentyEarthModel.STATIC, SDM_WGS84, .00405, 1e-12);
// Spherical Vincenty to Spherical Vincenty: with "geosphere" we have a high
// agreement.
testEarthModel(SphericalVincentyEarthModel.STATIC, GEOSPHERE_VINCENTY_SPHERE, 1.044e-14, 1e-12);
}
protected void testEarthModel(EarthModel model, final double[][] ref, final double relerror, final double abserror) {
double maxrel = 0.0, maxabs = 0.0;
for (int i = 0; i < DATA.length; i++) {
for (int j = i; j < DATA.length; j++) {
if (Double.isNaN(ref[i][j])) {
continue;
}
double d = model.distanceDeg(DATA[i][0], DATA[i][1], DATA[j][0], DATA[j][1]);
assertFalse("NaN in distance " + NAMES[i] + " to " + NAMES[j], Double.isNaN(d));
double test = (d > 0) ? (ref[i][j] / d - 1.) : (ref[i][j] - d);
if (Math.abs(d - ref[i][j]) > abserror) {
if (Math.abs(test) > relerror) {
assertEquals("Distances do not agree for " + NAMES[i] + " to " + NAMES[j] + " " + Math.abs(test), ref[i][j], d, abserror);
}
if (Math.abs(test) > maxrel) {
maxrel = Math.abs(test);
}
} else {
if (Math.abs(ref[i][j] - d) > maxabs) {
maxabs = Math.abs(ref[i][j] - d);
}
}
}
}
assertEquals("Relative error bound not tight.", maxrel, relerror, 1e-3 * relerror + 1e-12);
assertEquals("Absolute error bound not tight.", maxabs, abserror, 1e-3 * abserror + 1e-12);
}
}