diff --git a/Makefile b/Makefile
index 29b4aee4d993da53e295cd74df271e6b3fd48d13..5775cde6b5af0a464d029f14fff4a2dc58e316cb 100644
--- a/Makefile
+++ b/Makefile
@@ -1,4 +1,4 @@
-SOURCES=shakyground2 *.py
+SOURCES=shakyground2 tests *.py
 LENGTH=96
 
 check: $(SOURCES)
diff --git a/shakyground2/earthquake.py b/shakyground2/earthquake.py
index f4b293b607aa32a5ae987d1748ea489ac9ed367c..7dc90b89cd4e01b89748fae554fbb838d1324cf0 100644
--- a/shakyground2/earthquake.py
+++ b/shakyground2/earthquake.py
@@ -204,7 +204,7 @@ class Earthquake(object):
         width: float,
         lsd: float,
         usd: float = 0.0,
-        hypo_loc: Tuple[float, float] = (0.5, 0.5)
+        hypo_loc: Tuple[float, float] = (0.5, 0.5),
     ):
         """
         From a set of rupture properties returns a planar surface whose dimensions
@@ -266,12 +266,8 @@ class Earthquake(object):
         updip_dir = (dip - 90.0) % 360
         mid_left = centroid.point_at(left_length, 0.0, (strike + 180.0) % 360.0)
         mid_right = centroid.point_at(right_length, 0.0, strike)
-        top_left = mid_left.point_at(
-            updip_surface_length, -updip_depth_change, updip_dir
-        )
-        top_right = mid_right.point_at(
-            updip_surface_length, -updip_depth_change, updip_dir
-        )
+        top_left = mid_left.point_at(updip_surface_length, -updip_depth_change, updip_dir)
+        top_right = mid_right.point_at(updip_surface_length, -updip_depth_change, updip_dir)
         bottom_left = mid_left.point_at(
             downdip_surface_length, downdip_depth_change, downdip_dir
         )
@@ -288,7 +284,7 @@ class Earthquake(object):
             extended_error_message = [
                 "Rupture surface failed to build with the following properties:",
                 "Strike: %.2f, Dip: %.2f, Length: %.2f, Width: %.2f, Hypo. Pos: %s"
-                % (strike, dip, length, width, str(hypo_loc))
+                % (strike, dip, length, width, str(hypo_loc)),
             ]
             for pnt in [top_left, top_right, bottom_right, bottom_left]:
                 extended_error_message.append(str(pnt))
diff --git a/shakyground2/magnitude_scaling_relations.py b/shakyground2/magnitude_scaling_relations.py
index d44045255f415f07d4d1009c85c6390a37aec4b0..f7656262f7c45df953d3df0e03b2e6d65f9cfb69 100644
--- a/shakyground2/magnitude_scaling_relations.py
+++ b/shakyground2/magnitude_scaling_relations.py
@@ -12,10 +12,17 @@ class BaseScalingRelation(metaclass=abc.ABCMeta):
     Abstract base class representing an implementation of a magnitude scaling relation
 
     """
+
     @abc.abstractmethod
-    def get_rupture_dimensions(self, magnitude: float, rake: float = 0.0, dip: float = 90.0,
-                               aspect: float = 1.0, thickness: float = 20.0,
-                               epsilon: float = 0.0) -> Tuple[float, float, float]:
+    def get_rupture_dimensions(
+        self,
+        magnitude: float,
+        rake: float = 0.0,
+        dip: float = 90.0,
+        aspect: float = 1.0,
+        thickness: float = 20.0,
+        epsilon: float = 0.0,
+    ) -> Tuple[float, float, float]:
         """
         Args:
             magnitude:
@@ -41,9 +48,15 @@ class PEERScalingRelation(BaseScalingRelation):
     updated by the synthetic rupture generator
     """
 
-    def get_rupture_dimensions(self, magnitude: float, rake: float = 0.0, dip: float = 90.0,
-                               aspect: float = 1.0, thickness: float = 20.0,
-                               epsilon: float = 0.0) -> Tuple[float, float, float]:
+    def get_rupture_dimensions(
+        self,
+        magnitude: float,
+        rake: float = 0.0,
+        dip: float = 90.0,
+        aspect: float = 1.0,
+        thickness: float = 20.0,
+        epsilon: float = 0.0,
+    ) -> Tuple[float, float, float]:
         """
         Returns the area (km^2), length (km) and width (km) of the rupture using the PEER
         Magnitude-Area scaling relation and constrained by the crustal thickness
@@ -80,6 +93,7 @@ class StrasserEtAl2010Interface(PEERScalingRelation):
     Interface and Intraslab Subduction-zone Earthquakes with Moment Magnitude", Seismological
     Research Letters, 81: 941 - 950, doi:10.1785/gssrl.81.6.941
     """
+
     @staticmethod
     def get_area(magnitude: float) -> float:
         """
@@ -97,6 +111,7 @@ class StrasserEtAl2010Inslab(PEERScalingRelation):
     Interface and Intraslab Subduction-zone Earthquakes with Moment Magnitude", Seismological
     Research Letters, 81: 941 - 950, doi:10.1785/gssrl.81.6.941
     """
