diff --git a/CHANGES.rst b/CHANGES.rst
index 533fc56f..1bbdaad1 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -1,6 +1,11 @@
-1.7.1 (2025-11-2x)
+1.8.0 (2026-01-xx)
 ------------------
 
+New Features
+^^^^^^^^^^^^
+
+- Added uncertainty propagation to ``specreduce.extract.BoxcarExtract``.
+
 Other changes
 ^^^^^^^^^^^^^
 
diff --git a/specreduce/extract.py b/specreduce/extract.py
index 8553124f..61ccd713 100644
--- a/specreduce/extract.py
+++ b/specreduce/extract.py
@@ -178,7 +178,6 @@ class BoxcarExtract(SpecreduceOperation):
     width: float = 5
     disp_axis: int = 1
     crossdisp_axis: int = 0
-    # TODO: should disp_axis and crossdisp_axis be defined in the Trace object?
     mask_treatment: MaskingOption = "apply"
     _valid_mask_treatment_methods = (
         "apply",
@@ -246,19 +245,50 @@ def __call__(
         # the window multiplied by the window width.
         window_weights = _ap_weight_image(trace, width, disp_axis, cdisp_axis, self.image.shape)
 
+        # Extract variance for uncertainty propagation (if available)
+        if self.image.uncertainty is not None:
+            variance = self.image.uncertainty.represent_as(VarianceUncertainty).array
+            orig_uncty_type = type(self.image.uncertainty)
+        else:
+            variance = None
+            orig_uncty_type = None
+
+        extracted_variance = None
         if self.mask_treatment == "apply":
-            image_cleaned = np.where(~self.image.mask, self.image.data * window_weights, 0.0)
+            flux = np.where(~self.image.mask, self.image.data * window_weights, 0.0)
             weights = np.where(~self.image.mask, window_weights, 0.0)
-            spectrum = (
-                image_cleaned.sum(axis=cdisp_axis)
-                / weights.sum(axis=cdisp_axis)
-                * window_weights.sum(axis=cdisp_axis)
-            )
+            weights_sum = weights.sum(axis=cdisp_axis)
+            window_sum = window_weights.sum(axis=cdisp_axis)
+            extracted_flux = flux.sum(axis=cdisp_axis) / weights_sum * window_sum
+
+            if variance is not None:
+                extracted_variance = (
+                    np.where(~self.image.mask, variance * window_weights**2, 0.0).sum(
+                        axis=cdisp_axis
+                    )
+                    / weights_sum**2
+                    * window_sum**2
+                )
         else:
-            image_windowed = np.where(window_weights, self.image.data * window_weights, 0.0)
-            spectrum = np.sum(image_windowed, axis=cdisp_axis)
+            flux = np.where(window_weights, self.image.data * window_weights, 0.0)
+            extracted_flux = flux.sum(axis=cdisp_axis)
+
+            if variance is not None:
+                variance = np.where(window_weights, variance * window_weights**2, 0.0)
+                extracted_variance = np.sum(variance, axis=cdisp_axis)
 
-        return Spectrum(spectrum * self.image.unit, spectral_axis=self.image.spectral_axis)
+        if extracted_variance is not None:
+            spectrum_uncty = VarianceUncertainty(
+                extracted_variance * self.image.unit**2
+            ).represent_as(orig_uncty_type)
+        else:
+            spectrum_uncty = None
+
+        return Spectrum(
+            extracted_flux * self.image.unit,
+            spectral_axis=self.image.spectral_axis,
+            uncertainty=spectrum_uncty,
+        )
 
 
 @dataclass
diff --git a/specreduce/tests/test_extract.py b/specreduce/tests/test_extract.py
index 069996e3..9df36c6c 100644
--- a/specreduce/tests/test_extract.py
+++ b/specreduce/tests/test_extract.py
@@ -2,7 +2,13 @@
 import pytest
 from astropy import units as u
 from astropy.modeling import models
-from astropy.nddata import NDData, VarianceUncertainty, UnknownUncertainty
+from astropy.nddata import (
+    NDData,
+    VarianceUncertainty,
+    StdDevUncertainty,
+    InverseVariance,
+    UnknownUncertainty,
+)
 from astropy.tests.helper import assert_quantity_allclose
 
 from specreduce.background import Background
@@ -119,6 +125,154 @@ def test_boxcar_array_trace(mk_test_img):
     assert np.allclose(spectrum.flux.value, np.full_like(spectrum.flux.value, 75.0))
 
 
+def test_boxcar_uncertainty_propagation():
+    """Test that BoxcarExtract propagates uncertainties correctly."""
+    # Create an image with known uniform flux and uncertainty
+    nrows, ncols = 10, 20
+    flux = np.full((nrows, ncols), 100.0)
+    variance = np.full((nrows, ncols), 4.0)
+
+    img = Spectrum(
+        flux * u.DN,
+        uncertainty=VarianceUncertainty(variance),
+        spectral_axis=np.arange(ncols) * u.pix,
+    )
+
+    trace = FlatTrace(img, nrows // 2)
+    width = 3
+    extracted = BoxcarExtract(img, trace, width=width)()
+
+    # Check uncertainty exists and is the correct type
+    assert extracted.uncertainty is not None
+    assert isinstance(extracted.uncertainty, VarianceUncertainty)
+
+    # For width=3, summing 3 pixels each with variance=4 gives variance=12
+    # (weights are 1.0 for full pixels in aperture)
+    expected_variance = width * 4.0
+    np.testing.assert_allclose(extracted.uncertainty.array, expected_variance, rtol=0.01)
+
+    # Check units are correct (flux_unit^2)
+    assert extracted.uncertainty.unit == u.DN**2
+
+
+def test_boxcar_uncertainty_stddev_type():
+    """Test that StdDevUncertainty type is preserved."""
+    nrows, ncols = 10, 20
+    flux = np.full((nrows, ncols), 100.0)
+    stddev = np.full((nrows, ncols), 2.0)
+
+    img = Spectrum(
+        flux * u.DN, uncertainty=StdDevUncertainty(stddev), spectral_axis=np.arange(ncols) * u.pix
+    )
+
+    trace = FlatTrace(img, nrows // 2)
+    extracted = BoxcarExtract(img, trace, width=3)()
+
+    assert isinstance(extracted.uncertainty, StdDevUncertainty)
+    expected_stddev = np.sqrt(3 * 4.0)
+    np.testing.assert_allclose(extracted.uncertainty.array, expected_stddev, rtol=0.01)
+    assert extracted.uncertainty.unit == u.DN
+
+
+def test_boxcar_uncertainty_inverse_variance_type():
+    """Test that InverseVariance type is preserved."""
+    nrows, ncols = 10, 20
+    flux = np.full((nrows, ncols), 100.0)
+    ivar = np.full((nrows, ncols), 0.25)
+
+    img = Spectrum(
+        flux * u.DN, uncertainty=InverseVariance(ivar), spectral_axis=np.arange(ncols) * u.pix
+    )
+
+    trace = FlatTrace(img, nrows // 2)
+    extracted = BoxcarExtract(img, trace, width=3)()
+
+    assert isinstance(extracted.uncertainty, InverseVariance)
+    expected_ivar = 1.0 / (3 * 4.0)
+    np.testing.assert_allclose(extracted.uncertainty.array, expected_ivar, rtol=0.01)
+    assert extracted.uncertainty.unit == 1 / u.DN**2
+
+
+def test_boxcar_uncertainty_with_nonfinite():
+    """Test uncertainty propagation with nonfinite values in data."""
+    nrows, ncols = 10, 20
+    flux = np.full((nrows, ncols), 100.0)
+    variance = np.full((nrows, ncols), 4.0)
+    flux[5, 5] = np.nan
+
+    img = Spectrum(
+        flux * u.DN,
+        uncertainty=VarianceUncertainty(variance),
+        spectral_axis=np.arange(ncols) * u.pix,
+    )
+
+    trace = FlatTrace(img, 5)
+
+    # With mask_treatment="apply", NaN pixels excluded from both flux and variance
+    extracted = BoxcarExtract(img, trace, width=3, mask_treatment="apply")()
+    assert extracted.uncertainty is not None
+    # Column 5 has only 2 valid pixels instead of 3
+    # Variance formula: (sum w^2 sigma^2) / (sum w)^2 * W^2 = (2*4) / 4 * 9 = 18
+    np.testing.assert_allclose(extracted.uncertainty.array[5], 18.0, rtol=0.01)
+    # Other columns should have normal variance=12
+    np.testing.assert_allclose(extracted.uncertainty.array[0], 12.0, rtol=0.01)
+
+    # With mask_treatment="ignore", NaN in flux propagates to flux but variance
+    # is computed from the variance array which is finite
+    extracted_ignore = BoxcarExtract(img, trace, width=3, mask_treatment="ignore")()
+    assert np.isnan(extracted_ignore.flux.value[5])
+    # Variance is still computed from a finite variance array
+    np.testing.assert_allclose(extracted_ignore.uncertainty.array[5], 12.0, rtol=0.01)
+
+    # Test that NaN in a variance array does propagate
+    variance_with_nan = variance.copy()
+    variance_with_nan[5, 5] = np.nan
+    img_var_nan = Spectrum(
+        np.full((nrows, ncols), 100.0) * u.DN,
+        uncertainty=VarianceUncertainty(variance_with_nan),
+        spectral_axis=np.arange(ncols) * u.pix,
+    )
+    extracted_var_nan = BoxcarExtract(img_var_nan, trace, width=3, mask_treatment="ignore")()
+    assert np.isnan(extracted_var_nan.uncertainty.array[5])
+
+
+def test_boxcar_uncertainty_with_mask():
+    """Test uncertainty propagation with masked pixels."""
+    nrows, ncols = 10, 20
+    flux = np.full((nrows, ncols), 100.0)
+    variance = np.full((nrows, ncols), 4.0)
+    mask = np.zeros((nrows, ncols), dtype=bool)
+
+    # Mask one pixel in the extraction aperture at column 3
+    mask[5, 3] = True
+
+    img = Spectrum(
+        flux * u.DN,
+        uncertainty=VarianceUncertainty(variance),
+        mask=mask,
+        spectral_axis=np.arange(ncols) * u.pix,
+    )
+
+    trace = FlatTrace(img, 5)
+    extracted = BoxcarExtract(img, trace, width=3, mask_treatment="apply")()
+    assert extracted.uncertainty is not None
+    np.testing.assert_allclose(extracted.uncertainty.array[3], 18.0, rtol=0.01)
+    np.testing.assert_allclose(extracted.uncertainty.array[0], 12.0, rtol=0.01)
+
+
+def test_boxcar_no_input_uncertainty():
+    """Test that no uncertainty is returned when input has none."""
+    from astropy.nddata import CCDData
+
+    nrows, ncols = 10, 20
+    flux = np.full((nrows, ncols), 100.0)
+
+    img = CCDData(flux, unit=u.DN)
+    trace = FlatTrace(img, nrows // 2)
+    extracted = BoxcarExtract(img, trace, width=3)()
+    assert extracted.uncertainty is None
+
+
 def test_horne_image_validation(mk_test_img):
     #
     # Test HorneExtract scenarios specific to its use with an image of