Skip to content

Commit 9da2316

Browse files
kilianvolmercharlie0614mknaranjaHenrZureneSchm
authored
1342 make geographic location available to all of MEmilio (#1346)
- Moved GeographicalLocation to memilio/geography/locations.h - Added a distance calculation for two GeographicalLocations on the sphere - Added a wrapper for boosts r-Tree to work with GeographicalLocation - Added a distance struct to work with distances using the correct units Co-authored-by: Carlotta Gerstein <[email protected]> Co-authored-by: Martin J. Kühn <[email protected]> Co-authored-by: HenrZu <[email protected]> Co-authored-by: reneSchm <[email protected]>
1 parent c11273a commit 9da2316

12 files changed

Lines changed: 877 additions & 31 deletions

File tree

‎cpp/memilio/CMakeLists.txt‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -19,9 +19,12 @@ add_library(memilio
1919
epidemiology/lct_infection_state.h
2020
epidemiology/lct_populations.h
2121
epidemiology/adoption_rate.h
22+
epidemiology/simulation_day.h
23+
geography/distance.h
2224
geography/regions.h
2325
geography/regions.cpp
24-
epidemiology/simulation_day.h
26+
geography/geolocation.h
27+
geography/rtree.h
2528
geography/holiday_data.ipp
2629
compartments/compartmental_model.h
2730
compartments/flow_model.h
@@ -116,6 +119,7 @@ add_library(memilio
116119
utils/string_literal.h
117120
utils/type_list.h
118121
utils/base_dir.h
122+
utils/back_inserter_second_element.h
119123
ad/ad.h
120124
)
121125

‎cpp/memilio/geography/distance.h‎

Lines changed: 153 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,153 @@
1+
/*
2+
* Copyright (C) 2020-2025 MEmilio
3+
*
4+
* Authors: Kilian Volmer, Rene Schmiedling
5+
*
6+
* Contact: Martin J. Kuehn <[email protected]>
7+
*
8+
* Licensed under the Apache License, Version 2.0 (the "License");
9+
* you may not use this file except in compliance with the License.
10+
* You may obtain a copy of the License at
11+
*
12+
* http://www.apache.org/licenses/LICENSE-2.0
13+
*
14+
* Unless required by applicable law or agreed to in writing, software
15+
* distributed under the License is distributed on an "AS IS" BASIS,
16+
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17+
* See the License for the specific language governing permissions and
18+
* limitations under the License.
19+
*/
20+
#ifndef MIO_DISTANCE_H
21+
#define MIO_DISTANCE_H
22+
23+
#include "memilio/config.h"
24+
#include "memilio/io/default_serialize.h"
25+
26+
namespace mio
27+
{
28+
29+
namespace geo
30+
{
31+
32+
/**
33+
* @brief Represents a distance.
34+
* Internally, all distances are stored in meters.
35+
*/
36+
class Distance
37+
{
38+
public:
39+
/**
40+
* @brief Default ctor, unitinialized.
41+
*/
42+
Distance() = default;
43+
/**
44+
* @brief Creates a Distance from a specified distance in meters.
45+
* @param[in] meters The distance in meters.
46+
*/
47+
explicit constexpr Distance(ScalarType meters)
48+
: m_meters(meters)
49+
{
50+
}
51+
52+
/**
53+
* @brief return distance in kilometers.
54+
*/
55+
ScalarType kilometers() const
56+
{
57+
return ScalarType(m_meters) / 1000;
58+
}
59+
60+
/**
61+
* @brief return distance in meters.
62+
*/
63+
ScalarType meters() const
64+
{
65+
return m_meters;
66+
}
67+
68+
/**
69+
* @name Comparison operators.
70+
* @{
71+
*/
72+
bool operator==(const Distance& other) const
73+
{
74+
return m_meters == other.m_meters;
75+
}
76+
bool operator!=(const Distance& other) const
77+
{
78+
return !(*this == other);
79+
}
80+
bool operator<(const Distance& other) const
81+
{
82+
return m_meters < other.m_meters;
83+
}
84+
bool operator<=(const Distance& other) const
85+
{
86+
return m_meters <= other.m_meters;
87+
}
88+
bool operator>(const Distance& other) const
89+
{
90+
return m_meters > other.m_meters;
91+
}
92+
bool operator>=(const Distance& other) const
93+
{
94+
return m_meters >= other.m_meters;
95+
}
96+
/**@}*/
97+
98+
/**
99+
* @brief Add or subtract a Distance.
100+
* @{
101+
*/
102+
Distance operator+(const Distance& m) const
103+
{
104+
return Distance{m_meters + m.meters()};
105+
}
106+
Distance& operator+=(const Distance& m)
107+
{
108+
m_meters += m.meters();
109+
return *this;
110+
}
111+
Distance operator-(const Distance& m) const
112+
{
113+
return Distance{m_meters - m.meters()};
114+
}
115+
Distance& operator-=(const Distance& m)
116+
{
117+
m_meters -= m.meters();
118+
return *this;
119+
}
120+
/**@}*/
121+
122+
/// This method is used by the default serialization feature.
123+
auto default_serialize()
124+
{
125+
return Members("Distance").add("meters", m_meters);
126+
}
127+
128+
private:
129+
ScalarType m_meters; ///< The distance in meters.
130+
};
131+
132+
/**
133+
* @brief Create a Distance of a specified number of meters.
134+
* @param[in] meters distance in meters.
135+
*/
136+
constexpr inline Distance meters(ScalarType meters)
137+
{
138+
return Distance(meters);
139+
}
140+
141+
/**
142+
* @brief Create a Distance of a specified number of kilometers.
143+
* @param[in] kilometers distance in kilometers.
144+
*/
145+
constexpr inline Distance kilometers(ScalarType kilometers)
146+
{
147+
return Distance(kilometers * 1000);
148+
}
149+
150+
} // namespace geo
151+
} // namespace mio
152+
153+
#endif

‎cpp/memilio/geography/geolocation.h‎

Lines changed: 195 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,195 @@
1+
/*
2+
* Copyright (C) 2020-2025 MEmilio
3+
*
4+
* Authors: Kilian Volmer, Sascha Korf, Carlotta Gerstein, Daniel Abele, Elisabeth Kluth, Khoa Nguyen, David Kerkmann
5+
*
6+
* Contact: Martin J. Kuehn <[email protected]>
7+
*
8+
* Licensed under the Apache License, Version 2.0 (the "License");
9+
* you may not use this file except in compliance with the License.
10+
* You may obtain a copy of the License at
11+
*
12+
* http://www.apache.org/licenses/LICENSE-2.0
13+
*
14+
* Unless required by applicable law or agreed to in writing, software
15+
* distributed under the License is distributed on an "AS IS" BASIS,
16+
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17+
* See the License for the specific language governing permissions and
18+
* limitations under the License.
19+
*/
20+
#ifndef MIO_GEOGRAPHY_LOCATIONS_H
21+
#define MIO_GEOGRAPHY_LOCATIONS_H
22+
23+
#include "distance.h"
24+
#include "memilio/io/default_serialize.h"
25+
#include "memilio/geography/distance.h"
26+
#include "memilio/config.h"
27+
#include "memilio/utils/logging.h"
28+
#include <cassert>
29+
#include <cmath>
30+
#include <numbers>
31+
namespace mio
32+
{
33+
namespace geo
34+
{
35+
/**
36+
* @brief Class representing a geographical location on the Earth's surface.
37+
*
38+
* Provides latitude and longitude in degrees and a method to calculate the distance to another location.
39+
*/
40+
class GeographicalLocation
41+
{
42+
43+
public:
44+
/**
45+
* @brief Construct a new Geographical Location object.
46+
*
47+
* @param lat Latitude in degrees.
48+
* @param lon Longitude in degrees.
49+
*/
50+
GeographicalLocation(ScalarType lat, ScalarType lon)
51+
: m_latitude(lat)
52+
, m_longitude(lon)
53+
{
54+
check_input();
55+
}
56+
/**
57+
* @brief Construct a new Geographical Location object.
58+
*
59+
* @param coordinates Pair of latitude and longitude in degrees as ScalarTypes.
60+
*/
61+
GeographicalLocation(std::pair<ScalarType, ScalarType> coordinates)
62+
: m_latitude(coordinates.first)
63+
, m_longitude(coordinates.second)
64+
{
65+
check_input();
66+
}
67+
68+
/**
69+
* @brief Construct a new Geographical Location object using the default constructor to keep compatibility with ABM code.
70+
*
71+
*/
72+
GeographicalLocation() = default;
73+
74+
/**
75+
* @brief Compare two GeographicalLocation%s for equality
76+
*/
77+
bool operator==(const GeographicalLocation& other) const
78+
{
79+
return (m_latitude == other.m_latitude && m_longitude == other.m_longitude);
80+
}
81+
82+
/**
83+
* @brief Compare two GeographicalLocation%s for inequality
84+
*/
85+
bool operator!=(const GeographicalLocation& other) const
86+
{
87+
return !(*this == other);
88+
}
89+
90+
/**
91+
* @brief Check that this location is within a (small) distance of another.
92+
* @param[in] other Any location.
93+
* @param[in] tol The Absolute tolerance for considering two locations close.
94+
*/
95+
bool is_close(const GeographicalLocation& other,
96+
Distance tol = Distance(10 * Limits<ScalarType>::zero_tolerance())) const
97+
{
98+
return distance(other) < tol;
99+
}
100+
101+
/**
102+
* @brief Default serialize the GeographicalLocation object.
103+
*
104+
* This method is used by the default serialization feature.
105+
*/
106+
auto default_serialize()
107+
{
108+
return Members("GeographicalLocation").add("latitude", m_latitude).add("longitude", m_longitude);
109+
}
110+
111+
/*
112+
* @brief Calculate the distance between two GeographicalLocation%s.
113+
* @param other The other GeographicalLocation.
114+
* @return The distance between the two locations as @ref mio::geo::Distance .
115+
*
116+
* Uses the haversine formula (https://en.wikipedia.org/wiki/Haversine_formula) to calculate the distance between
117+
* two geographical locations. Uses an earth radius of 6371 km (https://en.wikipedia.org/wiki/Earth_radius).
118+
* The math functions use radians, whereas coordinates are given in degrees.
119+
*/
120+
Distance distance(const GeographicalLocation& other) const
121+
{
122+
const ScalarType delta_latitude = (m_latitude - other.m_latitude) * radians;
123+
const ScalarType delta_longitude = (m_longitude - other.m_longitude) * radians;
124+
const ScalarType sin_delta_latitude_half = sin(delta_latitude * 0.5);
125+
const ScalarType sin_delta_longitude_half = sin(delta_longitude * 0.5);
126+
const ScalarType first_part = sin_delta_latitude_half * sin_delta_latitude_half;
127+
const ScalarType second_part = cos(m_latitude * radians) * cos(other.m_latitude * radians) *
128+
sin_delta_longitude_half * sin_delta_longitude_half;
129+
const ScalarType distance = 2.0 * earth_radius.kilometers() * asin(sqrt(first_part + second_part));
130+
return kilometers(distance);
131+
}
132+
133+
/**
134+
* @brief Get the latitude object
135+
*
136+
* @return ScalarType latitude in degrees
137+
*/
138+
ScalarType get_latitude() const
139+
{
140+
return m_latitude;
141+
}
142+
143+
/**
144+
* @brief Get the longitude object
145+
*
146+
* @return ScalarType longitude in degrees
147+
*/
148+
ScalarType get_longitude() const
149+
{
150+
return m_longitude;
151+
}
152+
153+
/**
154+
* @brief Set the latitude object
155+
*
156+
* @param lat Latitude in degrees.
157+
*/
158+
void set_latitude(ScalarType lat)
159+
{
160+
m_latitude = lat;
161+
check_input();
162+
}
163+
164+
/**
165+
* @brief Set the longitude object
166+
*
167+
* @param lon Longitude in degrees.
168+
*/
169+
void set_longitude(ScalarType lon)
170+
{
171+
m_longitude = lon;
172+
check_input();
173+
}
174+
175+
private:
176+
/**
177+
* @brief Assert that the latitude and longitude are within valid ranges.
178+
*/
179+
void check_input() const
180+
{
181+
assert(m_latitude <= 90. && m_latitude >= -90. && "Latitude must be in [-90, 90]");
182+
assert(m_longitude <= 180. && m_longitude >= -180. && "Longitude must be in [-180, 180]");
183+
}
184+
185+
ScalarType m_latitude;
186+
ScalarType m_longitude;
187+
constexpr static mio::geo::Distance earth_radius = mio::geo::kilometers(6371.);
188+
constexpr static ScalarType radians = std::numbers::pi / 180.0;
189+
};
190+
191+
} // namespace geo
192+
193+
} // namespace mio
194+
195+
#endif // MIO_GEOGRAPHY_LOCATIONS_H

0 commit comments

Comments
 (0)