+
     @staticmethod
     def get_area(magnitude: float) -> float:
         """
@@ -117,27 +132,39 @@ class Stafford2014(BaseScalingRelation):
     # Model coefficients for the specific style of faulting
     COEFFS = {
         # Coefficients from Table 1
-        1: {"U": {"a0": -27.4922, "a1": 4.6656, "a2": -0.2033},
+        1: {
+            "U": {"a0": -27.4922, "a1": 4.6656, "a2": -0.2033},
             "SS": {"a0": -30.8395, "a1": 5.4184, "a2": -0.3044},
             "N": {"a0": -36.9770, "a1": 6.3070, "a2": -0.1696},
-            "R": {"a0": -35.8239, "a1": 5.0680, "a2": -0.0457}},
+            "R": {"a0": -35.8239, "a1": 5.0680, "a2": -0.0457},
+        },
         # Coefficients from Table 2
-        2: {"U": {"b0": -2.300, "b1": 0.7167, "sigma": 0.2337},
+        2: {
+            "U": {"b0": -2.300, "b1": 0.7167, "sigma": 0.2337},
             "SS": {"b0": -2.300, "b1": 0.7167, "sigma": 0.2337},
             "N": {"b0": -4.1055, "b1": 1.0370, "sigma": 0.2509},
-            "R": {"b0": -3.8300, "b1": 0.9982, "sigma": 0.2285}},
+            "R": {"b0": -3.8300, "b1": 0.9982, "sigma": 0.2285},
+        },
         # Coefficients from Table 3
-        3: {"U": {"gamma0": -9.3137, "sigma": 0.3138, "rho": 0.3104},
+        3: {
+            "U": {"gamma0": -9.3137, "sigma": 0.3138, "rho": 0.3104},
             "SS": {"gamma0": -9.3137, "sigma": 0.3138, "rho": 0.3104},
             "N": {"gamma0": -9.2483, "sigma": 0.3454, "rho": 0.4336},
-            "R": {"gamma0": -9.2749, "sigma": 0.2534, "rho": 0.1376}},
+            "R": {"gamma0": -9.2749, "sigma": 0.2534, "rho": 0.1376},
+        },
         # Coefficients from Table 4
-        4: {"U": 0.7574, "SS": 0.7574, "N": 0.8490, "R": 0.7496}
+        4: {"U": 0.7574, "SS": 0.7574, "N": 0.8490, "R": 0.7496},
     }
 
-    def get_rupture_dimensions(self, magnitude: float, rake: float = 0.0, dip: float = 90.0,
-                               aspect: float = 1.0, thickness: float = 20.0,
-                               epsilon: float = 0.0) -> Tuple[float, float, float]:
+    def get_rupture_dimensions(
+        self,
+        magnitude: float,
+        rake: float = 0.0,
+        dip: float = 90.0,
+        aspect: float = 1.0,
+        thickness: float = 20.0,
+        epsilon: float = 0.0,
+    ) -> Tuple[float, float, float]:
         """
         Gets the rupture dimensions from for the given magnitude subject to the physical
         constaints.
@@ -193,8 +220,9 @@ class Stafford2014(BaseScalingRelation):
             # normal
             return "N"
 
-    def get_rupture_width(self, magnitude: float, dip: float, sof: str, thickness: float) ->\
-            Tuple[float, float, float, float]:
+    def get_rupture_width(
+        self, magnitude: float, dip: float, sof: str, thickness: float
+    ) -> Tuple[float, float, float, float]:
         """
         Returns the parameters needed to define the censored rupture width
         distribution defined in equations 8 to 14
@@ -206,44 +234,45 @@ class Stafford2014(BaseScalingRelation):
         """
         # Gets the probability of a full width rupture
         rw_max = thickness / sin(radians(dip))
-        z_i = self.COEFFS[1][sof]["a0"] + \
-            self.COEFFS[1][sof]["a1"] * magnitude +\
-            self.COEFFS[1][sof]["a2"] * rw_max
+        z_i = (
+            self.COEFFS[1][sof]["a0"]
+            + self.COEFFS[1][sof]["a1"] * magnitude
+            + self.COEFFS[1][sof]["a2"] * rw_max
+        )
         # Probability of rupturing full seismogenic thickness from logistic regression
         p_i = 1.0 / (1.0 + exp(-z_i))
         # Median width from equation 16
-        ln_rw = self.COEFFS[2][sof]["b0"] + \
-            self.COEFFS[2][sof]["b1"] * magnitude
+        ln_rw = self.COEFFS[2][sof]["b0"] + self.COEFFS[2][sof]["b1"] * magnitude
         # Equations 18 - 20
         phi_rw = (log(rw_max) - ln_rw) / self.COEFFS[2][sof]["sigma"]
         phi_rw_ncdf = norm.cdf(phi_rw)
-        ln_rw_trunc = ln_rw - self.COEFFS[2][sof]["sigma"] *\
-            (norm.pdf(phi_rw) / phi_rw_ncdf)
+        ln_rw_trunc = ln_rw - self.COEFFS[2][sof]["sigma"] * (norm.pdf(phi_rw) / phi_rw_ncdf)
         mean_rw = p_i * log(rw_max) + (1.0 - p_i) * ln_rw_trunc
         # Equations 21 - 22
-        stddev_rw = self._get_rupture_width_sigma(self.COEFFS[2][sof]["sigma"],
-                                                  phi_rw,
-                                                  phi_rw_ncdf,
-                                                  p_i)
+        stddev_rw = self._get_rupture_width_sigma(
+            self.COEFFS[2][sof]["sigma"], phi_rw, phi_rw_ncdf, p_i
+        )
         return mean_rw, stddev_rw, rw_max, p_i
 
     @staticmethod
-    def _get_rupture_width_sigma(sigma: float, phi_rw: float, phi_rw_ncdf: float,
-                                 p_i: float) -> float:
+    def _get_rupture_width_sigma(
+        sigma: float, phi_rw: float, phi_rw_ncdf: float, p_i: float
+    ) -> float:
         """
         Returns the variabiliy in the rupture width described by equations 21 and 22
         """
         denom = sqrt(2.0 * pi) * phi_rw_ncdf
         if phi_rw_ncdf >= 0.0:
-            elem1 = sqrt(pi / 2.) * (1.0 + chi2.cdf(phi_rw, 3))
+            elem1 = sqrt(pi / 2.0) * (1.0 + chi2.cdf(phi_rw, 3))
         else:
