Skip to content

Commit

Permalink
Rewrite write raster method for hazard
Browse files Browse the repository at this point in the history
  • Loading branch information
Chahan Kropf committed Nov 22, 2023
1 parent 0cabeac commit 3729248
Show file tree
Hide file tree
Showing 2 changed files with 179 additions and 61 deletions.
93 changes: 73 additions & 20 deletions climada/hazard/base.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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):

Check warning on line 1654 in climada/hazard/base.py

View check run for this annotation

Jenkins - WCR / Pylint

too-many-locals 

LOW:
Too many local variables (19/15)
Raw output 
Used when a function or method has too many local variables.
"""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)

Check warning on line 1701 in climada/hazard/base.py

View check run for this annotation

Jenkins - WCR / Pylint

invalid-unary-operand-type 

HIGH:
bad operand type for unary -: NoneType
Raw output 
no description found
)
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):

Check warning on line 1734 in climada/hazard/base.py

View check run for this annotation

Jenkins - WCR / Pylint

too-complex 

LOW:
'write_hdf5' is too complex. The McCabe rating is 12
Raw output 
no description found
"""Write hazard in hdf5 format.
Expand Down
147 changes: 106 additions & 41 deletions climada/test/test_hazard.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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'
Expand All @@ -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):
Expand Down

0 comments on commit 3729248

Please sign in to comment.