Merge branch 'develop' into feature/restructure_fitfuncs_exceedance

# Conflicts:
#	CHANGELOG.md

Author: emanuel-schmid

.github/workflows/ci.yml

@@ -12,6 +12,7 @@ jobs:
   build-and-test:
     name: 'Core / Unit Test Pipeline'
     runs-on: ubuntu-latest
+    timeout-minutes: 20
     permissions:
       # For publishing results
       checks: write

CHANGELOG.md

@@ -12,9 +12,9 @@ Code freeze date: YYYY-MM-DD
 
 ### Added
 
-- `climada.util.interpolation` module for inter- and extrapolation util functions used in local exceedance intensity and return period functions [#930](https://github.com/CLIMADA-project/climada_python/pull/930)
 - `Hazard.local_exceedance_intensity`, `Hazard.local_return_period` and `Impact.local_exceedance_impact`, that all use the `climada.util.interpolation` module [#918](https://github.com/CLIMADA-project/climada_python/pull/918)
-  
+- `climada.util.interpolation` module for inter- and extrapolation util functions used in local exceedance intensity and return period functions [#930](https://github.com/CLIMADA-project/climada_python/pull/930)
+
 ### Changed
 
 - In `climada.util.plot.geo_im_from_array`, NaNs are plotted in gray while cells with no centroid are not plotted [#929](https://github.com/CLIMADA-project/climada_python/pull/929)
@@ -23,6 +23,9 @@ Code freeze date: YYYY-MM-DD
 
 ### Fixed
 
+- File handles are being closed after reading netcdf files with `climada.hazard` modules [#953](https://github.com/CLIMADA-project/climada_python/pull/953)
+- Avoids a ValueError in the impact calculation for cases with a single exposure point and MDR values of 0, by explicitly removing zeros in `climada.hazard.Hazard.get_mdr` [#933](https://github.com/CLIMADA-project/climada_python/pull/948)
+
 ### Deprecated
 
 ### Removed
@@ -471,4 +474,3 @@ updated:
 
 - `climada.enginge.impact.Impact.calc()` and `climada.enginge.impact.Impact.calc_impact_yearset()`
 [#436](https://github.com/CLIMADA-project/climada_python/pull/436).
-

climada/engine/test/test_impact_calc.py

@@ -28,8 +28,8 @@
 
 from climada import CONFIG
 from climada.entity.entity_def import Entity
-from climada.entity import Exposures, ImpactFuncSet, ImpactFunc
-from climada.hazard.base import Hazard
+from climada.entity import Exposures, ImpactFuncSet, ImpactFunc, ImpfTropCyclone
+from climada.hazard.base import Hazard, Centroids
 from climada.engine import ImpactCalc, Impact
 from climada.engine.impact_calc import LOGGER as ILOG
 from climada.util.constants import ENT_DEMO_TODAY, DEMO_DIR
@@ -471,6 +471,34 @@ def test_stitch_risk_metrics(self):
         np.testing.assert_array_equal(eai_exp, [2.25, 1.25, 4.5])
         self.assertEqual(aai_agg, 8.0)  # Sum of eai_exp
 
+    def test_single_exp_zero_mdr(self):
+        """Test for case where exposure has a single value and MDR or fraction contains zeros"""
+        centroids = Centroids.from_lat_lon([-26.16], [28.20])
+        haz = Hazard(
+            intensity=sparse.csr_matrix(np.array([[31.5], [19.0]])),
+            event_id=np.arange(2),
+            event_name=[0,1],
+            frequency=np.ones(2) / 2,
+            fraction=sparse.csr_matrix(np.zeros((2,1))),
+            date=np.array([0, 1]),
+            centroids=centroids,
+            haz_type='TC'
+        )
+        exp = Exposures({'value': [1.],
+                            'longitude': 28.22,
+                            'latitude': -26.17,
+                            'impf_TC': 1},
+                            crs="EPSG:4326")
+        imp_evt = 0.00250988804927603
+        aai_agg = imp_evt/2
+        eai_exp = np.array([aai_agg])
+        at_event = np.array([imp_evt, 0])
+        exp.set_geometry_points()
+        impf_tc = ImpfTropCyclone.from_emanuel_usa()
+        impf_set = ImpactFuncSet([impf_tc])
+        impf_set.check()
+        imp = ImpactCalc(exp, impf_set, haz).impact(save_mat=True)
+        check_impact(self, imp, haz, exp, aai_agg, eai_exp, at_event, at_event)
 
 class TestImpactMatrixCalc(unittest.TestCase):
     """Verify the computation of the impact matrix"""

climada/hazard/base.py

@@ -1061,7 +1061,9 @@ def get_mdr(self, cent_idx, impf):
             impf.id)
             mdr_array = impf.calc_mdr(mdr.toarray().ravel()).reshape(mdr.shape)
             mdr = sparse.csr_matrix(mdr_array)
-        return mdr[:, indices]
+        mdr_out = mdr[:, indices]
+        mdr_out.eliminate_zeros()
+        return mdr_out
 
     def get_paa(self, cent_idx, impf):
         """

climada/hazard/isimip_data.py

@@ -50,12 +50,12 @@ def _read_one_nc(file_name, bbox=None, years=None):
         Contains data in the specified bounding box and for the
         specified time period
     """
-    data = xr.open_dataset(file_name, decode_times=False)
-    if not bbox:
-        bbox = bbox_world
-    if not years:
-        return data.sel(lat=slice(bbox[3], bbox[1]), lon=slice(bbox[0], bbox[2]))
-
-    time_id = years - int(data['time'].units[12:16])
-    return data.sel(lat=slice(bbox[3], bbox[1]), lon=slice(bbox[0], bbox[2]),
-                    time=slice(time_id[0], time_id[1]))
+    with xr.open_dataset(file_name, decode_times=False) as data:
+        if not bbox:
+            bbox = bbox_world
+        if not years:
+            return data.sel(lat=slice(bbox[3], bbox[1]), lon=slice(bbox[0], bbox[2]))
+
+        time_id = years - int(data['time'].units[12:16])
+        return data.sel(lat=slice(bbox[3], bbox[1]), lon=slice(bbox[0], bbox[2]),
+                        time=slice(time_id[0], time_id[1]))

climada/hazard/storm_europe.py

@@ -228,37 +228,33 @@ def _read_one_nc(cls, file_name, centroids, intensity_thres):
         new_haz : StormEurope
             Hazard instance for one single storm.
         """
-        ncdf = xr.open_dataset(file_name)
-
-        if centroids.size != (ncdf.sizes['latitude'] * ncdf.sizes['longitude']):
-            ncdf.close()
-            LOGGER.warning(('Centroids size doesn\'t match NCDF dimensions. '
-                            'Omitting file %s.'), file_name)
-            return None
-
-        # xarray does not penalise repeated assignments, see
-        # http://xarray.pydata.org/en/stable/data-structures.html
-        stacked = ncdf['max_wind_gust'].stack(
-            intensity=('latitude', 'longitude', 'time')
-        )
-        stacked = stacked.where(stacked > intensity_thres)
-        stacked = stacked.fillna(0)
-
-        # fill in values from netCDF
-        ssi_wisc = np.array([float(ncdf.attrs['ssi'])])
-        intensity = sparse.csr_matrix(stacked)
-        new_haz = cls(ssi_wisc=ssi_wisc,
-                      intensity=intensity,
-                      event_name=[ncdf.attrs['storm_name']],
-                      date=np.array([datetime64_to_ordinal(ncdf['time'].data[0])]),
-                      # fill in default values
-                      centroids=centroids,
-                      event_id=np.array([1]),
-                      frequency=np.array([1]),
-                      orig=np.array([True]),)
-
-        ncdf.close()
-        return new_haz
+        with xr.open_dataset(file_name) as ncdf:
+            if centroids.size != (ncdf.sizes['latitude'] * ncdf.sizes['longitude']):
+                LOGGER.warning(('Centroids size doesn\'t match NCDF dimensions. '
+                                'Omitting file %s.'), file_name)
+                return None
+
+            # xarray does not penalise repeated assignments, see
+            # http://xarray.pydata.org/en/stable/data-structures.html
+            stacked = ncdf['max_wind_gust'].stack(
+                intensity=('latitude', 'longitude', 'time')
+            )
+            stacked = stacked.where(stacked > intensity_thres)
+            stacked = stacked.fillna(0)
+
+            # fill in values from netCDF
+            ssi_wisc = np.array([float(ncdf.attrs['ssi'])])
+            intensity = sparse.csr_matrix(stacked)
+            new_haz = cls(ssi_wisc=ssi_wisc,
+                          intensity=intensity,
+                          event_name=[ncdf.attrs['storm_name']],
+                          date=np.array([datetime64_to_ordinal(ncdf['time'].data[0])]),
+                          # fill in default values
+                          centroids=centroids,
+                          event_id=np.array([1]),
+                          frequency=np.array([1]),
+                          orig=np.array([True]),)
+            return new_haz
 
     def read_cosmoe_file(self, *args, **kwargs):
         """This function is deprecated, use StormEurope.from_cosmoe_file instead."""
@@ -309,69 +305,66 @@ def from_cosmoe_file(cls, fp_file, run_datetime, event_date=None,
         intensity_thres = cls.intensity_thres if intensity_thres is None else intensity_thres
 
         # read intensity from file
-        ncdf = xr.open_dataset(fp_file)
-        ncdf = ncdf.assign_coords(date=('time',ncdf["time"].dt.floor("D").values))
-
-        if event_date:
-            try:
-                stacked = ncdf.sel(
-                    time=event_date.strftime('%Y-%m-%d')
-                    ).groupby('date').max().stack(intensity=('y_1', 'x_1'))
-            except KeyError as ker:
-                raise ValueError('Extraction of date and coordinates failed. This is most likely '
-                                 'because the selected event_date '
-                                 f'{event_date.strftime("%Y-%m-%d")} is not contained in the '
-                                 'weather forecast selected by fp_file {fp_file}. Please adjust '
-                                 f'event_date or fp_file.') from ker
-            considered_dates = np.datetime64(event_date)
-        else:
-            time_covered_step = ncdf['time'].diff('time')
-            time_covered_day = time_covered_step.groupby('date').sum()
-            # forecast run should cover at least 18 hours of a day
-            considered_dates_bool = time_covered_day >= np.timedelta64(18,'h')
-            stacked = ncdf.groupby('date').max()\
-                          .sel(date=considered_dates_bool)\
-                          .stack(intensity=('y_1', 'x_1'))
-            considered_dates = stacked['date'].values
-        stacked = stacked.stack(date_ensemble=('date', 'epsd_1'))
-        stacked = stacked.where(stacked['VMAX_10M'] > intensity_thres)
-        stacked = stacked.fillna(0)
-
-        # fill in values from netCDF
-        intensity = sparse.csr_matrix(stacked['VMAX_10M'].T)
-        event_id = np.arange(stacked['date_ensemble'].size) + 1
-        date = np.repeat(
-            np.array(datetime64_to_ordinal(considered_dates)),
-            np.unique(ncdf['epsd_1']).size
-        )
-        orig = np.full_like(event_id, False)
-        orig[(stacked['epsd_1'] == 0).values] = True
-        if description is None:
-            description = (model_name +
-                           ' weather forecast windfield ' +
-                           'for run startet at ' +
-                           run_datetime.strftime('%Y%m%d%H'))
-
-        # Create Hazard
-        haz = cls(
-            intensity=intensity,
-            event_id=event_id,
-            centroids = cls._centroids_from_nc(fp_file),
-            # fill in default values
-            orig=orig,
-            date=date,
-            event_name=[date_i + '_ens' + str(ens_i)
-                        for date_i, ens_i
-                        in zip(date_to_str(date), stacked['epsd_1'].values + 1)],
-            frequency=np.divide(
-                np.ones_like(event_id),
-                np.unique(ncdf['epsd_1']).size),
-        )
+        with xr.open_dataset(fp_file) as ncdf:
+            ncdf = ncdf.assign_coords(date=('time',ncdf["time"].dt.floor("D").values))
+            if event_date:
+                try:
+                    stacked = ncdf.sel(
+                        time=event_date.strftime('%Y-%m-%d')
+                        ).groupby('date').max().stack(intensity=('y_1', 'x_1'))
+                except KeyError as ker:
+                    raise ValueError('Extraction of date and coordinates failed. This is most likely '
+                                     'because the selected event_date '
+                                     f'{event_date.strftime("%Y-%m-%d")} is not contained in the '
+                                     'weather forecast selected by fp_file {fp_file}. Please adjust '
+                                     f'event_date or fp_file.') from ker
+                considered_dates = np.datetime64(event_date)
+            else:
+                time_covered_step = ncdf['time'].diff('time')
+                time_covered_day = time_covered_step.groupby('date').sum()
+                # forecast run should cover at least 18 hours of a day
+                considered_dates_bool = time_covered_day >= np.timedelta64(18,'h')
+                stacked = ncdf.groupby('date').max()\
+                              .sel(date=considered_dates_bool)\
+                              .stack(intensity=('y_1', 'x_1'))
+                considered_dates = stacked['date'].values
+            stacked = stacked.stack(date_ensemble=('date', 'epsd_1'))
+            stacked = stacked.where(stacked['VMAX_10M'] > intensity_thres)
+            stacked = stacked.fillna(0)
+
+            # fill in values from netCDF
+            intensity = sparse.csr_matrix(stacked['VMAX_10M'].T)
+            event_id = np.arange(stacked['date_ensemble'].size) + 1
+            date = np.repeat(
+                np.array(datetime64_to_ordinal(considered_dates)),
+                np.unique(ncdf['epsd_1']).size
+            )
+            orig = np.full_like(event_id, False)
+            orig[(stacked['epsd_1'] == 0).values] = True
+            if description is None:
+                description = (model_name +
+                               ' weather forecast windfield ' +
+                               'for run startet at ' +
+                               run_datetime.strftime('%Y%m%d%H'))
+
+            # Create Hazard
+            haz = cls(
+                intensity=intensity,
+                event_id=event_id,
+                centroids = cls._centroids_from_nc(fp_file),
+                # fill in default values
+                orig=orig,
+                date=date,
+                event_name=[date_i + '_ens' + str(ens_i)
+                            for date_i, ens_i
+                            in zip(date_to_str(date), stacked['epsd_1'].values + 1)],
+                frequency=np.divide(
+                    np.ones_like(event_id),
+                    np.unique(ncdf['epsd_1']).size),
+            )
 
-        # close netcdf file
-        ncdf.close()
-        haz.check()
-        return haz
+            haz.check()
+            return haz
 
     def read_icon_grib(self, *args, **kwargs):
         """This function is deprecated, use StormEurope.from_icon_grib instead."""
@@ -444,11 +437,12 @@ def from_icon_grib(cls, run_datetime, event_date=None, model_name='icon-eu-eps',
             gripfile_path_i = Path(file_i[:-4])
             with open(file_i, 'rb') as source, open(gripfile_path_i, 'wb') as dest:
                 dest.write(bz2.decompress(source.read()))
-            ds_i = xr.open_dataset(gripfile_path_i, engine='cfgrib')
-            if ind_i == 0:
-                stacked = ds_i
-            else:
-                stacked = xr.concat([stacked,ds_i], 'valid_time')
+
+            with xr.open_dataset(gripfile_path_i, engine='cfgrib') as ds_i:
+                if ind_i == 0:
+                    stacked = ds_i
+                else:
+                    stacked = xr.concat([stacked,ds_i], 'valid_time')
 
         # create intensity matrix with max for each full day
         stacked = stacked.assign_coords(
@@ -524,35 +518,34 @@ def _centroids_from_nc(file_name):
         'longitude' variables in a netCDF file.
         """
         LOGGER.info('Constructing centroids from %s', file_name)
-        ncdf = xr.open_dataset(file_name)
-        create_meshgrid = True
-        if hasattr(ncdf, 'latitude'):
-            lats = ncdf['latitude'].data
-            lons = ncdf['longitude'].data
-        elif hasattr(ncdf, 'lat'):
-            lats = ncdf['lat'].data
-            lons = ncdf['lon'].data
-        elif hasattr(ncdf, 'lat_1'):
-            if len(ncdf['lon_1'].shape)>1 & \
-                (ncdf['lon_1'].shape == ncdf['lat_1'].shape) \
-                :
-                lats = ncdf['lat_1'].data.flatten()
-                lons = ncdf['lon_1'].data.flatten()
+        with xr.open_dataset(file_name) as ncdf:
+            create_meshgrid = True
+            if hasattr(ncdf, 'latitude'):
+                lats = ncdf['latitude'].data
+                lons = ncdf['longitude'].data
+            elif hasattr(ncdf, 'lat'):
+                lats = ncdf['lat'].data
+                lons = ncdf['lon'].data
+            elif hasattr(ncdf, 'lat_1'):
+                if len(ncdf['lon_1'].shape)>1 & \
+                    (ncdf['lon_1'].shape == ncdf['lat_1'].shape) \
+                    :
+                    lats = ncdf['lat_1'].data.flatten()
+                    lons = ncdf['lon_1'].data.flatten()
+                    create_meshgrid = False
+                else:
+                    lats = ncdf['lat_1'].data
+                    lons = ncdf['lon_1'].data
+            elif hasattr(ncdf, 'clat'):
+                lats = ncdf['clat'].data
+                lons = ncdf['clon'].data
+                if ncdf['clat'].attrs['units']=='radian':
+                    lats = np.rad2deg(lats)
+                    lons = np.rad2deg(lons)
                 create_meshgrid = False
             else:
-                lats = ncdf['lat_1'].data
-                lons = ncdf['lon_1'].data
-        elif hasattr(ncdf, 'clat'):
-            lats = ncdf['clat'].data
-            lons = ncdf['clon'].data
-            if ncdf['clat'].attrs['units']=='radian':
-                lats = np.rad2deg(lats)
-                lons = np.rad2deg(lons)
-            create_meshgrid = False
-        else:
-            raise AttributeError('netcdf file has no field named latitude or '
-                                 'other know abrivation for coordinates.')
-        ncdf.close()
+                raise AttributeError('netcdf file has no field named latitude or '
+                                     'other know abrivation for coordinates.')
 
         if create_meshgrid:
             lats, lons = np.array([np.repeat(lats, len(lons)),