Repository: ecmwf/earthkit-data (python)

Issue:
CDS source does not work with zipped vector data
### What happened?

The CDS source does not work with zipped vector data.

I think something is not right in the way earthkit handles this product.
I just tried, and geopandas is able to read the downloaded file: `geopandas.read_file("rgi6_vector.zip")`.

### What are the steps to reproduce the bug?

```python
import earthkit.data

collection_id = "insitu-glaciers-extent"
request = {
    "variable": "glacier_area",
    "format": "zip",
    "version": "6_0",
    "product_type": "vector",
}
ds = earthkit.data.from_source("cds", collection_id, request)
assert not isinstance(ds, earthkit.data.sources.empty.EmptySource)
```

### Version

0.5.1

### Platform (OS and architecture)

Darwin MacBook-Pro-3.local 22.6.0 Darwin Kernel Version 22.6.0: Wed Jul  5 22:21:56 PDT 2023; root:xnu-8796.141.3~6/RELEASE_X86_64 x86_64

### Relevant log output

_No response_

### Accompanying data

_No response_

### Organisation

B-Open / CADS-EQC

Pull request description:
Add shapefile support
Fixes #299 

Patch:
diff --git a/earthkit/data/readers/shapefile.py b/earthkit/data/readers/shapefile.py
new file mode 100644
index 0000000..de24ca4
--- /dev/null
+++ b/earthkit/data/readers/shapefile.py
@@ -0,0 +1,116 @@
+# (C) Copyright 2020 ECMWF.
+#
+# This software is licensed under the terms of the Apache Licence Version 2.0
+# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+# In applying this licence, ECMWF does not waive the privileges and immunities
+# granted to it by virtue of its status as an intergovernmental organisation
+# nor does it submit to any jurisdiction.
+#
+
+import os
+
+from earthkit.data.readers import Reader
+from earthkit.data.utils.bbox import BoundingBox
+
+
+class ShapeFileReader(Reader):
+    def __init__(self, source, path):
+        super().__init__(source, path)
+        self.__gdf = None
+        assert path.endswith(".shp")
+
+    def __repr__(self):
+        return f"{self.__class__.__name__}({self.path})"
+
+    def _repr_html_(self):
+        html_repr = (
+            f"<h4>{self.__class__.__name__}(represented as a geopandas object):</h3>"
+            f"{self.to_pandas()._repr_html_()}"
+        )
+
+        return html_repr
+
+    @property
+    def _gdf(self):
+        if self.__gdf is None:
+            self.__gdf = self.to_pandas()
+        return self.__gdf
+
+    def __len__(self):
+        return len(self._gdf)
+
+    def __getitem__(self, n):
+        return self._gdf.iloc[n]
+
+    def mutate_source(self):
+        # A ShapeFileReader is a source itself
+        return self
+
+    def bounding_box(self):
+        return BoundingBox.from_geopandas(self._gdf)
+
+    def ls(self, **kwargs):
+        return self._gdf
+
+    def head(self, **kwargs):
+        return self._gdf.head(**kwargs)
+
+    def tail(self, **kwargs):
+        return self._gdf.tail(**kwargs)
+
+    def describe(self, **kwargs):
+        return self._gdf
+
+    def to_numpy(self, flatten=False, **kwargs):
+        v = self.to_pandas().to_numpy(**kwargs)
+        if flatten:
+            import math
+
+            v = v.reshape(
+                math.prod(v.shape),
+            )
+        return v
+
+    def to_pandas(self, **kwargs):
+        if not kwargs and self.__gdf is not None:
+            return self.__gdf
+        else:
+            import geopandas as gpd
+            import numpy as np
+
+            geo_df = gpd.read_file(self.path, **kwargs)
+            return geo_df.set_index(np.arange(len(geo_df)))
+
+    def to_geopandas(self, **kwargs):
+        # TODO: handle multiple paths
+        return self.to_pandas(**kwargs)
+
+    def to_xarray(self, **kwargs):
+        return self.to_pandas(**kwargs).to_xarray()
+
+
+SHAPE_EXT = ".shp"
+MANDATORY = (".shp", ".shx", ".dbf")
+NON_MANDATORY = (".sbn", ".sbx", ".shp.xml", ".prj", ".CPG")
+DOUBLE_DOT_EXT = tuple([e for e in NON_MANDATORY if e.count(".") == 2])
+
+
+def reader(source, path, *, magic=None, deeper_check=False, **kwargs):
+    root, extension = os.path.splitext(path)
+    for e in DOUBLE_DOT_EXT:
+        if path.endswith(e):
+            root = path[: -len(e)]
+            extension = e
+
+    path = root + ".shp"
+
+    # a shapefile consists of multiple files, but we only create
+    # a reader for the .shp file
+    if extension == SHAPE_EXT:
+        if all(os.path.exists(root + e) for e in MANDATORY):
+            return ShapeFileReader(source, path)
+    else:
+        if extension in MANDATORY or extension in NON_MANDATORY:
+            from .unknown import UnknownReader
+
+            return UnknownReader(source, "", skip_warning=True)
diff --git a/earthkit/data/readers/unknown.py b/earthkit/data/readers/unknown.py
index 824a3fe..70da87a 100644
--- a/earthkit/data/readers/unknown.py
+++ b/earthkit/data/readers/unknown.py
@@ -15,10 +15,19 @@ LOG = logging.getLogger(__name__)
 
 
 class UnknownReaderBase(Reader):
