diff --git a/CHANGELOG.md b/CHANGELOG.md
index 6fbdeb2..08da392 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,11 @@
 # `cdshealpix-python` Change Log
 
+## Unreleased
+
+### Fix
+
+* bilinear_interpolation now accepts longitudes and latitudes with nan values (will be reflected by a masked value in the output) [#21]
+
 ## 0.6.5
 
 Released 2023-11-28
diff --git a/python/cdshealpix/nested/healpix.py b/python/cdshealpix/nested/healpix.py
index f47da09..7355ea0 100644
--- a/python/cdshealpix/nested/healpix.py
+++ b/python/cdshealpix/nested/healpix.py
@@ -1004,31 +1004,28 @@ def bilinear_interpolation(lon, lat, depth, num_threads=0):
     if not (isinstance(lat, Latitude)):
         raise ValueError("`lat` must be of type `astropy.coordinates.Latitude`")
 
-    # We could have continued to use `.to_value(u.rad)` instead of `.rad`.
-    # Although `to_value` is more generical (method of Quantity),
-    # Longitude/Latitude ensure that the values given to the contain are in the correct ranges.
-    lon = np.atleast_1d(lon.rad)
-    lat = np.atleast_1d(lat.rad)
-
     if depth < 0 or depth > 29:
         raise ValueError("Depth must be in the [0, 29] closed range")
 
+    lon = np.atleast_1d(lon.rad)
+    lat = np.atleast_1d(lat.rad)
+
     if lon.shape != lat.shape:
         raise ValueError(
             "The number of longitudes does not match with the number of latitudes given"
         )
 
-    # Useless since we test isinstance at the beginning of the function
-    # if ((lat < np.pi/2.0) | (lat > np.pi/2.0)).any():
-    #    raise ValueError("Lat must be in [-pi/2, 2pi/2]")
-
     num_coords = lon.shape
 
     # Retrieve nan and infinite values
-    mask_lon_invalid = np.isnan(lon) | ~np.isfinite(lon)
-    mask_lat_invalid = np.isnan(lat) | ~np.isfinite(lat)
+    mask_lon_invalid = np.ma.masked_invalid(lon).mask
+    mask_lat_invalid = np.ma.masked_invalid(lat).mask
+
     mask_invalid = mask_lon_invalid | mask_lat_invalid
 
+    lon = np.ma.fix_invalid(lon, mask_invalid, fill_value=0 * u.deg)
+    lat = np.ma.fix_invalid(lat, mask_invalid, fill_value=0 * u.deg)
+
     mask_invalid = np.repeat(mask_invalid[:, np.newaxis], 4, axis=mask_invalid.ndim)
 
     ipix = np.empty(shape=num_coords + (4,), dtype=np.uint64)
diff --git a/python/cdshealpix/tests/test_nested_healpix.py b/python/cdshealpix/tests/test_nested_healpix.py
index 068809b..a6b439e 100644
--- a/python/cdshealpix/tests/test_nested_healpix.py
+++ b/python/cdshealpix/tests/test_nested_healpix.py
@@ -415,10 +415,11 @@ def test_bilinear_interpolation(depth):
     assert ((ipix >= 0) & (ipix < 12 * (4**depth))).all()
 
 
-@pytest.mark.parametrize("depth", [5, 0, 7, 12, 20, 29])
-def test_bilinear_interpolation2(depth):
-    lon = Longitude([10, 25, 0], u.deg)
-    lat = Latitude([5, 10, 45], u.deg)
+def test_bilinear_interpolation_accepts_nan():
+    lon = Longitude([10, np.nan], unit="deg")
+    lat = Latitude([5, np.nan], unit="deg")
     depth = 5
 
-    ipix, weights = bilinear_interpolation(lon, lat, depth)
+    ipix, _ = bilinear_interpolation(lon, lat, depth)
+
+    assert np.all(ipix.mask == [[False, False, False, False], [True, True, True, True]])