-            elem1 = sqrt(pi / 2.) * (1.0 - chi2.cdf(phi_rw, 3))
+            elem1 = sqrt(pi / 2.0) * (1.0 - chi2.cdf(phi_rw, 3))
         elem2 = exp(-(phi_rw ** 2)) / denom
-        sigma_truncated = sqrt(((sigma ** 2.) / denom) * (elem1 - elem2))
+        sigma_truncated = sqrt(((sigma ** 2.0) / denom) * (elem1 - elem2))
         return (1.0 - p_i) * sigma_truncated
 
-    def get_rupture_area(self, magnitude: float, sof: str, rw_max: float,
-                         sigma_lnrw: float) -> Tuple[float, float]:
+    def get_rupture_area(
+        self, magnitude: float, sof: str, rw_max: float, sigma_lnrw: float
+    ) -> Tuple[float, float]:
         """
         Returns the rupture area conditioned upon the maximum rupture width and style of
         faulting
@@ -258,15 +287,16 @@ class Stafford2014(BaseScalingRelation):
             ln_ra: Natural logarithm of rupture area
             sigma_ra: Standard deviation of the natural logarithm of rupture area
         """
-        mw_crit = (log(rw_max) - self.COEFFS[2][sof]["b0"]) /\
-            self.COEFFS[2][sof]["b1"]
-        ln_ra = self.COEFFS[3][sof]["gamma0"] + log(10.) * magnitude
+        mw_crit = (log(rw_max) - self.COEFFS[2][sof]["b0"]) / self.COEFFS[2][sof]["b1"]
+        ln_ra = self.COEFFS[3][sof]["gamma0"] + log(10.0) * magnitude
         if magnitude > mw_crit:
             # Equation 23
-            ln_ra -= ((log(10.) / 4.) * (magnitude - mw_crit))
+            ln_ra -= (log(10.0) / 4.0) * (magnitude - mw_crit)
         # Get the sigma log rupture area (equation 28)
         sigma = self.COEFFS[3][sof]["sigma"]
         sigma_ra = sqrt(
-            (sigma ** 2.) + (sigma_lnrw ** 2.) +
-            (2.0 * self.COEFFS[3][sof]["rho"] * sigma_lnrw * sigma))
+            (sigma ** 2.0)
+            + (sigma_lnrw ** 2.0)
+            + (2.0 * self.COEFFS[3][sof]["rho"] * sigma_lnrw * sigma)
+        )
         return ln_ra, sigma_ra
diff --git a/shakyground2/rupture_mechanism.py b/shakyground2/rupture_mechanism.py
index 6ab453ce5075c884350128d7295250aea2c58bbc..03c43e73133ed57e6aae978985141e2f51ba9c07 100644
--- a/shakyground2/rupture_mechanism.py
+++ b/shakyground2/rupture_mechanism.py
@@ -17,6 +17,7 @@ class DIP_RANGES(Enum):
     for each style-of-faulting is provided and they are assumed to be uniformly
     distributed.
     """
+
     R = (20.0, 45.0)
     SS = (75.0, 91.0)
     N = (45.0, 75.0)
@@ -40,6 +41,7 @@ class RuptureMechanism(object):
                    of :class:`openquake.hazardlib.pmf.PMF`
 
     """
+
     def __init__(self, mechanism: Optional[PMF] = None):
         if mechanism:
             assert isinstance(mechanism, PMF)
@@ -76,8 +78,12 @@ class RuptureMechanism(object):
         return strikes, dips, rakes
 
     @classmethod
-    def from_strike_dip_rake(cls, strike: Optional[float] = None, dip: Optional[float] = None,
-                             rake: Optional[float] = None) -> RuptureMechanism:
+    def from_strike_dip_rake(
+        cls,
+        strike: Optional[float] = None,
+        dip: Optional[float] = None,
+        rake: Optional[float] = None,
+    ) -> RuptureMechanism:
         """
         Constructs the rupture mechanism distribution from a simple strike, dip
         and rake combination, permitting None values if undefined
@@ -87,13 +93,14 @@ class RuptureMechanism(object):
             dip: Dip of fault (in decimal degrees between 0 and 90)
             rake: Rake of fault (in decimal degrees between -180 and 180)
         """
-        return cls(cls.build_mechanism_distribution(valid.strike(strike),
-                                                    valid.dip(dip),
-                                                    valid.rake(rake)))
+        return cls(
+            cls.build_mechanism_distribution(
+                valid.strike(strike), valid.dip(dip), valid.rake(rake)
+            )
+        )
 
     @classmethod
-    def from_nodal_planes(cls, nodal_plane_1: Dict, nodal_plane_2: Dict)\
-            -> RuptureMechanism:
+    def from_nodal_planes(cls, nodal_plane_1: Dict, nodal_plane_2: Dict) -> RuptureMechanism:
         """
         Constructs the rupture mechanism distribution from either one single
         valid rupture plane or two valud nodal planes
@@ -108,13 +115,21 @@ class RuptureMechanism(object):
         strike_2 = valid.strike(nodal_plane_2["strike"])
         dip_2 = valid.dip(nodal_plane_2["dip"])
         rake_2 = valid.rake(nodal_plane_2["rake"])
-        return cls(PMF([(0.5, NodalPlane(strike_1, dip_1, rake_1)),
-                        (0.5, NodalPlane(strike_2, dip_2, rake_2))]))
+        return cls(
+            PMF(
+                [
+                    (0.5, NodalPlane(strike_1, dip_1, rake_1)),
+                    (0.5, NodalPlane(strike_2, dip_2, rake_2)),
+                ]
+            )
+        )
 
     @staticmethod
-    def build_mechanism_distribution(strike: Optional[float] = None,
-                                     dip: Optional[float] = None,
-                                     rake: Optional[float] = None) -> PMF:
+    def build_mechanism_distribution(
+        strike: Optional[float] = None,
+        dip: Optional[float] = None,
+        rake: Optional[float] = None,
+    ) -> PMF:
         """
         Builds a mechanism distribution from a partial (or complete) characterisation of the
         rupture mechanism
