diff --git a/python/src/exactextract/writer.py b/python/src/exactextract/writer.py
index 9bd3442..5cb1b2f 100644
--- a/python/src/exactextract/writer.py
+++ b/python/src/exactextract/writer.py
@@ -44,7 +44,6 @@ def __init__(
         self.array_type = array_type
         self.feature_list = []
         self.map_fields = map_fields or {}
-        self.op_fields = {}
         self.ops = []
         self.remove_temporary_fields = False
 
@@ -56,6 +55,9 @@ def add_operation(self, op):
     def write(self, feature):
         f = JSONFeature()
         feature.copy_to(f)
+        for op in self.ops:
+            if op.name not in f.feature["properties"]:
+                f.feature["properties"][op.name] = None
 
         self._create_map_fields(f)
         self._convert_arrays(f)
@@ -140,16 +142,21 @@ def write(self, feature):
         f = JSONFeature()
         feature.copy_to(f)
 
-        for field_name, value in f.feature["properties"].items():
-            self.fields[field_name].append(value)
-        if "id" in f.feature:
-            self.fields["id"].append(f.feature["id"])
-        if "geometry" in self.fields and "geometry" in f.feature:
-            import shapely
+        props = f.feature["properties"]
 
-            self.fields["geometry"].append(
-                shapely.geometry.shape(f.feature["geometry"])
-            )
+        for field in self.fields:
+            if field == "geometry" and "geometry" in f.feature:
+                import shapely
+
+                self.fields["geometry"].append(
+                    shapely.geometry.shape(f.feature["geometry"])
+                )
+            elif field == "id" and "id" in f.feature:
+                self.fields["id"].append(f.feature["id"])
+            elif field in props:
+                self.fields[field].append(props[field])
+            else:
+                self.fields[field].append(None)
 
     def features(self):
         if "geometry" in self.fields:
@@ -282,7 +289,10 @@ def _set_fields_types(fields_list: list[str], feature: JSONFeature) -> list:
 
         mod_fields_list = []
         for field_name in fields_list:
-            value = feature.get(field_name)
+            try:
+                value = feature.get(field_name)
+            except KeyError:
+                value = None
 
             if type(value) is str:
                 field_type = QVariant.String
@@ -383,9 +393,7 @@ def _collect_fields(feature):
 
                 value = feature.get(field_name)
 
-                if type(value) is str:
-                    field_type = ogr.OFTString
-                elif type(value) is float:
+                if type(value) is float:
                     field_type = ogr.OFTReal
                 elif type(value) is int:
                     field_type = ogr.OFTInteger
@@ -396,6 +404,8 @@ def _collect_fields(feature):
                         field_type = ogr.OFTIntegerList
                     elif value.dtype == np.float64:
                         field_type = ogr.OFTRealList
+                else:
+                    field_type = ogr.OFTString
 
                 ogr_fields[field_name] = ogr.FieldDefn(field_name, field_type)
 
diff --git a/python/tests/test_exact_extract.py b/python/tests/test_exact_extract.py
index 67e3298..d14a360 100644
--- a/python/tests/test_exact_extract.py
+++ b/python/tests/test_exact_extract.py
@@ -510,6 +510,79 @@ def test_gdal_mask_band(tmp_path, libname):
     np.testing.assert_array_equal(values, [1, 2, 3, 4, 5, 6])
 
 
+def test_all_nodata_geojson():
+    data = np.full((3, 3), -999, dtype=np.int32)
+    rast = NumPyRasterSource(data, nodata=-999)
+
+    square = make_rect(0.5, 0.5, 2.5, 2.5)
+    results = exact_extract(rast, square, ["mean", "mode", "variety"], output="geojson")
+
+    props = results[0]["properties"]
+
+    assert math.isnan(props["mean"])
+    assert props["variety"] == 0
+    assert props["mode"] is None
+
+
+def test_all_nodata_pandas():
+    pytest.importorskip("pandas")
+
+    data = np.full((3, 3), -999, dtype=np.int32)
+    rast = NumPyRasterSource(data, nodata=-999)
+
+    square = make_rect(0.5, 0.5, 2.5, 2.5)
+    results = exact_extract(rast, square, ["mean", "mode", "variety"], output="pandas")
+
+    assert math.isnan(results["mean"][0])
+    assert results["variety"][0] == 0
+    assert results["mode"][0] is None
+
+
+def test_all_nodata_qgis():
+
+    pytest.importorskip("qgis.core")
+
+    data = np.full((3, 3), -999, dtype=np.int32)
+    rast = NumPyRasterSource(data, nodata=-999)
+
+    square = make_rect(0.5, 0.5, 2.5, 2.5)
+    results = exact_extract(
+        rast, square, ["mean", "mode", "variety"], output="qgis", include_geom=True
+    )
+
+    props = next(results.getFeatures()).attributeMap()
+
+    assert math.isnan(props["mean"])
+    assert props["variety"] == 0
+    assert props["mode"] is None
+
+
+def test_all_nodata_gdal():
+
+    ogr = pytest.importorskip("osgeo.ogr")
+
+    data = np.array([[1, 1, 1], [-999, -999, -999], [-999, -999, -999]], dtype=np.int32)
+    rast = NumPyRasterSource(data, nodata=-999)
+
+    ds = ogr.GetDriverByName("Memory").CreateDataSource("")
+
+    squares = [make_rect(0.5, 0.5, 2.5, 2.5), make_rect(0.5, 0.5, 1.5, 1.5)]
+    exact_extract(
+        rast,
+        squares,
+        ["mean", "mode", "variety"],
+        output="gdal",
+        include_geom=True,
+        output_options={"dataset": ds},
+    )
+
+    features = [f for f in ds.GetLayer(0)]
+
+    assert math.isnan(features[1]["mean"])
+    assert features[1]["variety"] == 0
+    assert features[1]["mode"] is None
+
+
 def test_default_value():
     data = np.array([[1, 2, 3], [4, -99, -99], [-99, 8, 9]])
     rast = NumPyRasterSource(data, nodata=-99)
diff --git a/python/tests/test_writer.py b/python/tests/test_writer.py
index 236175a..292e97d 100644
--- a/python/tests/test_writer.py
+++ b/python/tests/test_writer.py
@@ -131,15 +131,6 @@ def test_pandas_writer(np_raster_source, point_features):
 
     for f in point_features:
         f.feature["properties"]["mean_result"] = f.feature["id"] * f.feature["id"]
-
-        # we are explicitly declaring columns (add_operation has been called)
-        # but we have an unexpected column
-        with pytest.raises(KeyError, match="type"):
-            w.write(f)
-
-    for f in point_features:
-        del f.feature["properties"]["type"]
-
         w.write(f)
 
     df = w.features()
diff --git a/src/operation.cpp b/src/operation.cpp
index 162ef29..0a65ec0 100644
--- a/src/operation.cpp
+++ b/src/operation.cpp
@@ -1,4 +1,3 @@
-// Copyright (c) 2018-2024 ISciences, LLC.
 // All rights reserved.
 //
 // This software is licensed under the Apache License, Version 2.0 (the "License").
@@ -213,10 +212,11 @@ class OperationImpl : public Operation
             auto&& value = static_cast<const Derived*>(this)->get(s);
 
             if constexpr (is_optional<decltype(value)>) {
-                std::visit([this, &f_out, &value](const auto& m) {
-                    f_out.set(name, value.value_or(m));
-                },
-                           m_missing);
+                if (value.has_value()) {
+                    f_out.set(name, value.value());
+                }
+                // TODO: set result to NaN if this is a floating point type?
+
             } else {
                 f_out.set(name, value);
             }
@@ -410,7 +410,6 @@ Operation::
   , name{ std::move(p_name) }
   , values{ p_values }
   , weights{ p_weights }
-  , m_missing{ get_missing_value() }
   , m_options{ options }
 {
     m_min_coverage = static_cast<float>(extract_arg<double>(options, "min_coverage_frac", 0));
@@ -596,20 +595,6 @@ Operation::create(std::string stat,
 #undef CONSTRUCT
 }
 
-Operation::missing_value_t
-Operation::get_missing_value()
-{
-    const auto& empty_rast = values->read_empty();
-
-    return std::visit([](const auto& r) -> missing_value_t {
-        if (r->has_nodata()) {
-            return r->nodata();
-        }
-        return std::numeric_limits<double>::quiet_NaN();
-    },
-                      empty_rast);
-}
-
 const StatsRegistry::RasterStatsVariant&
 Operation::empty_stats() const
 {
diff --git a/src/operation.h b/src/operation.h
index 968a9e7..bda88d8 100644
--- a/src/operation.h
+++ b/src/operation.h
@@ -177,15 +177,9 @@ class Operation
 
     std::string m_key;
 
-    using missing_value_t = std::variant<float, double, std::int8_t, std::uint8_t, std::int16_t, std::uint16_t, std::int32_t, std::uint32_t, std::int64_t, std::uint64_t>;
-
-    missing_value_t m_missing;
-
     std::optional<std::string> m_default_value;
     std::optional<double> m_default_weight;
 
-    missing_value_t get_missing_value();
-
     const StatsRegistry::RasterStatsVariant& empty_stats() const;
 
     float m_min_coverage;
diff --git a/src/output_writer.h b/src/output_writer.h
index fcd0b74..0ff2328 100644
--- a/src/output_writer.h
+++ b/src/output_writer.h
@@ -31,7 +31,8 @@ class OutputWriter
     /// by the `write` method.
     virtual std::unique_ptr<Feature> create_feature();
 
-    /// Write the provided feature
+    /// Write the provided feature. The feature may not contain
+    /// fields for the results of all Operations.
     virtual void write(const Feature&) = 0;
 
     /// Method to be called for each `Operation` whose results will
diff --git a/test/test_cli.py b/test/test_cli.py
index 56b6ec8..1588b03 100644
--- a/test/test_cli.py
+++ b/test/test_cli.py
@@ -435,7 +435,7 @@ def test_feature_intersecting_nodata(
         "id": "1",
         "metric_count": "0",
         "metric_mean": "nan",
-        "metric_mode": str(nodata) if nodata else "nan",
+        "metric_mode": "",
     }
 
 
diff --git a/test/test_operation.cpp b/test/test_operation.cpp
index 0a643ed..9418992 100644
--- a/test/test_operation.cpp
+++ b/test/test_operation.cpp
@@ -291,7 +291,7 @@ TEMPLATE_TEST_CASE("no error if feature does not intersect raster", "[processor]
 
     const MapFeature& f = writer.m_feature;
     CHECK(f.get_double("count") == 0);
-    CHECK(std::isnan(f.get_double("median")));
+    CHECK_THROWS(f.get_double("median"));
 }
 
 TEST_CASE("progress callback is called once for each feature", "[processor]")