Dec 28th, 2017

The primary unit in which longitude and latitude are given is **degrees** (°). There are 360° of longitude (180° E ↔ 180° W) and 180° of latitude (90° N ↔ 90° S). Each degree can be broken into 60 **minutes** (’). Each minute can be divided into 60 **seconds** (”). For finer accuracy, fractions of seconds given by a decimal point are used. A base-sixty notation is called a **sexagesimal** notation.

1° = 60’ = 3600”

For example, a spot of ground in upstate New York can be designated by 43°2’27” N, 77°14’30.60” W. Sometimes instead of using minutes and seconds to measure the fraction of a degree, a decimal value is used. With such a convention the coordinates above are 43.040833° N, 77.241833° W. The first number was converted by taking the minutes divided by 60 and the seconds divided by 3600 and adding them together. That is: 43.040833° = 43° + 2’ × (1°/60’) + 27” × (1°/3600”).

The distance between (xc,yc) and (xp,yp) is given by the Pythagorean theorem as

$\displaystyle d=\sqrt{(x_p - x_c)^2 + (y_p - y_c)^2} \lt r$

The point (xp,yp) is inside the circle if d<r, on the circle if d=r, and outside the circle if d>r. You can save yourself a little work by comparing d^2 with r^2 instead: the point is inside the circle if d^2<r^2 , on the circle if d^2=r^2, and outside the circle if d^2>r^2. Thus, you want to compare the number (xp-xc)^2+(yp-yc)^2 with r^2.

JPA Criteria API implementation using **Spherical Law of Cosines** is as below:

```
Expression<Double> earthRadiusKilometers = cb.literal(6371.0);
Expression<Double> LAT1 = pointofinterest.get(PointofinterestEntity_.latitudeRadian);
Expression<Double> LON1 = pointofinterest.get(PointofinterestEntity_.longitudeRadian);
Expression<Double> SIN_LAT1 = cb.function("SIN", Double.class, LAT1);
Expression<Double> COS_LAT1 = cb.function("COS", Double.class, LAT1);
Expression<Double> LAT2 = cb.literal(rfp.getPointofinterest().getLatitudeRadian());
Expression<Double> LON2 = cb.literal(rfp.getPointofinterest().getLongitudeRadian());
Expression<Double> SIN_LAT2 = cb.function("SIN", Double.class, LAT2);
Expression<Double> COS_LAT2 = cb.function("COS", Double.class, LAT2);
Expression<Double> LON2_MINUS_LON1 = cb.diff(LON2, LON1);
Expression<Double> COS_LON2_MINUS_LON1 = cb.function("COS", Double.class, LON2_MINUS_LON1);
Expression<Double> SIN_LAT1_PROD_SIN_LAT2 = cb.prod(SIN_LAT1, SIN_LAT2);
Expression<Double> COS_LAT1_PROD_COS_LAT2 = cb.prod(COS_LAT1, COS_LAT2);
Expression<Double> COS_LAT1_PROD_COS_LAT2_PROD_COS_LON2_MINUS_LON1 = cb.prod(COS_LAT1_PROD_COS_LAT2, COS_LON2_MINUS_LON1);
Expression<Double> SIN_LAT1_PROD_SIN_LAT2_SUM_COS_LAT1_PROD_COS_LAT2_PROD_COS_LON2_MINUS_LON1 = cb.sum(SIN_LAT1_PROD_SIN_LAT2, COS_LAT1_PROD_COS_LAT2_PROD_COS_LON2_MINUS_LON1);
Expression<Double> ACOS_SIN_LAT1_PROD_SIN_LAT2_SUM_COS_LAT1_PROD_COS_LAT2_PROD_COS_LON2_MINUS_LON1 = cb.function("ACOS", Double.class, SIN_LAT1_PROD_SIN_LAT2_SUM_COS_LAT1_PROD_COS_LAT2_PROD_COS_LON2_MINUS_LON1);
Expression<Double> d = cb.prod(ACOS_SIN_LAT1_PROD_SIN_LAT2_SUM_COS_LAT1_PROD_COS_LAT2_PROD_COS_LON2_MINUS_LON1, earthRadiusKilometers);
Predicate predicate = cb.equal(root.get(UserEntity_.userGroup), UserGroup.Sellers);
```

