III. Background
3.3 Geospatial Analysis Background
3.3.1 Calculating distance between two given points
3.3.1 Calculating distance between two given points Great-‐circle distance
Great-‐circle distance formula is based on a spherical model of the earth and uses average great-‐circle radius equal 6371 kilometers.
To calculate the distance between two given points Point1=(latitute1, longitude1) and Point2=(latitude2, longitude2) great-‐circle distance formula was used:
where R is the average great-‐circle radius.
For example, to calculate distance between NCTU (24.789345,120.999867) and NTHU (24.794463,120.990142), the following steps are required:
• Convert latitude and longitude measured in degrees to radians using formula:
𝑅𝑎𝑑𝑖𝑎𝑛𝑠 = 𝐷𝑒𝑔𝑟𝑒𝑒𝑠 ∗ 𝑃𝐼 / 180
NCTU (24.789345, 120.999867) = (0.4326556896627937 rad, 2.1118460736252334 rad) NTHU (24.794463, 120.990142) = (0.43274501561391077 rad, 2.1116763403554772 rad)
• Calculate distance using great-‐circle distance formula mentioned above 𝑑𝑖𝑠𝑡𝑎𝑛𝑐𝑒
= 𝑎𝑟𝑐𝑜𝑠 𝑠𝑖𝑛 0.4326556896627937 · 𝑠𝑖𝑛 0.43274501561391077 + 𝑐𝑜𝑠 0.4326556896627937 · 𝑐𝑜𝑠 0.43274501561391077
· 𝑐𝑜𝑠 2.1118460736252334 − 2.1116763403554772 · 6371
= 1.1347336565854425 𝑘𝑚
Great-‐circle distance in Python could be calculated using math module:
>>> from math import cos, sin, acos
>>> acos(sin(0.4326556896627937) * sin(0.43274501561391077) + cos (0.4326556896627937)
* cos (0.43274501561391077) * cos (2.1118460736252334 -‐ 2.1116763403554772)) * 6371 1.1347336565854425
This distance calculation algorithm is based on based on spherical trigonometry and assumption that the Earth is a sphere. However, because the earth is only nearly spherical, this formula might result small errors in calculation around 0.3%.
Vincenty distance formula
Vincenty distance formula is based on a more accurate ellipsoidal model of the Earth and considered to be accurate to within 0.5mm on ellipsoids.
The main idea is that the Earth is not perfectly spherical; it has an irregular shape usually called “geoid”. The geoid is represented as a reference (theoretical) ellipsoid, which is still spherical but flattened at the poles.
There are several accepted models for the earth ellipsoid. The most accurate and globally accepted model is so-‐called WGS-‐84 ellipsoid. WGS-‐84 is used as the default model in most of the distance calculation scripts and libraries.
As mentioned in [34] Vincenty formula could be used in a program/function in the following way in order to calculate distance between two points (with given latitude/longitude):
Notations:
a, b = major & minor semiaxes of the ellipsoid
f = flattening (a−b)/a
φ1, φ2 = geodetic latitude
L = difference in longitude
U1 = atan((1−f) ∙ tanφ1) (U is ‘reduced latitude’)
U2 = atan((1−f) ∙ tanφ2)
λ = L (first approximation)
Authors [34] propose to iterate until change in λ is negligible (e.g. 10-‐12 ≈ 0.006mm) {
Δσ = B∙sinσ∙{cos2σm+B/4∙[cosσ∙(-‐1+2∙cos²2σm) -‐ B/6∙cos2σm∙(-‐3+4∙sin²σ)∙(-‐
3+4∙cos²2σm)]}
Python libraries to calculate distance
Geopy (Python Geocoding Toolbox) is an open source Python library that allows performing geocoding tasks (currently 6 third-‐party geocoders are supported including Google Maps, Bing and GeoNames) together with geographical computations.
Geopy includes several core modules, including distance.py and point.py used for this project.
distance.py
This module allows calculating distance between two points A and B. Calculations can be performed based on either Great-‐circle distance formula, or the Vincenty Distance formula.
As the latter one considered to be more accurate, it is used by default when calculating distances using geopy. WGS-‐84 is a default ellipsoid model used by distance.py module with the Vincenty Distance formula.
Calculate distance using non-‐relational (NoSQL) database and spatial indexes MongoDb
MongoDB is a document oriented non-‐relational database, which uses JSON format to store objects (documents). Moreover, MongoDB has native support for geospatial indexes and provides easy access to geospatial documents.
A standard document in MongoDB has the following structure:
{ "_id" : ObjectId("5188f7309ada8d42fa826788"), "MarkerID" : "03", "title" : "Bench NCTU",
"loc" : [ 121.0011879, 24.7873803 ], "category" : "bench", "description" : "Bench in the NCTU" }
In the sample above information about a bench in NCTU is stored. This information includes geographical coordinates (lat/lon). In this way, MongoDB can be populated with wide range of POIs – other benches, traffic lights, drug stores, etc.
As have been mentioned before, MondoDB allows creating Geospatial Indexes and provides several commands for geospatial queries. There are two types of Geospatial Indexes to choose from:
• a 2dsphere index and
• a 2d index
2dsphere indexes support calculations on a sphere, whereas 2d indexes compute based on the flat geometry. The default ellipsoid model for 2dsphere indexes is WGS84.
The following commands can be used in MongoDB shell to create Indexes for a collection:
> db.PoiPoi.ensureIndex({ loc : "2dsphere" } ) (to create 2dsphere index)
> db.PoiPoi.ensureIndex({ loc : "2d" } ) (to create 2d index)
If cases there were no indexes created (at least one), the following error will be generated when trying to perform geospatial query:
error: {
"$err" : "can't find any special indices: 2d (needs index), 2dsphere (needs index), for:
{ loc: { $near: [ 121.001228, 24.788367 ], $maxDistance: 0.0001 } }", "code" : 13038
}
To check whether 2d and/or 2dsphere indexes were generated, the following command should be executed in MongoDB shell:
> db.PoiPoi.getIndexes()
In this case no geospatial indexes were created. The sample below illustrates .getIndexes() command output after both indexes were generated:
> db.PoiPoi.getIndexes()
"name" : "loc_2dsphere"
},
{
"v" : 1,
"key" : {
"loc" : "2d"
},
To find a bench within five meters from point(lon: 121.001228, lat: 24.788367) the following query should be performed:
> db.PoiPoi.find({loc: { $near : [ 121.001228, 24.788367 ], $maxDistance : 0.005 } })
{ "_id" : ObjectId("51b75630aea18129d71e6e0a"), "MarkerID" : "01", "title" : "Bench NCTU", "loc" : [ 121.001199, 24.78832 ], "category" : "bench", "description" : "Bench in the NCTU" }
3.3.2 Calculate bearing between two given points
Bearing is an angle between point A forward direction to point B, measured in degrees (°) from the north line in a clockwise direction.
To calculate the bearing for a great-‐circle route from point A to point B, the following
This formula gives only initial bearing, which will be changing along the path. The bearing should be compared to a mobile phone user heading in order to form a navigation guidance notification.
As an example, we calculate the initial bearing from Point A (NCTU 24.789345,120.999867) to Point B (NTHU 24.794463,120.990142) in Python by using the formula mentioned above:
>>> from math import cos, sin, atan2, degrees, radians
(importing required names from math module)
>>> lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
(converting latitude and longitude measured in degrees to radians)
>>> lon1, lat1, lon2, lat2
(0.036858667280074905, 0.007551266312102546, 0.036855704875667486,
0.007552825344057015)
>>> dlon = lon2 -‐ lon1