diff --git a/examples/advection_correction.py b/examples/advection_correction.py
index 5bab1d9bd..567fbc5f7 100644
--- a/examples/advection_correction.py
+++ b/examples/advection_correction.py
@@ -111,7 +111,7 @@ def advection_correction(R, T=5, t=1):
 
         Rd += (T - i) * R1 + i * R2
 
-    return t / T ** 2 * Rd
+    return t / T**2 * Rd
 
 
 ###############################################################################
diff --git a/examples/plot_cascade_decomposition.py b/examples/plot_cascade_decomposition.py
index 2149716d0..a9f027993 100644
--- a/examples/plot_cascade_decomposition.py
+++ b/examples/plot_cascade_decomposition.py
@@ -61,7 +61,7 @@
 M, N = F.shape
 fig, ax = plt.subplots()
 im = ax.imshow(
-    np.log(F ** 2), vmin=4, vmax=24, cmap=cm.jet, extent=(-N / 2, N / 2, -M / 2, M / 2)
+    np.log(F**2), vmin=4, vmax=24, cmap=cm.jet, extent=(-N / 2, N / 2, -M / 2, M / 2)
 )
 cb = fig.colorbar(im)
 ax.set_xlabel("Wavenumber $k_x$")
diff --git a/pysteps/cascade/bandpass_filters.py b/pysteps/cascade/bandpass_filters.py
index 8698f20f1..74d2591c7 100644
--- a/pysteps/cascade/bandpass_filters.py
+++ b/pysteps/cascade/bandpass_filters.py
@@ -208,7 +208,7 @@ def __init__(self, c, s):
 
         def __call__(self, x):
             x = log_e(x) - self.c
-            return np.exp(-(x ** 2.0) / (2.0 * self.s ** 2.0))
+            return np.exp(-(x**2.0) / (2.0 * self.s**2.0))
 
     weight_funcs = []
     central_wavenumbers = [0.0]
diff --git a/pysteps/datasets.py b/pysteps/datasets.py
index 26b70d6db..c75bdf02c 100644
--- a/pysteps/datasets.py
+++ b/pysteps/datasets.py
@@ -117,10 +117,10 @@ def __call__(self, count, block_size, total_size, exact=True):
 
         self._clear_line()
 
-        downloaded_size = count * block_size / (1024 ** 2)
+        downloaded_size = count * block_size / (1024**2)
 
         if self.total_size is None and total_size > 0:
-            self.total_size = total_size / (1024 ** 2)
+            self.total_size = total_size / (1024**2)
 
         if count == 0:
             self.init_time = time.time()
