add HALO flight circle dataset decomposition

eurec4a_environment/decomposition/space.py

@@ -0,0 +1,70 @@
+import cartopy.crs as ccrs
+import cartopy.geodesic as cgeod
+import numpy as np
+import xarray as xr
+
+
+from .. import nomenclature as nom, get_field
+
+
+class Circle:
+    def __init__(self, lat, lon, r_degrees):
+        self.lat = lat
+        self.lon = lon
+        self.r_degrees = r_degrees
+
+        self.radius = self._compute_radius()
+        self.geod = cgeod.Geodesic()  # Earth geodesic
+
+    def _compute_radius(self):
+        # Define the projection used to display the circle:
+        ortho = ccrs.Orthographic(central_longitude=self.lon, central_latitude=self.lat)
+
+        # Compute the required radius in projection native coordinates:
+        phi1 = self.lat + self.r_degrees if self.lat <= 0 else self.lat - self.r_degrees
+        _, y1 = ortho.transform_point(self.lon, phi1, ccrs.PlateCarree())
+        return abs(y1)
+
+    def distance_to_center(self, lat, lon):
+        @np.vectorize
+        def _dist(lat1, lon1, lat2, lon2):
+            geometry = np.array(
+                [
+                    [lon1, lat1],
+                    [lon2, lat2],
+                ]
+            )
+            return self.geod.geometry_length(geometry)
+
+        return _dist(lat1=self.lat, lon1=self.lon, lat2=lat, lon2=lon)
+
+    def contains(self, lat, lon):
+        """
+        Return whether a point is contained in the circle or not
+        The points are assumed to be in spherical (lon, lat) coordinates
+        """
+        p_dist = self.distance_to_center(lat=lat, lon=lon)
+        return p_dist < self.radius
+
+
+class HALOCircle(Circle):
+    # circle location and radius in degrees:
+    lat = 13 + (18 / 60)
+    lon = -(57 + (43 / 60))
+    r_degrees = 1
+
+    def __init__(self):
+        super().__init__(lat=self.lat, lon=self.lon, r_degrees=self.r_degrees)
+
+
+def inside_halo_circle(ds, lat=nom.LATITUDE, lon=nom.LONGITUDE):
+    circle = HALOCircle()
+    da_lat = get_field(ds=ds, name=lat, units="degrees")
+    da_lon = get_field(ds=ds, name=lon, units="degrees")
+    is_inside = circle.contains(lat=da_lat, lon=da_lon)
+    return xr.DataArray(
+        is_inside,
+        coords=da_lat.coords,
+        dims=da_lat.dims,
+        attrs=dict(long_name="Inside HALO flight circle", units="1"),
+    )

eurec4a_environment/nomenclature.py

@@ -46,6 +46,10 @@
 MERIDIONAL_WIND = "v"
 
 
+LATITUDE = "lat"
+LONGITUDE = "lon"
+
+
 CF_STANDARD_NAMES = {
     TEMPERATURE: "air_temperature",
     ALTITUDE: "geopotential_height",
@@ -56,6 +60,8 @@
     WIND_SPEED: "wind_speed",
     ZONAL_WIND: "eastward_wind",
     MERIDIONAL_WIND: "northward_wind",
+    LATITUDE: "latitude",
+    LONGITUDE: "longitude",
 }
 
 

setup.py

@@ -17,6 +17,8 @@
     "aiohttp",
     "cfunits",
     "intake-xarray @ git+https://github.com/leifdenby/intake-xarray#egg=intake-xarray",
+    "numpy",
+    "cartopy",
 ]
 setup(
     name="eurec4a-environment",

tests/test_decomposition.py

@@ -1,4 +1,7 @@
+import numpy as np
+
 from eurec4a_environment.decomposition import time as time_decomposition
+from eurec4a_environment.decomposition import space as space_decomposition
 from eurec4a_environment.source_data import open_joanne_dataset
 
 
@@ -10,3 +13,32 @@ def test_time_decomposition():
 
     ds_mean_per_segment = ds.groupby("halo_segment_id").mean()
     assert ds_mean_per_segment.halo_segment_id.count() > 0
+
+
+def test_halo_circle_space_decomposition():
+    circle = space_decomposition.HALOCircle()
+
+    # check for single set point
+    assert circle.contains(lat=circle.lat, lon=circle.lon)
+    # check for array of points
+    assert np.all(
+        circle.contains(
+            lat=np.array([circle.lat + circle.r_degrees, circle.lat]),
+            lon=np.array([circle.lon, circle.lon + circle.r_degrees]),
+        )
+    )
+    assert not circle.contains(
+        lat=circle.lat + circle.r_degrees, lon=circle.lon + circle.r_degrees
+    )
+
+    ds_joanne = open_joanne_dataset()
+    # take mean per-sounding so we find if the mean lat/lon position is inside
+    # the cirlce or not (just to make things a bit faster)
+    ds_per_sounding = ds_joanne.mean(dim="height")
+    ds_joanne["inside_circle"] = space_decomposition.inside_halo_circle(
+        ds=ds_per_sounding
+    )
+
+    # there should be points both inside and outside the circle for the dropsondes
+    assert ds_joanne.where(ds_joanne.inside_circle, drop=True).count() > 0
+    assert ds_joanne.where(~ds_joanne.inside_circle, drop=True).count() > 0