diff --git a/shakyground2/site_model.py b/shakyground2/site_model.py
index e435c65de09db12169165e498fa4814f78a9aad4..e835e0a71964c9f8c84d06b607d0be3ebe13f147 100644
--- a/shakyground2/site_model.py
+++ b/shakyground2/site_model.py
@@ -23,7 +23,7 @@ SITE_PROPERTIES = {
     "xvf": np.float64,  # Distance (km) to the volcanic front
     "backarc": np.bool,  # Site is in the subduction backarc (True) or else
     "region": np.int32,  # Region to which the site belongs
-    "geology": (np.string_, 20)  # Geological classification for the site
+    "geology": (np.string_, 20),  # Geological classification for the site
 }
 
 
@@ -218,15 +218,11 @@ class SiteModel(object):
             if key not in dframe.columns:
                 # Use the default
                 if SITE_PROPERTIES[key] in ((np.bytes_, 20),):
-                    site_array[key] = cls._get_string_array(
-                        nsites, value, SITE_PROPERTIES[key]
-                    )
+                    site_array[key] = cls._get_string_array(nsites, value, SITE_PROPERTIES[key])
                 elif SITE_PROPERTIES[key] == np.bool:
                     site_array[key] = cls._get_boolean_array(nsites, value)
                 else:
-                    site_array[key] = value * np.ones(
-                        nsites, dtype=SITE_PROPERTIES[key]
-                    )
+                    site_array[key] = value * np.ones(nsites, dtype=SITE_PROPERTIES[key])
         return cls(site_array)
 
     @staticmethod
@@ -277,9 +273,7 @@ class SiteModel(object):
         if japan:
             c1 = 412.0 ** 2.0
             c2 = 1360.0 ** 2.0
-            return np.exp(
-                (-5.23 / 2.0) * np.log((np.power(vs30, 2.0) + c1) / (c2 + c1))
-            )
+            return np.exp((-5.23 / 2.0) * np.log((np.power(vs30, 2.0) + c1) / (c2 + c1)))
         else:
             c1 = 571.0 ** 4.0
             c2 = 1360.0 ** 4.0
diff --git a/shakyground2/valid.py b/shakyground2/valid.py
index ca90c64e050671e3fea78a4a14631a1c0e70921e..07c769fce0ece8be144f5eb5bd71d4b544c836d4 100644
--- a/shakyground2/valid.py
+++ b/shakyground2/valid.py
@@ -4,8 +4,7 @@ consistent quantities
 """
 from typing import Tuple, Optional, Union
 from openquake.hazardlib.geo.nodalplane import NodalPlane
-from shakyground2.magnitude_scaling_relations import (BaseScalingRelation,
-                                                      PEERScalingRelation)
+from shakyground2.magnitude_scaling_relations import BaseScalingRelation, PEERScalingRelation
 
 
 def longitude(lon: float) -> float:
@@ -93,8 +92,10 @@ def seismogenic_thickness(
     usd = positive_float(upper_seismo_depth, "upper seismogenic depth")
     lsd = positive_float(lower_seismo_depth, "lower seismogenic depth")
     if lsd < usd:
-        raise ValueError("Lower seismogenic depth %.2f km shallower than upper seismogenic "
-                         "depth %.2f km" % (lsd, usd))
+        raise ValueError(
+            "Lower seismogenic depth %.2f km shallower than upper seismogenic "
+            "depth %.2f km" % (lsd, usd)
+        )
     return usd, lsd
 
 
@@ -105,8 +106,9 @@ def hypocenter_position(hypo_pos: Tuple[float, float]) -> Tuple[float, float]:
     """
     along_strike, down_dip = hypo_pos
     if along_strike < 0.0 or along_strike > 1.0:
-        raise ValueError("Along strike position %.3f should be in the range 0 to 1"
-                         % along_strike)
+        raise ValueError(
+            "Along strike position %.3f should be in the range 0 to 1" % along_strike
+        )
     if down_dip < 0.0 or down_dip > 1.0:
         raise ValueError("Down dip position %.3f should be in the range 0 to 1" % down_dip)
     return along_strike, down_dip
@@ -121,6 +123,7 @@ def scaling_relation(msr: Optional[BaseScalingRelation]):
         # If no relation is defined then use the default
         return PEERScalingRelation()
     if not isinstance(msr, BaseScalingRelation):
-        raise TypeError("Magnitude Scaling Relation %s not instance of BaseScalingRelation" %
-                        str(msr))
+        raise TypeError(
+            "Magnitude Scaling Relation %s not instance of BaseScalingRelation" % str(msr)
+        )
     return msr
diff --git a/tests/earthquake_test.py b/tests/earthquake_test.py
index 1aaeda4b5e192e7ac73c3f022241ac2e4f99fccc..038303d589f955f606a480d9ee66c2e156d6927f 100644
--- a/tests/earthquake_test.py
+++ b/tests/earthquake_test.py
@@ -11,7 +11,6 @@ from shakyground2.magnitude_scaling_relations import PEERScalingRelation
 
 
 class EarthquakeTestCase(unittest.TestCase):
-
     def setUp(self):
         self.lon = 35.0
         self.lat = 40.0
@@ -27,10 +26,18 @@ class EarthquakeTestCase(unittest.TestCase):
         self.assertEqual(str(eq0), "XYZ (35.00000E, 40.00000N, 10.00 km) M 6.00")
 
         # Minimal instantiation - no datetime
