adding datamesh readers

opendrift/readers/reader_datamesh_remote.py

@@ -0,0 +1,97 @@
+import os
+import datetime
+import pyproj
+import pandas as pd
+import numpy as np
+from oceanum.datamesh import Connector
+from opendrift.readers.basereader import BaseReader, ContinuousReader, StructuredReader
+
+
+datamesh = Connector()
+
+
+class DataException(Exception):
+    pass
+
+
+class DatameshReader(BaseReader, StructuredReader):
+    name = "datamesh"
+
+    def __init__(self, datasource_id, mapping={}):
+        self.proj4 = "+proj=lonlat +ellps=WGS84"  # Only working for WGS84 grids
+        self.proj = pyproj.Proj(self.proj4)
+        self.mapping = mapping
+        try:
+            self.dset = datamesh.load_datasource(datasource_id)
+            self._vars = list(self.dset.variables.keys())
+            self.variables = [v for v in mapping if mapping[v] in self._vars]
+            self.latax, self.lonax = self._get_grid_axes()
+            self.lon = self.dset[self.lonax].values
+            self.lat = self.dset[self.latax].values
+            self.time = self.dset["time"].values
+            self.xmin = self.lon.min()
+            self.xmax = self.lon.max()
+            self.ymin = self.lat.min()
+            self.ymax = self.lat.max()
+            self.delta_x = self.lon[1] - self.lon[0]
+            self.delta_y = self.lat[1] - self.lat[0]
+            self.start_time = (
+                pd.to_datetime(self.time[0]).tz_localize("UTC").to_pydatetime()
+            )
+            self.end_time = (
+                pd.to_datetime(self.time[-1]).tz_localize("UTC").to_pydatetime()
+            )
+            self.time_step = (
+                pd.to_datetime(self.time[1]) - pd.to_datetime(self.time[0])
+            ).to_pytimedelta()
+
+        except Exception as e:
+            raise DataException(e)
+
+        # Run constructor of parent Reader class
+        super(DatameshReader, self).__init__()
+        # Do this again to reset UTC
+        self.start_time = (
+            pd.to_datetime(self.time[0]).tz_localize("UTC").to_pydatetime()
+        )
+        self.end_time = pd.to_datetime(self.time[-1]).tz_localize("UTC").to_pydatetime()
+
+    def _get_grid_axes(self):
+        da_tmp = self.dset.get(self.mapping[self.variables[0]])
+        lonax = da_tmp.dims[-1]
+        latax = da_tmp.dims[-2]
+        return latax, lonax
+
+    def get_variables(self, requested_variables, time=None, x=None, y=None, z=None):
+
+        requested_variables, time, x, y, z, outside = self.check_arguments(
+            requested_variables, time, x, y, z
+        )
+
+        nearestTime, dummy1, dummy2, indxTime, dummy3, dummy4 = self.nearest_time(time)
+
+        variables = {}
+        delta = self.buffer * self.delta_x
+        lonmin = np.maximum(x.min() - delta, self.xmin)
+        lonmax = np.minimum(x.max() + delta, self.xmax)
+        latmin = np.maximum(y.min() - delta, self.ymin)
+        latmax = np.minimum(y.max() + delta, self.ymax)
+
+        if self.delta_y > 0:
+            latslice = slice(latmin, latmax)
+        else:
+            latslice = slice(latmax, latmin)
+        lonslice = slice(lonmin, lonmax)
+
+        for var in requested_variables:
+            subset = self.dset.get(self.mapping[var])
+            subset = subset.isel(time=indxTime)
+            variables[var] = subset.sel(
+                {self.latax: latslice, self.lonax: lonslice}
+            ).values
+        variables["x"] = self.dset[self.lonax].sel({self.lonax: lonslice}).values
+        variables["y"] = self.dset[self.latax].sel({self.latax: latslice}).values
+        variables["z"] = None
+        variables["time"] = nearestTime
+
+        return variables

opendrift/readers/reader_datamesh_shoreline.py

@@ -0,0 +1,70 @@
+import os
+import pygeos
+import geopandas
+from collections import OrderedDict
+from sqlalchemy import create_engine
+from opendrift.readers.basereader import BaseReader, ContinuousReader, StructuredReader
+
+from oceanum.datamesh import Connector
+
+datamesh = Connector()
+
+
+class ShorelineException(Exception):
+    pass
+
+
+# This reader utilises the shoreline database
+class ShorelineReader(BaseReader, ContinuousReader):
+    name = "shoreline"
+    variables = ["land_binary_mask"]
+    proj4 = None
+    crs = None
+    skippoly = False
+    datasource_select = OrderedDict(
+        {
+            1.0: "osm-land-polygons",
+            5.0: "gshhs_f_l1",
+            20.0: "gshhs_h_l1",
+            50.0: "gshhs_i_l1",
+            180.0: "gshhs_c_l1",
+            1000.0: "gshhs_l_l1",
+        }
+    )
+
+    def __init__(self, extent=None, skippoly=False):
+        self.proj4 = "+proj=lonlat +ellps=WGS84"
+        self.skippoly = skippoly
+
+        super(ShorelineReader, self).__init__()
+        self.z = None
+        if extent is not None:
+            self.xmin, self.ymin, self.xmax, self.ymax = extent
+        else:
+            self.xmin, self.ymin = -180, -90
+            self.xmax, self.ymax = 180, 90
+
+        domain_size = 1.0 * max(self.xmax - self.xmin, self.ymax - self.ymin)
+        for res in self.datasource_select:
+            if domain_size <= res:
+                datasource = self.datasource_select[res]
+                break
+
+        query = {
+            "datasource": datasource,
+            "geofilter": {
+                "type": "bbox",
+                "geom": [self.xmin, self.ymin, self.xmax, self.ymax],
+            },
+        }
+
+        self.shorelines = datamesh.query(query)
+
+    def __on_land__(self, x, y):
+        points = geopandas.GeoDataFrame({"geometry": pygeos.creation.points(x, y)})
+        test = geopandas.sjoin(self.shorelines, points, how="right")
+        return test.index_left >= 0
+
+    def get_variables(self, requestedVariables, time=None, x=None, y=None, z=None):
+        self.check_arguments(requestedVariables, time, x, y, z)
+        return {"land_binary_mask": self.__on_land__(x, y)}