diff --git a/docs/whats-new.rst b/docs/whats-new.rst
index 68114c6..024a22d 100644
--- a/docs/whats-new.rst
+++ b/docs/whats-new.rst
@@ -133,7 +133,7 @@ Enhancements
 
 - :py:func:`~transform_geopandas` can now handle grid to proj transformations.
 - Thw WRF vertical interpolation methods now accept a `fill_value` kwarg
-  which is `np.NaN` per default but can be set to `'extrapolate'` if
+  which is `np.nan` per default but can be set to `'extrapolate'` if
   needed.
 - New :py:func:`~reduce` function, useful to aggregate structured high-res
   grids to lower-res grids (:ref:`sphx_glr_auto_examples_plot_subgrid_mask.py`)
diff --git a/salem/gis.py b/salem/gis.py
index 16b32c3..27ff56d 100644
--- a/salem/gis.py
+++ b/salem/gis.py
@@ -835,10 +835,11 @@ def lookup_transform(self, data, grid=None, method=np.mean, lut=None,
         out_shape = list(in_shape)
         out_shape[-2:] = [self.ny, self.nx]
 
-        if data.dtype.kind == 'i':
-            out_data = np.zeros(out_shape, dtype=float) * np.NaN
-        else:
-            out_data = np.zeros(out_shape, dtype=data.dtype) * np.NaN
+        out_data = np.full(
+            out_shape,
+            np.nan,
+            dtype=float if data.dtype.kind == 'i' else data.dtype,
+        )
 
         def _2d_trafo(ind, outd):
             for ji, l in lut.items():
diff --git a/salem/graphics.py b/salem/graphics.py
index fa512ad..9f61582 100644
--- a/salem/graphics.py
+++ b/salem/graphics.py
@@ -515,14 +515,15 @@ def _check_data(self, data=None, crs=None, interp='nearest',
                     warnings.filterwarnings("ignore", message=mess)
                     mess = "Passing `np.nan` to mean no clipping in np.clip"
                     warnings.filterwarnings("ignore", message=mess)
+                    nans = data.filled(np.nan)
                     try:
-                        data = imresize(data.filled(np.NaN),
+                        data = imresize(nans,
                                         (self.grid.ny, self.grid.nx),
                                         order=interp, mode='edge',
                                         anti_aliasing=True)
                     except RuntimeError:
                         # For some order anti_aliasing doesnt work with 'edge'
-                        data = imresize(data.filled(np.NaN),
+                        data = imresize(nans,
                                         (self.grid.ny, self.grid.nx),
                                         order=interp, mode='edge',
                                         anti_aliasing=False)
diff --git a/salem/sio.py b/salem/sio.py
index a6ca87a..a8fda32 100644
--- a/salem/sio.py
+++ b/salem/sio.py
@@ -747,7 +747,7 @@ def deacc(self, as_rate=True):
 
         return out
 
-    def interpz(self, zcoord, levels, dim_name='', fill_value=np.NaN,
+    def interpz(self, zcoord, levels, dim_name='', fill_value=np.nan,
                 use_multiprocessing=True):
         """Interpolates the array along the vertical dimension
 
@@ -759,7 +759,7 @@ def interpz(self, zcoord, levels, dim_name='', fill_value=np.NaN,
           the levels at which to interpolate
         dim_name: str
           the name of the new dimension
-        fill_value : np.NaN or 'extrapolate', optional
+        fill_value : np.nan or 'extrapolate', optional
           how to handle levels below the topography. Default is to mark them
           as invalid, but you might want the have them extrapolated.
         use_multiprocessing: bool
@@ -835,7 +835,7 @@ def transform_and_add(self, other, grid=None, interp='nearest', ks=3,
                     new_name = v
                 self._obj[new_name] = out[v]
 
-    def wrf_zlevel(self, varname, levels=None, fill_value=np.NaN,
+    def wrf_zlevel(self, varname, levels=None, fill_value=np.nan,
                    use_multiprocessing=True):
         """Interpolates to a specified height above sea level.
 
@@ -845,7 +845,7 @@ def wrf_zlevel(self, varname, levels=None, fill_value=np.NaN,
           the name of the variable to interpolate
         levels: 1d array
           levels at which to interpolate (default: some levels I thought of)
-        fill_value : np.NaN or 'extrapolate', optional
+        fill_value : np.nan or 'extrapolate', optional
           how to handle levels below the topography. Default is to mark them
           as invalid, but you might want the have them extrapolated.
         use_multiprocessing: bool
@@ -867,7 +867,7 @@ def wrf_zlevel(self, varname, levels=None, fill_value=np.NaN,
         out['z'].attrs['units'] = 'm'
         return out
 
-    def wrf_plevel(self, varname, levels=None, fill_value=np.NaN,
+    def wrf_plevel(self, varname, levels=None, fill_value=np.nan,
                    use_multiprocessing=True):
         """Interpolates to a specified pressure level (hPa).
 
@@ -877,7 +877,7 @@ def wrf_plevel(self, varname, levels=None, fill_value=np.NaN,
           the name of the variable to interpolate
         levels: 1d array
           levels at which to interpolate (default: some levels I thought of)
-        fill_value : np.NaN or 'extrapolate', optional
+        fill_value : np.nan or 'extrapolate', optional
           how to handle levels below the topography. Default is to mark them
           as invalid, but you might want the have them extrapolated.
         use_multiprocessing: bool
diff --git a/salem/tests/test_gis.py b/salem/tests/test_gis.py
index 5f7eb25..f73c956 100644
--- a/salem/tests/test_gis.py
+++ b/salem/tests/test_gis.py
@@ -706,7 +706,7 @@ def test_map_gridded_data_over(self):
             nx, ny = (4, 5)
             data = np.arange(nx * ny).reshape((ny, nx)).astype(float)
 
-            in_data = data * np.NaN
+            in_data = data * np.nan
             in_data[0, :] = 78
 
             # Nearest Neighbor
@@ -823,13 +823,13 @@ def test_map_real_data(self):
         odata = grid_to.map_gridded_data(data, grid_from, interp='linear')
         # At the borders IDL and Python take other decision on wether it
         # should be a NaN or not (Python seems to be more conservative)
-        ref_data[np.where(odata.mask)] = np.NaN
+        ref_data[np.where(odata.mask)] = np.nan
         assert np.sum(np.isfinite(ref_data)) != 0
-        assert_allclose(ref_data, odata.filled(np.NaN), atol=1e-3)
+        assert_allclose(ref_data, odata.filled(np.nan), atol=1e-3)
 
         odata = grid_to.map_gridded_data(data, grid_from, interp='spline')
-        odata[np.where(~ np.isfinite(ref_data))] = np.NaN
-        ref_data[np.where(~ np.isfinite(odata))] = np.NaN
+        odata[np.where(~ np.isfinite(ref_data))] = np.nan
+        ref_data[np.where(~ np.isfinite(odata))] = np.nan
         assert np.sum(np.isfinite(ref_data)) != 0
         assert_allclose(ref_data, odata, rtol=0.2, atol=3)
 
@@ -837,14 +837,14 @@ def test_map_real_data(self):
         data = np.array([data, data])
         ref_data = np.array([ref_data, ref_data])
         odata = grid_to.map_gridded_data(data, grid_from, interp='linear')
-        odata = odata.filled(np.NaN)
-        ref_data[np.where(~ np.isfinite(odata))] = np.NaN
+        odata = odata.filled(np.nan)
+        ref_data[np.where(~ np.isfinite(odata))] = np.nan
         assert np.sum(np.isfinite(ref_data)) != 0
         assert_allclose(ref_data, odata, atol=1e-3)
 
         odata = grid_to.map_gridded_data(data, grid_from, interp='spline')
-        odata[np.where(~ np.isfinite(ref_data))] = np.NaN
-        ref_data[np.where(~ np.isfinite(odata))] = np.NaN
+        odata[np.where(~ np.isfinite(ref_data))] = np.nan
+        ref_data[np.where(~ np.isfinite(odata))] = np.nan
         assert np.sum(np.isfinite(ref_data)) != 0
         assert_allclose(ref_data, odata, rtol=0.2, atol=3)
 
diff --git a/salem/tests/test_graphics.py b/salem/tests/test_graphics.py
index 03cb193..68c186c 100644
--- a/salem/tests/test_graphics.py
+++ b/salem/tests/test_graphics.py
@@ -327,7 +327,7 @@ def test_map(self):
         cmap.set_bad('pink')
 
         # Add masked arrays
-        a[1, 1] = np.NaN
+        a[1, 1] = np.nan
         c.set_data(a)
         rgb1 = c.to_rgb()
         c.set_data(a, crs=g)
@@ -436,7 +436,7 @@ def test_datalevels():
     c.visualize(next(ax), title="Auto levels with oob data")
 
     # Missing data
-    a[3, 0] = np.NaN
+    a[3, 0] = np.nan
     c = DataLevels(nlevels=127, vmin=0, vmax=3, data=a, cmap=cm)
     c.visualize(next(ax), title="missing data")
 
@@ -801,7 +801,7 @@ def test_hef_topo_withnan():
     c.set_cmap(get_cmap('topo'))
     c.set_plot_params(nlevels=256)
     # Try with nan data
-    h[-100:, -100:] = np.NaN
+    h[-100:, -100:] = np.nan
     c.set_data(h)
     fig, ax = plt.subplots(1, 1)
     c.visualize(ax=ax, title='color with NaN')
diff --git a/salem/tests/test_misc.py b/salem/tests/test_misc.py
index 3ae5721..1cc437c 100644
--- a/salem/tests/test_misc.py
+++ b/salem/tests/test_misc.py
@@ -274,8 +274,8 @@ def test_roi(self):
         import xarray as xr
         # Check that all attrs are preserved
         with sio.open_xr_dataset(get_demo_file('era_interim_tibet.nc')) as ds:
-            ds.encoding = {'_FillValue': np.NaN}
-            ds['t2m'].encoding = {'_FillValue': np.NaN}
+            ds.encoding = {'_FillValue': np.nan}
+            ds['t2m'].encoding = {'_FillValue': np.nan}
             ds_ = ds.salem.roi(roi=np.ones_like(ds.t2m.values[0, ...]))
             xr.testing.assert_identical(ds, ds_)
             assert ds.encoding == ds_.encoding
@@ -524,13 +524,13 @@ def test_prcp(self):
         uns = nc['RAINNC'].isel(Time=slice(1, len(nc.bottom_top_stag))).values - \
               nc['RAINNC'].isel(Time=slice(0, -1)).values
         nc['REF_PRCP_NC'].values[1:, ...] = uns * 60 / 180.  # for three hours
-        nc['REF_PRCP_NC'].values[0, ...] = np.NaN
+        nc['REF_PRCP_NC'].values[0, ...] = np.nan
 
         nc['REF_PRCP_C'] = nc['RAINC']*0.
         uns = nc['RAINC'].isel(Time=slice(1, len(nc.bottom_top_stag))).values - \
               nc['RAINC'].isel(Time=slice(0, -1)).values
         nc['REF_PRCP_C'].values[1:, ...] = uns * 60 / 180.  # for three hours
-        nc['REF_PRCP_C'].values[0, ...] = np.NaN
+        nc['REF_PRCP_C'].values[0, ...] = np.nan
 
         nc['REF_PRCP'] = nc['REF_PRCP_C'] + nc['REF_PRCP_NC']
 
diff --git a/salem/wrftools.py b/salem/wrftools.py
index f2749cc..8b95d32 100644
--- a/salem/wrftools.py
+++ b/salem/wrftools.py
@@ -286,7 +286,7 @@ def __getitem__(self, item):
                 except:
                     pass
                 out[1:, ...] -= var[item]
-                out[0, ...] = np.NaN
+                out[0, ...] = np.nan
             else:
                 out = var[itemr]
                 out -= var[item]
@@ -508,7 +508,7 @@ def _interp1d(args):
     return f(args[2])
 
 
-def interp3d(data, zcoord, levels, fill_value=np.NaN,
+def interp3d(data, zcoord, levels, fill_value=np.nan,
              use_multiprocessing=True):
     """Interpolate on the first dimension of a 3d var
 
@@ -522,7 +522,7 @@ def interp3d(data, zcoord, levels, fill_value=np.NaN,
       same dims as data, the z coordinates of the data points
     levels: 1darray
       the levels at which to interpolate
-    fill_value : np.NaN or 'extrapolate', optional
+    fill_value : np.nan or 'extrapolate', optional
       how to handle levels below the topography. Default is to mark them
       as invalid, but you might want the have them extrapolated.
     use_multiprocessing: bool