diff --git a/ra2ce/analyses/direct/analyses_direct.py b/ra2ce/analyses/direct/analyses_direct.py
index ecdad1654..6f0ec6592 100644
--- a/ra2ce/analyses/direct/analyses_direct.py
+++ b/ra2ce/analyses/direct/analyses_direct.py
@@ -35,7 +35,7 @@
 
 class DirectAnalyses:  ### THIS SHOULD ONLY DO COORDINATION
     """
-    Coordination classs for all direct damage analysis
+    Coordination class for all direct damage analysis
 
     Methods of this class are independent modules to do:
      - direct damage analysis
@@ -91,7 +91,7 @@ def execute(self):
 
     def road_damage(self, analysis: dict) -> gpd.GeoDataFrame:
         """
-        ### CONTROLER FOR CALCULATING THE ROAD DAMAGE
+        ### CONTROLLER FOR CALCULATING THE ROAD DAMAGE
 
         Arguments:
             *analysis* (dict) : contains part of the settings from the analysis ini
@@ -247,7 +247,7 @@ def rename_road_gdf_to_conventions(road_gdf_columns):
     ### Handle return period columns
     new_cols = []
     for c in cs:
-        if c.startswith("RP"):
+        if c.startswith("RP") or c.startswith("EV"):
             new_cols.append("F_" + c)
         else:
             new_cols.append(c)
diff --git a/ra2ce/analyses/direct/damage_calculation/damage_network_base.py b/ra2ce/analyses/direct/damage_calculation/damage_network_base.py
index 709cf9b72..2ae330892 100644
--- a/ra2ce/analyses/direct/damage_calculation/damage_network_base.py
+++ b/ra2ce/analyses/direct/damage_calculation/damage_network_base.py
@@ -27,6 +27,7 @@
 import pandas as pd
 
 from ra2ce.analyses.direct.direct_lookup import LookUp as lookup
+from ra2ce.analyses.direct.direct_lookup import dataframe_lookup
 from ra2ce.analyses.direct.direct_utils import (
     clean_lane_data,
     create_summary_statistics,
@@ -229,7 +230,7 @@ def calculate_damage_HZ(self, events):
 
         df = self._gdf_mask
         df["lanes"] = df["lanes"].astype(int)
-        df["max_dam_hz"] = df_max_damages_huizinga.lookup(df["lanes"], df["road_type"])
+        df["max_dam_hz"] = df.apply(dataframe_lookup, args=(df_max_damages_huizinga, ['lanes', 'road_type']), axis=1)
 
         for event in events:
             df["dam_{}_{}".format(event, curve_name)] = round(
@@ -249,7 +250,7 @@ def calculate_damage_HZ(self, events):
 
         self.gdf[dam_cols] = df[dam_cols]
         logging.info(
-            "calculate_damage_HZ(): Damage calculation with the Huizinga damage functions was succesfull"
+            "calculate_damage_HZ(): Damage calculation with the Huizinga damage functions was successful"
         )
 
     def calculate_damage_OSdaMage(self, events):
diff --git a/ra2ce/analyses/direct/direct_lookup.py b/ra2ce/analyses/direct/direct_lookup.py
index 3beda843a..8dc5427c5 100644
--- a/ra2ce/analyses/direct/direct_lookup.py
+++ b/ra2ce/analyses/direct/direct_lookup.py
@@ -23,12 +23,18 @@
 import os
 from collections import OrderedDict
 from pathlib import Path
+from typing import Any
 
 import numpy as np
 import pandas as pd
+from pandas import DataFrame
 from scipy.interpolate import interp1d
 
 
+def dataframe_lookup(row: pd.Series, lookup_df: DataFrame, columns: list) -> Any:
+    row_values = [row[column] for column in columns]
+    return lookup_df.loc[tuple(row_values)]
+
 class LookUp:
     """ " This namespace contains several lookup tables, used e.g. for road damage calculation."""
 
diff --git a/ra2ce/graph/networks_utils.py b/ra2ce/graph/networks_utils.py
index 8049902ef..584132598 100644
--- a/ra2ce/graph/networks_utils.py
+++ b/ra2ce/graph/networks_utils.py
@@ -38,6 +38,7 @@
 from tqdm import tqdm
 from geopy import distance
 from networkx import Graph, set_edge_attributes
+from numpy.ma import MaskedArray
 from osgeo import gdal
 from osmnx import graph_to_gdfs, simplify_graph
 from rasterio.features import shapes
@@ -1831,10 +1832,14 @@ def clean_memory(list_delete: list) -> None:
         del to_delete
 
 
-def get_valid_mean(x_value: float) -> Optional[float]:
-    if not isinstance(x_value, float):
+def get_valid_mean(x_value: MaskedArray, **kwargs) -> Optional[float]:
+    # **kwargs should not be removed. properties var is passed in zonal_stats. So properties should not be removed,
+    #  else this will not be activated
+    if not isinstance(x_value, MaskedArray):
         return np.nan
-    return x_value.mean()  # You know it's a valid type, so return the mean.
+    if x_value.mask.all():
+        return np.nan
+    return np.mean(x_value)  # You know it's a valid type, so return the mean.
 
 
 def buffer_geometry(