-        eq0 = Earthquake("XYZ", self.lon, self.lat, self.depth, self.mag,
-                         date=date(2010, 10, 6), time=time(11, 35, 45))
-        self.assertEqual(str(eq0),
-                         "XYZ 2010-10-06 11:35:45 (35.00000E, 40.00000N, 10.00 km) M 6.00")
+        eq0 = Earthquake(
+            "XYZ",
+            self.lon,
+            self.lat,
+            self.depth,
+            self.mag,
+            date=date(2010, 10, 6),
+            time=time(11, 35, 45),
+        )
+        self.assertEqual(
+            str(eq0), "XYZ 2010-10-06 11:35:45 (35.00000E, 40.00000N, 10.00 km) M 6.00"
+        )
 
     def test_create_planar_surface_simple(self):
         # Test the creation of a simple planar surface that requires no
@@ -41,10 +48,10 @@ class EarthquakeTestCase(unittest.TestCase):
         dip = 90.0
         usd = 0.0
         lsd = 20.0
-        area, length, width = PEERScalingRelation().get_rupture_dimensions(mag, dip=dip,
-                                                                           thickness=lsd - usd)
-        plane0 = Earthquake.build_planar_surface(centroid, strike, dip,
-                                                 length, width, lsd, usd)
+        area, length, width = PEERScalingRelation().get_rupture_dimensions(
+            mag, dip=dip, thickness=lsd - usd
+        )
+        plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width, lsd, usd)
         # Length should be approximately 10 km, but allow a slight precision difference
         self.assertAlmostEqual(plane0.length, 10.0, 1)
         self.assertAlmostEqual(plane0.width, 10.0, 5)
@@ -66,10 +73,10 @@ class EarthquakeTestCase(unittest.TestCase):
         dip = 90.0
         usd = 0.0
         lsd = 20.0
-        area, length, width = PEERScalingRelation().get_rupture_dimensions(mag, dip=dip,
-                                                                           thickness=lsd - usd)
-        plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width,
-                                                 lsd, usd)
+        area, length, width = PEERScalingRelation().get_rupture_dimensions(
+            mag, dip=dip, thickness=lsd - usd
+        )
+        plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width, lsd, usd)
 
         # Verify correct plane from top left and bottom right points
         self.assertAlmostEqual(plane0.top_left.longitude, -self.target_width / 2.0, 7)
@@ -88,10 +95,10 @@ class EarthquakeTestCase(unittest.TestCase):
         dip = 90.0
         usd = 0.0
         lsd = 10.0
-        area, length, width = PEERScalingRelation().get_rupture_dimensions(mag, dip=dip,
-                                                                           thickness=lsd - usd)
-        plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width,
-                                                 lsd, usd)
+        area, length, width = PEERScalingRelation().get_rupture_dimensions(
+            mag, dip=dip, thickness=lsd - usd
+        )
+        plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width, lsd, usd)
 
         # Verify correct plane from top left and bottom right points
         self.assertAlmostEqual(plane0.top_left.longitude, -self.target_width / 2.0, 7)
@@ -110,17 +117,16 @@ class EarthquakeTestCase(unittest.TestCase):
         dip = 90.0
         usd = 0.0
         lsd = 10.0
-        area, length, width = PEERScalingRelation().get_rupture_dimensions(mag, dip=dip,
-                                                                           thickness=lsd - usd)
-        plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width,
-                                                 lsd, usd)
+        area, length, width = PEERScalingRelation().get_rupture_dimensions(
+            mag, dip=dip, thickness=lsd - usd
+        )
+        plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width, lsd, usd)
 
         # Verify correct plane from top left and bottom right points
         self.assertAlmostEqual(plane0.top_left.longitude, -10.0 * self.target_width / 2.0, 7)
         self.assertAlmostEqual(plane0.top_left.latitude, 0.0, 7)
         self.assertAlmostEqual(plane0.top_left.depth, 0.0, 7)
-        self.assertAlmostEqual(plane0.bottom_right.longitude,
-                               10.0 * self.target_width / 2.0, 7)
+        self.assertAlmostEqual(plane0.bottom_right.longitude, 10.0 * self.target_width / 2.0, 7)
         self.assertAlmostEqual(plane0.bottom_right.latitude, 0.0, 7)
         self.assertAlmostEqual(plane0.bottom_right.depth, 10.0, 7)
 
@@ -133,10 +139,10 @@ class EarthquakeTestCase(unittest.TestCase):
         dip = 60.0
         usd = 0.0
         lsd = 10.0
-        area, length, width = PEERScalingRelation().get_rupture_dimensions(mag, dip=dip,
-                                                                           thickness=lsd - usd)
-        plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width,
-                                                 lsd, usd)
+        area, length, width = PEERScalingRelation().get_rupture_dimensions(
+            mag, dip=dip, thickness=lsd - usd
+        )
+        plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width, lsd, usd)
         self.assertAlmostEqual(width, 10.0 / sin(radians(dip)))
 
         # Verify correct plane from top left and bottom right points - calculated by hand
@@ -154,34 +160,58 @@ class EarthquakeTestCase(unittest.TestCase):
         strike = 0.0
         dip = 90.0
         # Create a simple surface
-        area, length, width = PEERScalingRelation().get_rupture_dimensions(self.mag, dip=dip,
-                                                                           thickness=lsd - usd)
+        area, length, width = PEERScalingRelation().get_rupture_dimensions(
+            self.mag, dip=dip, thickness=lsd - usd
+        )
 
-        plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width,
-                                                 lsd, usd)
+        plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width, lsd, usd)
         # Instantiate earthquake object with known surface