Pi radians are equal to 180 degrees:

```
π rad = 180 degree
```

The angle α in degrees is equal to the angle α in radians times 180 degrees divided by pi constant:

```
α(degrees) = α(radians) × 180° / π
```

Calculate distance, bearing and more between Latitude/Longitude points

**Haversine Formula**

This uses the ‘haversine’ formula to calculate the great-circle distance between two points – that is, the shortest distance over the earth’s surface – giving an ‘as-the-crow-flies’ distance between the points (ignoring any hills they fly over, of course!).

```
a = sin²(Δφ/2) + cos φ1 ⋅ cos φ2 ⋅ sin²(Δλ/2)
c = 2 ⋅ atan2( √a, √(1−a) )
d = R ⋅ c
where φ is latitude, λ is longitude, R is earth’s radius (mean radius = 6,371km);
note that angles need to be in radians to pass to trig functions!
```

**Spherical Law of Cosines**

```
d = acos( sin φ1 ⋅ sin φ2 + cos φ1 ⋅ cos φ2 ⋅ cos Δλ ) ⋅ R
# Excel calculation in radians
=ACOS( SIN(lat1)*SIN(lat2) + COS(lat1)*COS(lat2)*COS(lon2-lon1) ) * 6371000
# Excel calculation in degrees
=ACOS( SIN(lat1*PI()/180)*SIN(lat2*PI()/180) + COS(lat1*PI()/180)*COS(lat2*PI()/180)*COS(lon2*PI()/180-lon1*PI()/180) ) * 6371000
```

```
/**
* Calculation of the geographic midpoint (also known as the geographic center,
* or center of gravity) for two or more points on the earth's surface Ref:
* http://www.geomidpoint.com/calculation.html
*
* @param points
* @return
*/
public static Point2D.Double midpoint(List<Point2D.Double> points) {
double Totweight = 0;
double xt = 0;
double yt = 0;
double zt = 0;
for (Point2D.Double point : points) {
Double latitude = point.x;
Double longitude = point.y;
/**
* Convert Lat and Lon from degrees to radians.
*/
double latn = latitude * Math.PI / 180;
double lonn = longitude * Math.PI / 180;
/**
* Convert lat/lon to Cartesian coordinates
*/
double xn = Math.cos(latn) * Math.cos(lonn);
double yn = Math.cos(latn) * Math.sin(lonn);
double zn = Math.sin(latn);
/**
* Compute weight (by time) If locations are to be weighted equally,
* set wn to 1
*/
double years = 0;
double months = 0;
double days = 0;
double wn = true ? 1 : (years * 365.25) + (months * 30.4375) + days;
/**
* Compute combined total weight for all locations.
*/
Totweight = Totweight + wn;
xt += xn * wn;
yt += yn * wn;
zt += zn * wn;
}
/**
* Compute weighted average x, y and z coordinates.
*/
double x = xt / Totweight;
double y = yt / Totweight;
double z = zt / Totweight;
/**
* If abs(x) < 10-9 and abs(y) < 10-9 and abs(z) < 10-9 then the
* geographic midpoint is the center of the earth.
*/
double lat = -0.001944;
double lon = -78.455833;
if (Math.abs(x) < Math.pow(10, -9) && Math.abs(y) < Math.pow(10, -9) && Math.abs(z) < Math.pow(10, -9)) {
} else {
/**
* Convert average x, y, z coordinate to latitude and longitude.
* Note that in Excel and possibly some other applications, the
* parameters need to be reversed in the atan2 function, for
* example, use atan2(X,Y) instead of atan2(Y,X).
*/
lon = Math.atan2(y, x);
double hyp = Math.sqrt(x * x + y * y);
lat = Math.atan2(z, hyp);
/**
* Convert lat and lon to degrees.
*/
lat = lat * 180 / Math.PI;
lon = lon * 180 / Math.PI;
}
LOG.log(Level.INFO, "MidPoint: {0}, {1}", new Object[]{lat, lon});
return new Point2D.Double(lat, lon);
}
```