diff --git a/climada/hazard/base.py b/climada/hazard/base.py index 963cbe288..12c590039 100644 --- a/climada/hazard/base.py +++ b/climada/hazard/base.py @@ -1651,32 +1651,85 @@ def size(self): """Return number of events.""" return self.event_id.size - def write_raster(self, file_name, intensity=True): - """Write intensity or fraction as GeoTIFF file. Each band is an event + def write_raster(self, file_name, variable='intensity', output_resolution=None): + """Write intensity or fraction as GeoTIFF file. Each band is an event. + Output raster is always a regular grid (same resolution for lat/lon). + + Note that is output_resolution is not None, the data is rasterized to that + resolution. This is an expensive operation. For hazards that are already + a raster, output_resolution=None saves on this raster which is efficient. + + If you want to save both fraction and intensity, create two separate files. + These can then be read together with the from_raster method. Parameters ---------- file_name: str file name to write in tif format - intensity: bool - if True, write intensity, otherwise write fraction + variable: str + if 'intensity', write intensity, if 'fraction' write fraction. + Default is 'intensity' + output_resolution: int + If not None, the data is rasterized to this resolution. + Default is None (resolution is estimated from the data). + + See also + -------- + from_raster: + method to read intensity and fraction raster files. """ - variable = self.intensity - if not intensity: - variable = self.fraction - pixel_geom = self.centroids.calc_pixels_polygons() - profile = self.centroids.meta - profile.update(driver='GTiff', dtype=rasterio.float32, count=self.size) - with rasterio.open(file_name, 'w', **profile) as dst: - LOGGER.info('Writing %s', file_name) - for i_ev in range(variable.shape[0]): - raster = rasterize( - [(x, val) for (x, val) in - zip(pixel_geom, np.array(variable[i_ev, :].toarray()).reshape(-1))], - out_shape=(profile['height'], profile['width']), - transform=profile['transform'], fill=0, - all_touched=True, dtype=profile['dtype'], ) - dst.write(raster.astype(profile['dtype']), i_ev + 1) + + if variable == 'intensity': + var_to_write = self.intensity + elif variable =='fraction': + var_to_write = self.fraction + else: + raise ValueError( + f"The variable {variable} is not valid. Please use 'intensity' or 'fraction'." + ) + + centroids = self.centroids + + res = np.abs(u_coord.get_resolution(centroids.lat, centroids.lon)).min() + xmin, ymin, xmax, ymax = centroids.gdf.total_bounds + + if output_resolution: + res = output_resolution + + # construct raster + rows, cols, ras_trans = u_coord.pts_to_raster_meta( + (xmin, ymin, xmax, ymax), (res, -res) + ) + meta = { + 'crs': self.centroids.crs, + 'height': rows, + 'width': cols, + 'transform': ras_trans, + } + meta.update(driver='GTiff', dtype=rasterio.float32, count=self.size) + + if rows*cols == centroids.shape[0]: + u_coord.write_raster(file_name, var_to_write.toarray(), meta) + else: + geometry = centroids.gdf.geometry.buffer(distance=res/2, resolution=1, cap_style=3) + #resolution=1, cap_style=3: squared buffers + #https://shapely.readthedocs.io/en/latest/manual.html#object.buffer + with rasterio.open(file_name, 'w', **meta) as dst: + LOGGER.info('Writing %s', file_name) + for i_ev in range(self.size): + raster = rasterio.features.rasterize( + ( + (geom, value) + for geom, value + in zip(geometry, var_to_write[i_ev].toarray().flatten()) + ), + out_shape=(meta['height'], meta['width']), + transform=meta['transform'], + fill=0, + all_touched=True, + dtype=meta['dtype'] + ) + dst.write(raster.astype(meta['dtype']), i_ev + 1) def write_hdf5(self, file_name, todense=False): """Write hazard in hdf5 format. diff --git a/climada/test/test_hazard.py b/climada/test/test_hazard.py index 200752678..c55dac5d9 100644 --- a/climada/test/test_hazard.py +++ b/climada/test/test_hazard.py @@ -55,7 +55,7 @@ def test_read_write_raster_pass(self): self.assertEqual(haz_fl.intensity.min(), -9999.0) self.assertAlmostEqual(haz_fl.intensity.max(), 4.662774085998535) - haz_fl.write_raster(DATA_DIR.joinpath('test_write_hazard.tif')) + haz_fl.write_raster(DATA_DIR.joinpath('test_write_hazard.tif'), variable='intensity') haz_read = Hazard.from_raster([DATA_DIR.joinpath('test_write_hazard.tif')]) haz_fl.haz_type = 'FL' @@ -77,53 +77,118 @@ def test_read_raster_pool_pass(self): pool.join() def test_read_write_vector_pass(self): - """Test write_raster: Hazard from vector data""" - haz_fl = Hazard('FL', - event_id=np.array([1]), - date=np.array([1]), - frequency=np.array([1]), - orig=np.array([1]), - event_name=['1'], - intensity=sparse.csr_matrix(np.array([0.5, 0.2, 0.1])), - fraction=sparse.csr_matrix(np.array([0.5, 0.2, 0.1]) / 2), - centroids=Centroids( - latitude=np.array([1, 2, 3]), longitude=np.array([1, 2, 3]), crs=DEF_CRS) - ) - haz_fl.check() + """Test write_raster: Rasterize intensity from vector data""" + haz_fl = Hazard( + 'FL', + event_id=np.array([1]), + date=np.array([1]), + frequency=np.array([1]), + orig=np.array([1]), + event_name=['1'], + intensity=sparse.csr_matrix(np.array([0.11, 0.22, 0.33, 0.31])), + fraction=sparse.csr_matrix(np.array([0, 1, 2, 3]) ), + centroids=Centroids( + longitude=np.array([1, 2, 3, 3]), latitude=np.array([1, 2, 3, 1]), crs=DEF_CRS + ) + ) - haz_fl.write_raster(DATA_DIR.joinpath('test_write_hazard.tif')) + haz_fl.write_raster(DATA_DIR.joinpath('test_write_hazard.tif'), variable='intensity') haz_read = Hazard.from_raster([DATA_DIR.joinpath('test_write_hazard.tif')], haz_type='FL') self.assertEqual(haz_read.intensity.shape, (1, 9)) - self.assertTrue(np.allclose(np.unique(np.array(haz_read.intensity.toarray())), - np.array([0.0, 0.1, 0.2, 0.5]))) - - def test_write_fraction_pass(self): - """Test write_raster with fraction""" - haz_fl = Hazard('FL', - event_id=np.array([1]), - date=np.array([1]), - frequency=np.array([1]), - orig=np.array([1]), - event_name=['1'], - intensity=sparse.csr_matrix(np.array([0.5, 0.2, 0.1])), - fraction=sparse.csr_matrix(np.array([0.5, 0.2, 0.1]) / 2), - centroids=Centroids( - latitude=np.array([1, 2, 3]),longitude=np.array([1, 2, 3]), crs=DEF_CRS) - ) - haz_fl.check() - haz_fl.write_raster(DATA_DIR.joinpath('test_write_hazard.tif'), intensity=False) + output_raster = np.array([ + [1, 3], [2, 3], [3, 3], + [1, 2], [2, 2], [3, 2], + [1, 1], [2, 1], [3, 1] + ]) + output_instensity = np.array([ + 0, 0, 0.33, + 0, 0.22, 0, + 0.11, 0, 0.31 + ]) + + np.testing.assert_array_equal( + haz_read.centroids.lon, + output_raster[:, 0] + ) + np.testing.assert_array_equal( + haz_read.centroids.lat, + output_raster[:, 1] + ) + np.testing.assert_array_almost_equal( + haz_read.intensity.toarray().flatten(), + output_instensity + ) - haz_read = Hazard.from_raster([DATA_DIR.joinpath('test_write_hazard.tif')], - files_fraction=[DATA_DIR.joinpath('test_write_hazard.tif')], - haz_type='FL') - self.assertEqual(haz_read.intensity.shape, (1, 9)) + DATA_DIR.joinpath('test_write_hazard.tif').unlink() + + def test_read_write_vector_fraction_pass(self): + """Test write_raster: Rasterize fraction from vector data""" + haz_fl = Hazard( + 'FL', + event_id=np.array([1]), + date=np.array([1]), + frequency=np.array([1]), + orig=np.array([1]), + event_name=['1'], + intensity=sparse.csr_matrix(np.array([-0.11, -0.22, -0.33, -0.31])), + fraction=sparse.csr_matrix(np.array([0.11, 0.22, 0.33, 0.31])), + centroids=Centroids( + longitude=np.array([1, 2, 3, 3]), latitude=np.array([1, 2, 3, 1]), crs=DEF_CRS + ) + ) + + intensity_file = DATA_DIR.joinpath('test_write_hazard_intensity.tif') + fraction_file = DATA_DIR.joinpath('test_write_hazard_fraction.tif') + + haz_fl.write_raster(fraction_file, variable='fraction') + haz_fl.write_raster(intensity_file, variable='intensity') + + haz_read = Hazard.from_raster( + [intensity_file], [fraction_file], haz_type='FL' + ) self.assertEqual(haz_read.fraction.shape, (1, 9)) - self.assertTrue(np.allclose(np.unique(np.array(haz_read.fraction.toarray())), - np.array([0.0, 0.05, 0.1, 0.25]))) - self.assertTrue(np.allclose(np.unique(np.array(haz_read.intensity.toarray())), - np.array([0.0, 0.05, 0.1, 0.25]))) + self.assertEqual(haz_read.intensity.shape, (1, 9)) + + + output_raster = np.array([ + [1, 3], [2, 3], [3, 3], + [1, 2], [2, 2], [3, 2], + [1, 1], [2, 1], [3, 1] + ]) + output_fraction = np.array([ + 0, 0, 0.33, + 0, 0.22, 0, + 0.11, 0, 0.31 + ]) + + output_intensity = np.array([ + 0, 0, -0.33, + 0, -0.22, 0, + -0.11, 0, -0.31 + ]) + + np.testing.assert_array_equal( + haz_read.centroids.lon, + output_raster[:, 0] + ) + np.testing.assert_array_equal( + haz_read.centroids.lat, + output_raster[:, 1] + ) + np.testing.assert_array_almost_equal( + haz_read.fraction.toarray().flatten(), + output_fraction + ) + np.testing.assert_array_almost_equal( + haz_read.intensity.toarray().flatten(), + output_intensity + ) + + DATA_DIR.joinpath(intensity_file).unlink() + DATA_DIR.joinpath(fraction_file).unlink() + class TestStormEurope(unittest.TestCase):