diff --git a/pysteps/downscaling/rainfarm.py b/pysteps/downscaling/rainfarm.py
index f33b489ad..c173748a8 100644
--- a/pysteps/downscaling/rainfarm.py
+++ b/pysteps/downscaling/rainfarm.py
@@ -98,7 +98,7 @@ def downscale(precip, alpha=None, ds_factor=16, threshold=None, return_alpha=Fal
     if alpha is None:
         fp = np.fft.fft2(precip)
         fp_abs = abs(fp)
-        log_power_spectrum = np.log(fp_abs ** 2)
+        log_power_spectrum = np.log(fp_abs**2)
         valid = (k != 0) & np.isfinite(log_power_spectrum)
         alpha = _log_slope(np.log(k[valid]), log_power_spectrum[valid])
 
@@ -114,7 +114,7 @@ def downscale(precip, alpha=None, ds_factor=16, threshold=None, return_alpha=Fal
     P_u = np.repeat(np.repeat(precip, ds_factor, axis=0), ds_factor, axis=1)
     rad = int(round(ds_factor / np.sqrt(np.pi)))
     (mx, my) = np.mgrid[-rad : rad + 0.01, -rad : rad + 0.01]
-    tophat = ((mx ** 2 + my ** 2) <= rad ** 2).astype(float)
+    tophat = ((mx**2 + my**2) <= rad**2).astype(float)
     tophat /= tophat.sum()
 
     P_agg = _balanced_spatial_average(P_u, tophat)
diff --git a/pysteps/noise/fftgenerators.py b/pysteps/noise/fftgenerators.py
index 370370324..aa423b45e 100644
--- a/pysteps/noise/fftgenerators.py
+++ b/pysteps/noise/fftgenerators.py
@@ -637,21 +637,21 @@ def initialize_nonparam_2d_nested_filter(field, gridres=1.0, **kwargs):
     # prepare indices
     Idxi = np.array([[0, dim_y]])
     Idxj = np.array([[0, dim_x]])
-    Idxipsd = np.array([[0, 2 ** max_level]])
-    Idxjpsd = np.array([[0, 2 ** max_level]])
+    Idxipsd = np.array([[0, 2**max_level]])
+    Idxjpsd = np.array([[0, 2**max_level]])
 
     # generate the FFT sample frequencies
     freqx = fft.fftfreq(dim_x, gridres)
     freqy = fft.fftfreq(dim_y, gridres)
     fx, fy = np.meshgrid(freqx, freqy)
-    freq_grid = np.sqrt(fx ** 2 + fy ** 2)
+    freq_grid = np.sqrt(fx**2 + fy**2)
 
     # domain fourier filter
     F0 = initialize_nonparam_2d_fft_filter(
         field, win_fun=win_fun, donorm=True, use_full_fft=True, fft_method=fft
     )["field"]
     # and allocate it to the final grid
-    F = np.zeros((2 ** max_level, 2 ** max_level, F0.shape[0], F0.shape[1]))
+    F = np.zeros((2**max_level, 2**max_level, F0.shape[0], F0.shape[1]))
     F += F0[np.newaxis, np.newaxis, :, :]
 
     # now loop levels and build composite spectra
@@ -710,9 +710,9 @@ def initialize_nonparam_2d_nested_filter(field, gridres=1.0, **kwargs):
 
         # update indices
         level += 1
-        Idxi, Idxj = _split_field((0, dim[0]), (0, dim[1]), 2 ** level)
+        Idxi, Idxj = _split_field((0, dim[0]), (0, dim[1]), 2**level)
         Idxipsd, Idxjpsd = _split_field(
-            (0, 2 ** max_level), (0, 2 ** max_level), 2 ** level
+            (0, 2**max_level), (0, 2**max_level), 2**level
         )
 
     return {"field": F, "input_shape": field.shape[1:], "use_full_fft": True}
@@ -842,8 +842,8 @@ def _split_field(idxi, idxj, Segments):
     winsizei = int(sizei / Segments)
     winsizej = int(sizej / Segments)
 
-    Idxi = np.zeros((Segments ** 2, 2))
-    Idxj = np.zeros((Segments ** 2, 2))
+    Idxi = np.zeros((Segments**2, 2))
+    Idxj = np.zeros((Segments**2, 2))
 
     count = -1
     for i in range(Segments):
diff --git a/pysteps/nowcasts/anvil.py b/pysteps/nowcasts/anvil.py
index 76ac45e3e..32a2478a5 100644
--- a/pysteps/nowcasts/anvil.py
+++ b/pysteps/nowcasts/anvil.py
@@ -475,14 +475,14 @@ def _moving_window_corrcoef(x, y, window_radius):
     if window_radius is not None:
         n = gaussian_filter(mask, window_radius, mode="constant")
 
-        ssx = gaussian_filter(x ** 2, window_radius, mode="constant")
-        ssy = gaussian_filter(y ** 2, window_radius, mode="constant")
+        ssx = gaussian_filter(x**2, window_radius, mode="constant")
+        ssy = gaussian_filter(y**2, window_radius, mode="constant")
         sxy = gaussian_filter(x * y, window_radius, mode="constant")
     else:
         n = np.mean(mask)
 
-        ssx = np.mean(x ** 2)
-        ssy = np.mean(y ** 2)
+        ssx = np.mean(x**2)
+        ssy = np.mean(y**2)
         sxy = np.mean(x * y)
 
     stdx = np.sqrt(ssx / n)
diff --git a/pysteps/nowcasts/lagrangian_probability.py b/pysteps/nowcasts/lagrangian_probability.py
index 2bc8a5a08..ed1534125 100644
--- a/pysteps/nowcasts/lagrangian_probability.py
+++ b/pysteps/nowcasts/lagrangian_probability.py
@@ -128,4 +128,4 @@ def _get_kernel(size):
     else:
         xx, yy = np.mgrid[:size, :size]
         circle = (xx - middle) ** 2 + (yy - middle) ** 2
-        return np.asarray(circle <= (middle ** 2), dtype=np.float32)
+        return np.asarray(circle <= (middle**2), dtype=np.float32)
diff --git a/pysteps/nowcasts/linda.py b/pysteps/nowcasts/linda.py
index 8f824d55a..52861116c 100644
--- a/pysteps/nowcasts/linda.py
+++ b/pysteps/nowcasts/linda.py
@@ -430,8 +430,8 @@ def _compute_inverse_acf_mapping(target_dist, target_dist_params, n_intervals=10
     """Compute the inverse ACF mapping between two distributions."""
     phi = (
         lambda x1, x2, rho: 1.0
-        / (2 * np.pi * np.sqrt(1 - rho ** 2))
-        * np.exp(-(x1 ** 2 + x2 ** 2 - 2 * rho * x1 * x2) / (2 * (1 - rho ** 2)))
+        / (2 * np.pi * np.sqrt(1 - rho**2))
+        * np.exp(-(x1**2 + x2**2 - 2 * rho * x1 * x2) / (2 * (1 - rho**2)))
     )
 
     rho_1 = np.linspace(-0.9, 0.9, n_intervals)
@@ -479,7 +479,7 @@ def _compute_kernel_anisotropic(params, cutoff=6.0):
 
     x2 = xy_grid[0, :] * xy_grid[0, :]
     y2 = xy_grid[1, :] * xy_grid[1, :]
-    result = np.exp(-(x2 / sigma1 ** 2 + y2 / sigma2 ** 2))
+    result = np.exp(-(x2 / sigma1**2 + y2 / sigma2**2))
 
     return np.reshape(result / np.sum(result), x_grid.shape)
 
@@ -576,8 +576,8 @@ def _compute_window_weights(coords, grid_height, grid_width, window_radius):
             dx = c[1] - grid_x
 
             w[i, :] = np.exp(
-                -dy * dy / (2 * window_radius_1 ** 2)
-                - dx * dx / (2 * window_radius_2 ** 2)
+                -dy * dy / (2 * window_radius_1**2)
+                - dx * dx / (2 * window_radius_2**2)
             )
     else:
         w[0, :] = np.ones((grid_height, grid_width))
@@ -873,10 +873,10 @@ def _fit_dist(err, dist, wf, mask):
     """
     Fit a lognormal distribution by maximizing the log-likelihood function
     with the constraint that the mean value is one."""
-    f = lambda p: -np.sum(np.log(stats.lognorm.pdf(err[mask], p, -0.5 * p ** 2)))
+    f = lambda p: -np.sum(np.log(stats.lognorm.pdf(err[mask], p, -0.5 * p**2)))
     p_opt = opt.minimize_scalar(f, bounds=(1e-3, 20.0), method="Bounded")
 
-    return (p_opt.x, -0.5 * p_opt.x ** 2)
+    return (p_opt.x, -0.5 * p_opt.x**2)
 
 
 # TODO: restrict the perturbation generation inside the radar mask
diff --git a/pysteps/scripts/fit_vel_pert_params.py b/pysteps/scripts/fit_vel_pert_params.py
index e18bd1cae..cfc685f31 100644
--- a/pysteps/scripts/fit_vel_pert_params.py
+++ b/pysteps/scripts/fit_vel_pert_params.py
@@ -42,7 +42,7 @@
     mu = dp_par_sum / dp_par_n
 
     std_par.append(
-        np.sqrt((dp_par_sq_sum - 2 * mu * dp_par_sum + dp_par_n * mu ** 2) / dp_par_n)
+        np.sqrt((dp_par_sq_sum - 2 * mu * dp_par_sum + dp_par_n * mu**2) / dp_par_n)
     )
 
     dp_perp_sum = results[lt]["dp_perp_sum"]
@@ -52,7 +52,7 @@
 
     std_perp.append(
         np.sqrt(
-            (dp_perp_sq_sum - 2 * mu * dp_perp_sum + dp_perp_n * mu ** 2) / dp_perp_n
+            (dp_perp_sq_sum - 2 * mu * dp_perp_sum + dp_perp_n * mu**2) / dp_perp_n
         )
     )
 
diff --git a/pysteps/scripts/run_vel_pert_analysis.py b/pysteps/scripts/run_vel_pert_analysis.py
index 8ae241553..7990b6298 100644
--- a/pysteps/scripts/run_vel_pert_analysis.py
+++ b/pysteps/scripts/run_vel_pert_analysis.py
@@ -170,9 +170,9 @@
             n_samples = DP_par.size
 
         results[lt]["dp_par_sum"] += np.sum(DP_par)
-        results[lt]["dp_par_sq_sum"] += np.sum(DP_par ** 2)
+        results[lt]["dp_par_sq_sum"] += np.sum(DP_par**2)
         results[lt]["dp_perp_sum"] += np.sum(DP_perp)
-        results[lt]["dp_perp_sq_sum"] += np.sum(DP_perp ** 2)
+        results[lt]["dp_perp_sq_sum"] += np.sum(DP_perp**2)
         results[lt]["n_samples"] += n_samples
 
 with open("%s" % args.outfile, "wb") as f:
diff --git a/pysteps/timeseries/autoregression.py b/pysteps/timeseries/autoregression.py
index 9ad00c2b8..729525480 100644
--- a/pysteps/timeseries/autoregression.py
+++ b/pysteps/timeseries/autoregression.py
@@ -72,7 +72,7 @@ def adjust_lag2_corrcoef2(gamma_1, gamma_2):
     """
     gamma_2 = np.maximum(gamma_2, 2 * gamma_1 * gamma_2 - 1)
     gamma_2 = np.maximum(
-        gamma_2, (3 * gamma_1 ** 2 - 2 + 2 * (1 - gamma_1 ** 2) ** 1.5) / gamma_1 ** 2
+        gamma_2, (3 * gamma_1**2 - 2 + 2 * (1 - gamma_1**2) ** 1.5) / gamma_1**2
     )
 
     return gamma_2
diff --git a/pysteps/timeseries/correlation.py b/pysteps/timeseries/correlation.py
index b7f0d31b9..a43b41a7e 100644
--- a/pysteps/timeseries/correlation.py
+++ b/pysteps/timeseries/correlation.py
@@ -245,20 +245,20 @@ def _moving_window_corrcoef(x, y, window_radius, window="gaussian", mask=None):
     else:
         window_size = window_radius
 
-    n = convol_filter(mask, window_size, mode="constant") * window_size ** 2
+    n = convol_filter(mask, window_size, mode="constant") * window_size**2
 
-    sx = convol_filter(x, window_size, mode="constant") * window_size ** 2
-    sy = convol_filter(y, window_size, mode="constant") * window_size ** 2
+    sx = convol_filter(x, window_size, mode="constant") * window_size**2
+    sy = convol_filter(y, window_size, mode="constant") * window_size**2
 
-    ssx = convol_filter(x ** 2, window_size, mode="constant") * window_size ** 2
-    ssy = convol_filter(y ** 2, window_size, mode="constant") * window_size ** 2
-    sxy = convol_filter(x * y, window_size, mode="constant") * window_size ** 2
+    ssx = convol_filter(x**2, window_size, mode="constant") * window_size**2
+    ssy = convol_filter(y**2, window_size, mode="constant") * window_size**2
+    sxy = convol_filter(x * y, window_size, mode="constant") * window_size**2
 
     mux = sx / n
     muy = sy / n
 
-    stdx = np.sqrt(ssx - 2 * mux * sx + n * mux ** 2)
-    stdy = np.sqrt(ssy - 2 * muy * sy + n * muy ** 2)
+    stdx = np.sqrt(ssx - 2 * mux * sx + n * mux**2)
+    stdy = np.sqrt(ssy - 2 * muy * sy + n * muy**2)
     cov = sxy - muy * sx - mux * sy + n * mux * muy
 
     mask = np.logical_and(stdx > 1e-8, stdy > 1e-8)
diff --git a/pysteps/utils/conversion.py b/pysteps/utils/conversion.py
index f68e736e3..d5587c880 100644
--- a/pysteps/utils/conversion.py
+++ b/pysteps/utils/conversion.py
@@ -282,7 +282,7 @@ def to_reflectivity(R, metadata, zr_a=None, zr_b=None):
         if zr_b is None:
             zr_b = metadata.get("zr_b", 1.6)  # Default to Marshall–Palmer
 
-        R = zr_a * R ** zr_b
+        R = zr_a * R**zr_b
         metadata["threshold"] = zr_a * metadata["threshold"] ** zr_b
         metadata["zerovalue"] = zr_a * metadata["zerovalue"] ** zr_b
         metadata["zr_a"] = zr_a
@@ -300,7 +300,7 @@ def to_reflectivity(R, metadata, zr_a=None, zr_b=None):
             zr_a = metadata.get("zr_a", 200.0)  # Default to Marshall-Palmer
         if zr_b is None:
             zr_b = metadata.get("zr_b", 1.6)  # Default to Marshall-Palmer
-        R = zr_a * R ** zr_b
+        R = zr_a * R**zr_b
         metadata["threshold"] = zr_a * metadata["threshold"] ** zr_b
         metadata["zerovalue"] = zr_a * metadata["zerovalue"] ** zr_b
         metadata["zr_a"] = zr_a
diff --git a/pysteps/utils/transformation.py b/pysteps/utils/transformation.py
index 945151505..46be85ab3 100644
--- a/pysteps/utils/transformation.py
+++ b/pysteps/utils/transformation.py
@@ -107,7 +107,7 @@ def boxcox_transform(
 
         else:
             R[~zeros] = (R[~zeros] ** Lambda - 1) / Lambda
-            threshold = (threshold ** Lambda - 1) / Lambda
+            threshold = (threshold**Lambda - 1) / Lambda
 
         # Set value for zeros
         if zerovalue is None:
@@ -375,7 +375,7 @@ def sqrt_transform(R, metadata=None, inverse=False, **kwargs):
         metadata["threshold"] = np.sqrt(metadata["threshold"])
     else:
         # inverse sqrt transform
-        R = R ** 2
+        R = R**2
 
         metadata["transform"] = None
         metadata["zerovalue"] = metadata["zerovalue"] ** 2
diff --git a/pysteps/verification/detcontscores.py b/pysteps/verification/detcontscores.py
index 494d45a06..c5a11dcf0 100644
--- a/pysteps/verification/detcontscores.py
+++ b/pysteps/verification/detcontscores.py
@@ -356,8 +356,8 @@ def det_cont_fct_accum(err, pred, obs):
     mobs = np.nanmean(obs, axis=axis)
     mpred = np.nanmean(pred, axis=axis)
     me = np.nanmean(res, axis=axis)
-    mse = np.nanmean(res ** 2, axis=axis)
-    mss = np.nanmean(sum ** 2, axis=axis)
+    mse = np.nanmean(res**2, axis=axis)
+    mss = np.nanmean(sum**2, axis=axis)
     mae = np.nanmean(np.abs(res), axis=axis)
 
     # expand axes for broadcasting
diff --git a/pysteps/verification/probscores.py b/pysteps/verification/probscores.py
index 51247c5cf..f1eb398ba 100644
--- a/pysteps/verification/probscores.py
+++ b/pysteps/verification/probscores.py
@@ -125,7 +125,7 @@ def CRPS_accum(CRPS, X_f, X_o):
     beta[mask, -1] = 0.0
 
     p = 1.0 * np.arange(m + 1) / m
-    res = np.sum(alpha * p ** 2.0 + beta * (1.0 - p) ** 2.0, axis=1)
+    res = np.sum(alpha * p**2.0 + beta * (1.0 - p) ** 2.0, axis=1)
 
     CRPS["CRPS_sum"] += np.sum(res)
     CRPS["n"] += len(res)
diff --git a/pysteps/verification/spatialscores.py b/pysteps/verification/spatialscores.py
index 450791783..afce35ee4 100644
--- a/pysteps/verification/spatialscores.py
+++ b/pysteps/verification/spatialscores.py
@@ -609,9 +609,9 @@ def fss_accum(fss, X_f, X_o):
         S_f = I_f
         S_o = I_o
 
-    fss["sum_obs_sq"] += np.nansum(S_o ** 2)
+    fss["sum_obs_sq"] += np.nansum(S_o**2)
     fss["sum_fct_obs"] += np.nansum(S_f * S_o)
-    fss["sum_fct_sq"] += np.nansum(S_f ** 2)
+    fss["sum_fct_sq"] += np.nansum(S_f**2)
 
 
 def fss_merge(fss_1, fss_2):
diff --git a/pysteps/visualization/spectral.py b/pysteps/visualization/spectral.py
index 2582bf3ce..aac16248a 100644
--- a/pysteps/visualization/spectral.py
+++ b/pysteps/visualization/spectral.py
@@ -94,7 +94,7 @@ def plot_spectrum1d(
     # Y-axis
     if y_units is not None:
         # { -> {{ with f-strings
-        power_units = fr"$10log_{{ 10 }}(\frac{{ {y_units}^2 }}{{ {x_units} }})$"
+        power_units = rf"$10log_{{ 10 }}(\frac{{ {y_units}^2 }}{{ {x_units} }})$"
         ax.set_ylabel(f"Power {power_units}")
 
     return ax