-        eq0 = Earthquake("XYZ", self.lon, self.lat, self.depth, self.mag,
-                         strike=strike, dip=dip, rake=0.0,
-                         upper_seismogenic_depth=usd, lower_seismogenic_depth=lsd,
-                         surface=plane0)
+        eq0 = Earthquake(
+            "XYZ",
+            self.lon,
+            self.lat,
+            self.depth,
+            self.mag,
+            strike=strike,
+            dip=dip,
+            rake=0.0,
+            upper_seismogenic_depth=usd,
+            lower_seismogenic_depth=lsd,
+            surface=plane0,
+        )
         # Generate rupture
         self.assertAlmostEqual(eq0.rupture.mag, self.mag)
         self.assertAlmostEqual(eq0.rupture.rake, 0.0)
-        np.testing.assert_array_almost_equal(plane0.corner_lons,
-                                             eq0.rupture.surface.corner_lons)
-        np.testing.assert_array_almost_equal(plane0.corner_lats,
-                                             eq0.rupture.surface.corner_lats)
-        np.testing.assert_array_almost_equal(plane0.corner_depths,
-                                             eq0.rupture.surface.corner_depths)
+        np.testing.assert_array_almost_equal(
+            plane0.corner_lons, eq0.rupture.surface.corner_lons
+        )
+        np.testing.assert_array_almost_equal(
+            plane0.corner_lats, eq0.rupture.surface.corner_lats
+        )
+        np.testing.assert_array_almost_equal(
+            plane0.corner_depths, eq0.rupture.surface.corner_depths
+        )
         # Instantiate earthquake object without surface
-        eq1 = Earthquake("XYZ", self.lon, self.lat, self.depth, self.mag,
-                         strike=strike, dip=dip, rake=0.0,
-                         upper_seismogenic_depth=usd, lower_seismogenic_depth=lsd)
+        eq1 = Earthquake(
+            "XYZ",
+            self.lon,
+            self.lat,
+            self.depth,
+            self.mag,
+            strike=strike,
+            dip=dip,
+            rake=0.0,
+            upper_seismogenic_depth=usd,
+            lower_seismogenic_depth=lsd,
+        )
         self.assertAlmostEqual(eq1.rupture.mag, self.mag)
         self.assertAlmostEqual(eq1.rupture.rake, 0.0)
-        np.testing.assert_array_almost_equal(plane0.corner_lons,
-                                             eq1.rupture.surface.corner_lons)
-        np.testing.assert_array_almost_equal(plane0.corner_lats,
-                                             eq1.rupture.surface.corner_lats)
-        np.testing.assert_array_almost_equal(plane0.corner_depths,
-                                             eq1.rupture.surface.corner_depths)
+        np.testing.assert_array_almost_equal(
+            plane0.corner_lons, eq1.rupture.surface.corner_lons
+        )
+        np.testing.assert_array_almost_equal(
+            plane0.corner_lats, eq1.rupture.surface.corner_lats
+        )
+        np.testing.assert_array_almost_equal(
+            plane0.corner_depths, eq1.rupture.surface.corner_depths
+        )
diff --git a/tests/magnitude_scaling_relations_test.py b/tests/magnitude_scaling_relations_test.py
index f05c9a0a434f2e879d54ca369211139815ddf5af..2486cbd52dd16ca6b604a3554ff06c3979e9cb9d 100644
--- a/tests/magnitude_scaling_relations_test.py
+++ b/tests/magnitude_scaling_relations_test.py
@@ -1,9 +1,11 @@
 import unittest
 from math import radians, sin
-from shakyground2.magnitude_scaling_relations import (PEERScalingRelation,
-                                                      StrasserEtAl2010Interface,
-                                                      StrasserEtAl2010Inslab,
-                                                      Stafford2014)
+from shakyground2.magnitude_scaling_relations import (
+    PEERScalingRelation,
+    StrasserEtAl2010Interface,
+    StrasserEtAl2010Inslab,
+    Stafford2014,
+)
 
 
 class PEERScalingRelationTestCase(unittest.TestCase):
@@ -11,6 +13,7 @@ class PEERScalingRelationTestCase(unittest.TestCase):
     Simple tests for the PEER scaling relation, including rescaling according to the
     seismogenic thickness
     """
+
     def setUp(self):
         self.msr = PEERScalingRelation()
 
@@ -64,6 +67,7 @@ class StrasserEtAl2010InslabTestCase(unittest.TestCase):
 
 # Expected values of the Stafford (2014) scaling relation generated by
 # code supplied by the author
+# fmt: off
 STAFFORD_2014_TEST_TABLE = [
     [4.0, -90.0, 60.0, 10.0,    0.9627518,   0.9226897,  1.0434189],
     [4.0, -90.0, 60.0, 20.0,    0.9627518,   0.9226917,  1.0434166],
diff --git a/tests/rupture_mechanism_test.py b/tests/rupture_mechanism_test.py
index 0118b6ede06383cbfddd412c80bc65233ed76bc4..79a7f3df9cdbb7751c536ecd19a4a65a2f676f4a 100644
--- a/tests/rupture_mechanism_test.py
+++ b/tests/rupture_mechanism_test.py
@@ -3,7 +3,6 @@ Unit tests for rupture mechanism
 """
 import unittest
 import numpy as np
-from openquake.hazardlib.pmf import PMF
 from shakyground2.rupture_mechanism import RuptureMechanism
 
 
@@ -11,6 +10,7 @@ class RuptureMechanismTestCase(unittest.TestCase):
     """
     Test the instatiation and operation of the rupture mechanism class
     """
