diff --git a/reproject/mosaicking/coadd.py b/reproject/mosaicking/coadd.py
index be6573c5f..5a08cf7a8 100644
--- a/reproject/mosaicking/coadd.py
+++ b/reproject/mosaicking/coadd.py
@@ -1,5 +1,7 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 
+import warnings
+
 import numpy as np
 
 from ..utils import parse_input_data, parse_input_weights, parse_output_projection
@@ -71,7 +73,7 @@ def reproject_and_coadd(
         `~astropy.io.fits.HDUList` instance, specifies the HDU to use.
     reproject_function : callable
         The function to use for the reprojection
-    combine_function : { 'mean', 'sum', 'median' }
+    combine_function : { 'mean', 'sum', 'median', 'min' }
         The type of function to use for combining the values into the final
         image.
     match_background : bool
@@ -92,7 +94,7 @@ def reproject_and_coadd(
 
     # Validate inputs
 
-    if combine_function not in ("mean", "sum", "median"):
+    if combine_function not in ("mean", "sum", "median", "min"):
         raise ValueError("combine_function should be one of mean/sum/median")
 
     if reproject_function is None:
@@ -228,12 +230,72 @@ def reproject_and_coadd(
         if combine_function == "mean":
             with np.errstate(invalid="ignore"):
                 final_array /= final_footprint
-
+    elif combine_function == "min":
+        final_array = np.ones(shape_out) * np.inf
+        for array in arrays:
+            array.array[array.footprint == 0] = np.nan
+            final_array[array.view_in_original_array] = np.fmin(
+                final_array[array.view_in_original_array], array.array * array.footprint
+            )
+            final_footprint[array.view_in_original_array] += array.footprint
     elif combine_function == "median":
+        BLOCK_SIZE = 100
+        for i in np.arange(0, shape_out[0], BLOCK_SIZE):
+            for j in np.arange(0, shape_out[1], BLOCK_SIZE):
+                blocks = []
+                foots = []
+                for array in arrays:
+                    part, foot = _get_block(array, i, j, BLOCK_SIZE)
+                    part[foot == 0] = np.nan
+                    foots.append(foot)
+                    blocks.append(part)
+                with warnings.catch_warnings():
+                    warnings.filterwarnings(
+                        "ignore", category=RuntimeWarning, message="All-NaN slice encountered"
+                    )
+                    patch_shape = final_array[i : i + BLOCK_SIZE, j : j + BLOCK_SIZE].shape
+                    final_array[i : i + BLOCK_SIZE, j : j + BLOCK_SIZE] = np.nanmedian(
+                        blocks, axis=0
+                    )[: patch_shape[0], : patch_shape[1]]
+                    final_footprint[i : i + BLOCK_SIZE, j : j + BLOCK_SIZE] = np.nansum(
+                        foots, axis=0
+                    )[: patch_shape[0], : patch_shape[1]]
+    return final_array, final_footprint
 
-        # Here we need to operate in chunks since we could otherwise run
-        # into memory issues
 
-        raise NotImplementedError("combine_function='median' is not yet implemented")
+def _get_block(array, block_index_x, block_index_y, block_size=100):
+    """Get square block of proper dimensions from an array.
 
-    return final_array, final_footprint
+    Args:
+        array (ReprojectedArraySubset): Input array to select data from.
+        block_index_x (int): X of the lower left corner
+        block_index_y (int): Y of the lower left corner
+        block_size (int, optional): block size in X and Y. Defaults to 100.
+
+    Returns:
+        np.ndarray, np.ndarray: block and its integer-valued footprint.
+    """
+    result = np.ones((block_size, block_size)) * np.nan
+    result_footprint = np.zeros((block_size, block_size))
+    jlower = block_index_x
+    jupper = block_index_x + block_size
+    ilower = block_index_y
+    iupper = block_index_y + block_size
+    if (
+        jlower <= array.jmax
+        and jupper >= array.jmin
+        and ilower <= array.imax
+        and iupper >= array.imin
+    ):
+        jlower1 = max(0, array.jmin - jlower)
+        jupper1 = min(jupper - jlower, array.jmax - jlower)
+        ilower1 = max(0, array.imin - ilower)
+        iupper1 = min(iupper - ilower, array.imax - ilower)
+        local = (slice(jlower1, jupper1), slice(ilower1, iupper1))
+        external = (
+            slice(jlower1 + jlower - array.jmin, jupper1 + jlower - array.jmin),
+            slice(ilower1 + ilower - array.imin, iupper1 + ilower - array.imin),
+        )
+        result[local] = array.array[external]
+        result_footprint[local] = array.footprint[external]
+    return result, result_footprint
diff --git a/reproject/mosaicking/tests/test_coadd.py b/reproject/mosaicking/tests/test_coadd.py
index 95e552fc6..79f36f212 100644
--- a/reproject/mosaicking/tests/test_coadd.py
+++ b/reproject/mosaicking/tests/test_coadd.py
@@ -78,7 +78,7 @@ def _overlapping_views(self):
 
         return views
 
-    @pytest.mark.parametrize("combine_function", ["mean", "sum"])
+    @pytest.mark.parametrize("combine_function", ["mean", "sum", "min", "median"])
     def test_coadd_no_overlap(self, combine_function, reproject_function):
 
         # Make sure that if all tiles are exactly non-overlapping, and