-    def __init__(self, source, path="", magic=None, content_type=None, **kwargs):
+    def __init__(
+        self,
+        source,
+        path="",
+        magic=None,
+        content_type=None,
+        skip_warning=False,
+        **kwargs,
+    ):
         super().__init__(source, path)
         self.magic = magic
         self.content_type = content_type
+        self.skip_warning = skip_warning
 
     def ignore(self):
         # Used by multi-source
@@ -31,31 +40,34 @@ class UnknownReaderBase(Reader):
 class UnknownReader(UnknownReaderBase):
     def __init__(self, source, path, **kwargs):
         super().__init__(source, path=path, **kwargs)
-        LOG.warning(
-            (
-                f"Unknown file type, no reader available. "
-                f"path={path} magic={self.magic} content_type={self.content_type}"
+        if not self.skip_warning:
+            LOG.warning(
+                (
+                    f"Unknown file type, no reader available. "
+                    f"path={path} magic={self.magic} content_type={self.content_type}"
+                )
             )
-        )
 
 
 class UnknownStreamReader(UnknownReaderBase):
     def __init__(self, source, data, **kwargs):
         super().__init__(source, **kwargs)
-        LOG.warning(
-            (
-                f"Unknown stream data type, no reader available. "
-                f"magic={self.magic} content_type={self.content_type}"
+        if not self.skip_warning:
+            LOG.warning(
+                (
+                    f"Unknown stream data type, no reader available. "
+                    f"magic={self.magic} content_type={self.content_type}"
+                )
             )
-        )
 
 
 class UnknownMemoryReader(UnknownReaderBase):
     def __init__(self, source, data, **kwargs):
         super().__init__(source, **kwargs)
-        LOG.warning(
-            (
-                f"Unknown memory data type, no reader available. "
-                f"magic={self.magic} content_type={self.content_type}"
+        if not self.skip_warning:
+            LOG.warning(
+                (
+                    f"Unknown memory data type, no reader available. "
+                    f"magic={self.magic} content_type={self.content_type}"
+                )
             )
-        )
diff --git a/earthkit/data/utils/bbox.py b/earthkit/data/utils/bbox.py
index 0952048..1166bd9 100644
--- a/earthkit/data/utils/bbox.py
+++ b/earthkit/data/utils/bbox.py
@@ -77,6 +77,12 @@ class BoundingBox:
 
         return self.as_tuple() == other.as_tuple()
 
+    @classmethod
+    def make_invalid(cls):
+        import numpy as np
+
+        return cls(north=np.nan, west=np.nan, south=np.nan, east=np.nan)
+
     @property
     def width(self):
         """number: Returns the East-West size (degrees)"""
@@ -237,6 +243,14 @@ class BoundingBox:
     def as_dict(self):
         return dict(north=self.north, west=self.west, south=self.south, east=self.east)
 
+    @classmethod
+    def from_geopandas(cls, gdf):
+        try:
+            bb = gdf.crs.area_of_use.bounds
+            return BoundingBox(east=bb[0], south=bb[1], west=bb[2], north=bb[3])
+        except AttributeError:
+            return BoundingBox.make_invalid()
+
 
 def bounding_box(obj):
     if isinstance(obj, BoundingBox):


Test patch:
diff --git a/tests/readers/test_shapefile_reader.py b/tests/readers/test_shapefile_reader.py
new file mode 100644
index 0000000..bda80fd
--- /dev/null
+++ b/tests/readers/test_shapefile_reader.py
@@ -0,0 +1,69 @@
+#!/usr/bin/env python3
+
+# (C) Copyright 2020 ECMWF.
+#
+# This software is licensed under the terms of the Apache Licence Version 2.0
+# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+# In applying this licence, ECMWF does not waive the privileges and immunities
+# granted to it by virtue of its status as an intergovernmental organisation
+# nor does it submit to any jurisdiction.
+#
+import numpy as np
+
+from earthkit.data import from_source
+from earthkit.data.testing import earthkit_test_data_file
+from earthkit.data.utils.bbox import BoundingBox
+
+
+def test_shapefile():
+    ds = from_source("file", earthkit_test_data_file("NUTS_RG_20M_2021_3035.shp.zip"))
+    assert len(ds) == 2010
+
+    # slicing
+    r = ds[1]
+    assert r.NUTS_ID == "HR"
+    assert r.NAME_LATN == "Hrvatska"
+    assert len(r.geometry.geoms) == 16
+
+    # iteration
+    cnt = 0
+    for f in ds:
+        assert isinstance(f.NUTS_ID, str)
+        cnt += 1
+
+    assert cnt == 2010
+
+    # conversion
+    gpf = ds.to_pandas()
+    assert len(gpf) == 2010
+    r = gpf.iloc[1]
+    assert r.NUTS_ID == "HR"
+    assert r.NAME_LATN == "Hrvatska"
+    assert len(r.geometry.geoms) == 16
+
+    gpf = ds.to_geopandas()
+    assert len(gpf) == 2010
+    r = gpf.iloc[1]
+    assert r.NUTS_ID == "HR"
+    assert r.NAME_LATN == "Hrvatska"
+    assert len(r.geometry.geoms) == 16
+
+    v = ds.to_numpy()
+    assert v.shape == (2010, 10)
+
+    v = ds.to_numpy(flatten=True)
+    assert v.shape == (2010 * 10,)
+
+
+def test_shapefile_bounding_box():
+    ds = from_source("file", earthkit_test_data_file("NUTS_RG_20M_2021_3035.shp.zip"))
+
+    ref = BoundingBox(north=84.73, west=44.83, south=24.6, east=324.42)
+    bb = ds.bounding_box()
+    assert np.allclose(np.array(bb.as_list()), np.array(ref.as_list()))
+
+
+if __name__ == "__main__":
+    from earthkit.data.testing import main
+
+    main(__file__)
diff --git a/tests/utils/test_bbox.py b/tests/utils/test_bbox.py
index 6b28c2d..ca2aad8 100644
--- a/tests/utils/test_bbox.py
+++ b/tests/utils/test_bbox.py
@@ -250,6 +250,13 @@ def test_bbox_union_with():
     assert b == BoundingBox(north=60, west=-20, south=-40, east=55)
 
 
+def test_bbox_make_invalid():
+    import numpy as np
+
+    b = BoundingBox.make_invalid()
+    assert b.as_tuple() == (np.nan, np.nan, np.nan, np.nan)
+
+
 if __name__ == "__main__":
     from earthkit.data.testing import main