GIS
December 28, 2017

## Degrees, Minutes, Seconds

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”)$.

## How to know if a point is inside a circle?

The distance between $(x_c,y_c)$ and $(x_p,y_p)$ is given by the Pythagorean theorem as

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

The point $(x_p,y_p)$ 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 $(x_p-x_c)^2+(y_p-y_c)^2$ with $r^2$.

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

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23  Expression earthRadiusKilometers = cb.literal(6371.0); Expression LAT1 = pointofinterest.get(PointofinterestEntity_.latitudeRadian); Expression LON1 = pointofinterest.get(PointofinterestEntity_.longitudeRadian); Expression SIN_LAT1 = cb.function("SIN", Double.class, LAT1); Expression COS_LAT1 = cb.function("COS", Double.class, LAT1); Expression LAT2 = cb.literal(rfp.getPointofinterest().getLatitudeRadian()); Expression LON2 = cb.literal(rfp.getPointofinterest().getLongitudeRadian()); Expression SIN_LAT2 = cb.function("SIN", Double.class, LAT2); Expression COS_LAT2 = cb.function("COS", Double.class, LAT2); Expression LON2_MINUS_LON1 = cb.diff(LON2, LON1); Expression COS_LON2_MINUS_LON1 = cb.function("COS", Double.class, LON2_MINUS_LON1); Expression SIN_LAT1_PROD_SIN_LAT2 = cb.prod(SIN_LAT1, SIN_LAT2); Expression COS_LAT1_PROD_COS_LAT2 = cb.prod(COS_LAT1, COS_LAT2); Expression COS_LAT1_PROD_COS_LAT2_PROD_COS_LON2_MINUS_LON1 = cb.prod(COS_LAT1_PROD_COS_LAT2, COS_LON2_MINUS_LON1); Expression 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 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 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);

## How to convert radians to degrees

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

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!).

 1 2 3 4 5 6  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

 1 2 3 4 5 6 7  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

## Java calculation of the geographic midpoint for two or more points on the earth’s surface

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83  /** * 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 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); }