+
     @staticmethod
     def _mechanism_pmf_to_arrays(mech):
         """
@@ -35,16 +35,13 @@ class RuptureMechanismTestCase(unittest.TestCase):
         # is defined - i.e. global distribution of mechanisms
         mech = RuptureMechanism()
         probs, strikes, dips, rakes = self._mechanism_pmf_to_arrays(mech)
-        np.testing.assert_array_almost_equal(np.unique(strikes),
-                                             np.arange(0.0, 359.0, 15.0))
-        np.testing.assert_array_almost_equal(np.unique(dips),
-                                             np.arange(20.0, 95.0, 5.0))
-        np.testing.assert_array_almost_equal(np.unique(rakes),
-                                             np.arange(-165.0, 181.0, 15.0))
+        np.testing.assert_array_almost_equal(np.unique(strikes), np.arange(0.0, 359.0, 15.0))
+        np.testing.assert_array_almost_equal(np.unique(dips), np.arange(20.0, 95.0, 5.0))
+        np.testing.assert_array_almost_equal(np.unique(rakes), np.arange(-165.0, 181.0, 15.0))
 
     def test_rupture_mechanism_single_plane(self):
         # Checks instantiation when a single combination of strike, dip and rake are input
-        mech = RuptureMechanism.from_strike_dip_rake(strike=20., dip=90.0, rake=10.0)
+        mech = RuptureMechanism.from_strike_dip_rake(strike=20.0, dip=90.0, rake=10.0)
         probs, strikes, dips, rakes = self._mechanism_pmf_to_arrays(mech)
         np.testing.assert_array_almost_equal(probs, np.array([1.0]))
         np.testing.assert_array_almost_equal(strikes, np.array([20.0]))
@@ -53,7 +50,7 @@ class RuptureMechanismTestCase(unittest.TestCase):
 
     def test_rupture_mechanism_from_nodal_planes(self):
         # Checks the creation of the class from a pair of nodal planes
-        npl1 = {"strike": 45.0, "dip": 80.0, "rake": 10.}
+        npl1 = {"strike": 45.0, "dip": 80.0, "rake": 10.0}
         npl2 = {"strike": 135.0, "dip": 90.0, "rake": -170.0}
         mech = RuptureMechanism.from_nodal_planes(npl1, npl2)
         probs, strikes, dips, rakes = self._mechanism_pmf_to_arrays(mech)
@@ -64,14 +61,15 @@ class RuptureMechanismTestCase(unittest.TestCase):
 
     def test_sample_mechanism(self):
         # Test the random sampling
-        npl1 = {"strike": 45.0, "dip": 80.0, "rake": 10.}
+        npl1 = {"strike": 45.0, "dip": 80.0, "rake": 10.0}
         npl2 = {"strike": 135.0, "dip": 90.0, "rake": -170.0}
         mech = RuptureMechanism.from_nodal_planes(npl1, npl2)
         np.random.seed(1000)
         strikes, dips, rakes = mech.sample(5)
-        np.testing.assert_array_almost_equal(strikes,
-                                             np.array([135.0, 45.0, 135.0, 45.0, 135.0]))
-        np.testing.assert_array_almost_equal(dips,
-                                             np.array([90.0, 80.0, 90.0, 80.0, 90.0]))
-        np.testing.assert_array_almost_equal(rakes,
-                                             np.array([-170.0, 10.0, -170.0, 10.0, -170.0]))
+        np.testing.assert_array_almost_equal(
+            strikes, np.array([135.0, 45.0, 135.0, 45.0, 135.0])
+        )
+        np.testing.assert_array_almost_equal(dips, np.array([90.0, 80.0, 90.0, 80.0, 90.0]))
+        np.testing.assert_array_almost_equal(
+            rakes, np.array([-170.0, 10.0, -170.0, 10.0, -170.0])
+        )
diff --git a/tests/site_model_test.py b/tests/site_model_test.py
index 490338481a70a09a63db00c6723b6ba1ba363c1d..9e7869d3529009489063b922846b6ee88c4f7fe5 100644
--- a/tests/site_model_test.py
+++ b/tests/site_model_test.py
@@ -177,9 +177,7 @@ class SiteModelTestCase(unittest.TestCase):
             if site_array1[name].dtype not in (np.float,):
                 self.assertTrue(np.array_equal(site_array1[name], site_array2[name]))
             else:
-                np.testing.assert_array_almost_equal(
-                    site_array1[name], site_array2[name]
-                )
+                np.testing.assert_array_almost_equal(site_array1[name], site_array2[name])
 
     def test_build_site_collection(self):
         # Checks that the returned site collection is accurately reproduced in
diff --git a/tests/valid_test.py b/tests/valid_test.py
index b311d27d656eb45e6f72b38d09c59468a89b3d99..0de176771bfcdb638a83405a3273162a70e54b8f 100644
--- a/tests/valid_test.py
+++ b/tests/valid_test.py
@@ -13,6 +13,7 @@ class AtomicValidationsTestCase(unittest.TestCase):
     Basic test set for "atomic" validation checks - i.e. validations on
     a single value
     """
+
     def test_valid_longitude(self):
         self.assertAlmostEqual(valid.longitude(10.0), 10.0)
         self.assertAlmostEqual(valid.longitude(-90.0), -90.0)
@@ -20,8 +21,7 @@ class AtomicValidationsTestCase(unittest.TestCase):
     def test_invalid_longitude(self):
         with self.assertRaises(ValueError) as ve:
             valid.longitude(200.0)
-        self.assertEqual("Longitude 200.0000 is not between -180 and 180",
-                         str(ve.exception))
+        self.assertEqual("Longitude 200.0000 is not between -180 and 180", str(ve.exception))
 
     def test_valid_latitude(self):
         self.assertAlmostEqual(valid.latitude(30.0), 30.0)
@@ -30,8 +30,7 @@ class AtomicValidationsTestCase(unittest.TestCase):
     def test_invalid_latitude(self):
         with self.assertRaises(ValueError) as ve:
             valid.latitude(100.0)
