Repository: natcap/invest
Pull request #1789 (2025-02-18 16:04:32)

auto merge main into release/*
Merge main into release/3.15.0

This PR was automatically generated in response to a push to `main`,
and is a chance to review any changes that will be included in the release
branch before merging.  Under most circumstances, this PR will probably be
a formality.  However, there are a few cases where we may need to do some
extra work to make sure `release/3.15.0` contains what it should after the
merge:

## There is a merge conflict in this PR

1. Leave a comment on this PR about the merge conflict and close the PR.
2. In your fork, `git pull` the upstream `main` and `release/3.15.0`.
3. Make a new `pr-resolution` branch off of `release/3.15.0`:
   ```shell
   $ git checkout main
   $ git pull upstream main  # Include the latest changes on the upstream main
   $ git checkout release/3.15.0
   $ git pull upstream release/3.15.0  # Include the latest changes on the upstream release
   $ git checkout -b pr-resolution
   $ git merge main
   ```
4. Resolve the conflicts locally
5. Commit the changes to `pr-resolution`.
6. Create a PR from `pr-resolution` into `release/3.15.0`, and include a link
   to the origin PR in the description.
7. When the PR is complete, delete the `pr-resolution` branch.  That will
   help us avoid confusion and extra work down the road when we do this again.

## This PR contains content that should not be in `release/3.15.0`

1. Leave a comment on this PR about the content that should not be included
   and close the PR.
2. In your fork, `git pull` the upstream `main` and `release/3.15.0`.
3. Make a new `pr-resolution` branch off of `release/3.15.0`:
   ```shell
   $ git checkout main
   $ git pull upstream main  # Include the latest changes on the upstream main
   $ git checkout release/3.15.0
   $ git pull upstream release/3.15.0  # Include the latest changes on the upstream release
   $ git checkout -b pr-resolution
   $ git merge main
   ```
4. Handle the content that should not end up in `release/3.15.0` however it
   needs to be handled.
5. Commit the updated content to `pr-resolution`.
6. Create a PR from `pr-resolution` into `release/3.15.0`, and include a link
   to the origin PR in the description.
7. When the PR is complete, delete the `pr-resolution` branch.  That will
   help us avoid confusion and extra work down the road when we do this again.

## What happens if we accidentally merge something we shouldn't?

There are several possibilities for recovery if we get to such a state.

1. A merge can be undone through the github interface if the error is caught
   directly after the PR is merged.
2. If we're commits in past the erroneous merge, create a branch off of
   `release/3.15.0`, back out of the changes or edit files needed to resolve
   the issue, and PR the branch back into `release/3.15.0`.

### Why was this PR created?

The workflow defining this PR is located at
`.github/workflows/auto-pr-from-main-into-releases.yml`.  In short, this PR
was created because there was a push to `main` that triggered this
workflow.  Some events that can trigger this include:

* Other pull requests being merged into `main`
* Automated releases on `main`
* Any manual push to `main`, if ever that happens (which shouldn't be the
  case given our branch protections)




Patch:
diff --git a/.github/workflows/release-part-1.yml b/.github/workflows/release-part-1.yml
index b106bcec2..9656dc3d7 100644
--- a/.github/workflows/release-part-1.yml
+++ b/.github/workflows/release-part-1.yml
@@ -22,6 +22,7 @@ env:
 
 jobs:
   start-release:
+    if: github.repository_owner == 'natcap'
     runs-on: ubuntu-latest
     steps:
       # skip scheduled runs that aren't the first tuesday of the month
@@ -44,11 +45,11 @@ jobs:
         run: pip install rst2html5 setuptools_scm
 
       - name: Get current version number (scheduled runs only)
-        if: github.event_name == 'schedule'
+        if: env.RUN == 'true' && github.event_name == 'schedule'
         run: echo "CURRENT_VERSION=$(python -m setuptools_scm)" > $GITHUB_ENV
 
       - name: Get next bugfix version number (scheduled runs only)
-        if: github.event_name == 'schedule'
+        if: env.RUN == 'true' && github.event_name == 'schedule'
         shell: python
         run: |
           import os
diff --git a/HISTORY.rst b/HISTORY.rst
index bf685cdb8..8cb6256d8 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -9,7 +9,7 @@
   - Crop Pollination
   - Crop Production
   - DelineateIt
-  - Forest Carbon Edge Effect
+  - Forest Carbon Edge Effects
   - Globio
   - Habitat Quality
   - HRA
@@ -48,10 +48,6 @@ Unreleased Changes
     * InVEST's windows binaries are now distributed once again with a valid
       signature, signed by Stanford University.
       https://github.com/natcap/invest/issues/1580
-* Workbench
-    * Auto-scrolling of log output is halted on user-initiated scrolling,
-      enabling easier inspection of log output while a model is running
-      (`InVEST #1533 <https://github.com/natcap/invest/issues/1533>`_).
 * Annual Water Yield
     * Fixed an issue where the model would crash if the valuation table was
       provided, but the demand table was not.  Validation will now warn about
@@ -62,42 +58,7 @@ Unreleased Changes
     * Updated styling of the HTML report generated by the carbon model, for
       visual consistency with the Workbench (`InVEST #1732
       <https://github.com/natcap/invest/issues/1732>`_).
-    * Raster outputs that previously contained per-pixel values (e.g., Mg/px)
-      now contain per-hectare values (e.g., Mg/ha). (`InVEST #1270
-      <https://github.com/natcap/invest/issues/1270>`_).
-* Coastal Blue Carbon
-    * The ``code`` column in the model's biophysical table input, as well as
-      the ``code`` column in the preprocessor's LULC lookup table input and
-      ``carbon_pool_transient_template`` output, have been renamed ``lucode``,
-      for consistency with other InVEST models (`InVEST #1249
-      <https://github.com/natcap/invest/issues/1249>`_).
-* Crop Production
-    * Raster outputs that previously contained per-pixel values (e.g., Mg/px)
-      now contain per-hectare values (e.g., Mg/ha). This change affects both
-      the Percentile and Regression models (`InVEST #1270
-      <https://github.com/natcap/invest/issues/1270>`_).
-* Forest Carbon Edge Effect
-    * Raster outputs that previously contained per-pixel values (e.g., Mg/px)
-      now contain per-hectare values (e.g., Mg/ha). (`InVEST #1270
-      <https://github.com/natcap/invest/issues/1270>`_).
-* Habitat Quality
-    * The ``lulc`` column in the sensitivity table input, and the ``lulc_code``
-      column in the rarity table outputs, have been renamed ``lucode``, for
-      consistency with other InVEST models (`InVEST #1249
-      <https://github.com/natcap/invest/issues/1249>`_).
-* NDR
-    * Align rasters to the grid of the DEM raster
-      (`#1488 <https://github.com/natcap/invest/issues/1488>`_).
-    * Raster outputs that previously contained per-pixel values (e.g., kg/px)
-      now contain per-hectare values (e.g., kg/ha). (`InVEST #1270
-      <https://github.com/natcap/invest/issues/1270>`_).
-* SDR
-    * Raster outputs that previously contained per-pixel values (e.g., Mg/px)
-      now contain per-hectare values (e.g., Mg/ha). (`InVEST #1270
-      <https://github.com/natcap/invest/issues/1270>`_).
 * Urban Cooling
-    * Align rasters to the grid of the LULC raster, rather than the ET0 raster
-      (`#1488 <https://github.com/natcap/invest/issues/1488>`_).
     * Updated the documentation for the ``mean_t_air`` attribute of the
       ``buildings_with_stats.shp`` output to clarify how the value is
       calculated.  https://github.com/natcap/invest/issues/1746
diff --git a/Makefile b/Makefile
index 472e2b363..a400384fe 100644
--- a/Makefile
+++ b/Makefile
@@ -2,15 +2,15 @@
 DATA_DIR := data
 GIT_SAMPLE_DATA_REPO        := https://bitbucket.org/natcap/invest-sample-data.git
 GIT_SAMPLE_DATA_REPO_PATH   := $(DATA_DIR)/invest-sample-data
-GIT_SAMPLE_DATA_REPO_REV    := dfcb820269fa50f65eed55b830314c853c2eb943
+GIT_SAMPLE_DATA_REPO_REV    := 0f8b41557753dad3670ba8220f41650b51435a93
 
 GIT_TEST_DATA_REPO          := https://bitbucket.org/natcap/invest-test-data.git
 GIT_TEST_DATA_REPO_PATH     := $(DATA_DIR)/invest-test-data
-GIT_TEST_DATA_REPO_REV      := 4eb5a2fac4818a37663780c4cbdff1b41e7019be
+GIT_TEST_DATA_REPO_REV      := 324abde73e1d770ad75921466ecafd1ec6297752
 
 GIT_UG_REPO                 := https://github.com/natcap/invest.users-guide
 GIT_UG_REPO_PATH            := doc/users-guide
-GIT_UG_REPO_REV             := 79cc80687052c870bbba67e9ee2e68e67dcec145
+GIT_UG_REPO_REV             := fd3194f35bc1652a93cf1c0241f98be1e7c9d43d
 
 ENV = "./env"
 ifeq ($(OS),Windows_NT)
diff --git a/src/natcap/invest/carbon.py b/src/natcap/invest/carbon.py
index 944cda977..9a899558a 100644
--- a/src/natcap/invest/carbon.py
+++ b/src/natcap/invest/carbon.py
@@ -4,6 +4,7 @@ import codecs
 import logging
 import os
 import time
+from functools import reduce
 
 from osgeo import gdal
 import numpy
@@ -26,7 +27,7 @@ CARBON_OUTPUTS = {
             "scenario, mapped from the Carbon Pools table to the LULC."),
         "bands": {1: {
             "type": "number",
-            "units": u.metric_ton/u.hectare
+            "units": u.metric_ton/u.pixel
         }}
     } for pool, pool_name in [
         ('above', 'aboveground'),
@@ -153,7 +154,7 @@ MODEL_SPEC = {
             "about": gettext(
                 "The calendar year of the future scenario depicted in the "
                 "future LULC map. Required if Run Valuation model is selected."),
-            "name": gettext("future LULC year")
+            "name": f"future LULC year"
         },
         "do_valuation": {
             "type": "boolean",
@@ -200,14 +201,14 @@ MODEL_SPEC = {
             "about": "Raster showing the amount of carbon stored in each pixel for the current scenario. It is a sum of all of the carbon pools provided by the biophysical table.",
             "bands": {1: {
                 "type": "number",
-                "units": u.metric_ton/u.hectare
+                "units": u.metric_ton/u.pixel
             }}
         },
         "tot_c_fut.tif": {
             "about": "Raster showing the amount of carbon stored in each pixel for the future scenario. It is a sum of all of the carbon pools provided by the biophysical table.",
             "bands": {1: {
                 "type": "number",
-                "units": u.metric_ton/u.hectare
+                "units": u.metric_ton/u.pixel
             }},
             "created_if": "lulc_fut_path"
         },
@@ -215,7 +216,7 @@ MODEL_SPEC = {
             "about": "Raster showing the amount of carbon stored in each pixel for the REDD scenario. It is a sum of all of the carbon pools provided by the biophysical table.",
             "bands": {1: {
                 "type": "number",
-                "units": u.metric_ton/u.hectare
+                "units": u.metric_ton/u.pixel
             }},
             "created_if": "lulc_redd_path"
         },
@@ -223,7 +224,7 @@ MODEL_SPEC = {
             "about": "Raster showing the difference in carbon stored between the future landscape and the current landscape. In this map some values may be negative and some positive. Positive values indicate sequestered carbon, negative values indicate carbon that was lost.",
             "bands": {1: {
                 "type": "number",
-                "units": u.metric_ton/u.hectare
+                "units": u.metric_ton/u.pixel
             }},
             "created_if": "lulc_fut_path"
         },
@@ -231,7 +232,7 @@ MODEL_SPEC = {
             "about": "Raster showing the difference in carbon stored between the REDD landscape and the current landscape. In this map some values may be negative and some positive. Positive values indicate sequestered carbon, negative values indicate carbon that was lost.",
             "bands": {1: {
                 "type": "number",
-                "units": u.metric_ton/u.hectare
+                "units": u.metric_ton/u.pixel
             }},
             "created_if": "lulc_redd_path"
         },
@@ -239,7 +240,7 @@ MODEL_SPEC = {
             "about": "Rasters showing the economic value of carbon sequestered between the current and the future landscape dates.",
             "bands": {1: {
                 "type": "number",
-                "units": u.currency/u.hectare
+                "units": u.currency/u.pixel
             }},
             "created_if": "lulc_fut_path"
         },
@@ -247,7 +248,7 @@ MODEL_SPEC = {
             "about": "Rasters showing the economic value of carbon sequestered between the current and the REDD landscape dates.",
             "bands": {1: {
                 "type": "number",
-                "units": u.currency/u.hectare
+                "units": u.currency/u.pixel
             }},
             "created_if": "lulc_redd_path"
         },
@@ -296,7 +297,6 @@ _TMP_BASE_FILES = {
 # -1.0 since carbon stocks are 0 or greater
 _CARBON_NODATA = -1.0
 
-
 def execute(args):
     """Carbon.
 
@@ -541,15 +541,17 @@ def _generate_carbon_map(
     Args:
         lulc_path (string): landcover raster with integer pixels.
         out_carbon_stock_path (string): path to output raster that will have
-            pixels with carbon storage values in them with units of Mg/ha.
+            pixels with carbon storage values in them with units of Mg*C
         carbon_pool_by_type (dict): a dictionary that maps landcover values
             to carbon storage densities per area (Mg C/Ha).
 
     Returns:
         None.
     """
+    lulc_info = pygeoprocessing.get_raster_info(lulc_path)
+    pixel_area = abs(numpy.prod(lulc_info['pixel_size']))
     carbon_stock_by_type = dict([
-        (lulcid, stock)
+        (lulcid, stock * pixel_area / 10**4)
         for lulcid, stock in carbon_pool_by_type.items()])
 
     reclass_error_details = {
@@ -702,17 +704,14 @@ def _generate_report(raster_file_set, model_args, file_registry):
             '<table><thead><tr><th>Description</th><th>Value</th><th>Units'
             '</th><th>Raw File</th></tr></thead><tbody>')
 
-        carbon_units = 'metric tons'
-
         # value lists are [sort priority, description, statistic, units]
         report = [
-            (file_registry['tot_c_cur'], 'Total cur', carbon_units),
-            (file_registry['tot_c_fut'], 'Total fut', carbon_units),
-            (file_registry['tot_c_redd'], 'Total redd', carbon_units),
-            (file_registry['delta_cur_fut'], 'Change in C for fut',
-             carbon_units),
+            (file_registry['tot_c_cur'], 'Total cur', 'Mg of C'),
+            (file_registry['tot_c_fut'], 'Total fut', 'Mg of C'),
+            (file_registry['tot_c_redd'], 'Total redd', 'Mg of C'),
+            (file_registry['delta_cur_fut'], 'Change in C for fut', 'Mg of C'),
             (file_registry['delta_cur_redd'],
-             'Change in C for redd', carbon_units),
+             'Change in C for redd', 'Mg of C'),
             (file_registry['npv_fut'],
              'Net present value from cur to fut', 'currency units'),
             (file_registry['npv_redd'],
@@ -721,17 +720,11 @@ def _generate_report(raster_file_set, model_args, file_registry):
 
         for raster_uri, description, units in report:
             if raster_uri in raster_file_set:
-                total = _accumulate_totals(raster_uri)
-                raster_info = pygeoprocessing.get_raster_info(raster_uri)
-                pixel_area = abs(numpy.prod(raster_info['pixel_size']))
-                # Since each pixel value is in Mg/ha, ``total`` is in (Mg/ha * px) = Mg•px/ha.
-                # Adjusted sum = ([total] Mg•px/ha) * ([pixel_area] m^2 / 1 px) * (1 ha / 10000 m^2) = Mg.
-                summary_stat = total * pixel_area / 10000
+                summary_stat = _accumulate_totals(raster_uri)
                 report_doc.write(
-                    '<tr><td>%s</td><td class="number" data-summary-stat="%s">'
-                    '%.2f</td><td>%s</td><td>%s</td></tr>' % (
-                        description, description, summary_stat, units,
-                        raster_uri))
+                    '<tr><td>%s</td><td class="number">%.2f</td><td>%s</td>'
+                    '<td>%s</td></tr>' % (
+                        description, summary_stat, units, raster_uri))
         report_doc.write('</tbody></table></body></html>')
 
 
diff --git a/src/natcap/invest/coastal_blue_carbon/coastal_blue_carbon.py b/src/natcap/invest/coastal_blue_carbon/coastal_blue_carbon.py
index 5423908ec..7e118cdd5 100644
--- a/src/natcap/invest/coastal_blue_carbon/coastal_blue_carbon.py
+++ b/src/natcap/invest/coastal_blue_carbon/coastal_blue_carbon.py
@@ -207,9 +207,9 @@ MODEL_SPEC = {
         "biophysical_table_path": {
             "name": gettext("biophysical table"),
             "type": "csv",
-            "index_col": "lucode",
+            "index_col": "code",
             "columns": {
-                "lucode": {
+                "code": {
                     "type": "integer",
                     "about": gettext(
                         "The LULC code that represents this LULC "
@@ -1679,7 +1679,7 @@ def _track_disturbance(
         disturbed_carbon_volume[:] = NODATA_FLOAT32_MIN
         disturbed_carbon_volume[
             ~pygeoprocessing.array_equals_nodata(disturbance_magnitude_matrix,
-                                                 NODATA_FLOAT32_MIN)] = 0.0
+                           NODATA_FLOAT32_MIN)] = 0.0
 
         if year_of_disturbance_band:
             known_transition_years_matrix = (
diff --git a/src/natcap/invest/coastal_blue_carbon/preprocessor.py b/src/natcap/invest/coastal_blue_carbon/preprocessor.py
index d78191711..03fda7714 100644
--- a/src/natcap/invest/coastal_blue_carbon/preprocessor.py
+++ b/src/natcap/invest/coastal_blue_carbon/preprocessor.py
@@ -36,9 +36,9 @@ MODEL_SPEC = {
                 "A table mapping LULC codes from the snapshot rasters to the "
                 "corresponding LULC class names, and whether or not the "
                 "class is a coastal blue carbon habitat."),
-            "index_col": "lucode",
+            "index_col": "code",
             "columns": {
-                "lucode": {
+                "code": {
                     "type": "integer",
                     "about": gettext(
                         "LULC code. Every value in the "
@@ -86,9 +86,9 @@ MODEL_SPEC = {
             "index_col": "lulc-class",
             "columns": {
                 "lulc-class": {
-                    "type": "freestyle_string",
+                    "type": "integer",
                     "about": gettext(
-                        "LULC class names matching those in the biophysical "
+                        "LULC codes matching the codes in the biophysical "
                         "table.")},
                 "[LULC]": {
                     "type": "option_string",
@@ -114,7 +114,7 @@ MODEL_SPEC = {
                 "Table mapping each LULC type to impact and accumulation "
                 "information. This is a template that you will fill out to "
                 "create the biophysical table input to the main model."),
-            "index_col": "lucode",
+            "index_col": "code",
             "columns": {
                 **BIOPHYSICAL_COLUMNS_SPEC,
                 # remove "expression" property which doesn't go in output spec
@@ -394,7 +394,7 @@ def _create_biophysical_table(landcover_df, target_biophysical_table_path):
             # have commas between fields.
             row = []
             for colname in target_column_names:
-                if colname == 'lucode':
+                if colname == 'code':
                     row.append(str(lulc_code))
                 else:
                     try:
diff --git a/src/natcap/invest/crop_production_percentile.py b/src/natcap/invest/crop_production_percentile.py
index fefc3de14..ea4dc057f 100644
--- a/src/natcap/invest/crop_production_percentile.py
+++ b/src/natcap/invest/crop_production_percentile.py
@@ -761,6 +761,7 @@ def execute(args):
                     args['landcover_raster_path'],
                     interpolated_yield_percentile_raster_path,
                     crop_lucode,
+                    pixel_area_ha,
                     percentile_crop_production_raster_path),
                 target_path_list=[percentile_crop_production_raster_path],
                 dependent_task_list=[
@@ -835,7 +836,7 @@ def execute(args):
             args=([(args['landcover_raster_path'], 1),
                    (interpolated_observed_yield_raster_path, 1),
                    (observed_yield_nodata, 'raw'), (landcover_nodata, 'raw'),
-                   (crop_lucode, 'raw')],
+                   (crop_lucode, 'raw'), (pixel_area_ha, 'raw')],
                   _mask_observed_yield_op, observed_production_raster_path,
                   gdal.GDT_Float32, observed_yield_nodata),
             target_path_list=[observed_production_raster_path],
@@ -852,7 +853,7 @@ def execute(args):
         output_dir, 'result_table%s.csv' % file_suffix)
 
     crop_names = crop_to_landcover_df.index.to_list()
-    _ = task_graph.add_task(
+    tabulate_results_task = task_graph.add_task(
         func=tabulate_results,
         args=(nutrient_df, yield_percentile_headers,
               crop_names, pixel_area_ha,
@@ -869,14 +870,14 @@ def execute(args):
             output_dir, _AGGREGATE_VECTOR_FILE_PATTERN % (file_suffix))
         aggregate_results_table_path = os.path.join(
             output_dir, _AGGREGATE_TABLE_FILE_PATTERN % file_suffix)
-        _ = task_graph.add_task(
+        aggregate_results_task = task_graph.add_task(
             func=aggregate_to_polygons,
             args=(args['aggregate_polygon_path'],
                   target_aggregate_vector_path,
                   landcover_raster_info['projection_wkt'],
                   crop_names, nutrient_df,
-                  yield_percentile_headers, pixel_area_ha, output_dir,
-                  file_suffix, aggregate_results_table_path),
+                  yield_percentile_headers, output_dir, file_suffix,
+                  aggregate_results_table_path),
             target_path_list=[target_aggregate_vector_path,
                               aggregate_results_table_path],
             dependent_task_list=dependent_task_list,
@@ -887,14 +888,14 @@ def execute(args):
 
 
 def calculate_crop_production(lulc_path, yield_path, crop_lucode,
-                              target_path):
+                              pixel_area_ha, target_path):
     """Calculate crop production for a particular crop.
 
     The resulting production value is:
 
     - nodata, where either the LULC or yield input has nodata
     - 0, where the LULC does not match the given LULC code
-    - yield (in Mg/ha), where the given LULC code exists
+    - yield * pixel area, where the given LULC code exists
 
     Args:
         lulc_path (str): path to a raster of LULC codes
@@ -902,6 +903,7 @@ def calculate_crop_production(lulc_path, yield_path, crop_lucode,
             by ``crop_lucode``, in units per hectare
         crop_lucode (int): LULC code that identifies the crop of interest in
             the ``lulc_path`` raster.
+        pixel_area_ha (number): Pixel area in hectares for both input rasters
         target_path (str): Path to write the output crop production raster
 
     Returns:
@@ -909,7 +911,7 @@ def calculate_crop_production(lulc_path, yield_path, crop_lucode,
     """
     pygeoprocessing.raster_map(
         op=lambda lulc, _yield: numpy.where(
-            lulc == crop_lucode, _yield, 0),
+            lulc == crop_lucode, _yield * pixel_area_ha, 0),
         rasters=[lulc_path, yield_path],
         target_path=target_path,
         target_nodata=_NODATA_YIELD)
@@ -939,7 +941,7 @@ def _zero_observed_yield_op(observed_yield_array, observed_yield_nodata):
 
 def _mask_observed_yield_op(
         lulc_array, observed_yield_array, observed_yield_nodata,
-        landcover_nodata, crop_lucode):
+        landcover_nodata, crop_lucode, pixel_area_ha):
     """Mask total observed yield to crop lulc type.
 
     Args:
@@ -948,6 +950,7 @@ def _mask_observed_yield_op(
         observed_yield_nodata (float): yield raster nodata value
         landcover_nodata (float): landcover raster nodata value
         crop_lucode (int): code used to mask in the current crop
+        pixel_area_ha (float): area of lulc raster cells (hectares)
 
     Returns:
         numpy.ndarray with float values of yields masked to crop_lucode
@@ -956,13 +959,13 @@ def _mask_observed_yield_op(
     result = numpy.empty(lulc_array.shape, dtype=numpy.float32)
     if landcover_nodata is not None:
         result[:] = observed_yield_nodata
-        valid_mask = ~pygeoprocessing.array_equals_nodata(lulc_array,
-                                                          landcover_nodata)
+        valid_mask = ~pygeoprocessing.array_equals_nodata(lulc_array, landcover_nodata)
         result[valid_mask] = 0
     else:
         result[:] = 0
     lulc_mask = lulc_array == crop_lucode
-    result[lulc_mask] = observed_yield_array[lulc_mask]
+    result[lulc_mask] = (
+        observed_yield_array[lulc_mask] * pixel_area_ha)
     return result
 
 
@@ -983,7 +986,7 @@ def tabulate_results(
         landcover_raster_path (string): path to landcover raster
         landcover_nodata (float): landcover raster nodata value
         output_dir (string): the file path to the output workspace.
-        file_suffix (string): string to append to any output filenames.
+        file_suffix (string): string to appened to any output filenames.
         target_table_path (string): path to 'result_table.csv' in the output
             workspace
 
@@ -1004,11 +1007,6 @@ def tabulate_results(
         for nutrient_id in _EXPECTED_NUTRIENT_TABLE_HEADERS
         for yield_percentile_id in sorted(yield_percentile_headers) + [
             'yield_observed']]
-
-    # Since pixel values in observed and percentile rasters are Mg/(ha•yr),
-    # raster sums are (Mg•px)/(ha•yr). Before recording sums in
-    # production_lookup dictionary, convert to Mg/yr by multiplying by ha/px.
-
     with open(target_table_path, 'w') as result_table:
         result_table.write(
             'crop,area (ha),' + 'production_observed,' +
@@ -1040,7 +1038,6 @@ def tabulate_results(
                 production_pixel_count += numpy.count_nonzero(
                     valid_mask & (yield_block > 0))
                 yield_sum += numpy.sum(yield_block[valid_mask])
-            yield_sum *= pixel_area_ha
             production_area = production_pixel_count * pixel_area_ha
             production_lookup['observed'] = yield_sum
             result_table.write(',%f' % production_area)
@@ -1058,7 +1055,6 @@ def tabulate_results(
                     yield_sum += numpy.sum(
                         yield_block[~pygeoprocessing.array_equals_nodata(
                             yield_block, _NODATA_YIELD)])
-                yield_sum *= pixel_area_ha
                 production_lookup[yield_percentile_id] = yield_sum
                 result_table.write(",%f" % yield_sum)
 
@@ -1084,8 +1080,7 @@ def tabulate_results(
                 (landcover_raster_path, 1)):
             if landcover_nodata is not None:
                 total_area += numpy.count_nonzero(
-                    ~pygeoprocessing.array_equals_nodata(band_values,
-                                                         landcover_nodata))
+                    ~pygeoprocessing.array_equals_nodata(band_values, landcover_nodata))
             else:
                 total_area += band_values.size
         result_table.write(
@@ -1095,8 +1090,8 @@ def tabulate_results(
 
 def aggregate_to_polygons(
         base_aggregate_vector_path, target_aggregate_vector_path,
-        landcover_raster_projection, crop_names, nutrient_df,
-        yield_percentile_headers, pixel_area_ha, output_dir, file_suffix,
+        landcover_raster_projection, crop_names,
+        nutrient_df, yield_percentile_headers, output_dir, file_suffix,
         target_aggregate_table_path):
     """Write table with aggregate results of yield and nutrient values.
 
@@ -1113,9 +1108,8 @@ def aggregate_to_polygons(
         nutrient_df (pandas.DataFrame): a table of nutrient values by crop
         yield_percentile_headers (list): list of strings indicating percentiles
             at which yield was calculated.
-        pixel_area_ha (float): area of lulc raster cells (hectares)
         output_dir (string): the file path to the output workspace.
-        file_suffix (string): string to append to any output filenames.
+        file_suffix (string): string to appened to any output filenames.
         target_aggregate_table_path (string): path to 'aggregate_results.csv'
             in the output workspace
 
@@ -1130,10 +1124,6 @@ def aggregate_to_polygons(
         target_aggregate_vector_path,
         driver_name='ESRI Shapefile')
 
-    # Since pixel values are Mg/(ha•yr), zonal stats sum is (Mg•px)/(ha•yr).
-    # Before writing sum to results tables or when using sum to calculate
-    # nutrient yields, convert to Mg/yr by multiplying by ha/px.
-
     # loop over every crop and query with pgp function
     total_yield_lookup = {}
     total_nutrient_table = collections.defaultdict(
@@ -1163,12 +1153,11 @@ def aggregate_to_polygons(
                         crop_name, yield_percentile_id)]:
                     total_nutrient_table[nutrient_id][
                         yield_percentile_id][id_index] += (
-                            nutrient_factor
-                            * total_yield_lookup[
-                                    '%s_%s' % (crop_name, yield_percentile_id)
-                                ][id_index]['sum']
-                            * pixel_area_ha
-                            * nutrient_df[nutrient_id][crop_name])
+                            nutrient_factor *
+                            total_yield_lookup['%s_%s' % (
+                                crop_name, yield_percentile_id)][
+                                    id_index]['sum'] *
+                            nutrient_df[nutrient_id][crop_name])
 
         # process observed
         observed_yield_path = os.path.join(
@@ -1182,11 +1171,10 @@ def aggregate_to_polygons(
             for id_index in total_yield_lookup[f'{crop_name}_observed']:
                 total_nutrient_table[
                     nutrient_id]['observed'][id_index] += (
-                        nutrient_factor
-                        * total_yield_lookup[
-                            f'{crop_name}_observed'][id_index]['sum']
-                        * pixel_area_ha
-                        * nutrient_df[nutrient_id][crop_name])
+                        nutrient_factor *
+                        total_yield_lookup[
+                            f'{crop_name}_observed'][id_index]['sum'] *
+                        nutrient_df[nutrient_id][crop_name])
 
     # report everything to a table
     with open(target_aggregate_table_path, 'w') as aggregate_table:
@@ -1205,8 +1193,7 @@ def aggregate_to_polygons(
         for id_index in list(total_yield_lookup.values())[0]:
             aggregate_table.write('%s,' % id_index)
             aggregate_table.write(','.join([
-                str(total_yield_lookup[yield_header][id_index]['sum']
-                    * pixel_area_ha)
+                str(total_yield_lookup[yield_header][id_index]['sum'])
                 for yield_header in sorted(total_yield_lookup)]))
 
             for nutrient_id in _EXPECTED_NUTRIENT_TABLE_HEADERS:
diff --git a/src/natcap/invest/crop_production_regression.py b/src/natcap/invest/crop_production_regression.py
index 3c18c221f..b795ee0d8 100644
--- a/src/natcap/invest/crop_production_regression.py
+++ b/src/natcap/invest/crop_production_regression.py
@@ -1,4 +1,4 @@
-"""InVEST Crop Production Regression Model."""
+"""InVEST Crop Production Percentile Model."""
 import collections
 import logging
 import os
@@ -218,8 +218,8 @@ MODEL_SPEC = {
                         "about": f"{x} {name} production within the polygon",
                         "type": "number",
                         "units": units
-                    } for (nutrient, name, units) in NUTRIENTS
-                    for x in ["modeled", "observed"]
+                    } for nutrient, name, units in NUTRIENTS
+                      for x in ["modeled", "observed"]
                 }
             }
         },
@@ -251,8 +251,8 @@ MODEL_SPEC = {
                         "about": f"{x} {name} production from the crop",
                         "type": "number",
                         "units": units
-                    } for (nutrient, name, units) in NUTRIENTS
-                    for x in ["modeled", "observed"]
+                    } for nutrient, name, units in NUTRIENTS
+                      for x in ["modeled", "observed"]
                 }
             }
         },
@@ -660,7 +660,8 @@ def execute(args):
                    (regression_parameter_raster_path_lookup['c_n'], 1),
                    (args['landcover_raster_path'], 1),
                    (crop_to_fertilization_rate_df['nitrogen_rate'][crop_name],
-                    'raw'), (crop_lucode, 'raw')],
+                    'raw'),
+                   (crop_lucode, 'raw'), (pixel_area_ha, 'raw')],
                   _x_yield_op,
                   nitrogen_yield_raster_path, gdal.GDT_Float32, _NODATA_YIELD),
             target_path_list=[nitrogen_yield_raster_path],
@@ -678,7 +679,8 @@ def execute(args):
                    (regression_parameter_raster_path_lookup['c_p2o5'], 1),
                    (args['landcover_raster_path'], 1),
                    (crop_to_fertilization_rate_df['phosphorus_rate'][crop_name],
-                    'raw'), (crop_lucode, 'raw')],
+                    'raw'),
+                   (crop_lucode, 'raw'), (pixel_area_ha, 'raw')],
                   _x_yield_op,
                   phosphorus_yield_raster_path, gdal.GDT_Float32, _NODATA_YIELD),
             target_path_list=[phosphorus_yield_raster_path],
@@ -696,7 +698,8 @@ def execute(args):
                    (regression_parameter_raster_path_lookup['c_k2o'], 1),
                    (args['landcover_raster_path'], 1),
                    (crop_to_fertilization_rate_df['potassium_rate'][crop_name],
-                    'raw'), (crop_lucode, 'raw')],
+                    'raw'),
+                   (crop_lucode, 'raw'), (pixel_area_ha, 'raw')],
                   _x_yield_op,
                   potassium_yield_raster_path, gdal.GDT_Float32, _NODATA_YIELD),
             target_path_list=[potassium_yield_raster_path],
@@ -791,7 +794,7 @@ def execute(args):
             args=([(args['landcover_raster_path'], 1),
                    (interpolated_observed_yield_raster_path, 1),
                    (observed_yield_nodata, 'raw'), (landcover_nodata, 'raw'),
-                   (crop_lucode, 'raw')],
+                   (crop_lucode, 'raw'), (pixel_area_ha, 'raw')],
                   _mask_observed_yield_op, observed_production_raster_path,
                   gdal.GDT_Float32, observed_yield_nodata),
             target_path_list=[observed_production_raster_path],
@@ -831,10 +834,10 @@ def execute(args):
             func=aggregate_regression_results_to_polygons,
             args=(args['aggregate_polygon_path'],
                   target_aggregate_vector_path,
-                  aggregate_results_table_path,
                   landcover_raster_info['projection_wkt'],
-                  crop_names, nutrient_df, pixel_area_ha,
-                  output_dir, file_suffix),
+                  crop_names, nutrient_df,
+                  output_dir, file_suffix,
+                  aggregate_results_table_path),
             target_path_list=[target_aggregate_vector_path,
                               aggregate_results_table_path],
             dependent_task_list=dependent_task_list,
@@ -845,12 +848,12 @@ def execute(args):
 
 
 def _x_yield_op(
-        y_max, b_x, c_x, lulc_array, fert_rate, crop_lucode):
+        y_max, b_x, c_x, lulc_array, fert_rate, crop_lucode, pixel_area_ha):
     """Calc generalized yield op, Ymax*(1-b_NP*exp(-cN * N_GC)).
 
     The regression model has identical mathematical equations for
-    the nitrogen, phosphorus, and potassium. The only difference is
-    the scalar in the equation (fertilization rate).
+    the nitrogen, phosphorus, and potassium.  The only difference is
+    the scalars in the equation (fertilization rate and pixel area).
     """
     result = numpy.empty(b_x.shape, dtype=numpy.float32)
     result[:] = _NODATA_YIELD
@@ -859,7 +862,7 @@ def _x_yield_op(
         ~pygeoprocessing.array_equals_nodata(b_x, _NODATA_YIELD) &
         ~pygeoprocessing.array_equals_nodata(c_x, _NODATA_YIELD) &
         (lulc_array == crop_lucode))
-    result[valid_mask] = y_max[valid_mask] * (
+    result[valid_mask] = pixel_area_ha * y_max[valid_mask] * (
         1 - b_x[valid_mask] * numpy.exp(
             -c_x[valid_mask] * fert_rate))
 
@@ -894,7 +897,7 @@ def _zero_observed_yield_op(observed_yield_array, observed_yield_nodata):
 
 def _mask_observed_yield_op(
         lulc_array, observed_yield_array, observed_yield_nodata,
-        landcover_nodata, crop_lucode):
+        landcover_nodata, crop_lucode, pixel_area_ha):
     """Mask total observed yield to crop lulc type.
 
     Args:
@@ -903,6 +906,7 @@ def _mask_observed_yield_op(
         observed_yield_nodata (float): yield raster nodata value
         landcover_nodata (float): landcover raster nodata value
         crop_lucode (int): code used to mask in the current crop
+        pixel_area_ha (float): area of lulc raster cells (hectares)
 
     Returns:
         numpy.ndarray with float values of yields masked to crop_lucode
@@ -911,13 +915,13 @@ def _mask_observed_yield_op(
     result = numpy.empty(lulc_array.shape, dtype=numpy.float32)
     if landcover_nodata is not None:
         result[:] = observed_yield_nodata
-        valid_mask = ~pygeoprocessing.array_equals_nodata(
-            lulc_array, landcover_nodata)
+        valid_mask = ~pygeoprocessing.array_equals_nodata(lulc_array, landcover_nodata)
         result[valid_mask] = 0
     else:
         result[:] = 0
     lulc_mask = lulc_array == crop_lucode
-    result[lulc_mask] = observed_yield_array[lulc_mask]
+    result[lulc_mask] = (
+        observed_yield_array[lulc_mask] * pixel_area_ha)
     return result
 
 
@@ -948,11 +952,6 @@ def tabulate_regression_results(
         nutrient_id + '_' + mode
         for nutrient_id in _EXPECTED_NUTRIENT_TABLE_HEADERS
         for mode in ['modeled', 'observed']]
-
-    # Since pixel values in observed and percentile rasters are Mg/(ha•yr),
-    # raster sums are (Mg•px)/(ha•yr). Before recording sums in
-    # production_lookup dictionary, convert to Mg/yr by multiplying by ha/px.
-
     with open(target_table_path, 'w') as result_table:
         result_table.write(
             'crop,area (ha),' + 'production_observed,production_modeled,' +
@@ -983,7 +982,6 @@ def tabulate_regression_results(
                 production_pixel_count += numpy.count_nonzero(
                     valid_mask & (yield_block > 0.0))
                 yield_sum += numpy.sum(yield_block[valid_mask])
-            yield_sum *= pixel_area_ha
             production_area = production_pixel_count * pixel_area_ha
             production_lookup['observed'] = yield_sum
             result_table.write(',%f' % production_area)
@@ -999,7 +997,6 @@ def tabulate_regression_results(
                     # _NODATA_YIELD will always have a value (defined above)
                     yield_block[~pygeoprocessing.array_equals_nodata(
                         yield_block, _NODATA_YIELD)])
-            yield_sum *= pixel_area_ha
             production_lookup['modeled'] = yield_sum
             result_table.write(",%f" % yield_sum)
 
@@ -1034,8 +1031,9 @@ def tabulate_regression_results(
 
 def aggregate_regression_results_to_polygons(
         base_aggregate_vector_path, target_aggregate_vector_path,
-        aggregate_results_table_path, landcover_raster_projection,
-        crop_names, nutrient_df, pixel_area_ha, output_dir, file_suffix):
+        landcover_raster_projection, crop_names,
+        nutrient_df, output_dir, file_suffix,
+        target_aggregate_table_path):
     """Write table with aggregate results of yield and nutrient values.
 
     Use zonal statistics to summarize total observed and interpolated
@@ -1044,16 +1042,15 @@ def aggregate_regression_results_to_polygons(
 
     Args:
         base_aggregate_vector_path (string): path to polygon vector
-        target_aggregate_vector_path (string): path to re-projected copy of
-            polygon vector
-        aggregate_results_table_path (string): path to CSV file where aggregate
-            results will be reported.
+        target_aggregate_vector_path (string):
+            path to re-projected copy of polygon vector
         landcover_raster_projection (string): a WKT projection string
         crop_names (list): list of crop names
         nutrient_df (pandas.DataFrame): a table of nutrient values by crop
-        pixel_area_ha (float): area of lulc raster cells (hectares)
         output_dir (string): the file path to the output workspace.
         file_suffix (string): string to append to any output filenames.
+        target_aggregate_table_path (string): path to 'aggregate_results.csv'
+            in the output workspace
 
     Returns:
         None
@@ -1065,10 +1062,6 @@ def aggregate_regression_results_to_polygons(
         target_aggregate_vector_path,
         driver_name='ESRI Shapefile')
 
-    # Since pixel values are Mg/(ha•yr), zonal stats sum is (Mg•px)/(ha•yr).
-    # Before writing sum to results tables or when using sum to calculate
-    # nutrient yields, convert to Mg/yr by multiplying by ha/px.
-
     # loop over every crop and query with pgp function
     total_yield_lookup = {}
     total_nutrient_table = collections.defaultdict(
@@ -1092,11 +1085,10 @@ def aggregate_regression_results_to_polygons(
             for fid_index in total_yield_lookup['%s_modeled' % crop_name]:
                 total_nutrient_table[nutrient_id][
                     'modeled'][fid_index] += (
-                        nutrient_factor
-                        * total_yield_lookup['%s_modeled' % crop_name][
-                            fid_index]['sum']
-                        * pixel_area_ha
-                        * nutrient_df[nutrient_id][crop_name])
+                        nutrient_factor *
+                        total_yield_lookup['%s_modeled' % crop_name][
+                            fid_index]['sum'] *
+                        nutrient_df[nutrient_id][crop_name])
 
         # process observed
         observed_yield_path = os.path.join(
@@ -1111,14 +1103,15 @@ def aggregate_regression_results_to_polygons(
                     '%s_observed' % crop_name]:
                 total_nutrient_table[
                     nutrient_id]['observed'][fid_index] += (
-                        nutrient_factor  # percent crop used * 1000 [100g per Mg]
-                        * total_yield_lookup[
-                            '%s_observed' % crop_name][fid_index]['sum']
-                        * pixel_area_ha
-                        * nutrient_df[nutrient_id][crop_name])  # nutrient unit per 100g crop
+                        nutrient_factor * # percent crop used * 1000 [100g per Mg]
+                        total_yield_lookup[
+                            '%s_observed' % crop_name][fid_index]['sum'] *
+                        nutrient_df[nutrient_id][crop_name])  # nutrient unit per 100g crop
 
     # report everything to a table
-    with open(aggregate_results_table_path, 'w') as aggregate_table:
+    aggregate_table_path = os.path.join(
+        output_dir, _AGGREGATE_TABLE_FILE_PATTERN % file_suffix)
+    with open(aggregate_table_path, 'w') as aggregate_table:
         # write header
         aggregate_table.write('FID,')
         aggregate_table.write(','.join(sorted(total_yield_lookup)) + ',')
@@ -1134,8 +1127,7 @@ def aggregate_regression_results_to_polygons(
         for id_index in list(total_yield_lookup.values())[0]:
             aggregate_table.write('%s,' % id_index)
             aggregate_table.write(','.join([
-                str(total_yield_lookup[yield_header][id_index]['sum']
-                    * pixel_area_ha)
+                str(total_yield_lookup[yield_header][id_index]['sum'])
                 for yield_header in sorted(total_yield_lookup)]))
 
             for nutrient_id in _EXPECTED_NUTRIENT_TABLE_HEADERS:
diff --git a/src/natcap/invest/forest_carbon_edge_effect.py b/src/natcap/invest/forest_carbon_edge_effect.py
index d506d8342..9d18eb26a 100644
--- a/src/natcap/invest/forest_carbon_edge_effect.py
+++ b/src/natcap/invest/forest_carbon_edge_effect.py
@@ -7,6 +7,7 @@ import logging
 import os
 import pickle
 import time
+import uuid
 
 import numpy
 import pandas
@@ -186,13 +187,13 @@ MODEL_SPEC = {
     "outputs": {
         "carbon_map.tif": {
             "about": (
-                "A map of carbon stock per hectare, with the amount in forest "
-                "derived from the regression based on distance to forest "
-                "edge, and the amount in non-forest classes according to the "
-                "biophysical table. "),
+                "A map of carbon stock per pixel, with the amount in forest derived from the regression based on "
+                "distance to forest edge, and the amount in non-forest classes according to the biophysical table. "
+                "Note that because the map displays carbon per pixel, coarser resolution maps should have higher "
+                "values for carbon, because the pixel areas are larger."),
             "bands": {1: {
                 "type": "number",
-                "units": u.metric_ton/u.hectare
+                "units": u.metric_ton/u.pixel
             }}
         },
         "aggregated_carbon_stocks.shp": {
@@ -204,7 +205,7 @@ MODEL_SPEC = {
                     "units": u.metric_ton,
                     "about": "Total carbon in the area."
                 },
-                "c_ha_mean": {
+                "c_ha_mean":{
                     "type": "number",
                     "units": u.metric_ton/u.hectare,
                     "about": "Mean carbon density in the area."
@@ -579,24 +580,17 @@ def _aggregate_carbon_map(
     target_aggregate_layer.ResetReading()
     target_aggregate_layer.StartTransaction()
 
-    # Since pixel values are Mg/ha, raster sum is (Mg•px)/ha.
-    # To convert to Mg, multiply by ha/px.
-    raster_info = pygeoprocessing.get_raster_info(carbon_map_path)
-    pixel_area = abs(numpy.prod(raster_info['pixel_size']))
-    ha_per_px = pixel_area / 10000
-
     for poly_feat in target_aggregate_layer:
         poly_fid = poly_feat.GetFID()
         poly_feat.SetField(
-            'c_sum', float(serviceshed_stats[poly_fid]['sum'] * ha_per_px))
+            'c_sum', float(serviceshed_stats[poly_fid]['sum']))
         # calculates mean pixel value per ha in for each feature in AOI
         poly_geom = poly_feat.GetGeometryRef()
         poly_area_ha = poly_geom.GetArea() / 1e4  # converts m^2 to hectare
         poly_geom = None
         poly_feat.SetField(
             'c_ha_mean',
-            float(serviceshed_stats[poly_fid]['sum'] / poly_area_ha
-                  * ha_per_px))
+            float(serviceshed_stats[poly_fid]['sum'] / poly_area_ha))
 
         target_aggregate_layer.SetFeature(poly_feat)
     target_aggregate_layer.CommitTransaction()
@@ -635,6 +629,9 @@ def _calculate_lulc_carbon_map(
         biophysical_table_path, **MODEL_SPEC['args']['biophysical_table_path'])
 
     lucode_to_per_cell_carbon = {}
+    cell_size = pygeoprocessing.get_raster_info(
+        lulc_raster_path)['pixel_size']  # in meters
+    cell_area_ha = abs(cell_size[0]) * abs(cell_size[1]) / 10000
 
     # Build a lookup table
     for lucode, row in biophysical_df.iterrows():
@@ -651,7 +648,8 @@ def _calculate_lulc_carbon_map(
                     "Could not interpret carbon pool value as a number. "
                     f"lucode: {lucode}, pool_type: {carbon_pool_type}, "
                     f"value: {row[carbon_pool_type]}")
-            lucode_to_per_cell_carbon[lucode] = row[carbon_pool_type]
+            lucode_to_per_cell_carbon[lucode] = row[carbon_pool_type] * cell_area_ha
+
 
     # map aboveground carbon from table to lulc that is not forest
     reclass_error_details = {
@@ -875,6 +873,7 @@ def _calculate_tropical_forest_edge_carbon_map(
     cell_xsize, cell_ysize = pygeoprocessing.get_raster_info(
         edge_distance_path)['pixel_size']
     cell_size_km = (abs(cell_xsize) + abs(cell_ysize))/2 / 1000
+    cell_area_ha = (abs(cell_xsize) * abs(cell_ysize)) / 10000
 
     # Loop memory block by memory block, calculating the forest edge carbon
     # for every forest pixel.
@@ -958,19 +957,19 @@ def _calculate_tropical_forest_edge_carbon_map(
         biomass[mask_1] = (
             thetas[mask_1][:, 0] - thetas[mask_1][:, 1] * numpy.exp(
                 -thetas[mask_1][:, 2] * valid_edge_distances_km[mask_1])
-        )
+        ) * cell_area_ha
 
         # logarithmic model
         # biomass_2 = t1 + t2 * numpy.log(edge_dist_km)
         biomass[mask_2] = (
             thetas[mask_2][:, 0] + thetas[mask_2][:, 1] * numpy.log(
-                valid_edge_distances_km[mask_2]))
+                valid_edge_distances_km[mask_2])) * cell_area_ha
 
         # linear regression
         # biomass_3 = t1 + t2 * edge_dist_km
         biomass[mask_3] = (
             thetas[mask_3][:, 0] + thetas[mask_3][:, 1] *
-            valid_edge_distances_km[mask_3])
+            valid_edge_distances_km[mask_3]) * cell_area_ha
 
         # reshape the array so that each set of points is in a separate
         # dimension, here distances are distances to each valid model
diff --git a/src/natcap/invest/habitat_quality.py b/src/natcap/invest/habitat_quality.py
index 08eacd90c..c167bca01 100644
--- a/src/natcap/invest/habitat_quality.py
+++ b/src/natcap/invest/habitat_quality.py
@@ -172,9 +172,9 @@ MODEL_SPEC = {
         },
         "sensitivity_table_path": {
             "type": "csv",
-            "index_col": "lucode",
+            "index_col": "lulc",
             "columns": {
-                "lucode": spec_utils.LULC_TABLE_COLUMN,
+                "lulc": spec_utils.LULC_TABLE_COLUMN,
                 "name": {
                     "type": "freestyle_string",
                     "required": False
@@ -282,9 +282,9 @@ MODEL_SPEC = {
                 "rarity_c.csv": {
                     "about": ("Table of rarity values by LULC code for the "
                               "current landscape."),
-                    "index_col": "lucode",
+                    "index_col": "lulc_code",
                     "columns": {
-                        "lucode": {
+                        "lulc_code": {
                             "type": "number",
                             "units": u.none,
                             "about": "LULC class",
@@ -314,9 +314,9 @@ MODEL_SPEC = {
                 "rarity_f.csv": {
                     "about": ("Table of rarity values by LULC code for the "
                               "future landscape."),
-                    "index_col": "lucode",
+                    "index_col": "lulc_code",
                     "columns": {
-                        "lucode": {
+                        "lulc_code": {
                             "type": "number",
                             "units": u.none,
                             "about": "LULC class",
@@ -986,7 +986,7 @@ def _generate_rarity_csv(rarity_dict, target_csv_path):
     lulc_codes = sorted(rarity_dict)
     with open(target_csv_path, 'w', newline='') as csvfile:
         writer = csv.writer(csvfile, delimiter=',')
-        writer.writerow(['lucode', 'rarity_value'])
+        writer.writerow(['lulc_code', 'rarity_value'])
         for lulc_code in lulc_codes:
             writer.writerow([lulc_code, rarity_dict[lulc_code]])
 
diff --git a/src/natcap/invest/ndr/ndr.py b/src/natcap/invest/ndr/ndr.py
index b552c4ad1..4958a6252 100644
--- a/src/natcap/invest/ndr/ndr.py
+++ b/src/natcap/invest/ndr/ndr.py
@@ -1,5 +1,6 @@
 """InVEST Nutrient Delivery Ratio (NDR) module."""
 import copy
+import itertools
 import logging
 import os
 import pickle
@@ -93,13 +94,16 @@ MODEL_SPEC = {
                     "about": gettext(
                         "The distance after which it is assumed that this "
                         "LULC type retains the nutrient at its maximum "
-                        "capacity.")},
+                        "capacity. If nutrients travel a shorter distance "
+                        "that this, the retention "
+                        "efficiency will be less than the maximum value "
+                        "eff_x, following an exponential decay.")},
                 "proportion_subsurface_n": {
                     "type": "ratio",
                     "required": "calc_n",
                     "about": gettext(
                         "The proportion of the total amount of nitrogen that "
-                        "is dissolved into the subsurface. By default, this "
+                        "are dissolved into the subsurface. By default, this "
                         "value should be set to 0, indicating that all "
                         "nutrients are delivered via surface flow. There is "
                         "no equivalent of this for phosphorus.")}
@@ -206,28 +210,28 @@ MODEL_SPEC = {
             "about": "A pixel level map showing how much phosphorus from each pixel eventually reaches the stream by surface flow.",
             "bands": {1: {
                 "type": "number",
-                "units": u.kilogram/u.hectare
+                "units": u.kilogram/u.pixel
             }}
         },
         "n_surface_export.tif": {
             "about": "A pixel level map showing how much nitrogen from each pixel eventually reaches the stream by surface flow.",
             "bands": {1: {
                 "type": "number",
-                "units": u.kilogram/u.hectare
+                "units": u.kilogram/u.pixel
             }}
         },
         "n_subsurface_export.tif": {
             "about": "A pixel level map showing how much nitrogen from each pixel eventually reaches the stream by subsurface flow.",
             "bands": {1: {
                 "type": "number",
-                "units": u.kilogram/u.hectare
+                "units": u.kilogram/u.pixel
             }}
         },
         "n_total_export.tif": {
             "about": "A pixel level map showing how much nitrogen from each pixel eventually reaches the stream by either flow.",
             "bands": {1: {
                 "type": "number",
-                "units": u.kilogram/u.hectare
+                "units": u.kilogram/u.pixel
             }}
         },
         "intermediate_outputs": {
@@ -597,10 +601,7 @@ def execute(args):
             base_raster_list, aligned_raster_list,
             ['near']*len(base_raster_list), dem_info['pixel_size'],
             'intersection'),
-        kwargs={
-            'base_vector_path_list': [args['watersheds_path']],
-            'raster_align_index': 0  # align to the grid of the DEM
-        },
+        kwargs={'base_vector_path_list': [args['watersheds_path']]},
         target_path_list=aligned_raster_list,
         task_name='align rasters')
 
@@ -1028,15 +1029,12 @@ def execute(args):
 # raster_map equation: Multiply a series of arrays element-wise
 def _mult_op(*array_list): return numpy.prod(numpy.stack(array_list), axis=0)
 
-
 # raster_map equation: Sum a list of arrays element-wise
 def _sum_op(*array_list): return numpy.sum(array_list, axis=0)
 
-
 # raster_map equation: calculate inverse of S factor
 def _inverse_op(base_val): return numpy.where(base_val == 0, 0, 1 / base_val)
 
-
 # raster_map equation: rescale and threshold slope between 0.005 and 1
 def _slope_proportion_and_threshold_op(slope):
     slope_fraction = slope / 100
@@ -1137,13 +1135,8 @@ def _add_fields_to_shapefile(field_pickle_map, target_vector_path):
     for feature in target_layer:
         fid = feature.GetFID()
         for field_name in field_pickle_map:
-            # Since pixel values are kg/(ha•yr), raster sum is (kg•px)/(ha•yr).
-            # To convert to kg/yr, multiply by ha/px.
-            pixel_area = field_summaries[field_name]['pixel_area']
-            ha_per_px = pixel_area / 10000
             feature.SetField(
-                field_name, float(
-                    field_summaries[field_name][fid]['sum']) * ha_per_px)
+                field_name, float(field_summaries[field_name][fid]['sum']))
         # Save back to datasource
         target_layer.SetFeature(feature)
     target_layer = None
@@ -1211,8 +1204,8 @@ def _normalize_raster(base_raster_path_band, target_normalized_raster_path):
 
     """
     value_sum, value_count = pygeoprocessing.raster_reduce(
-        function=lambda sum_count, block: (  # calculate both in one pass
-            sum_count[0] + numpy.sum(block), sum_count[1] + block.size),
+        function=lambda sum_count, block:  # calculate both in one pass
+            (sum_count[0] + numpy.sum(block), sum_count[1] + block.size),
         raster_path_band=base_raster_path_band,
         initializer=(0, 0))
 
@@ -1228,25 +1221,29 @@ def _normalize_raster(base_raster_path_band, target_normalized_raster_path):
 
 
 def _calculate_load(lulc_raster_path, lucode_to_load, target_load_raster):
-    """Calculate load raster by mapping landcover.
+    """Calculate load raster by mapping landcover and multiplying by area.
 
     Args:
         lulc_raster_path (string): path to integer landcover raster.
         lucode_to_load (dict): a mapping of landcover IDs to per-area
             nutrient load.
         target_load_raster (string): path to target raster that will have
-            load values (kg/ha) mapped to pixels based on LULC.
+            total load per pixel.
 
     Returns:
         None.
 
     """
+    cell_area_ha = abs(numpy.prod(pygeoprocessing.get_raster_info(
+        lulc_raster_path)['pixel_size'])) * 0.0001
+
     def _map_load_op(lucode_array):
         """Convert unit load to total load & handle nodata."""
         result = numpy.empty(lucode_array.shape)
         for lucode in numpy.unique(lucode_array):
             try:
-                result[lucode_array == lucode] = (lucode_to_load[lucode])
+                result[lucode_array == lucode] = (
+                    lucode_to_load[lucode] * cell_area_ha)
             except KeyError:
                 raise KeyError(
                     'lucode: %d is present in the landuse raster but '
@@ -1442,12 +1439,6 @@ def _aggregate_and_pickle_total(
         base_raster_path_band, aggregate_vector_path,
         working_dir=os.path.dirname(target_pickle_path))
 
-    # Write pixel area to pickle file so that _add_fields_to_shapefile
-    # can adjust totals as needed.
-    raster_info = pygeoprocessing.get_raster_info(base_raster_path_band[0])
-    pixel_area = abs(numpy.prod(raster_info['pixel_size']))
-    result['pixel_area'] = pixel_area
-
     with open(target_pickle_path, 'wb') as target_pickle_file:
         pickle.dump(result, target_pickle_file)
 
diff --git a/src/natcap/invest/sdr/sdr.py b/src/natcap/invest/sdr/sdr.py
index a6fb96c01..f1540982b 100644
--- a/src/natcap/invest/sdr/sdr.py
+++ b/src/natcap/invest/sdr/sdr.py
@@ -147,20 +147,20 @@ MODEL_SPEC = {
             "about": "The contribution of vegetation to keeping soil from eroding from each pixel. (Eq. (82))",
             "bands": {1: {
                 "type": "number",
-                "units": u.metric_ton/u.hectare
+                "units": u.metric_ton/u.pixel
             }}
         },
         "avoided_export.tif": {
             "about": "The contribution of vegetation to keeping erosion from entering a stream. This combines local/on-pixel sediment retention with trapping of erosion from upslope of the pixel. (Eq. (83))",
             "bands": {1: {
                 "type": "number",
-                "units": u.metric_ton/u.hectare
+                "units": u.metric_ton/u.pixel
             }}
         },
         "rkls.tif": {
             "bands": {1: {
                 "type": "number",
-                "units": u.metric_ton/u.hectare
+                "units": u.metric_ton/u.pixel
             }},
             "about": "Total potential soil loss per pixel in the original land cover from the RKLS equation. Equivalent to the soil loss for bare soil. (Eq. (68), without applying the C or P factors)."
         },
@@ -168,14 +168,14 @@ MODEL_SPEC = {
             "about": "The total amount of sediment deposited on the pixel from upslope sources as a result of trapping. (Eq. (80))",
             "bands": {1: {
                 "type": "number",
-                "units": u.metric_ton/u.hectare
+                "units": u.metric_ton/u.pixel
             }}
         },
         "sed_export.tif": {
             "about": "The total amount of sediment exported from each pixel that reaches the stream. (Eq. (76))",
             "bands": {1: {
                 "type": "number",
-                "units": u.metric_ton/u.hectare
+                "units": u.metric_ton/u.pixel
             }}
         },
         "stream.tif": spec_utils.STREAM,
@@ -185,10 +185,10 @@ MODEL_SPEC = {
             "bands": {1: {"type": "integer"}}
         },
         "usle.tif": {
-            "about": "Total potential soil loss per hectare in the original land cover calculated from the USLE equation. (Eq. (68))",
+            "about": "Total potential soil loss per pixel in the original land cover calculated from the USLE equation. (Eq. (68))",
             "bands": {1: {
                 "type": "number",
-                "units": u.metric_ton/u.hectare
+                "units": u.metric_ton/u.pixel
             }}
         },
         "watershed_results_sdr.shp": {
@@ -960,7 +960,6 @@ def _avoided_export_op(avoided_erosion, sdr, sed_deposition):
     # known as a modified USLE)
     return avoided_erosion * sdr + sed_deposition
 
-
 def add_drainage_op(stream, drainage):
     """raster_map equation: add drainage mask to stream layer.
 
@@ -975,11 +974,9 @@ def add_drainage_op(stream, drainage):
     """
     return numpy.where(drainage == 1, 1, stream)
 
-
 # raster_map equation: calculate USLE
 def usle_op(rkls, cp_factor): return rkls * cp_factor
 
-
 # raster_map equation: calculate the inverse ws factor
 def inverse_ws_op(w_factor, s_factor): return 1 / (w_factor * s_factor)
 
@@ -1188,7 +1185,7 @@ def _calculate_ls_factor(
 def _calculate_rkls(
         ls_factor_path, erosivity_path, erodibility_path, stream_path,
         rkls_path):
-    """Calculate potential soil loss (tons / (ha * year)) using RKLS.
+    """Calculate potential soil loss (tons / (pixel * year)) using RKLS.
 
     (revised universal soil loss equation with no C or P).
 
@@ -1214,6 +1211,10 @@ def _calculate_rkls(
     stream_nodata = pygeoprocessing.get_raster_info(
         stream_path)['nodata'][0]
 
+    cell_size = abs(
+        pygeoprocessing.get_raster_info(ls_factor_path)['pixel_size'][0])
+    cell_area_ha = cell_size**2 / 10000.0  # hectares per pixel
+
     def rkls_function(ls_factor, erosivity, erodibility, stream):
         """Calculate the RKLS equation.
 
@@ -1226,7 +1227,7 @@ def _calculate_rkls(
             stream (numpy.ndarray): stream mask (1 stream, 0 no stream)
 
         Returns:
-            numpy.ndarray of RKLS values in tons / (ha * year))
+            numpy.ndarray of RKLS values in tons / (pixel * year))
         """
         rkls = numpy.empty(ls_factor.shape, dtype=numpy.float32)
         nodata_mask = (
@@ -1242,10 +1243,11 @@ def _calculate_rkls(
         valid_mask = nodata_mask & (stream == 0)
         rkls[:] = _TARGET_NODATA
 
-        rkls[valid_mask] = (           # rkls units are tons / (ha * year)
+        rkls[valid_mask] = (           # rkls units are tons / (pixel * year)
             ls_factor[valid_mask] *    # unitless
             erosivity[valid_mask] *    # MJ * mm / (ha * hr * yr)
-            erodibility[valid_mask])   # t * ha * hr / (MJ * ha * mm)
+            erodibility[valid_mask] *  # t * ha * hr / (MJ * ha * mm)
+            cell_area_ha)              # ha / pixel
         return rkls
 
     # aligning with index 3 that's the stream and the most likely to be
@@ -1384,6 +1386,18 @@ def _calculate_d_up(
         target_path=out_d_up_path)
 
 
+def _calculate_d_up_bare(
+        s_bar_path, flow_accumulation_path, out_d_up_bare_path):
+    """Calculate s_bar * sqrt(flow accumulation * cell area)."""
+    cell_area = abs(
+        pygeoprocessing.get_raster_info(s_bar_path)['pixel_size'][0])**2
+    pygeoprocessing.raster_map(
+        op=lambda s_bar, flow_accum: (
+            numpy.sqrt(flow_accum * cell_area) * s_bar),
+        rasters=[s_bar_path, flow_accumulation_path],
+        target_path=out_d_up_bare_path)
+
+
 def _calculate_ic(d_up_path, d_dn_path, out_ic_factor_path):
     """Calculate log10(d_up/d_dn)."""
     # ic can be positive or negative, so float.min is a reasonable nodata value
@@ -1507,20 +1521,13 @@ def _generate_report(
         field_def.SetPrecision(11)
         target_layer.CreateField(field_def)
 
-    # Since pixel values are t/(ha•yr), raster sum is (t•px)/(ha•yr).
-    # To convert to t/yr, multiply by ha/px.
-    raster_info = pygeoprocessing.get_raster_info(usle_path)
-    pixel_area = abs(numpy.prod(raster_info['pixel_size']))
-    ha_per_px = pixel_area / 10000
-
     target_layer.ResetReading()
     for feature in target_layer:
         feature_id = feature.GetFID()
         for field_name in field_summaries:
             feature.SetField(
                 field_name,
-                float(field_summaries[field_name][feature_id]['sum']
-                      * ha_per_px))
+                float(field_summaries[field_name][feature_id]['sum']))
         target_layer.SetFeature(feature)
     target_vector = None
     target_layer = None
diff --git a/src/natcap/invest/urban_cooling_model.py b/src/natcap/invest/urban_cooling_model.py
index d9b8cd4b1..931cff3aa 100644
--- a/src/natcap/invest/urban_cooling_model.py
+++ b/src/natcap/invest/urban_cooling_model.py
@@ -487,7 +487,7 @@ def execute(args):
               'intersection'),
         kwargs={
             'base_vector_path_list': [args['aoi_vector_path']],
-            'raster_align_index': 0,
+            'raster_align_index': 1,
             'target_projection_wkt': lulc_raster_info['projection_wkt']},
         target_path_list=aligned_raster_path_list,
         task_name='align rasters')
diff --git a/workbench/src/renderer/components/LogTab/index.jsx b/workbench/src/renderer/components/LogTab/index.jsx
index e4fd8ddef..ea4a7417e 100644
--- a/workbench/src/renderer/components/LogTab/index.jsx
+++ b/workbench/src/renderer/components/LogTab/index.jsx
@@ -1,4 +1,4 @@
-import React, { useEffect, useRef, useState } from 'react';
+import React, { useEffect, useRef } from 'react';
 import ReactDom from 'react-dom';
 import PropTypes from 'prop-types';
 
@@ -13,65 +13,15 @@ const { ipcRenderer } = window.Workbench.electron;
 function LogDisplay(props) {
   const ref = useRef();
 
-  const [autoScroll, setAutoScroll] = useState(true);
-  const [prevScrollTop, setPrevScrollTop] = useState(0);
-
-  // A scroll event doesn't tell us whether it was initiated by a user or by code,
-  // so we assume all scroll events are user-initiated unless otherwise specified.
-  const [userInitiatedScroll, setUserInitiatedScroll] = useState(true);
-
-  let scrollHandlerTimer;
-
-  // `scrollOffsetThreshold` is used to determine when user has scrolled to bottom of window.
-  // It includes a buffer to account for imprecision stemming from rounding errors and/or
-  // scroll event throttling. 24px (the height of one line of log output) seems to be
-  // large enough to detect scroll events we want, and small enough to avoid false positives.
-  const scrollOffsetThreshold = 24;
-
   useEffect(() => {
-    if (autoScroll) {
-      // Setting `ref.current.scrollTop` will fire a scroll event, which will
-      // result in a call to `handleScroll`. To avoid unnecessary operations in
-      // `handleScroll`, we flag the next scroll event as _not_ user-initiated.
-      setUserInitiatedScroll(false);
-      ref.current.scrollTop = ref.current.scrollHeight;
-    }
+    ref.current.scrollTop = ref.current.scrollHeight;
   }, [props.logdata]);
 
-  // Check scroll direction or position IFF scroll event was user-initiated.
-  // Always update `prevScrollTop` and reset `userInitiatedScroll` to `true`.
-  const handleScroll = () => {
-    if (scrollHandlerTimer) {
-      clearTimeout(scrollHandlerTimer);
-    }
-    scrollHandlerTimer = setTimeout(() => {
-      const currentScrollTop = ref.current.scrollTop;
-      if (userInitiatedScroll) {
-        if (autoScroll) {
-          // If user has scrolled up, halt auto-scrolling.
-          const scrollingUp = (currentScrollTop < prevScrollTop);
-          if (scrollingUp) {
-            setAutoScroll(false);
-          }
-        } else {
-          // If user has scrolled back to the bottom, resume auto-scrolling.
-          const currentScrollOffset = ref.current.scrollHeight - currentScrollTop;
-          if (Math.abs(ref.current.offsetHeight - currentScrollOffset) <= scrollOffsetThreshold) {
-            setAutoScroll(true);
-          }
-        }
-      }
-      setPrevScrollTop(currentScrollTop);
-      setUserInitiatedScroll(true);
-    }, 10);
-  };
-
   return (
     <Col
       className="text-break"
       id="log-display"
       ref={ref}
-      onScroll={handleScroll}
     >
       {
         props.logdata.map(([line, cls], idx) => (


Test patch:
diff --git a/tests/test_carbon.py b/tests/test_carbon.py
index bee15ff1c..ed0e11109 100644
--- a/tests/test_carbon.py
+++ b/tests/test_carbon.py
@@ -4,19 +4,16 @@ import unittest
 import tempfile
 import shutil
 import os
-import re
 
 from osgeo import gdal
 from osgeo import osr
 import numpy
 import numpy.random
 import numpy.testing
-
 import pygeoprocessing
 
 gdal.UseExceptions()
 
-
 def make_simple_raster(base_raster_path, fill_val, nodata_val):
     """Create a 10x10 raster on designated path with fill value.
 
@@ -71,8 +68,7 @@ def assert_raster_equal_value(base_raster_path, val_to_compare):
 
     array_to_compare = numpy.empty(base_array.shape)
     array_to_compare.fill(val_to_compare)
-    numpy.testing.assert_allclose(base_array, array_to_compare,
-                                  rtol=0, atol=1e-3)
+    numpy.testing.assert_allclose(base_array, array_to_compare, rtol=0, atol=1e-6)
 
 
 def make_pools_csv(pools_csv_path):
@@ -94,28 +90,6 @@ def make_pools_csv(pools_csv_path):
         open_table.write('2,1,5,0,3,"lulc code 3"\n')
 
 
-def assert_aggregate_result_equal(html_report_path, stat_name, val_to_compare):
-    """Assert that the given stat in the HTML report has a specific value.
-
-    Args:
-        html_report_path (str): path to the HTML report generated by the model.
-        stat_name (str): name of the stat to find. Must match the name listed
-            in the HTML.
-        val_to_compare (float): the value to check against.
-
-    Returns:
-        None.
-    """
-    with open(html_report_path) as file:
-        report = file.read()
-        pattern = (r'data-summary-stat="'
-                   + stat_name
-                   + r'">(\-?\d+\.\d{2})</td>')
-        match = re.search(pattern, report)
-        stat_str = match.groups()[0]
-        assert float(stat_str) == val_to_compare
-
-
 class CarbonTests(unittest.TestCase):
     """Tests for the Carbon Model."""
 
@@ -157,46 +131,12 @@ class CarbonTests(unittest.TestCase):
 
         carbon.execute(args)
 
-        # Ensure every pixel has the correct total C value.
-        # Current: 15 + 10 + 60 + 1 = 86 Mg/ha
-        assert_raster_equal_value(
-            os.path.join(args['workspace_dir'], 'tot_c_cur.tif'), 86)
-        # Future: 5 + 3 + 20 + 0 = 28 Mg/ha
-        assert_raster_equal_value(
-            os.path.join(args['workspace_dir'], 'tot_c_fut.tif'), 28)
-        # REDD: 2 + 1 + 5 + 0 = 8 Mg/ha
-        assert_raster_equal_value(
-            os.path.join(args['workspace_dir'], 'tot_c_redd.tif'), 8)
-
-        # Ensure deltas are correct.
-        assert_raster_equal_value(
-            os.path.join(args['workspace_dir'], 'delta_cur_fut.tif'), -58)
-        assert_raster_equal_value(
-            os.path.join(args['workspace_dir'], 'delta_cur_redd.tif'), -78)
-
-        # Ensure NPV calculations are correct.
-        # Valuation constant based on provided args is 59.00136.
-        # Future: 59.00136 * -58 Mg/ha = -3422.079 Mg/ha
+        # Add assertions for npv for future and REDD scenarios.
+        # The npv was calculated based on _calculate_npv in carbon.py.
         assert_raster_equal_value(
-            os.path.join(args['workspace_dir'], 'npv_fut.tif'), -3422.079)
-        # REDD: 59.00136 * -78 Mg/ha = -4602.106 Mg/ha
+            os.path.join(args['workspace_dir'], 'npv_fut.tif'), -0.3422078)
         assert_raster_equal_value(
-            os.path.join(args['workspace_dir'], 'npv_redd.tif'), -4602.106)
-
-        # Ensure aggregate results are correct.
-        report_path = os.path.join(args['workspace_dir'], 'report.html')
-        # Raster size is 100 m^2; therefore, raster total is as follows:
-        # (x Mg / 1 ha) * (1 ha / 10000 m^2) * (100 m^2) = (x / 100) Mg
-        for (stat, expected_value) in [
-                ('Total cur', 0.86),
-                ('Total fut', 0.28),
-                ('Total redd', 0.08),
-                ('Change in C for fut', -0.58),
-                ('Change in C for redd', -0.78),
-                ('Net present value from cur to fut', -34.22),
-                ('Net present value from cur to redd', -46.02),
-                ]:
-            assert_aggregate_result_equal(report_path, stat, expected_value)
+            os.path.join(args['workspace_dir'], 'npv_redd.tif'), -0.4602106)
 
     def test_carbon_zero_rates(self):
         """Carbon: test with 0 discount and rate change."""
@@ -228,16 +168,16 @@ class CarbonTests(unittest.TestCase):
 
         # Add assertions for npv for future and REDD scenarios.
         # carbon change from cur to fut:
-        # -58 Mg/ha * 43 $/Mg = -2494 $/ha
+        # -58 Mg/ha * .0001 ha/pixel * 43 $/Mg = -0.2494 $/pixel
         assert_raster_equal_value(
-            os.path.join(args['workspace_dir'], 'npv_fut.tif'), -2494)
+            os.path.join(args['workspace_dir'], 'npv_fut.tif'), -0.2494)
         # carbon change from cur to redd:
-        # -78 Mg/ha * 43 $/Mg = -3354 $/ha
+        # -78 Mg/ha * .0001 ha/pixel * 43 $/Mg = -0.3354 $/pixel
         assert_raster_equal_value(
-            os.path.join(args['workspace_dir'], 'npv_redd.tif'), -3354)
+            os.path.join(args['workspace_dir'], 'npv_redd.tif'), -0.3354)
 
-    def test_carbon_without_redd(self):
-        """Carbon: regression testing for future scenario only (no REDD)."""
+    def test_carbon_future(self):
+        """Carbon: regression testing future scenario."""
         from natcap.invest import carbon
         args = {
             'workspace_dir': self.workspace_dir,
@@ -261,12 +201,10 @@ class CarbonTests(unittest.TestCase):
         make_pools_csv(args['carbon_pools_path'])
 
         carbon.execute(args)
-
-        # Ensure NPV calculations are correct.
-        # Valuation constant based on provided args is 59.00136.
-        # Future: 59.00136 * -58 Mg/ha = -3422.079 Mg/ha
+        # Add assertions for npv for the future scenario.
+        # The npv was calculated based on _calculate_npv in carbon.py.
         assert_raster_equal_value(
-            os.path.join(args['workspace_dir'], 'npv_fut.tif'), -3422.079)
+            os.path.join(args['workspace_dir'], 'npv_fut.tif'), -0.3422078)
 
     def test_carbon_missing_landcover_values(self):
         """Carbon: testing expected exception on missing LULC codes."""
@@ -323,14 +261,12 @@ class CarbonTests(unittest.TestCase):
 
         carbon.execute(args)
 
-        # Ensure NPV calculations are correct.
-        # Valuation constant based on provided args is 59.00136.
-        # Future: 59.00136 * -58 Mg/ha = -3422.079 Mg/ha
+        # Add assertions for npv for future and REDD scenarios.
+        # The npv was calculated based on _calculate_npv in carbon.py.
         assert_raster_equal_value(
-            os.path.join(args['workspace_dir'], 'npv_fut.tif'), -3422.079)
-        # REDD: 59.00136 * -78 Mg/ha = -4602.106 Mg/ha
+            os.path.join(args['workspace_dir'], 'npv_fut.tif'), -0.3422078)
         assert_raster_equal_value(
-            os.path.join(args['workspace_dir'], 'npv_redd.tif'), -4602.106)
+            os.path.join(args['workspace_dir'], 'npv_redd.tif'), -0.4602106)
 
     def test_generate_carbon_map(self):
         """Test `_generate_carbon_map`"""
@@ -377,7 +313,7 @@ class CarbonTests(unittest.TestCase):
         band = actual_carbon_stock.GetRasterBand(1)
         actual_carbon_stock = band.ReadAsArray()
 
-        expected_carbon_stock = numpy.array([[5000, 5000], [60, 120]],
+        expected_carbon_stock = numpy.array([[0.5, 0.5], [0.006, 0.012]],
                                             dtype=numpy.float32)
 
         numpy.testing.assert_array_equal(actual_carbon_stock,
diff --git a/tests/test_coastal_blue_carbon.py b/tests/test_coastal_blue_carbon.py
index 93df53aaa..e5e2b6d2c 100644
--- a/tests/test_coastal_blue_carbon.py
+++ b/tests/test_coastal_blue_carbon.py
@@ -98,7 +98,7 @@ class TestPreprocessor(unittest.TestCase):
         with open(lulc_attributes_path, 'w') as attributes_table:
             attributes_table.write(textwrap.dedent(
                 """\
-                lulc-class,lucode,is_coastal_blue_carbon_habitat
+                lulc-class,code,is_coastal_blue_carbon_habitat
                 mangrove,1,TRUE
                 parkinglot,2,FALSE
                 """))
@@ -177,7 +177,7 @@ class TestPreprocessor(unittest.TestCase):
         expected_landcover_codes = set(range(0, 24))
         found_landcover_codes = set(pandas.read_csv(
             os.path.join(outputs_dir, 'carbon_biophysical_table_template_150225.csv')
-        )['lucode'].values)
+        )['code'].values)
         self.assertEqual(expected_landcover_codes, found_landcover_codes)
 
     def test_transition_table(self):
@@ -207,7 +207,7 @@ class TestPreprocessor(unittest.TestCase):
         landcover_table_path = os.path.join(self.workspace_dir,
                                             'lulc_table.csv')
         with open(landcover_table_path, 'w') as lulc_csv:
-            lulc_csv.write('lucode,lulc-class,is_coastal_blue_carbon_habitat\n')
+            lulc_csv.write('code,lulc-class,is_coastal_blue_carbon_habitat\n')
             lulc_csv.write('0,mangrove,True\n')
             lulc_csv.write('1,parking lot,False\n')
 
@@ -451,7 +451,7 @@ class TestCBC2(unittest.TestCase):
         wkt = srs.ExportToWkt()
 
         biophysical_table = [
-            ['lucode', 'lulc-class', 'biomass-initial', 'soil-initial',
+            ['code', 'lulc-class', 'biomass-initial', 'soil-initial',
                 'litter-initial', 'biomass-half-life',
                 'biomass-low-impact-disturb', 'biomass-med-impact-disturb',
                 'biomass-high-impact-disturb', 'biomass-yearly-accumulation',
diff --git a/tests/test_crop_production.py b/tests/test_crop_production.py
index 3372b70f1..dc9b340c8 100644
--- a/tests/test_crop_production.py
+++ b/tests/test_crop_production.py
@@ -22,20 +22,6 @@ TEST_DATA_PATH = os.path.join(
     'crop_production_model')
 
 
-def _get_pixels_per_hectare(raster_path):
-    """Calculate number of pixels per hectare for a given raster.
-
-    Args:
-        raster_path (str): full path to the raster.
-
-    Returns:
-        A float representing the number of pixels per hectare.
-    """
-    raster_info = pygeoprocessing.get_raster_info(raster_path)
-    pixel_area = abs(numpy.prod(raster_info['pixel_size']))
-    return 10000 / pixel_area
-
-
 def make_aggregate_vector(path_to_shp):
     """
     Generate shapefile with two overlapping polygons
@@ -190,41 +176,6 @@ class CropProductionTests(unittest.TestCase):
         pandas.testing.assert_frame_equal(
             expected_result_table, result_table, check_dtype=False)
 
-        # Check raster outputs to make sure values are in Mg/ha.
-        # Raster sum is (Mg•px)/(ha•yr).
-        # Result table reports totals in Mg/yr.
-        # To convert from Mg/yr to (Mg•px)/(ha•yr), multiply by px/ha.
-        expected_raster_sums = {}
-        for (index, crop) in [(0, 'barley'), (1, 'soybean'), (2, 'wheat')]:
-            filename = crop + '_observed_production.tif'
-            pixels_per_hectare = _get_pixels_per_hectare(
-                os.path.join(args['workspace_dir'], filename))
-            expected_raster_sums[filename] = (
-                expected_result_table.loc[index]['production_observed']
-                * pixels_per_hectare)
-            for percentile in ['25', '50', '75', '95']:
-                filename = (
-                    crop + '_yield_' + percentile + 'th_production.tif')
-                col_name = 'production_' + percentile + 'th'
-                pixels_per_hectare = _get_pixels_per_hectare(
-                    os.path.join(args['workspace_dir'], filename))
-                expected_raster_sums[filename] = (
-                    expected_result_table.loc[index][col_name]
-                    * pixels_per_hectare)
-
-        for filename in expected_raster_sums:
-            raster_path = os.path.join(args['workspace_dir'], filename)
-            raster_info = pygeoprocessing.get_raster_info(raster_path)
-            nodata = raster_info['nodata'][0]
-            raster_sum = 0.0
-            for _, block in pygeoprocessing.iterblocks((raster_path, 1)):
-                raster_sum += numpy.sum(
-                    block[~pygeoprocessing.array_equals_nodata(
-                            block, nodata)], dtype=numpy.float32)
-            expected_sum = expected_raster_sums[filename]
-            numpy.testing.assert_allclose(raster_sum, expected_sum,
-                                          rtol=0, atol=0.1)
-
     def test_crop_production_percentile_no_nodata(self):
         """Crop Production: test percentile model with undefined nodata raster.
 
@@ -373,6 +324,9 @@ class CropProductionTests(unittest.TestCase):
             'model_data_path': MODEL_DATA_PATH,
             'fertilization_rate_table_path': os.path.join(
                 SAMPLE_DATA_PATH, 'crop_fertilization_rates.csv'),
+            'nitrogen_fertilization_rate': 29.6,
+            'phosphorus_fertilization_rate': 8.4,
+            'potassium_fertilization_rate': 14.2,
             'n_workers': '-1'
         }
 
@@ -401,6 +355,9 @@ class CropProductionTests(unittest.TestCase):
             'model_data_path': MODEL_DATA_PATH,
             'fertilization_rate_table_path': os.path.join(
                 SAMPLE_DATA_PATH, 'crop_fertilization_rates.csv'),
+            'nitrogen_fertilization_rate': 29.6,
+            'phosphorus_fertilization_rate': 8.4,
+            'potassium_fertilization_rate': 14.2,
             'n_workers': '-1'
         }
 
@@ -452,13 +409,15 @@ class CropProductionTests(unittest.TestCase):
             'model_data_path': MODEL_DATA_PATH,
             'fertilization_rate_table_path': os.path.join(
                 SAMPLE_DATA_PATH, 'crop_fertilization_rates.csv'),
+            'nitrogen_fertilization_rate': 29.6,
+            'phosphorus_fertilization_rate': 8.4,
+            'potassium_fertilization_rate': 14.2,
         }
 
         crop_production_regression.execute(args)
 
         expected_agg_result_table = pandas.read_csv(
-            os.path.join(TEST_DATA_PATH,
-                         'expected_regression_aggregate_results.csv'))
+            os.path.join(TEST_DATA_PATH, 'expected_regression_aggregate_results.csv'))
         agg_result_table = pandas.read_csv(
             os.path.join(args['workspace_dir'], 'aggregate_results.csv'))
         pandas.testing.assert_frame_equal(
@@ -476,38 +435,6 @@ class CropProductionTests(unittest.TestCase):
         pandas.testing.assert_frame_equal(
             expected_result_table, result_table, check_dtype=False)
 
-        # Check raster outputs to make sure values are in Mg/ha.
-        # Raster sum is (Mg•px)/(ha•yr).
-        # Result table reports totals in Mg/yr.
-        # To convert from Mg/yr to (Mg•px)/(ha•yr), multiply by px/ha.
-        expected_raster_sums = {}
-        for (index, crop) in [(0, 'barley'), (1, 'soybean'), (2, 'wheat')]:
-            filename = crop + '_observed_production.tif'
-            pixels_per_hectare = _get_pixels_per_hectare(
-                os.path.join(args['workspace_dir'], filename))
-            expected_raster_sums[filename] = (
-                expected_result_table.loc[index]['production_observed']
-                * pixels_per_hectare)
-            filename = crop + '_regression_production.tif'
-            pixels_per_hectare = _get_pixels_per_hectare(
-                os.path.join(args['workspace_dir'], filename))
-            expected_raster_sums[filename] = (
-                expected_result_table.loc[index]['production_modeled']
-                * pixels_per_hectare)
-
-        for filename in expected_raster_sums:
-            raster_path = os.path.join(args['workspace_dir'], filename)
-            raster_info = pygeoprocessing.get_raster_info(raster_path)
-            nodata = raster_info['nodata'][0]
-            raster_sum = 0.0
-            for _, block in pygeoprocessing.iterblocks((raster_path, 1)):
-                raster_sum += numpy.sum(
-                    block[~pygeoprocessing.array_equals_nodata(
-                            block, nodata)], dtype=numpy.float32)
-            expected_sum = expected_raster_sums[filename]
-            numpy.testing.assert_allclose(raster_sum, expected_sum,
-                                          rtol=0, atol=0.001)
-
     def test_crop_production_regression_no_nodata(self):
         """Crop Production: test regression model with undefined nodata raster.
 
@@ -526,6 +453,9 @@ class CropProductionTests(unittest.TestCase):
             'model_data_path': MODEL_DATA_PATH,
             'fertilization_rate_table_path': os.path.join(
                 SAMPLE_DATA_PATH, 'crop_fertilization_rates.csv'),
+            'nitrogen_fertilization_rate': 29.6,
+            'phosphorus_fertilization_rate': 8.4,
+            'potassium_fertilization_rate': 14.2,
         }
 
         # Create a raster based on the test data geotransform, but smaller and
@@ -574,11 +504,12 @@ class CropProductionTests(unittest.TestCase):
         lulc_array = numpy.array([[3, 3, 2], [3, -1, 3]])
         fert_rate = 0.6
         crop_lucode = 3
+        pixel_area_ha = 10
 
         actual_result = _x_yield_op(y_max, b_x, c_x, lulc_array, fert_rate,
-                                    crop_lucode)
-        expected_result = numpy.array([[-1, -1.9393047, -1],
-                                       [2.6776089, -1, 1.51231]])
+                                    crop_lucode, pixel_area_ha)
+        expected_result = numpy.array([[-1, -19.393047, -1],
+                                       [26.776089, -1, 15.1231]])
 
         numpy.testing.assert_allclose(actual_result, expected_result)
 
@@ -613,12 +544,13 @@ class CropProductionTests(unittest.TestCase):
 
         landcover_nodata = -9999
         crop_lucode = 3
+        pixel_area_ha = 10
 
         actual_result = _mask_observed_yield_op(
             lulc_array, observed_yield_array, observed_yield_nodata,
-            landcover_nodata, crop_lucode)
+            landcover_nodata, crop_lucode, pixel_area_ha)
 
-        expected_result = numpy.array([[-1, 0, -1], [8, -9999, 0]])
+        expected_result = numpy.array([[-10, 0, -1], [80, -99990, 0]])
 
         numpy.testing.assert_allclose(actual_result, expected_result)
 
@@ -631,75 +563,75 @@ class CropProductionTests(unittest.TestCase):
             """Creates the expected results DataFrame."""
             return pandas.DataFrame([
                 {'crop': 'corn', 'area (ha)': 20.0,
-                 'production_observed': 80.0, 'production_modeled': 40.0,
-                 'protein_modeled': 15624000.0, 'protein_observed': 31248000.0,
-                 'lipid_modeled': 2976000.0, 'lipid_observed': 5952000.0,
-                 'energy_modeled': 177072000.0, 'energy_observed': 354144000.0,
-                 'ca_modeled': 10044000.0, 'ca_observed': 20088000.0,
-                 'fe_modeled': 5840400.0, 'fe_observed': 11680800.0,
-                 'mg_modeled': 104160000.0, 'mg_observed': 208320000.0,
-                 'ph_modeled': 261888000.0, 'ph_observed': 523776000.0,
-                 'k_modeled': 642444000.0, 'k_observed': 1284888000.0,
-                 'na_modeled': 744000.0, 'na_observed': 1488000.0,
-                 'zn_modeled': 1822800.0, 'zn_observed': 3645600.0,
-                 'cu_modeled': 706800.0, 'cu_observed': 1413600.0,
-                 'fl_modeled': 2976000.0, 'fl_observed': 5952000.0,
-                 'mn_modeled': 1078800.0, 'mn_observed': 2157600.0,
-                 'se_modeled': 37200.0, 'se_observed': 74400.0,
-                 'vita_modeled': 1116000.0, 'vita_observed': 2232000.0,
-                 'betac_modeled': 5952000.0, 'betac_observed': 11904000.0,
-                 'alphac_modeled': 855600.0, 'alphac_observed': 1711200.0,
-                 'vite_modeled': 297600.0, 'vite_observed': 595200.0,
-                 'crypto_modeled': 595200.0, 'crypto_observed': 1190400.0,
-                 'lycopene_modeled': 133920.0, 'lycopene_observed': 267840.0,
-                 'lutein_modeled': 23436000.0, 'lutein_observed': 46872000.0,
-                 'betat_modeled': 186000.0, 'betat_observed': 372000.0,
-                 'gammat_modeled': 781200.0, 'gammat_observed': 1562400.0,
-                 'deltat_modeled': 706800.0, 'deltat_observed': 1413600.0,
-                 'vitc_modeled': 2529600.0, 'vitc_observed': 5059200.0,
-                 'thiamin_modeled': 148800.0, 'thiamin_observed': 297600.0,
-                 'riboflavin_modeled': 669600.0, 'riboflavin_observed': 1339200.0,
-                 'niacin_modeled': 3050400.0, 'niacin_observed': 6100800.0,
-                 'pantothenic_modeled': 334800.0, 'pantothenic_observed': 669600.0,
-                 'vitb6_modeled': 520800.0, 'vitb6_observed': 1041600.0,
-                 'folate_modeled': 143220000.0, 'folate_observed': 286440000.0,
-                 'vitb12_modeled': 744000.0, 'vitb12_observed': 1488000.0,
-                 'vitk_modeled': 15252000.0, 'vitk_observed': 30504000.0},
+                'production_observed': 8.0, 'production_modeled': 4.0,
+                'protein_modeled': 1562400.0, 'protein_observed': 3124800.0,
+                'lipid_modeled': 297600.0, 'lipid_observed': 595200.0,
+                'energy_modeled': 17707200.0, 'energy_observed': 35414400.0,
+                'ca_modeled': 1004400.0, 'ca_observed': 2008800.0,
+                'fe_modeled': 584040.0, 'fe_observed': 1168080.0,
+                'mg_modeled': 10416000.0, 'mg_observed': 20832000.0,
+                'ph_modeled': 26188800.0, 'ph_observed': 52377600.0,
+                'k_modeled': 64244400.0, 'k_observed': 128488800.0,
+                'na_modeled': 74400.0, 'na_observed': 148800.0,
+                'zn_modeled': 182280.0, 'zn_observed': 364560.0,
+                'cu_modeled': 70680.0, 'cu_observed': 141360.0,
+                'fl_modeled': 297600.0, 'fl_observed': 595200.0,
+                'mn_modeled': 107880.0, 'mn_observed': 215760.0,
+                'se_modeled': 3720.0, 'se_observed': 7440.0,
+                'vita_modeled': 111600.0, 'vita_observed': 223200.0,
+                'betac_modeled': 595200.0, 'betac_observed': 1190400.0,
+                'alphac_modeled': 85560.0, 'alphac_observed': 171120.0,
+                'vite_modeled': 29760.0, 'vite_observed': 59520.0,
+                'crypto_modeled': 59520.0, 'crypto_observed': 119040.0,
+                'lycopene_modeled': 13392.0, 'lycopene_observed': 26784.0,
+                'lutein_modeled': 2343600.0, 'lutein_observed': 4687200.0,
+                'betat_modeled': 18600.0, 'betat_observed': 37200.0,
+                'gammat_modeled': 78120.0, 'gammat_observed': 156240.0,
+                'deltat_modeled': 70680.0, 'deltat_observed': 141360.0,
+                'vitc_modeled': 252960.0, 'vitc_observed': 505920.0,
+                'thiamin_modeled': 14880.0, 'thiamin_observed': 29760.0,
+                'riboflavin_modeled': 66960.0, 'riboflavin_observed': 133920.0,
+                'niacin_modeled': 305040.0, 'niacin_observed': 610080.0,
+                'pantothenic_modeled': 33480.0, 'pantothenic_observed': 66960.0,
+                'vitb6_modeled': 52080.0, 'vitb6_observed': 104160.0,
+                'folate_modeled': 14322000.0, 'folate_observed': 28644000.0,
+                'vitb12_modeled': 74400.0, 'vitb12_observed': 148800.0,
+                'vitk_modeled': 1525200.0, 'vitk_observed': 3050400.0},
                 {'crop': 'soybean', 'area (ha)': 40.0,
-                 'production_observed': 120.0, 'production_modeled': 70.0,
-                 'protein_modeled': 21021000.0, 'protein_observed': 36036000.0,
-                 'lipid_modeled': 1274000.0, 'lipid_observed': 2184000.0,
-                 'energy_modeled': 63063000.0, 'energy_observed': 108108000.0,
-                 'ca_modeled': 163709000.0, 'ca_observed': 280644000.0,
-                 'fe_modeled': 10000900.0, 'fe_observed': 17144400.0,
-                 'mg_modeled': 178360000.0, 'mg_observed': 305760000.0,
-                 'ph_modeled': 448448000.0, 'ph_observed': 768768000.0,
-                 'k_modeled': 125489000.0, 'k_observed': 215124000.0,
-                 'na_modeled': 1274000.0, 'na_observed': 2184000.0,
-                 'zn_modeled': 3121300.0, 'zn_observed': 5350800.0,
-                 'cu_modeled': 1019200.0, 'cu_observed': 1747200.0,
-                 'fl_modeled': 1911000.0, 'fl_observed': 3276000.0,
-                 'mn_modeled': 3312400.0, 'mn_observed': 5678400.0,
-                 'se_modeled': 191100.0, 'se_observed': 327600.0,
-                 'vita_modeled': 1911000.0, 'vita_observed': 3276000.0,
-                 'betac_modeled': 10192000.0, 'betac_observed': 17472000.0,
-                 'alphac_modeled': 637000.0, 'alphac_observed': 1092000.0,
-                 'vite_modeled': 509600.0, 'vite_observed': 873600.0,
-                 'crypto_modeled': 382200.0, 'crypto_observed': 655200.0,
-                 'lycopene_modeled': 191100.0, 'lycopene_observed': 327600.0,
-                 'lutein_modeled': 38857000.0, 'lutein_observed': 66612000.0,
-                 'betat_modeled': 318500.0, 'betat_observed': 546000.0,
-                 'gammat_modeled': 1465100.0, 'gammat_observed': 2511600.0,
-                 'deltat_modeled': 764400.0, 'deltat_observed': 1310400.0,
-                 'vitc_modeled': 1911000.0, 'vitc_observed': 3276000.0,
-                 'thiamin_modeled': 267540.0, 'thiamin_observed': 458640.0,
-                 'riboflavin_modeled': 522340.0, 'riboflavin_observed': 895440.0,
-                 'niacin_modeled': 7771400.0, 'niacin_observed': 13322400.0,
-                 'pantothenic_modeled': 586040.0, 'pantothenic_observed': 1004640.0,
-                 'vitb6_modeled': 3439800.0, 'vitb6_observed': 5896800.0,
-                 'folate_modeled': 194285000.0, 'folate_observed': 333060000.0,
-                 'vitb12_modeled': 1911000.0, 'vitb12_observed': 3276000.0,
-                 'vitk_modeled': 26754000.0, 'vitk_observed': 45864000.0}])
+                'production_observed': 12.0, 'production_modeled': 7.0,
+                'protein_modeled': 2102100.0, 'protein_observed': 3603600.0,
+                'lipid_modeled': 127400.0, 'lipid_observed': 218400.0,
+                'energy_modeled': 6306300.0, 'energy_observed': 10810800.0,
+                'ca_modeled': 16370900.0, 'ca_observed': 28064400.0,
+                'fe_modeled': 1000090.0, 'fe_observed': 1714440.0,
+                'mg_modeled': 17836000.0, 'mg_observed': 30576000.0,
+                'ph_modeled': 44844800.0, 'ph_observed': 76876800.0,
+                'k_modeled': 12548900.0, 'k_observed': 21512400.0,
+                'na_modeled': 127400.0, 'na_observed': 218400.0,
+                'zn_modeled': 312130.0, 'zn_observed': 535080.0,
+                'cu_modeled': 101920.0, 'cu_observed': 174720.0,
+                'fl_modeled': 191100.0, 'fl_observed': 327600.0,
+                'mn_modeled': 331240.0, 'mn_observed': 567840.0,
+                'se_modeled': 19110.0, 'se_observed': 32760.0,
+                'vita_modeled': 191100.0, 'vita_observed': 327600.0,
+                'betac_modeled': 1019200.0, 'betac_observed': 1747200.0,
+                'alphac_modeled': 63700.0, 'alphac_observed': 109200.0,
+                'vite_modeled': 50960.0, 'vite_observed': 87360.0,
+                'crypto_modeled': 38220.0, 'crypto_observed': 65520.0,
+                'lycopene_modeled': 19110.0, 'lycopene_observed': 32760.0,
+                'lutein_modeled': 3885700.0, 'lutein_observed': 6661200.0,
+                'betat_modeled': 31850.0, 'betat_observed': 54600.0,
+                'gammat_modeled': 146510.0, 'gammat_observed': 251160.0,
+                'deltat_modeled': 76440.0, 'deltat_observed': 131040.0,
+                'vitc_modeled': 191100.0, 'vitc_observed': 327600.0,
+                'thiamin_modeled': 26754.0, 'thiamin_observed': 45864.0,
+                'riboflavin_modeled': 52234.0, 'riboflavin_observed': 89544.0,
+                'niacin_modeled': 777140.0, 'niacin_observed': 1332240.0,
+                'pantothenic_modeled': 58604.0, 'pantothenic_observed': 100464.0,
+                'vitb6_modeled': 343980.0, 'vitb6_observed': 589680.0,
+                'folate_modeled': 19428500.0, 'folate_observed': 33306000.0,
+                'vitb12_modeled': 191100.0, 'vitb12_observed': 327600.0,
+                'vitk_modeled': 2675400.0, 'vitk_observed': 4586400.0}])
 
         nutrient_df = create_nutrient_df()
 
@@ -743,76 +675,76 @@ class CropProductionTests(unittest.TestCase):
             """Create expected output results"""
             # Define the new values manually
             return pandas.DataFrame([
-                {"FID": 0, "corn_modeled": 10, "corn_observed": 40,
-                 "soybean_modeled": 20, "soybean_observed": 50,
-                 "protein_modeled": 9912000, "protein_observed": 30639000,
-                 "lipid_modeled": 1108000, "lipid_observed": 3886000,
-                 "energy_modeled": 62286000, "energy_observed": 222117000,
-                 "ca_modeled": 49285000, "ca_observed": 126979000,
-                 "fe_modeled": 4317500, "fe_observed": 12983900,
-                 "mg_modeled": 77000000, "mg_observed": 231560000,
-                 "ph_modeled": 193600000, "ph_observed": 582208000,
-                 "k_modeled": 196465000, "k_observed": 732079000,
-                 "na_modeled": 550000, "na_observed": 1654000,
-                 "zn_modeled": 1347500, "zn_observed": 4052300,
-                 "cu_modeled": 467900, "cu_observed": 1434800,
-                 "fl_modeled": 1290000, "fl_observed": 4341000,
-                 "mn_modeled": 1216100, "mn_observed": 3444800,
-                 "se_modeled": 63900, "se_observed": 173700,
-                 "vita_modeled": 825000, "vita_observed": 2481000,
-                 "betac_modeled": 4400000, "betac_observed": 13232000,
-                 "alphac_modeled": 395900, "alphac_observed": 1310600,
-                 "vite_modeled": 220000, "vite_observed": 661600,
-                 "crypto_modeled": 258000, "crypto_observed": 868200,
-                 "lycopene_modeled": 88080, "lycopene_observed": 270420,
-                 "lutein_modeled": 16961000, "lutein_observed": 51191000,
-                 "betat_modeled": 137500, "betat_observed": 413500,
-                 "gammat_modeled": 613900, "gammat_observed": 1827700,
-                 "deltat_modeled": 395100, "deltat_observed": 1252800,
-                 "vitc_modeled": 1178400, "vitc_observed": 3894600,
-                 "thiamin_modeled": 113640, "thiamin_observed": 339900,
-                 "riboflavin_modeled": 316640, "riboflavin_observed": 1042700,
-                 "niacin_modeled": 2983000, "niacin_observed": 8601400,
-                 "pantothenic_modeled": 251140, "pantothenic_observed": 753400,
-                 "vitb6_modeled": 1113000, "vitb6_observed": 2977800,
-                 "folate_modeled": 91315000, "folate_observed": 281995000,
-                 "vitb12_modeled": 732000, "vitb12_observed": 2109000,
-                 "vitk_modeled": 11457000, "vitk_observed": 34362000},
-                {"FID": 1, "corn_modeled": 40, "corn_observed": 80,
-                 "soybean_modeled": 70, "soybean_observed": 120,
-                 "protein_modeled": 36645000, "protein_observed": 67284000,
-                 "lipid_modeled": 4250000, "lipid_observed": 8136000,
-                 "energy_modeled": 240135000, "energy_observed": 462252000,
-                 "ca_modeled": 173753000, "ca_observed": 300732000,
-                 "fe_modeled": 15841300, "fe_observed": 28825200,
-                 "mg_modeled": 282520000, "mg_observed": 514080000,
-                 "ph_modeled": 710336000, "ph_observed": 1292544000,
-                 "k_modeled": 767933000, "k_observed": 1500012000,
-                 "na_modeled": 2018000, "na_observed": 3672000,
-                 "zn_modeled": 4944100, "zn_observed": 8996400,
-                 "cu_modeled": 1726000, "cu_observed": 3160800,
-                 "fl_modeled": 4887000, "fl_observed": 9228000,
-                 "mn_modeled": 4391200, "mn_observed": 7836000,
-                 "se_modeled": 228300, "se_observed": 402000,
-                 "vita_modeled": 3027000, "vita_observed": 5508000,
-                 "betac_modeled": 16144000, "betac_observed": 29376000,
-                 "alphac_modeled": 1492600, "alphac_observed": 2803200,
-                 "vite_modeled": 807200, "vite_observed": 1468800,
-                 "crypto_modeled": 977400, "crypto_observed": 1845600,
-                 "lycopene_modeled": 325020, "lycopene_observed": 595440,
-                 "lutein_modeled": 62293000, "lutein_observed": 113484000,
-                 "betat_modeled": 504500, "betat_observed": 918000,
-                 "gammat_modeled": 2246300, "gammat_observed": 4074000,
-                 "deltat_modeled": 1471200, "deltat_observed": 2724000,
-                 "vitc_modeled": 4440600, "vitc_observed": 8335200,
-                 "thiamin_modeled": 416340, "thiamin_observed": 756240,
-                 "riboflavin_modeled": 1191940, "riboflavin_observed": 2234640,
-                 "niacin_modeled": 10821800, "niacin_observed": 19423200,
-                 "pantothenic_modeled": 920840, "pantothenic_observed": 1674240,
-                 "vitb6_modeled": 3960600, "vitb6_observed": 6938400,
-                 "folate_modeled": 337505000, "folate_observed": 619500000,
-                 "vitb12_modeled": 2655000, "vitb12_observed": 4764000,
-                 "vitk_modeled": 42006000, "vitk_observed": 76368000}
+                {"FID": 0, "corn_modeled": 1, "corn_observed": 4,
+                 "soybean_modeled": 2, "soybean_observed": 5,
+                 "protein_modeled": 991200, "protein_observed": 3063900,
+                 "lipid_modeled": 110800, "lipid_observed": 388600,
+                 "energy_modeled": 6228600, "energy_observed": 22211700,
+                 "ca_modeled": 4928500, "ca_observed": 12697900,
+                 "fe_modeled": 431750, "fe_observed": 1298390,
+                 "mg_modeled": 7700000, "mg_observed": 23156000,
+                 "ph_modeled": 19360000, "ph_observed": 58220800,
+                 "k_modeled": 19646500, "k_observed": 73207900,
+                 "na_modeled": 55000, "na_observed": 165400,
+                 "zn_modeled": 134750, "zn_observed": 405230,
+                 "cu_modeled": 46790, "cu_observed": 143480,
+                 "fl_modeled": 129000, "fl_observed": 434100,
+                 "mn_modeled": 121610, "mn_observed": 344480,
+                 "se_modeled": 6390, "se_observed": 17370,
+                 "vita_modeled": 82500, "vita_observed": 248100,
+                 "betac_modeled": 440000, "betac_observed": 1323200,
+                 "alphac_modeled": 39590, "alphac_observed": 131060,
+                 "vite_modeled": 22000, "vite_observed": 66160,
+                 "crypto_modeled": 25800, "crypto_observed": 86820,
+                 "lycopene_modeled": 8808, "lycopene_observed": 27042,
+                 "lutein_modeled": 1696100, "lutein_observed": 5119100,
+                 "betat_modeled": 13750, "betat_observed": 41350,
+                 "gammat_modeled": 61390, "gammat_observed": 182770,
+                 "deltat_modeled": 39510, "deltat_observed": 125280,
+                 "vitc_modeled": 117840, "vitc_observed": 389460,
+                 "thiamin_modeled": 11364, "thiamin_observed": 33990,
+                 "riboflavin_modeled": 31664, "riboflavin_observed": 104270,
+                 "niacin_modeled": 298300, "niacin_observed": 860140,
+                 "pantothenic_modeled": 25114, "pantothenic_observed": 75340,
+                 "vitb6_modeled": 111300, "vitb6_observed": 297780,
+                 "folate_modeled": 9131500, "folate_observed": 28199500,
+                 "vitb12_modeled": 73200, "vitb12_observed": 210900,
+                 "vitk_modeled": 1145700, "vitk_observed": 3436200},
+                {"FID": 1, "corn_modeled": 4, "corn_observed": 8,
+                 "soybean_modeled": 7, "soybean_observed": 12,
+                 "protein_modeled": 3664500, "protein_observed": 6728400,
+                 "lipid_modeled": 425000, "lipid_observed": 813600,
+                 "energy_modeled": 24013500, "energy_observed": 46225200,
+                 "ca_modeled": 17375300, "ca_observed": 30073200,
+                 "fe_modeled": 1584130, "fe_observed": 2882520,
+                 "mg_modeled": 28252000, "mg_observed": 51408000,
+                 "ph_modeled": 71033600, "ph_observed": 129254400,
+                 "k_modeled": 76793300, "k_observed": 150001200,
+                 "na_modeled": 201800, "na_observed": 367200,
+                 "zn_modeled": 494410, "zn_observed": 899640,
+                 "cu_modeled": 172600, "cu_observed": 316080,
+                 "fl_modeled": 488700, "fl_observed": 922800,
+                 "mn_modeled": 439120, "mn_observed": 783600,
+                 "se_modeled": 22830, "se_observed": 40200,
+                 "vita_modeled": 302700, "vita_observed": 550800,
+                 "betac_modeled": 1614400, "betac_observed": 2937600,
+                 "alphac_modeled": 149260, "alphac_observed": 280320,
+                 "vite_modeled": 80720, "vite_observed": 146880,
+                 "crypto_modeled": 97740, "crypto_observed": 184560,
+                 "lycopene_modeled": 32502, "lycopene_observed": 59544,
+                 "lutein_modeled": 6229300, "lutein_observed": 11348400,
+                 "betat_modeled": 50450, "betat_observed": 91800,
+                 "gammat_modeled": 224630, "gammat_observed": 407400,
+                 "deltat_modeled": 147120, "deltat_observed": 272400,
+                 "vitc_modeled": 444060, "vitc_observed": 833520,
+                 "thiamin_modeled": 41634, "thiamin_observed": 75624,
+                 "riboflavin_modeled": 119194, "riboflavin_observed": 223464,
+                 "niacin_modeled": 1082180, "niacin_observed": 1942320,
+                 "pantothenic_modeled": 92084, "pantothenic_observed": 167424,
+                 "vitb6_modeled": 396060, "vitb6_observed": 693840,
+                 "folate_modeled": 33750500, "folate_observed": 61950000,
+                 "vitb12_modeled": 265500, "vitb12_observed": 476400,
+                 "vitk_modeled": 4200600, "vitk_observed": 7636800}
             ], dtype=float)
 
         workspace = self.workspace_dir
@@ -833,24 +765,25 @@ class CropProductionTests(unittest.TestCase):
         output_dir = os.path.join(workspace, "OUTPUT")
         os.makedirs(output_dir, exist_ok=True)
         file_suffix = 'test'
-
-        _AGGREGATE_TABLE_FILE_PATTERN = os.path.join(
-            '.', 'aggregate_results%s.csv')
-
-        aggregate_table_path = os.path.join(
-            output_dir, _AGGREGATE_TABLE_FILE_PATTERN % file_suffix)
-
-        pixel_area_ha = 10
+        target_aggregate_table_path = ''  # unused
 
         _create_crop_rasters(output_dir, crop_names, file_suffix)
 
         aggregate_regression_results_to_polygons(
             base_aggregate_vector_path, target_aggregate_vector_path,
-            aggregate_table_path, landcover_raster_projection, crop_names,
-            nutrient_df, pixel_area_ha, output_dir, file_suffix)
+            landcover_raster_projection, crop_names,
+            nutrient_df, output_dir, file_suffix,
+            target_aggregate_table_path)
+
+        _AGGREGATE_TABLE_FILE_PATTERN = os.path.join(
+            '.','aggregate_results%s.csv')
+
+        aggregate_table_path = os.path.join(
+            output_dir, _AGGREGATE_TABLE_FILE_PATTERN % file_suffix)
 
         actual_aggregate_table = pandas.read_csv(aggregate_table_path,
                                                  dtype=float)
+        print(actual_aggregate_table)
 
         expected_aggregate_table = _create_expected_agg_table()
 
@@ -858,6 +791,7 @@ class CropProductionTests(unittest.TestCase):
             actual_aggregate_table, expected_aggregate_table)
 
 
+
 class CropValidationTests(unittest.TestCase):
     """Tests for the Crop Productions' MODEL_SPEC and validation."""
 
diff --git a/tests/test_forest_carbon_edge.py b/tests/test_forest_carbon_edge.py
index dd8f12f16..cacbd13ae 100644
--- a/tests/test_forest_carbon_edge.py
+++ b/tests/test_forest_carbon_edge.py
@@ -7,8 +7,6 @@ import os
 from osgeo import gdal
 import numpy
 
-import pygeoprocessing
-
 gdal.UseExceptions()
 REGRESSION_DATA = os.path.join(
     os.path.dirname(__file__), '..', 'data', 'invest-test-data',
@@ -55,31 +53,21 @@ class ForestCarbonEdgeTests(unittest.TestCase):
             args['workspace_dir'])
 
         self._assert_vector_results_close(
-            'id',
-            ['c_sum', 'c_ha_mean'],
-            os.path.join(
+            args['workspace_dir'], 'id', ['c_sum', 'c_ha_mean'], os.path.join(
                 args['workspace_dir'], 'aggregated_carbon_stocks.shp'),
             os.path.join(REGRESSION_DATA, 'agg_results_base.shp'))
 
-        # Check raster output to make sure values are in Mg/ha.
-        raster_path = os.path.join(args['workspace_dir'], 'carbon_map.tif')
-        raster_info = pygeoprocessing.get_raster_info(raster_path)
-        nodata = raster_info['nodata'][0]
-        raster_sum = 0.0
-        for _, block in pygeoprocessing.iterblocks((raster_path, 1)):
-            raster_sum += numpy.sum(
-                block[~pygeoprocessing.array_equals_nodata(
-                        block, nodata)], dtype=numpy.float64)
-
-        # expected_sum_per_pixel_values is in Mg, calculated from raster
-        # generated by the model when each pixel value was in Mg/px.
-        # Since pixel values are now Mg/ha, raster sum is (Mg•px)/ha.
-        # To convert expected_sum_per_pixel_values from Mg, multiply by px/ha.
-        expected_sum_per_pixel_values = 21414391.997192383
-        pixel_area = abs(numpy.prod(raster_info['pixel_size']))
-        pixels_per_hectare = 10000 / pixel_area
-        expected_sum = expected_sum_per_pixel_values * pixels_per_hectare
-        numpy.testing.assert_allclose(raster_sum, expected_sum)
+        expected_carbon_raster = gdal.OpenEx(os.path.join(REGRESSION_DATA, 
+                                                          'carbon_map.tif'))
+        expected_carbon_band = expected_carbon_raster.GetRasterBand(1)
+        expected_carbon_array = expected_carbon_band.ReadAsArray()
+        actual_carbon_raster = gdal.OpenEx(os.path.join(REGRESSION_DATA, 
+                                                        'carbon_map.tif'))
+        actual_carbon_band = actual_carbon_raster.GetRasterBand(1)
+        actual_carbon_array = actual_carbon_band.ReadAsArray()
+        self.assertTrue(numpy.allclose(expected_carbon_array, 
+                                       actual_carbon_array))
+
 
     def test_carbon_dup_output(self):
         """Forest Carbon Edge: test for existing output overlap."""
@@ -140,8 +128,7 @@ class ForestCarbonEdgeTests(unittest.TestCase):
                 REGRESSION_DATA, 'file_list_no_edge_effect.txt'),
             args['workspace_dir'])
         self._assert_vector_results_close(
-            'id',
-            ['c_sum', 'c_ha_mean'],
+            args['workspace_dir'], 'id', ['c_sum', 'c_ha_mean'],
             os.path.join(
                 args['workspace_dir'],
                 'aggregated_carbon_stocks_small_no_edge_effect.shp'),
@@ -175,7 +162,7 @@ class ForestCarbonEdgeTests(unittest.TestCase):
         expected_message = 'Could not interpret carbon pool value'
         actual_message = str(cm.exception)
         self.assertTrue(expected_message in actual_message, actual_message)
-
+    
     def test_missing_lulc_value(self):
         """Forest Carbon Edge: test with missing LULC value."""
         from natcap.invest import forest_carbon_edge_effect
@@ -274,11 +261,12 @@ class ForestCarbonEdgeTests(unittest.TestCase):
                 '\n'.join(missing_files))
 
     def _assert_vector_results_close(
-            self, id_fieldname, field_list, result_vector_path,
+            self, workspace_dir, id_fieldname, field_list, result_vector_path,
             expected_vector_path):
         """Test workspace state against expected aggregate results.
 
         Args:
+            workspace_dir (string): path to the completed model workspace
             id_fieldname (string): fieldname of the unique ID.
             field_list (list of string): list of fields to check
                 near-equality.
diff --git a/tests/test_habitat_quality.py b/tests/test_habitat_quality.py
index fd3179797..db5d6b5d5 100644
--- a/tests/test_habitat_quality.py
+++ b/tests/test_habitat_quality.py
@@ -14,7 +14,6 @@ from shapely.geometry import Polygon
 
 gdal.UseExceptions()
 
-
 def make_raster_from_array(
         base_array, base_raster_path, nodata_val=-1, gdal_type=gdal.GDT_Int32):
     """Make a raster from an array on a designated path.
@@ -134,14 +133,14 @@ def make_sensitivity_samp_csv(
     """
     if include_threat:
         with open(csv_path, 'w') as open_table:
-            open_table.write('lucode,name,habitat,threat_1,threat_2\n')
+            open_table.write('LULC,NAME,HABITAT,threat_1,threat_2\n')
             open_table.write('1,"lulc 1",1,1,1\n')
             if not missing_lines:
                 open_table.write('2,"lulc 2",0.5,0.5,1\n')
                 open_table.write('3,"lulc 3",0,0.3,1\n')
     else:
         with open(csv_path, 'w') as open_table:
-            open_table.write('lucode,name,habitat\n')
+            open_table.write('LULC,NAME,HABITAT\n')
             open_table.write('1,"lulc 1",1\n')
             if not missing_lines:
                 open_table.write('2,"lulc 2",0.5\n')
@@ -267,11 +266,11 @@ class HabitatQualityTests(unittest.TestCase):
             with open(csv_path, newline='') as csvfile:
                 reader = csv.DictReader(csvfile, delimiter=',')
                 self.assertEqual(reader.fieldnames,
-                                 ['lucode', 'rarity_value'])
+                                 ['lulc_code', 'rarity_value'])
                 for (exp_lulc, exp_rarity,
                      places_to_round) in expected_csv_values[csv_filename]:
                     row = next(reader)
-                    self.assertEqual(int(row['lucode']), exp_lulc)
+                    self.assertEqual(int(row['lulc_code']), exp_lulc)
                     self.assertAlmostEqual(float(row['rarity_value']),
                                            exp_rarity, places_to_round)
 
@@ -600,7 +599,7 @@ class HabitatQualityTests(unittest.TestCase):
 
         with open(args['sensitivity_table_path'], 'w') as open_table:
             open_table.write(
-                'lucode,name,habitat,%s,%s\n' % tuple(threatnames))
+                'LULC,NAME,HABITAT,%s,%s\n' % tuple(threatnames))
             open_table.write('1,"lulc 1",1,1,1\n')
             open_table.write('2,"lulc 2",0.5,0.5,1\n')
             open_table.write('3,"lulc 3",0,0.3,1\n')
@@ -1047,7 +1046,7 @@ class HabitatQualityTests(unittest.TestCase):
             args['workspace_dir'], 'sensitivity_samp.csv')
 
         with open(args['sensitivity_table_path'], 'w') as open_table:
-            open_table.write('LuCode,Name,habitat,Threat_1,THREAT_2\n')
+            open_table.write('Lulc,Name,habitat,Threat_1,THREAT_2\n')
             open_table.write('1,"lulc 1",1,1,1\n')
             open_table.write('2,"lulc 2",0.5,0.5,1\n')
             open_table.write('3,"lulc 3",0,0.3,1\n')
@@ -2136,7 +2135,7 @@ class HabitatQualityTests(unittest.TestCase):
         self.assertEqual(validate_result, expected)
 
     def test_habitat_quality_missing_lulc_val_in_sens_table(self):
-        """Habitat Quality: test for empty value in lucode column of
+        """Habitat Quality: test for empty value in LULC column of
         sensitivity table. Expects TypeError"""
         from natcap.invest import habitat_quality
 
@@ -2155,9 +2154,9 @@ class HabitatQualityTests(unittest.TestCase):
         args['sensitivity_table_path'] = os.path.join(
             args['workspace_dir'], 'sensitivity_samp.csv')
         with open(args['sensitivity_table_path'], 'w') as open_table:
-            open_table.write('lucode,name,habitat,threat_1,threat_2\n')
+            open_table.write('LULC,NAME,HABITAT,threat_1,threat_2\n')
             open_table.write('1,"lulc 1",1,1,1\n')
-            open_table.write(',"lulc 2",0.5,0.5,1\n')  # missing lucode value
+            open_table.write(',"lulc 2",0.5,0.5,1\n')  # missing LULC value
             open_table.write('3,"lulc 3",0,0.3,1\n')
 
         make_threats_raster(
diff --git a/tests/test_ndr.py b/tests/test_ndr.py
index f70bdb9ca..a719c7dd5 100644
--- a/tests/test_ndr.py
+++ b/tests/test_ndr.py
@@ -1,4 +1,5 @@
 """InVEST NDR model tests."""
+import collections
 import os
 import shutil
 import tempfile
@@ -147,13 +148,13 @@ class NDRTests(unittest.TestCase):
         if not feature:
             raise AssertionError("No features were output.")
         for field, value in [
-                ('p_surface_load', 41.826904),
-                ('p_surface_export', 5.566120),
-                ('n_surface_load', 2977.551270),
-                ('n_surface_export', 274.020844),
-                ('n_subsurface_load', 28.558048),
-                ('n_subsurface_export', 15.578484),
-                ('n_total_export', 289.599314)]:
+                ('p_surface_load', 41.921),
+                ('p_surface_export', 5.59887886),
+                ('n_surface_load', 2978.520),
+                ('n_subsurface_load', 28.614),
+                ('n_surface_export', 289.0498),
+                ('n_subsurface_load', 28.614094),
+                ('n_total_export', 304.66061401)]:
             if not numpy.isclose(feature.GetField(field), value, atol=1e-2):
                 error_results[field] = (
                     'field', feature.GetField(field), value)
@@ -194,12 +195,11 @@ class NDRTests(unittest.TestCase):
         self.assertEqual(len(validation_messages), 1)
 
     def test_base_regression(self):
-        """NDR base regression test on test data.
+        """NDR base regression test on sample data.
 
-        Executes NDR with test data. Checks for accuracy of aggregate
-        values in summary vector, presence of drainage raster in
-        intermediate outputs, and accuracy of raster outputs (as
-        measured by the sum of their non-nodata pixel values).
+        Execute NDR with sample data and checks that the output files are
+        generated and that the aggregate shapefile fields are the same as the
+        regression case.
         """
         from natcap.invest.ndr import ndr
 
@@ -222,18 +222,14 @@ class NDRTests(unittest.TestCase):
         mismatch_list = []
         # these values were generated by manual inspection of regression
         # results
-        expected_watershed_totals = {
-            'p_surface_load': 41.826904,
-            'p_surface_export': 5.870544,
-            'n_surface_load': 2977.551270,
-            'n_surface_export': 274.020844,
-            'n_subsurface_load': 28.558048,
-            'n_subsurface_export': 15.578484,
-            'n_total_export': 289.599314
-        }
-
-        for field in expected_watershed_totals:
-            expected_value = expected_watershed_totals[field]
+        for field, expected_value in [
+                ('p_surface_load', 41.921860),
+                ('p_surface_export', 5.899117),
+                ('n_surface_load', 2978.519775),
+                ('n_surface_export', 289.0498),
+                ('n_subsurface_load', 28.614094),
+                ('n_subsurface_export', 15.61077),
+                ('n_total_export', 304.660614)]:
             val = result_feature.GetField(field)
             if not numpy.isclose(val, expected_value):
                 mismatch_list.append(
@@ -241,7 +237,7 @@ class NDRTests(unittest.TestCase):
                      'actual: %f' % val))
         result_feature = None
         if mismatch_list:
-            raise AssertionError("results not expected: %s" % mismatch_list)
+            raise RuntimeError("results not expected: %s" % mismatch_list)
 
         # We only need to test that the drainage mask exists.  Functionality
         # for that raster is tested in SDR.
@@ -251,28 +247,6 @@ class NDRTests(unittest.TestCase):
                     args['workspace_dir'], 'intermediate_outputs',
                     'what_drains_to_stream.tif')))
 
-        # Check raster outputs to make sure values are in kg/ha/yr.
-        raster_info = pygeoprocessing.get_raster_info(args['dem_path'])
-        pixel_area = abs(numpy.prod(raster_info['pixel_size']))
-        pixels_per_hectare = 10000 / pixel_area
-        for attr_name in ['p_surface_export',
-                          'n_surface_export',
-                          'n_subsurface_export',
-                          'n_total_export']:
-            # Since pixel values are kg/(ha•yr), raster sum is (kg•px)/(ha•yr),
-            # equal to the watershed total (kg/yr) * (pixels_per_hectare px/ha).
-            expected_sum = (expected_watershed_totals[attr_name]
-                            * pixels_per_hectare)
-            raster_name = attr_name + '.tif'
-            raster_path = os.path.join(args['workspace_dir'], raster_name)
-            nodata = pygeoprocessing.get_raster_info(raster_path)['nodata'][0]
-            raster_sum = 0.0
-            for _, block in pygeoprocessing.iterblocks((raster_path, 1)):
-                raster_sum += numpy.sum(
-                    block[~pygeoprocessing.array_equals_nodata(
-                            block, nodata)], dtype=numpy.float64)
-            numpy.testing.assert_allclose(raster_sum, expected_sum, rtol=1e-6)
-
     def test_regression_undefined_nodata(self):
         """NDR test when DEM, LULC and runoff proxy have undefined nodata."""
         from natcap.invest.ndr import ndr
@@ -304,13 +278,13 @@ class NDRTests(unittest.TestCase):
         # these values were generated by manual inspection of regression
         # results
         for field, expected_value in [
-                ('p_surface_load', 41.826904),
-                ('p_surface_export', 5.870544),
-                ('n_surface_load', 2977.551270),
-                ('n_surface_export', 274.020844),
-                ('n_subsurface_load', 28.558048),
-                ('n_subsurface_export', 15.578484),
-                ('n_total_export', 289.599314)]:
+                ('p_surface_load', 41.921860),
+                ('p_surface_export', 5.899117),
+                ('n_surface_load', 2978.519775),
+                ('n_surface_export', 289.0498),
+                ('n_subsurface_load', 28.614094),
+                ('n_subsurface_export', 15.61077),
+                ('n_total_export', 304.660614)]:
             val = result_feature.GetField(field)
             if not numpy.isclose(val, expected_value):
                 mismatch_list.append(
diff --git a/tests/test_sdr.py b/tests/test_sdr.py
index 0b517b395..341076a08 100644
--- a/tests/test_sdr.py
+++ b/tests/test_sdr.py
@@ -128,21 +128,20 @@ class SDRTests(unittest.TestCase):
             'expected a validation error but didn\'t get one')
 
     def test_base_regression(self):
-        """SDR base regression test on test data.
+        """SDR base regression test on sample data.
 
-        Executes SDR with test data. Checks for accuracy of aggregate
-        values in summary vector, presence of drainage raster in
-        intermediate outputs, absence of negative (non-nodata) values
-        in sed_deposition raster, and accuracy of raster outputs (as
-        measured by the sum of their non-nodata pixel values).
+        Execute SDR with sample data and checks that the output files are
+        generated and that the aggregate shapefile fields are the same as the
+        regression case.
         """
         from natcap.invest.sdr import sdr
 
         # use predefined directory so test can clean up files during teardown
         args = SDRTests.generate_base_args(self.workspace_dir)
+        # make args explicit that this is a base run of SWY
 
         sdr.execute(args)
-        expected_watershed_totals = {
+        expected_results = {
             'usle_tot': 2.62457418442,
             'sed_export': 0.09748090804,
             'sed_dep': 1.71672844887,
@@ -152,8 +151,7 @@ class SDRTests(unittest.TestCase):
 
         vector_path = os.path.join(
             args['workspace_dir'], 'watershed_results_sdr.shp')
-        assert_expected_results_in_vector(expected_watershed_totals,
-                                          vector_path)
+        assert_expected_results_in_vector(expected_results, vector_path)
 
         # We only need to test that the drainage mask exists.  Functionality
         # for that raster is tested elsewhere
@@ -177,29 +175,6 @@ class SDRTests(unittest.TestCase):
         self.assertEqual(
             numpy.count_nonzero(sed_dep_array[negative_non_nodata_mask]), 0)
 
-        # Check raster outputs to make sure values are in Mg/ha/yr.
-        raster_info = pygeoprocessing.get_raster_info(args['dem_path'])
-        pixel_area = abs(numpy.prod(raster_info['pixel_size']))
-        pixels_per_hectare = 10000 / pixel_area
-        for (raster_name,
-             attr_name) in [('usle.tif', 'usle_tot'),
-                            ('sed_export.tif', 'sed_export'),
-                            ('sed_deposition.tif', 'sed_dep'),
-                            ('avoided_export.tif', 'avoid_exp'),
-                            ('avoided_erosion.tif', 'avoid_eros')]:
-            # Since pixel values are Mg/(ha•yr), raster sum is (Mg•px)/(ha•yr),
-            # equal to the watershed total (Mg/yr) * (pixels_per_hectare px/ha).
-            expected_sum = (expected_watershed_totals[attr_name]
-                            * pixels_per_hectare)
-            raster_path = os.path.join(args['workspace_dir'], raster_name)
-            nodata = pygeoprocessing.get_raster_info(raster_path)['nodata'][0]
-            raster_sum = 0.0
-            for _, block in pygeoprocessing.iterblocks((raster_path, 1)):
-                raster_sum += numpy.sum(
-                    block[~pygeoprocessing.array_equals_nodata(
-                            block, nodata)], dtype=numpy.float64)
-            numpy.testing.assert_allclose(raster_sum, expected_sum)
-
     def test_regression_with_undefined_nodata(self):
         """SDR base regression test with undefined nodata values.
 
@@ -313,8 +288,8 @@ class SDRTests(unittest.TestCase):
         with self.assertRaises(ValueError) as context:
             sdr.execute(args)
         self.assertIn(
-            'A value in the biophysical table is not a number '
-            'within range 0..1.', str(context.exception))
+            f'A value in the biophysical table is not a number '
+            f'within range 0..1.', str(context.exception))
 
     def test_base_usle_p_nan(self):
         """SDR test expected exception for USLE_P not a number."""
@@ -329,7 +304,7 @@ class SDRTests(unittest.TestCase):
         with self.assertRaises(ValueError) as context:
             sdr.execute(args)
         self.assertIn(
-            'could not be interpreted as ratios', str(context.exception))
+            f'could not be interpreted as ratios', str(context.exception))
 
     def test_lucode_not_a_number(self):
         """SDR test expected exception for invalid data in lucode column."""
diff --git a/tests/test_ucm.py b/tests/test_ucm.py
index 21995d6a8..5fdeaed7e 100644
--- a/tests/test_ucm.py
+++ b/tests/test_ucm.py
@@ -67,11 +67,11 @@ class UCMTests(unittest.TestCase):
         results_feature = results_layer.GetFeature(1)
 
         expected_results = {
-            'avg_cc': 0.221991,
-            'avg_tmp_v': 37.321395,
-            'avg_tmp_an': 2.321395,
-            'avd_eng_cn': 3525589.694141,
-            'avg_wbgt_v': 32.600393,
+            'avg_cc': 0.222150472947109,
+            'avg_tmp_v': 37.325275675470998,
+            'avg_tmp_an': 2.325275675470998,
+            'avd_eng_cn': 3520217.313878,
+            'avg_wbgt_v': 32.60417266705069,
             'avg_ltls_v': 75.000000000000000,
             'avg_hvls_v': 75.000000000000000,
         }
@@ -89,7 +89,7 @@ class UCMTests(unittest.TestCase):
 
         # Assert that the decimal value of the energy savings value is what we
         # expect.
-        expected_energy_sav = 3562738.321307
+        expected_energy_sav = 3564038.678764
 
         energy_sav = 0.0
         n_nonetype = 0
@@ -110,7 +110,7 @@ class UCMTests(unittest.TestCase):
             # Expected energy savings is an accumulated value and may differ
             # past about 4 decimal places.
             numpy.testing.assert_allclose(energy_sav, expected_energy_sav, rtol=1e-4)
-            self.assertEqual(n_nonetype, 134)
+            self.assertEqual(n_nonetype, 119)
         finally:
             buildings_layer = None
             buildings_vector = None
@@ -149,7 +149,7 @@ class UCMTests(unittest.TestCase):
             # and may differ past about 4 decimal places.
             numpy.testing.assert_allclose(energy_sav, expected_energy_sav,
                                           rtol=1e-4)
-            self.assertEqual(n_nonetype, 134)
+            self.assertEqual(n_nonetype, 119)
         finally:
             buildings_layer = None
             buildings_vector = None
@@ -197,12 +197,12 @@ class UCMTests(unittest.TestCase):
         results_feature = results_layer.GetFeature(1)
 
         expected_results = {
-            'avg_cc': 0.422250,
-            'avg_tmp_v': 36.621779,
-            'avg_tmp_an': 1.621779,
-            'avd_eng_cn': 7148968.928616,
-            'avg_wbgt_v': 31.92365,
-            'avg_ltls_v': 29.380548,
+            'avg_cc': 0.428302583240327,
+            'avg_tmp_v': 36.60869797039769,
+            'avg_tmp_an': 1.608697970397692,
+            'avd_eng_cn': 7239992.744486,
+            'avg_wbgt_v': 31.91108630952381,
+            'avg_ltls_v': 28.73463901689708,
             'avg_hvls_v': 75.000000000000000,
         }
         try: