Generalized convex_hull calculation to use dynamical tlon and tlat instead of hardcoded lon, lat

Author: knutfrode

examples/example_nc_particles.py

@@ -7,6 +7,7 @@
 import xarray as xr
 import trajan as ta
 
+
 #%%
 # Demonstrating analysis of a trajecory_dataset from GNOME in the nc_particles format
 # https://noaa-orr-erd.github.io/nc_particles/nc_particle_standard.html
@@ -45,4 +46,4 @@
 
 #%%
 # Basic animation
-ds_ortho.traj.animate().set_title('nc_particles sample dataset').show()
+ds_ortho.traj.animate(land=None).set_title('nc_particles sample dataset').show()

tests/test_convex_hull.py

@@ -0,0 +1,12 @@
+import pytest
+import numpy as np
+import scipy
+import trajan as ta
+
+
+def test_convex_hull(openoil):
+    ch = openoil.isel(time=-1).traj.convex_hull()
+    assert isinstance(ch, scipy.spatial._qhull.ConvexHull)
+
+    cha = openoil.isel(time=-1).traj.get_area_convex_hull()
+    np.testing.assert_allclose(cha.values, 7.963e+09, rtol=1e-5)

trajan/traj/__init__.py

@@ -433,9 +433,9 @@ def crs(self) -> pyproj.crs.CRS:
         transform, set_crs
         """
         if len(self.ds.cf.grid_mapping_names) == 0:
-            logger.debug(
-                f'No grid-mapping specified, checking if coordinates are lon/lat..'
-            )
+            #logger.debug(
+            #    f'No grid-mapping specified, checking if coordinates are lon/lat..'
+            #)
             if self.tx.name == 'lon' or self.tx.name == 'longitude':
                 # assume this is in latlon projection
                 return self.__gcrs__
@@ -932,11 +932,8 @@ def convex_hull(self):
 
         from scipy.spatial import ConvexHull
 
-        lon = self.ds.lon
-        lat = self.ds.lat
-        if 'status' in self.ds.variables:
-            lon = lon.where(self.ds.status == 0)
-            lat = lat.where(self.ds.status == 0)
+        lon = self.ds.traj.tlon.values.ravel()
+        lat = self.ds.traj.tlat.values.ravel()
         fin = np.isfinite(lat + lon)
         if np.sum(fin) <= 3:
             return None
@@ -980,12 +977,10 @@ def get_area_convex_hull(self):
         """
         from scipy.spatial import ConvexHull
 
-        if 'status' in self.ds.variables:
-            lon = self.ds.lon.where(self.ds.status == 0)  # OpenDrift specific
-            lat = self.ds.lat.where(self.ds.status == 0)
-        else:
-            lon = self.ds.lon.where(np.isfinite(self.ds.lon) is True)
-            lat = self.ds.lat.where(np.isfinite(self.ds.lat) is True)
+        lon = self.ds.traj.tlon.values.ravel()
+        lat = self.ds.traj.tlat.values.ravel()
+        lon = lon[np.isfinite(lon)]
+        lat = lat[np.isfinite(lat)]
 
         fin = np.isfinite(lat + lon)
         if np.sum(fin) <= 3:
@@ -996,10 +991,10 @@ def get_area_convex_hull(self):
             return xr.DataArray(0,
                                 name="convex_hull_area",
                                 attrs={"units": "m2"})
-        lat = lat.where(fin)
-        lon = lon.where(fin)
+        lat = lat[fin]
+        lon = lon[fin]
         aea = pyproj.Proj(
-            f'+proj=aea +lat_0={lat.mean().values} +lat_1={lat.min().values} +lat_2={lat.max().values} +lon_0={lon.mean().values} +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs'
+            f'+proj=aea +lat_0={lat.mean()} +lat_1={lat.min()} +lat_2={lat.max()} +lon_0={lon.mean()} +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs'
         )
         x, y = aea(lat, lon, inverse=False)
         fin = np.isfinite(x + y)