-        self.assertEqual("Latitude 100.0000 is not between -90 and 90",
-                         str(ve.exception))
+        self.assertEqual("Latitude 100.0000 is not between -90 and 90", str(ve.exception))
 
     def test_valid_positive_float(self):
         self.assertEqual(valid.positive_float(1.0, "foo"), 1.0)
@@ -39,14 +38,16 @@ class AtomicValidationsTestCase(unittest.TestCase):
     def test_invalid_positive_float(self):
         with self.assertRaises(ValueError) as ve:
             valid.positive_float(-1.0, "foo")
-        self.assertEqual("foo must be positive float or zero, -1.000000 input",
-                         str(ve.exception))
+        self.assertEqual(
+            "foo must be positive float or zero, -1.000000 input", str(ve.exception)
+        )
 
 
 class MechanismValidationTestCase(unittest.TestCase):
     """
     Tests basic checks for the focal mechanims
     """
+
     def test_valid_strike(self):
         self.assertAlmostEqual(valid.strike(90.0), 90.0)
         self.assertAlmostEqual(valid.strike(0.0), 0.0)
@@ -54,8 +55,7 @@ class MechanismValidationTestCase(unittest.TestCase):
     def test_invalid_strike(self):
         with self.assertRaises(ValueError) as ve:
             valid.strike(400.0)
-        self.assertEqual("Strike 400.00 not in the range 0 to 360 degrees",
-                         str(ve.exception))
+        self.assertEqual("Strike 400.00 not in the range 0 to 360 degrees", str(ve.exception))
 
     def test_valid_dip(self):
         self.assertAlmostEqual(valid.dip(90.0), 90.0)
@@ -64,8 +64,7 @@ class MechanismValidationTestCase(unittest.TestCase):
     def test_invalid_dip(self):
         with self.assertRaises(ValueError) as ve:
             valid.dip(-20.0)
-        self.assertEqual("Dip -20.00 not in the range 0 to 90 degrees",
-                         str(ve.exception))
+        self.assertEqual("Dip -20.00 not in the range 0 to 90 degrees", str(ve.exception))
 
     def test_valid_rake(self):
         self.assertAlmostEqual(valid.rake(90.0), 90.0)
@@ -74,13 +73,14 @@ class MechanismValidationTestCase(unittest.TestCase):
     def test_invalid_rake(self):
         with self.assertRaises(ValueError) as ve:
             valid.rake(-200.0)
-        self.assertEqual("Rake -200.00 not in the range -180 to 180 degrees",
-                         str(ve.exception))
+        self.assertEqual("Rake -200.00 not in the range -180 to 180 degrees", str(ve.exception))
 
     def test_full_mechanism(self):
         strike = 45.0
         dip = 60.0
         rake = -90.0
+
+        # pylint: disable=E1101
         mechanism = valid.mechanism(strike, dip, rake)
         self.assertIsInstance(mechanism, NP)
         self.assertAlmostEqual(mechanism.strike, strike)
@@ -102,6 +102,7 @@ class EarthquakeDomainTestCase(unittest.TestCase):
     """
     Tests validations applied to the earthquake rupture or environment
     """
+
     def test_valid_seismogenic_thickness(self):
         usd, lsd = valid.seismogenic_thickness(0.0, 20.0)
         self.assertAlmostEqual(usd, 0.0)
@@ -110,8 +111,11 @@ class EarthquakeDomainTestCase(unittest.TestCase):
     def test_invalid_seismogenic_thickness(self):
         with self.assertRaises(ValueError) as ve:
             usd, lsd = valid.seismogenic_thickness(20.0, 0.0)
-        self.assertEqual("Lower seismogenic depth 0.00 km shallower than "
-                         "upper seismogenic depth 20.00 km", str(ve.exception))
+        self.assertEqual(
+            "Lower seismogenic depth 0.00 km shallower than "
+            "upper seismogenic depth 20.00 km",
+            str(ve.exception),
+        )
 
     def test_valid_hypocenter_position(self):
         hypo_pos = valid.hypocenter_position([0.3, 0.8])
@@ -122,25 +126,28 @@ class EarthquakeDomainTestCase(unittest.TestCase):
         # Invalid along-strike position
         with self.assertRaises(ValueError) as ve_as:
             valid.hypocenter_position([-0.1, 0.5])
-        self.assertEqual("Along strike position -0.100 should be in the range 0 to 1",
-                         str(ve_as.exception))
+        self.assertEqual(
+            "Along strike position -0.100 should be in the range 0 to 1", str(ve_as.exception)
+        )
         # Invalid down-dip position
         with self.assertRaises(ValueError) as ve_dd:
             valid.hypocenter_position([0.5, 1.2])
-        self.assertEqual("Down dip position 1.200 should be in the range 0 to 1",
-                         str(ve_dd.exception))
+        self.assertEqual(
+            "Down dip position 1.200 should be in the range 0 to 1", str(ve_dd.exception)
+        )
 
     def test_valid_scaling_relation(self):
         # Test the case that no scaling relation is input
         self.assertIsInstance(valid.scaling_relation(None), PEERScalingRelation)
         # Test the case that a valid MSR is passed
-        self.assertIsInstance(valid.scaling_relation(Stafford2014()),
-                              Stafford2014)
+        self.assertIsInstance(valid.scaling_relation(Stafford2014()), Stafford2014)
 
     def test_invalid_scaling_relation(self):
         # Verify that a correct TypeError is raised when an invalid scaling
         # relation type is passed
         with self.assertRaises(TypeError) as te:
             valid.scaling_relation("XYZ")
-        self.assertEqual("Magnitude Scaling Relation XYZ not instance of BaseScalingRelation",
-                         str(te.exception))
+        self.assertEqual(
+            "Magnitude Scaling Relation XYZ not instance of BaseScalingRelation",
+            str(te.exception),
+        )