Feature/fix rasterio read netcdf (CLIMADA-project#65)

* climada.util.coordinates: slightly increased readability of read_raster

* tests_runner: enable specific test selection

* tests: allow small discrepancy

* conda environment: update rasterio to 1.1.5 - and bear the consequences

* config: suppress FutureWarning

* conda environment: fix cfgrib dependency

* Reduce size of demo file for TCs in Florida

* Reduce size of MAT demo file for TCs in Florida

* Reduce size of MAT demo file for WS in ERA40

* litpop: fix several issues with updated SparseArray and SparseDataFrame

* conda: make sure matplotlib and matplotlib-base have a matching version

* open_street_map: turn FutureWarning and DeprecationWarning off temporarily, to save the logs from being flooded

* Update requirements/env_climada.yml

* Update requirements/env_climada.yml

* dependencies: sync setup.py + cosmetics in env_climada.yml: proj 7.0.0 from the default channel logs with batch echo on

Co-authored-by: schmide <[email protected]>
Co-authored-by: Thomas Vogt <[email protected]>

Author: 3 people

climada/entity/exposures/gpw_import.py

@@ -134,7 +134,7 @@ def get_box_gpw(**parameters):
 
     Returns
     -------
-    tile_temp : pandas SparseArray
+    tile_temp : pandas.arrays.SparseArray
         GPW data
     lon : list
         List with longitudinal infomation on the GPW data. Same
@@ -210,7 +210,7 @@ def get_box_gpw(**parameters):
             LOGGER.error('Error: Matrix has an invalid number of dimensions \
                          (more than 2). Could not continue operation.')
             raise TypeError
-        tile_temp = pd.SparseArray(
+        tile_temp = pd.arrays.SparseArray(
             tile_temp.reshape((tile_temp.size,), order='F'),
             fill_value=0)
         del arr1

climada/entity/exposures/litpop.py

@@ -13,6 +13,7 @@
 
 import os
 import time
+import logging
 from functools import partial
 
 #matplotlib.use('Qt5Agg', force=True)
@@ -562,10 +563,14 @@ def _get_midpoints(highValueArea):
         High_Value_Area_gdf.loc[index, "Midpoint"] = \
         High_Value_Area_gdf.loc[index, "geometry"].centroid.wkt
         s = shape(High_Value_Area_gdf.loc[index, "geometry"])
+        # turn warnings off, otherwise Future and Deprecation warnings are flooding the logs
+        logging.captureWarnings(True)
         proj = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'),
                        pyproj.Proj(init='epsg:3857'))
         High_Value_Area_gdf.loc[index, "projected_area"] = transform(proj, s).area
-
+        # turn warnings on again
+        logging.captureWarnings(False)
+        
     # change active geometry from polygons to midpoints
     from shapely.wkt import loads
     High_Value_Area_gdf = High_Value_Area_gdf.rename(columns={'geometry': 'geo_polys',

climada/entity/exposures/open_street_map.py

@@ -106,7 +106,7 @@ def test_mask_from_shape(self):
         curr_shp = lp._get_country_shape(curr_country, 0)
         mask = lp._mask_from_shape(curr_shp, resolution=60)
         self.assertEqual(mask.sp_index.indices.size, 5591)
-        self.assertTrue(mask.values.max())
+        self.assertTrue(mask.sp_values.max())
         self.assertIn(140 and 7663, mask.sp_index.indices)
 
     def test_litpop_box2coords(self):

climada/entity/exposures/test/test_litpop_unit.py

@@ -144,7 +144,7 @@ def haz_max_events(hazard, min_thresh=0):
     LOGGER.info("Condensing %d hazards to %d max events ...", inten.shape[0], exp_hazards.size)
     inten = inten[exp_hazards]
     inten_max_ids = np.asarray(inten.argmax(axis=1)).ravel()
-    inten_max = inten[range(inten.shape[0]), inten_max_ids].todense()
+    inten_max = inten[range(inten.shape[0]), inten_max_ids]
     dates = hazard.date[exp_hazards]
     dates = [datetime.date.fromordinal(d) for d in dates]
     return pd.DataFrame({

climada/hazard/emulator/stats.py

@@ -178,7 +178,7 @@ def test_gpw_import(self):
         self.assertEqual(len(gpw), 323)
         self.assertIn(np.around(max(gpw)), [103070.0, 137840.0])
         self.assertEqual(type(gpw),
-                         type(pd.SparseArray(data=1, fill_value=0)))
+                         type(pd.arrays.SparseArray(data=1, fill_value=0)))
         self.assertAlmostEqual(lat[0], -27.3164)
         self.assertAlmostEqual(lat[1], 0.083333333)
         self.assertAlmostEqual(lon[0], 30.78291)

climada/test/test_LitPop.py

@@ -55,8 +55,8 @@ def test_switzerland(self):
         self.assertEqual(hazard_set.size, 114)
         self.assertEqual(hazard_set.centroids.size, 130)
         self.assertEqual(exposure_agrar.latitude.values.size, 766 / 2)
-        self.assertEqual(exposure_agrar.value[3], 1720024.4)
-        self.assertEqual(damages_drought, 61995472.555223145)
+        self.assertAlmostEqual(exposure_agrar.value[3], 1720024.4)
+        self.assertAlmostEqual(damages_drought, 61995472.555223145)
 
 
 # Execute Tests

climada/test/test_drought_integr.py

@@ -33,6 +33,7 @@
 
 from climada.util.constants import SOURCE_DIR
 
+
 WORKING_DIR = os.getcwd()
 WINDOWS_END = WORKING_DIR[0:3]
 UNIX_END = '/'

climada/util/config.py

@@ -1093,40 +1093,48 @@ def read_raster(file_name, band=None, src_crs=None, window=None, geometry=None,
 
     with rasterio.Env():
         with rasterio.open(file_name, 'r') as src:
-            src_crs = src.crs if src_crs is None else src_crs
-            if not src_crs:
-                src_crs = rasterio.crs.CRS.from_dict(DEF_CRS)
             dst_meta = src.meta.copy()
 
             if dst_crs or transform:
                 LOGGER.debug('Reprojecting ...')
+
+                src_crs = src.crs if src_crs is None else src_crs
+                if not src_crs:
+                    src_crs = rasterio.crs.CRS.from_dict(DEF_CRS)
                 transform = (transform, width, height) if transform else None
                 inten = _read_raster_reproject(src, src_crs, dst_meta, band=band,
                                                geometry=geometry, dst_crs=dst_crs,
                                                transform=transform, resampling=resampling)
             else:
-                trans = dst_meta['transform']
                 if geometry:
                     inten, trans = rasterio.mask.mask(src, geometry, crop=True, indexes=band)
                     if dst_meta['nodata'] and np.isnan(dst_meta['nodata']):
                         inten[np.isnan(inten)] = 0
                     else:
                         inten[inten == dst_meta['nodata']] = 0
+
                 else:
                     masked_array = src.read(band, window=window, masked=True)
                     inten = masked_array.data
                     inten[masked_array.mask] = 0
+
                     if window:
                         trans = rasterio.windows.transform(window, src.transform)
+                    else:
+                        trans = dst_meta['transform']
+
                 dst_meta.update({
                     "height": inten.shape[1],
                     "width": inten.shape[2],
                     "transform": trans,
                 })
+
     if not dst_meta['crs']:
         dst_meta['crs'] = rasterio.crs.CRS.from_dict(DEF_CRS)
+
     intensity = inten[range(len(band)), :]
     dst_shape = (len(band), dst_meta['height'] * dst_meta['width'])
+
     return dst_meta, intensity.reshape(dst_shape)
 
 def read_raster_bounds(path, bounds, res=None, bands=None):

climada/util/coordinates.py

@@ -270,7 +270,7 @@ def test_get_land_geometry_all_pass(self):
         """get_land_geometry with all earth."""
         res = get_land_geometry(resolution=110)
         self.assertIsInstance(res, shapely.geometry.multipolygon.MultiPolygon)
-        self.assertEqual(res.area, 21496.99098799273)
+        self.assertAlmostEqual(res.area, 21496.99098799273)
 
     def test_on_land_pass(self):
         """check point on land with 1:50.000.000 resolution."""