diff --git a/pygeoapi/api/__init__.py b/pygeoapi/api/__init__.py
index b47541f23..12fae3ece 100644
--- a/pygeoapi/api/__init__.py
+++ b/pygeoapi/api/__init__.py
@@ -1138,6 +1138,18 @@ def describe_collections(self, request: Union[APIRequest, Any],
                     except ProviderTypeError:
                         pass
                     else:
+                        for name, mimetype in p.supported_formats.items():
+                            if name in (F_JSON, collection_data_format['name']):  # noqa
+                                continue
+                            title_ = l10n.translate('Coverage data as', request.locale)  # noqa
+                            title_ = f"{title_} {name}"
+                            collection['links'].append({
+                                'type': mimetype,
+                                'rel': f'{OGC_RELTYPES_BASE}/coverage',
+                                'title': title_,
+                                'href': f"{self.get_collections_url()}/{k}/coverage?f={name}"  # noqa
+                            })
+
                         collection['extent']['spatial']['grid'] = [{
                             'cellsCount': p._coverage_properties['width'],
                             'resolution': p._coverage_properties['resx']
diff --git a/pygeoapi/api/coverages.py b/pygeoapi/api/coverages.py
index c6687042c..c12e38675 100644
--- a/pygeoapi/api/coverages.py
+++ b/pygeoapi/api/coverages.py
@@ -189,6 +189,9 @@ def get_collection_coverage(
 
         headers['Content-Type'] = collection_def['format']['mimetype']
         return headers, HTTPStatus.OK, data
+    elif format_ in p.supported_formats:
+        headers['Content-Type'] = p.supported_formats[format_]
+        return headers, HTTPStatus.OK, data
     elif format_ == F_JSON:
         headers['Content-Type'] = 'application/prs.coverage+json'
         return headers, HTTPStatus.OK, to_json(data, api.pretty_print)
diff --git a/pygeoapi/provider/base.py b/pygeoapi/provider/base.py
index 5f9456870..e74ac9cfd 100644
--- a/pygeoapi/provider/base.py
+++ b/pygeoapi/provider/base.py
@@ -80,6 +80,7 @@ def __init__(self, provider_def):
         self.axes = []
         self.crs = None
         self.num_bands = None
+        self.supported_formats = {}  # e.g. {'netcdf':'application/x-netcdf'}
 
     def get_fields(self):
         """
diff --git a/pygeoapi/provider/xarray_.py b/pygeoapi/provider/xarray_.py
index ba835f033..2d81245c2 100644
--- a/pygeoapi/provider/xarray_.py
+++ b/pygeoapi/provider/xarray_.py
@@ -46,7 +46,7 @@
                                     ProviderConnectionError,
                                     ProviderNoDataError,
                                     ProviderQueryError)
-from pygeoapi.util import get_crs_from_uri, read_data
+from pygeoapi.util import get_crs_from_uri
 
 LOGGER = logging.getLogger(__name__)
 
@@ -63,6 +63,9 @@ def __init__(self, provider_def):
 
         super().__init__(provider_def)
 
+        self.supported_formats['netcdf'] = 'application/x-netcdf'
+        self.supported_formats['zarr'] = 'application/zip+zarr'
+
         try:
             if provider_def['data'].endswith('.zarr'):
                 open_func = xarray.open_zarr
@@ -136,13 +139,6 @@ def query(self, properties=[], subsets={}, bbox=[], bbox_crs=4326,
         :returns: coverage data as dict of CoverageJSON or native format
         """
 
-        if not properties and not subsets and format_ != 'json':
-            LOGGER.debug('No parameters specified, returning native data')
-            if format_ == 'zarr':
-                return _get_zarr_data(self._data)
-            else:
-                return read_data(self.data)
-
         if len(properties) < 1:
             properties = self.fields.keys()
 
@@ -229,26 +225,27 @@ def query(self, properties=[], subsets={}, bbox=[], bbox_crs=4326,
             # json does not support float32
             data = _convert_float32_to_float64(data)
 
-        out_meta = {
-            'bbox': [
-                data.coords[self.x_field].values[0],
-                data.coords[self.y_field].values[0],
-                data.coords[self.x_field].values[-1],
-                data.coords[self.y_field].values[-1]
-            ],
-            "driver": "xarray",
-            "height": data.sizes[self.y_field],
-            "width": data.sizes[self.x_field],
-            "variables": {var_name: var.attrs
-                          for var_name, var in data.variables.items()}
-        }
+            out_meta = {
+                'bbox': [
+                    data.coords[self.x_field].values[0],
+                    data.coords[self.y_field].values[0],
+                    data.coords[self.x_field].values[-1],
+                    data.coords[self.y_field].values[-1]
+                ],
+                "driver": "xarray",
+                "height": data.sizes[self.y_field],
+                "width": data.sizes[self.x_field],
+                "variables": {var_name: var.attrs
+                              for var_name, var in data.variables.items()}
+            }
 
-        if self.time_field is not None:
-            out_meta['time'] = [
-                _to_datetime_string(data.coords[self.time_field].values[0]),
-                _to_datetime_string(data.coords[self.time_field].values[-1]),
-            ]
-            out_meta["time_steps"] = data.sizes[self.time_field]
+            if self.time_field is not None:
+                time_values = data.coords[self.time_field].values
+                out_meta['time'] = [
+                    _to_datetime_string(time_values[0]),
+                    _to_datetime_string(time_values[-1]),
+                ]
+                out_meta["time_steps"] = data.sizes[self.time_field]
 
         LOGGER.debug('Serializing data in memory')
         if format_ == 'json':