Skip to content

Commit

Permalink
Change haversine location
Browse files Browse the repository at this point in the history
  • Loading branch information
thorbjoernl committed Nov 6, 2024
1 parent d7a6e68 commit 4387f64
Show file tree
Hide file tree
Showing 4 changed files with 30 additions and 26 deletions.
21 changes: 21 additions & 0 deletions src/pyaro/mathutils.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,21 @@
import numpy as np


EARTH_RADIUS = 6378137 # meters


def haversine(lon1, lat1, lon2, lat2):
"""
Calculate the great-circle distance between two points on the Earth (specified in decimal degrees).
returns:
Distance (in meters)
"""
# Convert decimal degrees to radians
lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])
# Haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = np.sin(dlat / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2) ** 2
c = 2 * np.arcsin(np.sqrt(a))
m = EARTH_RADIUS * c
return m
2 changes: 1 addition & 1 deletion src/pyaro/timeseries/AutoFilterReaderEngine.py
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@ def supported_filters(cls) -> list[Filter]:
:return: list of filters
"""
supported = "variables,stations,countries,bounding_boxes,duplicates,time_bounds,time_resolution,flags,time_variable_station,altitude,relaltitude,valleyfloorrelativealtitudefilter".split(
supported = "variables,stations,countries,bounding_boxes,duplicates,time_bounds,time_resolution,flags,time_variable_station,altitude,relaltitude,valleyfloor_relaltitude".split(
","
)
return [filters.get(name) for name in supported]
Expand Down
31 changes: 7 additions & 24 deletions src/pyaro/timeseries/Filter.py
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,9 @@
from .Data import Data, Flag
from .Station import Station

from ..mathutils import haversine


try:
# Optional dependencies required for relative altitude filter.
import xarray as xr
Expand Down Expand Up @@ -1188,8 +1191,6 @@ def filter_stations(self, stations: dict[str, Station]) -> dict[str, Station]:

@registered_filter
class ValleyFloorRelativeAltitudeFilter(StationFilter):
EARTH_RADIUS = 6378137 # meters

# https://cfconventions.org/Data/cf-conventions/cf-conventions-1.11/cf-conventions.html#latitude-coordinate
UNITS_LAT = set(
[
Expand Down Expand Up @@ -1227,11 +1228,11 @@ def __init__(
def topography(self):
if "cf_units" not in sys.modules:
raise ModuleNotFoundError(
"valleyfloorrelativealtitudefilter filter is missing required dependency 'cf-units'. Please install to use this filter."
"valleyfloor_relaltitude filter is missing required dependency 'cf-units'. Please install to use this filter."
)
if "xarray" not in sys.modules:
raise ModuleNotFoundError(
"valleyfloorrelativealtitudefilter filter is missing required dependency 'xarray'. Please install to use this filter."
"valleyfloor_relaltitude filter is missing required dependency 'xarray'. Please install to use this filter."
)

if self._topography is None:
Expand All @@ -1255,7 +1256,7 @@ def init_kwargs(self):
}

def name(self):
return "valleyfloorrelativealtitudefilter"
return "valleyfloor_relaltitude"

def filter_stations(self, stations: dict[str, Station]) -> dict[str, Station]:
filtered_stations = {}
Expand Down Expand Up @@ -1291,7 +1292,7 @@ def _calculate_relative_altitude(
lon=slice(lon - s, lon + s), lat=slice(lat - s, lat + s)
)

distances = self._haversine(topo["lon"], topo["lat"], lon, lat)
distances = haversine(topo["lon"], topo["lat"], lon, lat)
within_radius = distances <= radius

values_within_radius = topo[self._topo_var].where(within_radius, drop=True)
Expand Down Expand Up @@ -1325,21 +1326,3 @@ def _find_lat_lon_variables(self, topo_xr):
f"Required variable names for lat, lon dimensions could not be found in file '{self._topo_file}"
)
return lat, lon

def _haversine(self, lon1, lat1, lon2, lat2):
"""
Calculate the great-circle distance between two points on the Earth (specified in decimal degrees).
returns:
Distance (in meters)
"""
# Convert decimal degrees to radians
lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])

# Haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = np.sin(dlat / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2) ** 2
c = 2 * np.arcsin(np.sqrt(a))
m = ValleyFloorRelativeAltitudeFilter.EARTH_RADIUS * c
return m
2 changes: 1 addition & 1 deletion tests/test_CSVTimeSeriesReader.py
Original file line number Diff line number Diff line change
Expand Up @@ -674,7 +674,7 @@ def test_valley_floor_filter(self):
filename=self.elevation_file,
filters=[
pyaro.timeseries.filters.get(
"valleyfloorrelativealtitudefilter",
"valleyfloor_relaltitude",
topo_file="/lustre/storeB/project/aerocom/aerocom1/AEROCOM_OBSDATA/GTOPO30/merged/N.nc",
radius=5000,
upper=500,
Expand Down

0 comments on commit 4387f64

Please sign in to comment.