From 2b8ee589f061a50694b2698690e9730ada3b6d8a Mon Sep 17 00:00:00 2001 From: Carlos Velasco Date: Tue, 17 Aug 2021 22:57:25 +1000 Subject: [PATCH] Bom import xarray (#228) * Add import_bom_rf3 using xarray * Add tests to xarray version * Fix mrms importer tests * Pass **kwargs to internal functions * Add nwp_importers to read bom nwp sample data * Add bom nwp data to source file * Add tests for bom_nwp reader * Fix pystepsrc Co-authored-by: Andres Perez Hortal <16256571+aperezhortal@users.noreply.github.com> --- pysteps/io/__init__.py | 1 + pysteps/io/importers.py | 164 ++++++++++++++++- pysteps/io/nwp_importers.py | 277 +++++++++++++++++++++++++++++ pysteps/pystepsrc | 11 ++ pysteps/tests/test_io_bom_nwp.py | 82 +++++++++ pysteps/tests/test_io_bom_rf3.py | 71 +++++++- pysteps/tests/test_io_mrms_grib.py | 13 +- 7 files changed, 610 insertions(+), 9 deletions(-) create mode 100644 pysteps/io/nwp_importers.py create mode 100644 pysteps/tests/test_io_bom_nwp.py diff --git a/pysteps/io/__init__.py b/pysteps/io/__init__.py index 6f5fd1677..b37a2f5d4 100644 --- a/pysteps/io/__init__.py +++ b/pysteps/io/__init__.py @@ -8,3 +8,4 @@ from .importers import * from .nowcast_importers import * from .readers import * +from .nwp_importers import * diff --git a/pysteps/io/importers.py b/pysteps/io/importers.py index 081855c0b..40d865c77 100644 --- a/pysteps/io/importers.py +++ b/pysteps/io/importers.py @@ -457,6 +457,164 @@ def import_bom_rf3(filename, **kwargs): return precip, None, metadata +@postprocess_import() +def import_bom_rf3_xr(filename, **kwargs): + """Import a NetCDF radar rainfall product from the BoM Rainfields3 + using xarray. + + Parameters + ---------- + filename: str + Name of the file to import. + + {extra_kwargs_doc} + + Returns + ------- + out_da : xr.DataArray + A xarray DataArray containing the rainfall field in mm/h imported + from the Bureau RF3 netcdf, the quality field and the metadata. The + quality field is currently set to None. + """ + + if not NETCDF4_IMPORTED: + raise MissingOptionalDependency( + "netCDF4 package is required to import BoM Rainfields3 products " + "but it is not installed" + ) + + ds = _import_bom_rf3_data_xr(filename, **kwargs,) + ds_meta = _import_bom_rf3_geodata_xr(ds, **kwargs,) + + # rename valid_time to t if exists + if 'valid_time' in ds_meta: + ds_meta = ds_meta.rename({'valid_time': 't'}) + + return ds_meta.precipitation + + +def _import_bom_rf3_data_xr(filename, **kwargs,): + + varname_time = kwargs.get('varname_time', 'valid_time') + # Tested in python3.6 and chunks did not work properly + # commenting next line until find the reason + # chunks = kwargs.get('chunks', {varname_time: 1}) + + ds_rainfall = xr.open_mfdataset( + filename, + combine='nested', + concat_dim=varname_time, + # chunks=chunks, + lock=False, + parallel=True, + ) + + return ds_rainfall + + +def _import_bom_rf3_geodata_xr(ds_in, + **kwargs, + ): + + # select a default varname if none is passed + varname = kwargs.get('varname', 'precipitation') + + # extract useful information + # projection + projdef = None + if "proj" in ds_in: + projection = ds_in.proj + if projection.grid_mapping_name == "albers_conical_equal_area": + projdef = "+proj=aea " + lon_0 = projection.longitude_of_central_meridian + projdef += f" +lon_0={lon_0:.3f}" + lat_0 = projection.latitude_of_projection_origin + projdef += f" +lat_0={lat_0:.3f}" + standard_parallel = projection.standard_parallel + projdef += f" +lat_1={standard_parallel[0]:.3f}" + projdef += f" +lat_2={standard_parallel[1]:.3f}" + + # get the accumulation period + valid_time = None + if "valid_time" in ds_in: + valid_time = ds_in.valid_time + + start_time = None + if "start_time" in ds_in: + start_time = ds_in.start_time + + time_step = None + if start_time is not None: + if valid_time is not None: + time_step = (valid_time - start_time).isel(valid_time=0) + time_step = time_step.values.astype('timedelta64[m]') + + # get the units of precipitation + units = None + if "units" in ds_in[varname].attrs: + units = ds_in[varname].units + if units in ("kg m-2", "mm"): + units = "mm" + ds_in[varname].attrs.update({'units': units}) + # get spatial boundaries and pixelsize + # move to meters if coordiantes in kilometers + if "units" in ds_in.x.attrs: + if ds_in.x.units == "km": + ds_in['x'] = ds_in.x*1000. + ds_in.x.attrs.update({'units': 'm'}) + ds_in['y'] = ds_in.y*1000. + ds_in.y.attrs.update({'units': 'm'}) + + xmin = ds_in.x.min().values + xmax = ds_in.x.max().values + ymin = ds_in.y.min().values + ymax = ds_in.y.max().values + xpixelsize = abs(ds_in.x[1] - ds_in.x[0]) + ypixelsize = abs(ds_in.y[1] - ds_in.y[0]) + + cartesian_unit = ds_in.x.units + + # Add metadata needed by pySTEPS as attrs in X and Y variables + + ds_in.x.attrs.update({ + # TODO: Remove before final 2.0 version + "x1": xmin, + "x2": xmax, + "cartesian_unit": cartesian_unit, + } + ) + + ds_in.y.attrs.update({ + # TODO: Remove before final 2.0 version + "y1": ymin, + "y2": ymax, + "cartesian_unit": cartesian_unit, + } + ) + + # Add metadata needed by pySTEPS as attrs in rainfall variable + da_rainfall = ds_in[varname].isel(valid_time=0) + + ds_in[varname].attrs.update( + {"transform": None, + "unit": units, # copy 'units' in 'unit' for legacy reasons + "projection": projdef, + "accutime": time_step, + "zr_a": None, + "zr_b": None, + "zerovalue": np.nanmin(da_rainfall), + "institution": "Commonwealth of Australia, Bureau of Meteorology", + "threshold": _get_threshold_value(da_rainfall.values), + # TODO(_import_bom_rf3_geodata_xr): Remove before final 2.0 version + "yorigin": "upper", + "xpixelsize": xpixelsize.values, + "ypixelsize": ypixelsize.values, + } + ) + + return ds_in + + def _import_bom_rf3_data(filename): ds_rainfall = netCDF4.Dataset(filename) if "precipitation" in ds_rainfall.variables.keys(): @@ -523,7 +681,8 @@ def _import_bom_rf3_geodata(filename): calendar = "standard" if "calendar" in times.ncattrs(): calendar = times.calendar - valid_time = netCDF4.num2date(times[:], units=times.units, calendar=calendar) + valid_time = netCDF4.num2date(times[:], units=times.units, + calendar=calendar) start_time = None if "start_time" in ds_rainfall.variables.keys(): @@ -531,7 +690,8 @@ def _import_bom_rf3_geodata(filename): calendar = "standard" if "calendar" in times.ncattrs(): calendar = times.calendar - start_time = netCDF4.num2date(times[:], units=times.units, calendar=calendar) + start_time = netCDF4.num2date(times[:], units=times.units, + calendar=calendar) time_step = None diff --git a/pysteps/io/nwp_importers.py b/pysteps/io/nwp_importers.py new file mode 100644 index 000000000..680bea081 --- /dev/null +++ b/pysteps/io/nwp_importers.py @@ -0,0 +1,277 @@ +""" +pysteps.io.nwp_importers +==================== + +Methods for importing files containing two-dimensional NWP mosaics. + +The methods in this module implement the following interface:: + + nwp_import_xxx(filename, optional arguments) + +where **xxx** is the name (or abbreviation) of the file format and filename +is the name of the input file. + +The output of each method is a xarray DataArray containing +forecast rainfall fields in mm/timestep and the metadata as attributes. + +The metadata should contain the following recommended key-value pairs: + +.. tabularcolumns:: |p{2cm}|L| + ++------------------+----------------------------------------------------------+ +| Key | Value | ++==================+==========================================================+ +| projection | PROJ.4-compatible projection definition | ++------------------+----------------------------------------------------------+ +| x1 | x-coordinate of the lower-left corner of the data raster | ++------------------+----------------------------------------------------------+ +| y1 | y-coordinate of the lower-left corner of the data raster | ++------------------+----------------------------------------------------------+ +| x2 | x-coordinate of the upper-right corner of the data raster| ++------------------+----------------------------------------------------------+ +| y2 | y-coordinate of the upper-right corner of the data raster| ++------------------+----------------------------------------------------------+ +| xpixelsize | grid resolution in x-direction | ++------------------+----------------------------------------------------------+ +| ypixelsize | grid resolution in y-direction | ++------------------+----------------------------------------------------------+ +| cartesian_unit | the physical unit of the cartesian x- and y-coordinates: | +| | e.g. 'm' or 'km' | ++------------------+----------------------------------------------------------+ +| yorigin | a string specifying the location of the first element in | +| | the data raster w.r.t. y-axis: | +| | 'upper' = upper border | +| | 'lower' = lower border | ++------------------+----------------------------------------------------------+ +| institution | name of the institution who provides the data | ++------------------+----------------------------------------------------------+ +| unit | the physical unit of the data: 'mm/h', 'mm' or 'dBZ' | ++------------------+----------------------------------------------------------+ +| transform | the transformation of the data: None, 'dB', 'Box-Cox' or | +| | others | ++------------------+----------------------------------------------------------+ +| accutime | the accumulation time in minutes of the data, float | ++------------------+----------------------------------------------------------+ +| threshold | the rain/no rain threshold with the same unit, | +| | transformation and accutime of the data. | ++------------------+----------------------------------------------------------+ +| zerovalue | the value assigned to the no rain pixels with the same | +| | unit, transformation and accutime of the data. | ++------------------+----------------------------------------------------------+ +| zr_a | the Z-R constant a in Z = a*R**b | ++------------------+----------------------------------------------------------+ +| zr_b | the Z-R exponent b in Z = a*R**b | ++------------------+----------------------------------------------------------+ + +Available Importers +------------------- + +.. autosummary:: + :toctree: ../generated/ + + import_bom_nwp +""" + +import numpy as np +import xarray as xr + +from pysteps.decorators import postprocess_import +from pysteps.exceptions import MissingOptionalDependency + +try: + import netCDF4 + + NETCDF4_IMPORTED = True +except ImportError: + NETCDF4_IMPORTED = False + + +@postprocess_import() +def import_bom_nwp_xr(filename, **kwargs): + """Import a NetCDF with NWP rainfall forecasts regrided to a BoM Rainfields3 + using xarray. + + Parameters + ---------- + filename: str + Name of the file to import. + + {extra_kwargs_doc} + + Returns + ------- + out_da : xr.DataArray + A xarray DataArray containing forecast rainfall fields in mm/timestep + imported from a netcdf, and the metadata. + """ + + if not NETCDF4_IMPORTED: + raise MissingOptionalDependency( + "netCDF4 package is required to import BoM NWP regrided rainfall " + "products but it is not installed" + ) + + ds = _import_bom_nwp_data_xr(filename, **kwargs) + ds_meta = _import_bom_nwp_geodata_xr(ds, **kwargs) + + # rename varname_time (def: time) to t + varname_time = kwargs.get('varname_time', 'time') + ds_meta = ds_meta.rename({varname_time: 't'}) + varname_time = 't' + + # if data variable is named accum_prcp + # it is assumed that NWP rainfall data is accumulated + # so it needs to be disagregated by time step + varname = kwargs.get('varname', 'accum_prcp') + if varname == 'accum_prcp': + print("Rainfall values are accumulated. Disagreagating by time step") + accum_prcp = ds_meta[varname] + precipitation = (accum_prcp - accum_prcp.shift({varname_time: 1})) + precipitation = precipitation.dropna(varname_time, 'all') + # update/copy attributes + precipitation.name = 'precipitation' + # copy attributes + precipitation.attrs.update({**accum_prcp.attrs}) + # update attributes + precipitation.attrs.update({'standard_name': "precipitation_amount"}) + else: + precipitation = ds_meta[varname] + + return precipitation + + +def _import_bom_nwp_data_xr(filename, **kwargs): + + varname_time = kwargs.get('varname_time', 'time') + chunks = kwargs.get('chunks', {varname_time: 1}) + + ds_rainfall = xr.open_mfdataset( + filename, + combine='nested', + concat_dim=varname_time, + chunks=chunks, + lock=False, + parallel=True, + ) + + return ds_rainfall + + +def _import_bom_nwp_geodata_xr(ds_in, + **kwargs, + ): + + varname = kwargs.get('varname', 'accum_prcp') + varname_time = kwargs.get('varname_time', 'time') + + # extract useful information + # projection + projdef = None + if "proj" in ds_in: + projection = ds_in.proj + if projection.grid_mapping_name == "albers_conical_equal_area": + projdef = "+proj=aea " + lon_0 = projection.longitude_of_central_meridian + projdef += f" +lon_0={lon_0:.3f}" + lat_0 = projection.latitude_of_projection_origin + projdef += f" +lat_0={lat_0:.3f}" + standard_parallel = projection.standard_parallel + projdef += f" +lat_1={standard_parallel[0]:.3f}" + projdef += f" +lat_2={standard_parallel[1]:.3f}" + + # get the accumulation period + time = ds_in[varname_time] + # shift the array to calculate the time step + delta_time = time - time.shift({varname_time: 1}) + # assuming first valid delta_time is representative of all time steps + time_step = delta_time[1] + time_step = time_step.values.astype('timedelta64[m]') + + # get the units of precipitation + units = None + if "units" in ds_in[varname].attrs: + units = ds_in[varname].units + if units in ("kg m-2", "mm"): + units = "mm" + ds_in[varname].attrs.update({'units': units}) + # get spatial boundaries and pixelsize + # move to meters if coordiantes in kilometers + if "units" in ds_in.x.attrs: + if ds_in.x.units == "km": + ds_in['x'] = ds_in.x*1000. + ds_in.x.attrs.update({'units': 'm'}) + ds_in['y'] = ds_in.y*1000. + ds_in.y.attrs.update({'units': 'm'}) + + xmin = ds_in.x.min().values + xmax = ds_in.x.max().values + ymin = ds_in.y.min().values + ymax = ds_in.y.max().values + xpixelsize = abs(ds_in.x[1] - ds_in.x[0]) + ypixelsize = abs(ds_in.y[1] - ds_in.y[0]) + + cartesian_unit = ds_in.x.units + + # Add metadata needed by pySTEPS as attrs in X and Y variables + + ds_in.x.attrs.update({ + # TODO: Remove before final 2.0 version + "x1": xmin, + "x2": xmax, + "cartesian_unit": cartesian_unit, + } + ) + + ds_in.y.attrs.update({ + # TODO: Remove before final 2.0 version + "y1": ymin, + "y2": ymax, + "cartesian_unit": cartesian_unit, + } + ) + + # Add metadata needed by pySTEPS as attrs in rainfall variable + da_rainfall = ds_in[varname].isel({varname_time: 0}) + + ds_in[varname].attrs.update( + {"transform": None, + "unit": units, # copy 'units' in 'unit' for legacy reasons + "projection": projdef, + "accutime": time_step, + "zr_a": None, + "zr_b": None, + "zerovalue": np.nanmin(da_rainfall), + "institution": "Commonwealth of Australia, Bureau of Meteorology", + "threshold": _get_threshold_value(da_rainfall.values), + # TODO(_import_bom_rf3_geodata_xr): Remove before final 2.0 version + "yorigin": "upper", + "xpixelsize": xpixelsize.values, + "ypixelsize": ypixelsize.values, + } + ) + + return ds_in + + +def _get_threshold_value(precip): + """ + Get the the rain/no rain threshold with the same unit, transformation and + accutime of the data. + If all the values are NaNs, the returned value is `np.nan`. + Otherwise, np.min(precip[precip > precip.min()]) is returned. + + Returns + ------- + threshold: float + """ + valid_mask = np.isfinite(precip) + if valid_mask.any(): + _precip = precip[valid_mask] + min_precip = _precip.min() + above_min_mask = _precip > min_precip + if above_min_mask.any(): + return np.min(_precip[above_min_mask]) + else: + return min_precip + else: + return np.nan diff --git a/pysteps/pystepsrc b/pysteps/pystepsrc index 7e2b2c9d7..eb1bf8a37 100644 --- a/pysteps/pystepsrc +++ b/pysteps/pystepsrc @@ -89,6 +89,17 @@ "importer_kwargs": { "gzipped": true } + }, + "bom_nwp": { + "root_path": "./nwp/bom", + "path_fmt": "%Y/%m/%d", + "fn_pattern": "%Y%m%d_%H00_regrid_short", + "fn_ext": "nc", + "importer": "bom_nwp", + "timestep": 10, + "importer_kwargs": { + "gzipped": true + } } } } diff --git a/pysteps/tests/test_io_bom_nwp.py b/pysteps/tests/test_io_bom_nwp.py new file mode 100644 index 000000000..e2d5e5d1d --- /dev/null +++ b/pysteps/tests/test_io_bom_nwp.py @@ -0,0 +1,82 @@ +# -*- coding: utf-8 -*- + +import os +import numpy as np + +import xarray as xr +import pytest + +import pysteps +from pysteps.tests.helpers import smart_assert + +pytest.importorskip("netCDF4") + +root_path = pysteps.rcparams.data_sources["bom_nwp"]["root_path"] +rel_path = os.path.join("2020", "10", "31") +filename = os.path.join(root_path, rel_path, "20201031_0000_regrid_short.nc") +data_array_xr = pysteps.io.import_bom_nwp_xr(filename) + +expected_proj = ( + "+proj=aea +lon_0=153.240 +lat_0=-27.718 +lat_1=-26.200 +lat_2=-29.300" +) + + +def test_io_import_bom_nwp_xarray(): + """Test the importer Bom NWP.""" + assert isinstance(data_array_xr, xr.DataArray) + + +def test_io_import_bom_nwp_xarray_shape(): + """Test the importer Bom NWP shape.""" + assert isinstance(data_array_xr, xr.DataArray) + assert data_array_xr.shape == (144, 512, 512) + + +test_attrs_xr = [ + ("projection", expected_proj, None), + ("institution", "Commonwealth of Australia, Bureau of Meteorology", None), + ("transform", None, None), + ("zerovalue", 0.0, 0.1), + ("unit", "mm", None), + ("accutime", np.timedelta64(10, 'm'), None), + ("zr_a", None, None), + ("zr_b", None, None), + ("xpixelsize", 500.0, 0.1), + ("ypixelsize", 500.0, 0.1), + ("units", "mm", None), + ("yorigin", "upper", None), +] + + +@pytest.mark.parametrize("variable, expected, tolerance", test_attrs_xr) +def test_io_import_bom_nwp_xarray_attrs(variable, expected, tolerance): + """Test the importer Bom RF3.""" + smart_assert(data_array_xr.attrs[variable], expected, tolerance) + + +test_attrs_xr_coord_x = [ + ("units", "m", None), + ("x1", -127750.0, 0.1), + ("x2", 127750.0, 0.1), +] + + +@pytest.mark.parametrize("variable, expected, tolerance", + test_attrs_xr_coord_x) +def test_io_import_bom_nwp_xarray_attrs_coordx(variable, expected, tolerance): + """Test the importer Bom NWP.""" + smart_assert(data_array_xr.x.attrs[variable], expected, tolerance) + + +test_attrs_xr_coord_y = [ + ("units", "m", None), + ("y1", -127750.0, 0.1), + ("y2", 127750.0, 0.1), +] + + +@pytest.mark.parametrize("variable, expected, tolerance", + test_attrs_xr_coord_y) +def test_io_import_bom_nwp_xarray_attrs_coordy(variable, expected, tolerance): + """Test the importer Bom NWP.""" + smart_assert(data_array_xr.y.attrs[variable], expected, tolerance) diff --git a/pysteps/tests/test_io_bom_rf3.py b/pysteps/tests/test_io_bom_rf3.py index ef10a3624..b276c3a61 100644 --- a/pysteps/tests/test_io_bom_rf3.py +++ b/pysteps/tests/test_io_bom_rf3.py @@ -1,6 +1,7 @@ # -*- coding: utf-8 -*- import os +import numpy as np import xarray as xr import pytest @@ -37,7 +38,7 @@ def test_io_import_bom_rf3_shape(): @pytest.mark.parametrize("variable, expected, tolerance", test_attrs) -def test_io_import_mch_gif_dataset_attrs(variable, expected, tolerance): +def test_io_import_bom_rf3_dataset_attrs(variable, expected, tolerance): """Test the importer Bom RF3.""" smart_assert(data_array.attrs[variable], expected, tolerance) @@ -64,6 +65,72 @@ def test_io_import_bom_rf3_geodata(variable, expected, tolerance): """Test the importer Bom RF3.""" root_path = pysteps.rcparams.data_sources["bom"]["root_path"] rel_path = os.path.join("prcp-cscn", "2", "2018", "06", "16") - filename = os.path.join(root_path, rel_path, "2_20180616_100000.prcp-cscn.nc") + filename = os.path.join(root_path, rel_path, + "2_20180616_100000.prcp-cscn.nc") geodata = pysteps.io.importers._import_bom_rf3_geodata(filename) smart_assert(geodata[variable], expected, tolerance) + + +# TEST XARRAY IMPLEMENTATION +data_array_xr = pysteps.io.import_bom_rf3_xr(filename) + + +def test_io_import_bom_rf3_xarray(): + """Test the importer Bom RF3.""" + assert isinstance(data_array_xr, xr.DataArray) + + +def test_io_import_bom_rf3_xarray_shape(): + """Test the importer Bom RF3.""" + assert isinstance(data_array_xr, xr.DataArray) + assert data_array_xr.shape == (1, 512, 512) + + +test_attrs_xr = [ + ("projection", expected_proj, None), + ("institution", "Commonwealth of Australia, Bureau of Meteorology", None), + ("transform", None, None), + ("zerovalue", 0.0, 0.1), + ("unit", "mm", None), + ("accutime", np.timedelta64(6,'m'), None), + ("zr_a", None, None), + ("zr_b", None, None), + ("xpixelsize", 500.0, 0.1), + ("ypixelsize", 500.0, 0.1), + ("units", "mm", None), + ("yorigin", "upper", None), +] + + +@pytest.mark.parametrize("variable, expected, tolerance", test_attrs_xr) +def test_io_import_bom_rf3_xarray_attrs(variable, expected, tolerance): + """Test the importer Bom RF3.""" + smart_assert(data_array_xr.attrs[variable], expected, tolerance) + + +test_attrs_xr_coord_x = [ + ("units", "m", None), + ("x1", -128000.0, 0.1), + ("x2", 127500.0, 0.1), +] + + +@pytest.mark.parametrize("variable, expected, tolerance", + test_attrs_xr_coord_x) +def test_io_import_bom_rf3_xarray_attrs_coordx(variable, expected, tolerance): + """Test the importer Bom RF3.""" + smart_assert(data_array_xr.x.attrs[variable], expected, tolerance) + + +test_attrs_xr_coord_y = [ + ("units", "m", None), + ("y1", -127500.0, 0.1), + ("y2", 128000.0, 0.1), +] + + +@pytest.mark.parametrize("variable, expected, tolerance", + test_attrs_xr_coord_y) +def test_io_import_bom_rf3_xarray_attrs_coordy(variable, expected, tolerance): + """Test the importer Bom RF3.""" + smart_assert(data_array_xr.y.attrs[variable], expected, tolerance) diff --git a/pysteps/tests/test_io_mrms_grib.py b/pysteps/tests/test_io_mrms_grib.py index bb47670ad..afebb7b37 100644 --- a/pysteps/tests/test_io_mrms_grib.py +++ b/pysteps/tests/test_io_mrms_grib.py @@ -25,7 +25,7 @@ def test_io_import_mrms_grib(): assert data_Array.shape == (3500, 7000) assert data_Array.dtype == "single" - expected_attrs = { + expected_metadata = { "institution": "NOAA National Severe Storms Laboratory", "projection": "+proj=longlat +ellps=IAU76", "unit": "mm/h", @@ -37,17 +37,20 @@ def test_io_import_mrms_grib(): "x2": -60.00000199999991, "y1": 20.000001, "y2": 55.00000000000001, - "cartesian_unit": "degrees", + "units": "degrees", } metadata = data_Array.attrs - for key, value in expected_attrs.items(): + metadata.update(**data_Array.x.attrs) + metadata.update(**data_Array.y.attrs) + for key, value in expected_metadata.items(): if isinstance(value, float): - assert_array_almost_equal(metadata[key], expected_attrs[key]) + assert_array_almost_equal(metadata[key], expected_metadata[key]) else: - assert metadata[key] == expected_attrs[key] + assert metadata[key] == expected_metadata[key] x = np.arange(metadata["x1"], metadata["x2"], metadata["xpixelsize"]) y = np.arange(metadata["y1"], metadata["y2"], metadata["ypixelsize"]) + precip_full = data_Array.values assert y.size == precip_full.shape[0] assert x.size == precip_full.shape[1]