Fixes following pandas 2 (#327)

* Remove deprecated append method, use concat instead

* Division by zero returns `-inf`, `nan`, or `inf` depending on the numerator, instead of raising

* Fix DeprecationWarnings from scipy.ndimage

* Concatenate once

* Suppress pandas' SettingWithCopyWarning

Author: dnerini

pysteps/blending/steps.py

@@ -46,7 +46,7 @@
 import time
 
 import numpy as np
-import scipy.ndimage
+from scipy.ndimage import binary_dilation, generate_binary_structure, iterate_structure
 
 from pysteps import cascade
 from pysteps import extrapolation
@@ -1599,13 +1599,13 @@ def _compute_incremental_mask(Rbin, kr, r):
 
     # buffer observation mask
     Rbin = np.ndarray.astype(Rbin.copy(), "uint8")
-    Rd = scipy.ndimage.morphology.binary_dilation(Rbin, kr)
+    Rd = binary_dilation(Rbin, kr)
 
     # add grayscale rim
-    kr1 = scipy.ndimage.generate_binary_structure(2, 1)
+    kr1 = generate_binary_structure(2, 1)
     mask = Rd.astype(float)
     for n in range(r):
-        Rd = scipy.ndimage.morphology.binary_dilation(Rd, kr1)
+        Rd = binary_dilation(Rd, kr1)
         mask += Rd
     # normalize between 0 and 1
     return mask / mask.max()
@@ -1995,10 +1995,10 @@ def _prepare_forecast_loop(
         mask_rim = mask_kwargs.get("mask_rim", 10)
         mask_f = mask_kwargs.get("mask_f", 1.0)
         # initialize the structuring element
-        struct = scipy.ndimage.generate_binary_structure(2, 1)
+        struct = generate_binary_structure(2, 1)
         # iterate it to expand it nxn
         n = mask_f * timestep / kmperpixel
-        struct = scipy.ndimage.iterate_structure(struct, int((n - 1) / 2.0))
+        struct = iterate_structure(struct, int((n - 1) / 2.0))
     else:
         mask_rim, struct = None, None
 

pysteps/extrapolation/semilagrangian.py

@@ -15,7 +15,7 @@
 import warnings
 
 import numpy as np
-import scipy.ndimage.interpolation as ip
+from scipy.ndimage import map_coordinates
 
 
 def extrapolate(
@@ -182,10 +182,10 @@ def interpolate_motion(displacement, velocity_inc, td):
         coords_warped = xy_coords + displacement
         coords_warped = [coords_warped[1, :, :], coords_warped[0, :, :]]
 
-        velocity_inc_x = ip.map_coordinates(
+        velocity_inc_x = map_coordinates(
             velocity[0, :, :], coords_warped, mode="nearest", order=1, prefilter=False
         )
-        velocity_inc_y = ip.map_coordinates(
+        velocity_inc_y = map_coordinates(
             velocity[1, :, :], coords_warped, mode="nearest", order=1, prefilter=False
         )
 
@@ -222,7 +222,7 @@ def interpolate_motion(displacement, velocity_inc, td):
         coords_warped = [coords_warped[1, :, :], coords_warped[0, :, :]]
 
         if precip is not None:
-            precip_warped = ip.map_coordinates(
+            precip_warped = map_coordinates(
                 precip,
                 coords_warped,
                 mode=map_coordinates_mode,
@@ -232,7 +232,7 @@ def interpolate_motion(displacement, velocity_inc, td):
             )
 
             if interp_order > 1:
-                mask_warped = ip.map_coordinates(
+                mask_warped = map_coordinates(
                     mask_min,
                     coords_warped,
                     mode=map_coordinates_mode,
@@ -242,7 +242,7 @@ def interpolate_motion(displacement, velocity_inc, td):
                 )
                 precip_warped[mask_warped < 0.5] = minval
 
-                mask_warped = ip.map_coordinates(
+                mask_warped = map_coordinates(
                     mask_finite,
                     coords_warped,
                     mode=map_coordinates_mode,

pysteps/feature/tstorm.py

@@ -233,26 +233,33 @@ def get_profile(areas, binary, ref, loc_max, time, minref):
     cells = areas * binary
     cell_labels = cells[loc_max]
     labels = np.zeros(cells.shape)
+    cells_id = []
+    for n, cell_label in enumerate(cell_labels):
+        this_id = n + 1
+        x = np.where(cells == cell_label)[1]
+        y = np.where(cells == cell_label)[0]
+        cell_unique = np.zeros(cells.shape)
+        cell_unique[cells == cell_label] = 1
+        maxref = np.nanmax(ref[y, x])
+        contours = skime.find_contours(cell_unique, 0.8)
+        cells_id.append(
+            {
+                "ID": this_id,
+                "time": time,
+                "x": x,
+                "y": y,
+                "cen_x": np.round(np.nanmean(x)).astype(int),
+                "cen_y": np.round(np.nanmean(y)).astype(int),
+                "max_ref": maxref,
+                "cont": contours,
+                "area": len(x),
+            }
+        )
+        labels[cells == cell_labels[n]] = this_id
     cells_id = pd.DataFrame(
-        data=None,
+        data=cells_id,
         index=range(len(cell_labels)),
         columns=["ID", "time", "x", "y", "cen_x", "cen_y", "max_ref", "cont", "area"],
     )
-    cells_id.time = time
-    for n in range(len(cell_labels)):
-        ID = n + 1
-        cells_id.ID.iloc[n] = ID
-        cells_id.x.iloc[n] = np.where(cells == cell_labels[n])[1]
-        cells_id.y.iloc[n] = np.where(cells == cell_labels[n])[0]
-        cell_unique = np.zeros(cells.shape)
-        cell_unique[cells == cell_labels[n]] = 1
-        maxref = np.nanmax(ref[cells_id.y[n], cells_id.x[n]])
-        contours = skime.find_contours(cell_unique, 0.8)
-        cells_id.cont.iloc[n] = contours
-        cells_id.cen_x.iloc[n] = np.round(np.nanmean(cells_id.x[n])).astype(int)
-        cells_id.cen_y.iloc[n] = np.round(np.nanmean(cells_id.y[n])).astype(int)
-        cells_id.max_ref.iloc[n] = maxref
-        cells_id.area.iloc[n] = len(cells_id.x.iloc[n])
-        labels[cells == cell_labels[n]] = ID
 
     return cells_id, labels

pysteps/motion/constant.py

@@ -13,8 +13,8 @@
 """
 
 import numpy as np
-import scipy.ndimage.interpolation as ip
 import scipy.optimize as op
+from scipy.ndimage import map_coordinates
 
 
 def constant(R, **kwargs):
@@ -40,7 +40,7 @@ def constant(R, **kwargs):
 
     def f(v):
         XYW = [Y + v[1], X + v[0]]
-        R_w = ip.map_coordinates(
+        R_w = map_coordinates(
             R[-2, :, :], XYW, mode="constant", cval=np.nan, order=0, prefilter=False
         )
 

pysteps/motion/vet.py

@@ -36,7 +36,7 @@
 
 import numpy
 from numpy.ma.core import MaskedArray
-from scipy.ndimage.interpolation import zoom
+from scipy.ndimage import zoom
 from scipy.optimize import minimize
 
 from pysteps.decorators import check_input_frames

pysteps/nowcasts/sseps.py

@@ -19,8 +19,9 @@
 """
 
 import numpy as np
-import scipy.ndimage
 import time
+from scipy.ndimage import generate_binary_structure, iterate_structure
+
 
 from pysteps import cascade
 from pysteps import extrapolation
@@ -342,10 +343,10 @@ def forecast(
         mask_rim = mask_kwargs.get("mask_rim", 10)
         mask_f = mask_kwargs.get("mask_f", 1.0)
         # initialize the structuring element
-        struct = scipy.ndimage.generate_binary_structure(2, 1)
+        struct = generate_binary_structure(2, 1)
         # iterate it to expand it nxn
         n = mask_f * timestep / kmperpixel
-        struct = scipy.ndimage.iterate_structure(struct, int((n - 1) / 2.0))
+        struct = iterate_structure(struct, int((n - 1) / 2.0))
 
     noise_kwargs.update(
         {

pysteps/nowcasts/steps.py

@@ -12,7 +12,7 @@
 """
 
 import numpy as np
-import scipy.ndimage
+from scipy.ndimage import generate_binary_structure, iterate_structure
 import time
 
 from pysteps import cascade
@@ -598,10 +598,10 @@ def f(precip, i):
             mask_rim = mask_kwargs.get("mask_rim", 10)
             mask_f = mask_kwargs.get("mask_f", 1.0)
             # initialize the structuring element
-            struct = scipy.ndimage.generate_binary_structure(2, 1)
+            struct = generate_binary_structure(2, 1)
             # iterate it to expand it nxn
             n = mask_f * timestep / kmperpixel
-            struct = scipy.ndimage.iterate_structure(struct, int((n - 1) / 2.0))
+            struct = iterate_structure(struct, int((n - 1) / 2.0))
             # initialize precip mask for each member
             mask_prec = nowcast_utils.compute_dilated_mask(mask_prec, struct, mask_rim)
             mask_prec = [mask_prec.copy() for _ in range(n_ens_members)]

pysteps/nowcasts/utils.py

@@ -18,7 +18,8 @@
 
 import time
 import numpy as np
-import scipy.ndimage
+from scipy.ndimage import binary_dilation, generate_binary_structure
+
 from pysteps import extrapolation
 
 
@@ -85,13 +86,13 @@ def compute_dilated_mask(input_mask, kr, r):
     """
     # buffer the input mask
     input_mask = np.ndarray.astype(input_mask.copy(), "uint8")
-    mask_dilated = scipy.ndimage.morphology.binary_dilation(input_mask, kr)
+    mask_dilated = binary_dilation(input_mask, kr)
 
     # add grayscale rim
-    kr1 = scipy.ndimage.generate_binary_structure(2, 1)
+    kr1 = generate_binary_structure(2, 1)
     mask = mask_dilated.astype(float)
     for _ in range(r):
-        mask_dilated = scipy.ndimage.morphology.binary_dilation(mask_dilated, kr1)
+        mask_dilated = binary_dilation(mask_dilated, kr1)
         mask += mask_dilated
 
     # normalize between 0 and 1

pysteps/tracking/tdating.py

@@ -200,7 +200,7 @@ def tracking(cells_id, cells_id_prev, labels, V1, max_ID):
     cells_ad = advect(cells_id_prev, labels, V1)
     cells_ov, labels = match(cells_ad, labels)
     newlabels = np.zeros(labels.shape)
-    for ID, cell in cells_id_new.iterrows():
+    for index, cell in cells_id_new.iterrows():
         if cell.ID == 0 or np.isnan(cell.ID):
             continue
         new_ID = cells_ov[cells_ov.t_ID == cell.ID].ID.values
@@ -211,12 +211,12 @@ def tracking(cells_id, cells_id_prev, labels, V1, max_ID):
                 size.append(len(x))
             biggest = np.argmax(size)
             new_ID = new_ID[biggest]
-            cells_id_new.ID.iloc[ID] = new_ID
+            cells_id_new.loc[index, "ID"] = new_ID
         else:
             max_ID += 1
             new_ID = max_ID
-            cells_id_new.ID.iloc[ID] = new_ID
-        newlabels[labels == ID + 1] = new_ID
+            cells_id_new.loc[index, "ID"] = new_ID
+        newlabels[labels == index + 1] = new_ID
         del new_ID
     return cells_id_new, max_ID, newlabels
 
@@ -308,10 +308,11 @@ def couple_track(cell_list, max_ID, mintrack):
             index=None,
             columns=["ID", "time", "x", "y", "cen_x", "cen_y", "max_ref", "cont"],
         )
+        cell_track = []
         for t in range(len(cell_list)):
             mytime = cell_list[t]
-            mycell = mytime[mytime.ID == n]
-            cell_track = cell_track.append(mycell)
+            cell_track.append(mytime[mytime.ID == n])
+        cell_track = pd.concat(cell_track, axis=0)
 
         if len(cell_track) < mintrack:
             continue

pysteps/verification/salscores.py

@@ -16,7 +16,7 @@
 from math import sqrt, hypot
 
 import numpy as np
-from scipy.ndimage.measurements import center_of_mass
+from scipy.ndimage import center_of_mass
 
 from pysteps.exceptions import MissingOptionalDependency
 from pysteps.feature import tstorm as tstorm_detect
@@ -159,7 +159,7 @@ def sal_structure(
     observation_volume = _sal_scaled_volume(observation_objects).sum()
     nom = prediction_volume - observation_volume
     denom = prediction_volume + observation_volume
-    return nom / (0.5 * denom)
+    return np.divide(nom, (0.5 * denom))
 
 
 def sal_amplitude(prediction, observation):