diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml
index 6034709..2fa0a44 100644
--- a/.github/workflows/python-package.yml
+++ b/.github/workflows/python-package.yml
@@ -48,3 +48,4 @@ jobs:
         run: |
           shellcheck ./scripts/run.sh
           shellcheck ./scripts/slurm.sh
+          shellcheck ./scripts/generate_food_map.sh
diff --git a/prepare_layers/build_gaez_hyde.py b/prepare_layers/build_gaez_hyde.py
new file mode 100644
index 0000000..de0e632
--- /dev/null
+++ b/prepare_layers/build_gaez_hyde.py
@@ -0,0 +1,87 @@
+import argparse
+import math
+import os
+import tempfile
+from pathlib import Path
+
+import yirgacheffe as yg
+import yirgacheffe.operators as yo
+
+from make_area_map import make_area_map
+
+DISAGG_CUTOFF = yg.constant(0.95)
+
+def build_gaez_hyde(
+    gaez_path: Path,
+    hyde_path: Path,
+    output_dir_path: Path,
+) -> None:
+    os.makedirs(output_dir_path, exist_ok=True)
+
+    with yg.read_raster(gaez_path) as gaez:
+        with yg.read_raster(hyde_path) as hyde:
+            assert gaez.map_projection == hyde.map_projection
+            projection = gaez.map_projection
+
+            with tempfile.TemporaryDirectory() as tmpdir:
+                area_raster = Path(tmpdir) / "area.tif"
+                make_area_map(math.fabs(projection.ystep), area_raster)
+
+                with yg.read_narrow_raster(area_raster) as area:
+
+                    portional_hyde = (hyde.nan_to_num() * 1000000) / area
+                    portional_gaez = gaez / 100.0
+
+                    # where gaez and hyde disagree (sum greater than disagg cutoff), scale down
+                    uncapped_total = portional_gaez + portional_hyde
+                    # NaNs stop warnings about divide by zero
+                    uncapped_total_with_nan = yo.where(uncapped_total == 0.0, float("nan"), uncapped_total)
+
+                    # calculate ag-perc scalars
+                    total = yo.where(
+                        uncapped_total_with_nan >= DISAGG_CUTOFF,
+                        DISAGG_CUTOFF - (yg.constant(1) / yo.exp(uncapped_total_with_nan * 2)),
+                        uncapped_total_with_nan,
+                    )
+
+                    gaez_ratio = portional_gaez / uncapped_total_with_nan
+                    gaez_values = total * gaez_ratio
+                    gaez_values.to_geotiff(output_dir_path / "crop.tif")
+
+                    hyde_ratio = portional_hyde / uncapped_total_with_nan
+                    hyde_values = total * hyde_ratio
+                    hyde_values.to_geotiff(output_dir_path / "pasture.tif")
+
+def main() -> None:
+    parser = argparse.ArgumentParser(description="Generate a combined GAEZ and Hyde layers")
+    parser.add_argument(
+        '--gaez',
+        type=Path,
+        help='Gaez raster',
+        required=True,
+        dest='gaez_path',
+    )
+    parser.add_argument(
+        '--hyde',
+        type=Path,
+        help='Hyde raster',
+        required=True,
+        dest='hyde_path',
+    )
+    parser.add_argument(
+        '--output',
+        type=Path,
+        help='Path of directory for combined rasters raster',
+        required=True,
+        dest='output_dir_path',
+    )
+    args = parser.parse_args()
+
+    build_gaez_hyde(
+        args.gaez_path,
+        args.hyde_path,
+        args.output_dir_path,
+    )
+
+if __name__ == "__main__":
+    main()
diff --git a/prepare_layers/make_food_current_map.py b/prepare_layers/make_food_current_map.py
new file mode 100644
index 0000000..91ec634
--- /dev/null
+++ b/prepare_layers/make_food_current_map.py
@@ -0,0 +1,303 @@
+import argparse
+import math
+import os
+import resource
+import sys
+import time
+from pathlib import Path
+from multiprocessing import Manager, Process, cpu_count
+from queue import Queue
+from typing import List, NamedTuple, Optional, Tuple
+
+import numpy as np
+import yirgacheffe as yg
+from yirgacheffe.layers import RasterLayer
+
+NULL_CODE = 0
+CROP_CODE = 1401
+PASTURE_CODE = 1402
+# Codes not to touch. We're currently working at Level 1 except for artificial which is level 2
+PRESERVE_CODES = [600, 700, 900, 1000, 1100, 1200, 1300, 1405]
+
+class TileInfo(NamedTuple):
+    """Info about a tile to process"""
+    x_position : int
+    y_position : int
+    width : int
+    height : int
+    crop_diff : float
+    pasture_diff : float
+
+def process_tile(
+    current: yg.layers.RasterLayer,
+    pnv: yg.layers.RasterLayer,
+    tile: TileInfo,
+) -> np.ndarray:
+
+    data = current.read_array(tile.x_position, tile.y_position, tile.width, tile.height)
+
+    diffs = [
+        (tile.crop_diff, CROP_CODE),
+        (tile.pasture_diff, PASTURE_CODE),
+    ]
+    diffs.sort(key=lambda a: a[0])
+
+    # Ordered so removes will come first
+    for diff_value, habitat_code in diffs:
+        if np.isnan(diff_value):
+            continue
+        required_points = math.floor(data.shape[0] * data.shape[1] * math.fabs(diff_value))
+        if required_points == 0:
+            continue
+
+        if diff_value > 0:
+            valid_mask = ~np.isin(data, [CROP_CODE, PASTURE_CODE] + PRESERVE_CODES)
+        else:
+            valid_mask = data == habitat_code
+
+        valid_locations = valid_mask.nonzero()
+        possible_points = len(valid_locations[0])
+        if possible_points == 0:
+            continue
+        required_points = min(required_points, possible_points)
+
+        selected_locations = np.random.choice(
+            len(valid_locations[0]),
+            size=required_points,
+            replace=False
+        )
+        rows = valid_locations[0][selected_locations]
+        cols = valid_locations[1][selected_locations]
+        if diff_value > 0:
+            data[rows, cols] = habitat_code
+        else:
+            for y, x in zip(rows, cols):
+                absolute_x = x + tile.x_position
+                absolute_y = y + tile.y_position
+                lat, lng = current.latlng_for_pixel(absolute_x, absolute_y)
+                pnv_x, pnv_y = pnv.pixel_for_latlng(lat, lng)
+                val = pnv.read_array(pnv_x, pnv_y, 1, 1)[0][0]
+                data[y][x] = val
+
+    return data
+
+def process_tile_concurrently(
+    current_lvl1_path: Path,
+    pnv_path: Path,
+    input_queue: Queue,
+    result_queue: Queue,
+) -> None:
+    with yg.read_raster(current_lvl1_path) as current:
+        with yg.read_raster(pnv_path) as pnv:
+            while True:
+                tile: Optional[TileInfo] = input_queue.get()
+                if tile is None:
+                    break
+                if np.isnan(tile.crop_diff) and np.isnan(tile.pasture_diff):
+                    result_queue.put((tile, None))
+                else:
+                    data = process_tile(current, pnv, tile)
+                    result_queue.put((tile, data.tobytes()))
+
+    result_queue.put(None)
+
+def build_tile_list(
+    current_lvl1_path: Path,
+    crop_adjustment_path: Path,
+    pasture_adjustment_path: Path,
+) -> List[TileInfo]:
+    tiles = []
+    with yg.read_raster(current_lvl1_path) as current:
+        current_dimensions = current.window.xsize, current.window.ysize
+    with yg.read_raster(crop_adjustment_path) as crop_diff:
+        with yg.read_raster(pasture_adjustment_path) as pasture_diff:
+            assert crop_diff.window == pasture_diff.window
+            diff_dimensions = crop_diff.window.xsize, crop_diff.window.ysize
+
+            x_scale = current_dimensions[0] / diff_dimensions[0]
+            y_scale = current_dimensions[1] / diff_dimensions[1]
+
+            x_steps = [round(i * x_scale) for i in range(diff_dimensions[0])]
+            x_steps.append(current_dimensions[0])
+            y_steps = [round(i * y_scale) for i in range(diff_dimensions[1])]
+            y_steps.append(current_dimensions[1])
+
+            for y in range(crop_diff.window.ysize):
+                crop_row = crop_diff.read_array(0, y, crop_diff.window.xsize, 1)
+                pasture_row = pasture_diff.read_array(0, y, pasture_diff.window.xsize, 1)
+                for x in range(crop_diff.window.xsize):
+                    tiles.append(TileInfo(
+                        x_steps[x],
+                        y_steps[y],
+                        (x_steps[x+1] - x_steps[x]),
+                        (y_steps[y+1] - y_steps[y]),
+                        crop_row[0][x],
+                        pasture_row[0][x],
+                    ))
+    return tiles
+
+def assemble_map(
+    current_lvl1_path: Path,
+    output_path: Path,
+    result_queue: Queue,
+    sentinal_count: int,
+) -> None:
+
+    with yg.read_raster(current_lvl1_path) as current:
+        dtype = current.read_array(0, 0, 1, 1).dtype
+        with RasterLayer.empty_raster_layer_like(current, filename=output_path) as output:
+
+            # A cheat as we don't have a neat API for this on yirgacheffe yet
+            band = output._dataset.GetRasterBand(1) # pylint: disable=W0212
+
+            while True:
+                result : Optional[Tuple[TileInfo,Optional[bytearray]]] = result_queue.get()
+                if result is None:
+                    sentinal_count -= 1
+                    if sentinal_count == 0:
+                        break
+                    continue
+
+                tile, rawdata = result
+                if rawdata is None:
+                    data = current.read_array(tile.x_position, tile.y_position, tile.width, tile.height)
+                else:
+                    n = np.frombuffer(rawdata, dtype=dtype)
+                    data = np.reshape(n, (tile.height, tile.width))
+
+                band.WriteArray(data, tile.x_position, tile.y_position)
+
+
+def pipeline_source(
+    current_lvl1_path: Path,
+    crop_adjustment_path: Path,
+    pasture_adjustment_path: Path,
+    source_queue: Queue,
+    sentinal_count: int,
+) -> None:
+    tiles = build_tile_list(
+        current_lvl1_path,
+        crop_adjustment_path,
+        pasture_adjustment_path,
+    )
+    for tile in tiles:
+        source_queue.put(tile)
+    for _ in range(sentinal_count):
+        source_queue.put(None)
+
+def make_food_current_map(
+    current_lvl1_path: Path,
+    pnv_path: Path,
+    crop_adjustment_path: Path,
+    pasture_adjustment_path: Path,
+    output_path: Path,
+    processes_count: int,
+) -> None:
+    # We'll use a lot of processes which will talk back to the main process, so
+    # we need to adjust the ulimit, which is quite low by default
+    _, max_fd_limit = resource.getrlimit(resource.RLIMIT_NOFILE)
+    resource.setrlimit(resource.RLIMIT_NOFILE, (max_fd_limit, max_fd_limit))
+    print(f"Set fd limit to {max_fd_limit}")
+
+    os.makedirs(output_path.parent, exist_ok=True)
+
+    with Manager() as manager:
+        source_queue = manager.Queue()
+        result_queue = manager.Queue()
+
+        workers = [Process(target=process_tile_concurrently, args=(
+            current_lvl1_path,
+            pnv_path,
+            source_queue,
+            result_queue,
+        )) for index in range(processes_count)]
+        for worker_process in workers:
+            worker_process.start()
+
+        source_worker = Process(target=pipeline_source, args=(
+            current_lvl1_path,
+            crop_adjustment_path,
+            pasture_adjustment_path,
+            source_queue,
+            processes_count,
+        ))
+        source_worker.start()
+
+        assemble_map(
+            current_lvl1_path,
+            output_path,
+            result_queue,
+            processes_count,
+        )
+
+        processes = workers
+        processes.append(source_worker)
+        while processes:
+            candidates = [x for x in processes if not x.is_alive()]
+            for candidate in candidates:
+                candidate.join()
+                if candidate.exitcode:
+                    for victim in processes:
+                        victim.kill()
+                    sys.exit(candidate.exitcode)
+                processes.remove(candidate)
+            time.sleep(0.1)
+
+def main() -> None:
+    parser = argparse.ArgumentParser(description="Build the food current map")
+    parser.add_argument(
+        "--current_lvl1",
+        type=Path,
+        required=True,
+        help="Path of lvl1 current map",
+        dest="current_lvl1_path",
+    )
+    parser.add_argument(
+        '--pnv',
+        type=Path,
+        help='Path of PNV map',
+        required=True,
+        dest='pnv_path',
+    )
+    parser.add_argument(
+        "--crop_diff",
+        type=Path,
+        required=True,
+        help="Path of adjustment for crop diff",
+        dest="crop_adjustment_path",
+    )
+    parser.add_argument(
+        "--pasture_diff",
+        type=Path,
+        required=True,
+        help="Path of adjustment for pasture diff",
+        dest="pasture_adjustment_path",
+    )
+    parser.add_argument(
+        '--output',
+        type=Path,
+        help='Path of food current raster',
+        required=True,
+        dest='output_path',
+    )
+    parser.add_argument(
+        "-j",
+        type=int,
+        required=False,
+        default=round(cpu_count() / 1),
+        dest="processes_count",
+        help="Number of concurrent threads to use."
+    )
+    args = parser.parse_args()
+
+    make_food_current_map(
+        args.current_lvl1_path,
+        args.pnv_path,
+        args.crop_adjustment_path,
+        args.pasture_adjustment_path,
+        args.output_path,
+        args.processes_count,
+    )
+
+if __name__ == "__main__":
+    main()
diff --git a/scripts/generate_food_map.sh b/scripts/generate_food_map.sh
new file mode 100644
index 0000000..aa192c4
--- /dev/null
+++ b/scripts/generate_food_map.sh
@@ -0,0 +1,86 @@
+#!/bin/bash
+#
+# This script will generate the updated habitat map for the FOOD paper based on the original Jung habitat map.
+# Note that this is not derived directly from the Jung map, but rather from the simplified version used in LIFE,
+# which has all habitats at level 1 except anthropomorphic ones at level 2. As such this script assumes you have
+# downloaded and generated `current_raw.tif` from the original LIFE pipeline (see run.sh)
+
+set -e
+
+if [ -z "${DATADIR}" ]; then
+    echo "Please specify $DATADIR"
+    exit 1
+fi
+
+if [ -z "${VIRTUAL_ENV}" ]; then
+    echo "Please specify run in a virtualenv"
+    exit 1
+fi
+
+if [ ! -f "${DATADIR}"/habitat/current_raw.tif ]; then
+    echo "LIFE Level 1 current map required"
+    exit 1
+fi
+
+if [ ! -f "${DATADIR}"/habitat/pnv_raw.tif ]; then
+    echo "Jung PNV map required"
+    exit 1
+fi
+
+if [ ! -d "${DATADIR}"/food ]; then
+    mkdir -p "${DATADIR}"/food
+fi
+
+# Get GAEZ data
+if [ ! -f "${DATADIR}"/food/GLCSv11_02_5m.tif ]; then
+    if [ ! -f "${DATADIR}"/food/LR.zip ]; then
+        curl -o "${DATADIR}"/food/LR.zip https://s3.eu-west-1.amazonaws.com/data.gaezdev.aws.fao.org/LR.zip
+    fi
+    unzip -j "${DATADIR}"/food/LR.zip LR/lco/GLCSv11_02_5m.tif -d "${DATADIR}"/food/
+fi
+
+# Get HYDE 3.2 data
+if [ ! -f "${DATADIR}"/food/modified_grazing2017AD.asc ]; then
+    if [ ! -f "${DATADIR}"/food/baseline.zip ]; then
+        curl -o "${DATADIR}"/food/baseline.zip "https://geo.public.data.uu.nl/vault-hyde/HYDE%203.2%5B1710494848%5D/original/baseline.zip"
+    fi
+    if [ ! -f "${DATADIR}"/food/2017AD_lu.zip ]; then
+        unzip -j "${DATADIR}"/food/baseline.zip baseline/zip/2017AD_lu.zip -d "${DATADIR}"/food/
+    fi
+    if [ ! -f "${DATADIR}"/food/grazing2017AD.asc ]; then
+        unzip -j "${DATADIR}"/food/2017AD_lu.zip grazing2017AD.asc -d "${DATADIR}"/food/
+    fi
+    # The pixel scale in the two files have different rounding despite covering the same area
+    # and so this makes them align.
+    sed "s/0.0833333/0.08333333333333333/" "${DATADIR}"/food/grazing2017AD.asc > "${DATADIR}"/food/modified_grazing2017AD.asc
+fi
+
+if [ ! -f "${DATADIR}"/food/modified_grazing2017AD.prj ]; then
+    echo 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]' > "${DATADIR}"/food/modified_grazing2017AD.prj
+fi
+
+# We need rescaled versions of the current data
+if [ ! -d "${DATADIR}"/food/current_layers ]; then
+    python3 ./aoh-calculator/habitat_process.py --habitat "${DATADIR}"/habitat/current_raw.tif \
+                                                --scale 0.08333333333333333 \
+                                                --output "${DATADIR}"/food/current_layers/
+fi
+
+# Combine GAEZ and HYDE data
+python3 ./prepare_layers/build_gaez_hyde.py --gaez "${DATADIR}"/food/GLCSv11_02_5m.tif \
+                                            --hyde "${DATADIR}"/food/modified_grazing2017AD.asc \
+                                            --output "${DATADIR}"/food/
+
+python3 ./utils/raster_diff.py --raster_a "${DATADIR}"/food/crop.tif \
+                               --raster_b "${DATADIR}"/food/current_layers/lcc_1401.tif \
+                               --output "${DATADIR}"/food/crop_diff.tif
+
+python3 ./utils/raster_diff.py --raster_a "${DATADIR}"/food/pasture.tif \
+                              --raster_b "${DATADIR}"/food/current_layers/lcc_1402.tif \
+                              --output "${DATADIR}"/food/pasture_diff.tif
+
+python3 ./prepare_layers/make_food_current_map.py --current_lvl1 "${DATADIR}"/habitat/current_raw.tif \
+                                                  --pnv "${DATADIR}"/habitat/pnv_raw.tif \
+                                                  --crop_diff "${DATADIR}"/food/crop_diff.tif \
+                                                  --pasture_diff "${DATADIR}"/food/pasture_diff.tif \
+                                                  --output "${DATADIR}"/food/current_raw.tif
diff --git a/scripts/run.sh b/scripts/run.sh
index dcc75bb..4df375e 100755
--- a/scripts/run.sh
+++ b/scripts/run.sh
@@ -23,137 +23,194 @@ fi
 declare -a SCENARIOS=("arable" "restore" "restore_all" "urban" "pasture" "restore_agriculture")
 declare -a TAXAS=("AMPHIBIA" "AVES" "MAMMALIA" "REPTILIA")
 export CURVE=0.25
+export PIXEL_SCALE=0.016666666666667
 
-python3 ./prepare_layers/generate_crosswalk.py --output "${DATADIR}"/crosswalk.csv
+if [ ! -f "${DATADIR}"/crosswalk.csv ]; then
+    python3 ./prepare_layers/generate_crosswalk.py --output "${DATADIR}"/crosswalk.csv
+fi
 
 # Get habitat layer and prepare for use
-reclaimer zenodo --zenodo_id 4058819 \
-                 --filename iucn_habitatclassification_composite_lvl2_ver004.zip \
-                 --extract \
-                 --output "${DATADIR}"/habitat/jung_l2_raw.tif
-
-reclaimer zenodo --zenodo_id 4058819 \
-                 --filename lvl2_changemasks_ver004.zip \
-                 --extract \
-                 --output "${DATADIR}"/habitat/
-
-python3 ./prepare_layers/make_current_map.py --jung "${DATADIR}"/habitat/jung_l2_raw.tif \
-                  --update_masks "${DATADIR}"/habitat/lvl2_changemasks_ver004 \
-                  --crosswalk "${DATADIR}"/crosswalk.csv \
-                  --output "${DATADIR}"/habitat/current_raw.tif \
-                  -j 16
-
-python3 ./aoh-calculator/habitat_process.py --habitat "${DATADIR}"/habitat/current_raw.tif \
-                                            --scale 0.016666666666667 \
-                                            --output "${DATADIR}"/habitat_maps/current/
+if [ ! -f "${DATADIR}"/habitat/current_raw.tif ]; then
+    if [ ! -f "${DATADIR}"/habitat/jung_l2_raw.tif ]; then
+        reclaimer zenodo --zenodo_id 4058819 \
+                        --filename iucn_habitatclassification_composite_lvl2_ver004.zip \
+                        --extract \
+                        --output "${DATADIR}"/habitat/jung_l2_raw.tif
+    fi
+
+    if [ ! -d "${DATADIR}"/habitat/lvl2_changemasks_ver004 ]; then
+        reclaimer zenodo --zenodo_id 4058819 \
+                        --filename lvl2_changemasks_ver004.zip \
+                        --extract \
+                        --output "${DATADIR}"/habitat/
+    fi
+
+    python3 ./prepare_layers/make_current_map.py --jung "${DATADIR}"/habitat/jung_l2_raw.tif \
+                    --update_masks "${DATADIR}"/habitat/lvl2_changemasks_ver004 \
+                    --crosswalk "${DATADIR}"/crosswalk.csv \
+                    --output "${DATADIR}"/habitat/current_raw.tif \
+                    -j 16
+fi
 
-# Get PNV layer and prepare for use
-reclaimer zenodo --zenodo_id 4038749 \
-                 --filename pnv_lvl1_004.zip \
-                 --extract \
-                 --output "${DATADIR}"/habitat/pnv_raw.tif
 
-python3 ./aoh-calculator/habitat_process.py --habitat "${DATADIR}"/habitat/pnv_raw.tif \
-                                            --scale 0.016666666666667 \
-                                            --output "${DATADIR}"/habitat_maps/pnv/
+if [ ! -d "${DATADIR}"/habitat_maps/current ]; then
+    python3 ./aoh-calculator/habitat_process.py --habitat "${DATADIR}"/habitat/current_raw.tif \
+                                                --scale "${PIXEL_SCALE}" \
+                                                --output "${DATADIR}"/habitat_maps/current/
+fi
+
+# Get PNV layer and prepare for use
+if [ ! -f "${DATADIR}"/habitat/pnv_raw.tif ]; then
+    reclaimer zenodo --zenodo_id 4038749 \
+                    --filename pnv_lvl1_004.zip \
+                    --extract \
+                    --output "${DATADIR}"/habitat/pnv_raw.tif
+fi
 
+if [ ! -d "${DATADIR}"/habitat_maps/pnv ]; then
+    python3 ./aoh-calculator/habitat_process.py --habitat "${DATADIR}"/habitat/pnv_raw.tif \
+                                                --scale "${PIXEL_SCALE}" \
+                                                --output "${DATADIR}"/habitat_maps/pnv/
+fi
 
 # Generate an area scaling map
-python3 ./prepare_layers/make_area_map.py --scale 0.016666666666667 --output "${DATADIR}"/area-per-pixel.tif
+if [ ! -f "${DATADIR}"/area-per-pixel.tif ]; then
+    python3 ./prepare_layers/make_area_map.py --scale "${PIXEL_SCALE}" --output "${DATADIR}"/area-per-pixel.tif
+fi
 
 # Generate the arable scenario map
-python3 ./prepare_layers/make_arable_map.py --current "${DATADIR}"/habitat/current_raw.tif \
-                                  --output "${DATADIR}"/habitat/arable.tif
+if [ ! -f "${DATADIR}"/habitat/arable.tif ]; then
+    python3 ./prepare_layers/make_arable_map.py --current "${DATADIR}"/habitat/current_raw.tif \
+                                    --output "${DATADIR}"/habitat/arable.tif
+fi
 
-python3 ./aoh-calculator/habitat_process.py --habitat "${DATADIR}"/habitat/arable.tif \
-                                            --scale 0.016666666666667 \
-                                            --output "${DATADIR}"/habitat_maps/arable/
+if [ ! -d "${DATADIR}"/habitat_maps/arable ]; then
+    python3 ./aoh-calculator/habitat_process.py --habitat "${DATADIR}"/habitat/arable.tif \
+                                                --scale "${PIXEL_SCALE}" \
+                                                --output "${DATADIR}"/habitat_maps/arable/
+fi
 
-python3 ./prepare_layers/make_diff_map.py --current "${DATADIR}"/habitat/current_raw.tif \
-                                          --scenario "${DATADIR}"/habitat/arable.tif \
-                                          --area "${DATADIR}"/area-per-pixel.tif \
-                                          --scale 0.016666666666667 \
-                                          --output "${DATADIR}"/habitat/arable_diff_area.tif
+if [ ! -f "${DATADIR}"/habitat/arable_diff_area.tif ]; then
+    python3 ./prepare_layers/make_diff_map.py --current "${DATADIR}"/habitat/current_raw.tif \
+                                            --scenario "${DATADIR}"/habitat/arable.tif \
+                                            --area "${DATADIR}"/area-per-pixel.tif \
+                                            --scale "${PIXEL_SCALE}" \
+                                            --output "${DATADIR}"/habitat/arable_diff_area.tif
+fi
 
 # Generate the pasture scenario map
-python3 ./prepare_layers/make_pasture_map.py --current "${DATADIR}"/habitat/current_raw.tif \
-                                             --output "${DATADIR}"/habitat/pasture.tif
+if [ ! -f "${DATADIR}"/habitat/pasture.tif ]; then
+    python3 ./prepare_layers/make_pasture_map.py --current "${DATADIR}"/habitat/current_raw.tif \
+                                                --output "${DATADIR}"/habitat/pasture.tif
+fi
 
-python3 ./aoh-calculator/habitat_process.py --habitat "${DATADIR}"/habitat/pasture.tif \
-                                          --scale 0.016666666666667 \
-                                          --output "${DATADIR}"/habitat_maps/pasture/
+if [ ! -d "${DATADIR}"/habitat_maps/pasture ]; then
+    python3 ./aoh-calculator/habitat_process.py --habitat "${DATADIR}"/habitat/pasture.tif \
+                                            --scale "${PIXEL_SCALE}" \
+                                            --output "${DATADIR}"/habitat_maps/pasture/
+fi
 
-python3 ./prepare_layers/make_diff_map.py --current "${DATADIR}"/habitat/current_raw.tif \
-                                        --scenario "${DATADIR}"/habitat/pasture.tif \
-                                        --area "${DATADIR}"/area-per-pixel.tif \
-                                        --scale 0.016666666666667 \
-                                        --output "${DATADIR}"/habitat/pasture_diff_area.tif
+if [ ! -f "${DATADIR}"/habitat/pasture_diff_area.tif ]; then
+    python3 ./prepare_layers/make_diff_map.py --current "${DATADIR}"/habitat/current_raw.tif \
+                                            --scenario "${DATADIR}"/habitat/pasture.tif \
+                                            --area "${DATADIR}"/area-per-pixel.tif \
+                                            --scale "${PIXEL_SCALE}" \
+                                            --output "${DATADIR}"/habitat/pasture_diff_area.tif
+fi
 
 # Generate the restore map
-python3 ./prepare_layers/make_restore_map.py --pnv "${DATADIR}"/habitat/pnv_raw.tif \
-                                   --current "${DATADIR}"/habitat/current_raw.tif \
-                                   --crosswalk "${DATADIR}"/crosswalk.csv \
-                                   --output "${DATADIR}"/habitat/restore.tif
+if [ ! -f "${DATADIR}"/habitat/restore.tif ]; then
+    python3 ./prepare_layers/make_restore_map.py --pnv "${DATADIR}"/habitat/pnv_raw.tif \
+                                    --current "${DATADIR}"/habitat/current_raw.tif \
+                                    --crosswalk "${DATADIR}"/crosswalk.csv \
+                                    --output "${DATADIR}"/habitat/restore.tif
+fi
 
-python3 ./aoh-calculator/habitat_process.py --habitat "${DATADIR}"/habitat/restore.tif \
-                                             --scale 0.016666666666667 \
-                                             --output "${DATADIR}"/habitat_maps/restore/
+if [ ! -d "${DATADIR}"/habitat_maps/restore ]; then
+    python3 ./aoh-calculator/habitat_process.py --habitat "${DATADIR}"/habitat/restore.tif \
+                                                --scale "${PIXEL_SCALE}" \
+                                                --output "${DATADIR}"/habitat_maps/restore/
+fi
 
-python3 ./prepare_layers/make_diff_map.py --current "${DATADIR}"/habitat/current_raw.tif \
-                                          --scenario "${DATADIR}"/habitat/restore.tif \
-                                          --area "${DATADIR}"/area-per-pixel.tif \
-                                          --scale 0.016666666666667 \
-                                          --output "${DATADIR}"/habitat/restore_diff_area.tif
+if [ ! -f "${DATADIR}"/habitat/restore_diff_area.tif ]; then
+    python3 ./prepare_layers/make_diff_map.py --current "${DATADIR}"/habitat/current_raw.tif \
+                                            --scenario "${DATADIR}"/habitat/restore.tif \
+                                            --area "${DATADIR}"/area-per-pixel.tif \
+                                            --scale "${PIXEL_SCALE}" \
+                                            --output "${DATADIR}"/habitat/restore_diff_area.tif
+fi
 
 # Generate the restore map
-python3 ./prepare_layers/make_restore_agriculture_map.py --pnv "${DATADIR}"/habitat/pnv_raw.tif \
-                                                        --current "${DATADIR}"/habitat/current_raw.tif \
-                                                        --crosswalk "${DATADIR}"/crosswalk.csv \
-                                                        --output "${DATADIR}"/habitat/restore_agriculture.tif
+if [ ! -f "${DATADIR}"/habitat/restore_agriculture.tif ]; then
+    python3 ./prepare_layers/make_restore_agriculture_map.py --pnv "${DATADIR}"/habitat/pnv_raw.tif \
+                                                            --current "${DATADIR}"/habitat/current_raw.tif \
+                                                            --crosswalk "${DATADIR}"/crosswalk.csv \
+                                                            --output "${DATADIR}"/habitat/restore_agriculture.tif
+fi
 
-python3 ./aoh-calculator/habitat_process.py --habitat "${DATADIR}"/habitat/restore_agriculture.tif \
-                                           --scale 0.016666666666667 \
-                                           --output "${DATADIR}"/habitat_maps/restore_agriculture/
+if [ ! -d "${DATADIR}"/habitat_maps/restore_agriculture ]; then
+    python3 ./aoh-calculator/habitat_process.py --habitat "${DATADIR}"/habitat/restore_agriculture.tif \
+                                            --scale "${PIXEL_SCALE}" \
+                                            --output "${DATADIR}"/habitat_maps/restore_agriculture/
+fi
 
-python3 ./prepare_layers/make_diff_map.py --current "${DATADIR}"/habitat/current_raw.tif \
-                                        --scenario "${DATADIR}"/habitat/restore_agriculture.tif \
-                                        --area "${DATADIR}"/area-per-pixel.tif \
-                                        --scale 0.016666666666667 \
-                                        --output "${DATADIR}"/habitat/restore_agriculture_diff_area.tif
+if [ ! -f "${DATADIR}"/habitat/restore_agriculture_diff_area.tif ]; then
+    python3 ./prepare_layers/make_diff_map.py --current "${DATADIR}"/habitat/current_raw.tif \
+                                            --scenario "${DATADIR}"/habitat/restore_agriculture.tif \
+                                            --area "${DATADIR}"/area-per-pixel.tif \
+                                            --scale "${PIXEL_SCALE}" \
+                                            --output "${DATADIR}"/habitat/restore_agriculture_diff_area.tif
+fi
 
 # Generate the restore all map
-python3 ./prepare_layers/make_restore_all_map.py --pnv "${DATADIR}"/habitat/pnv_raw.tif \
-                                 --current "${DATADIR}"/habitat/current_raw.tif \
-                                 --crosswalk "${DATADIR}"/crosswalk.csv \
-                                 --output "${DATADIR}"/habitat/restore_all.tif
+if [ ! -f "${DATADIR}"/habitat/restore_all.tif ]; then
+    python3 ./prepare_layers/make_restore_all_map.py --pnv "${DATADIR}"/habitat/pnv_raw.tif \
+                                    --current "${DATADIR}"/habitat/current_raw.tif \
+                                    --crosswalk "${DATADIR}"/crosswalk.csv \
+                                    --output "${DATADIR}"/habitat/restore_all.tif
+fi
 
-python3 ./aoh-calculator/habitat_process.py --habitat "${DATADIR}"/habitat/restore_all.tif \
-                                           --scale 0.016666666666667 \
-                                           --output "${DATADIR}"/habitat_maps/restore_all/
+if [ ! -d "${DATADIR}"/habitat_maps/restore_all ]; then
+    python3 ./aoh-calculator/habitat_process.py --habitat "${DATADIR}"/habitat/restore_all.tif \
+                                            --scale "${PIXEL_SCALE}" \
+                                            --output "${DATADIR}"/habitat_maps/restore_all/
+fi
 
-python3 ./prepare_layers/make_diff_map.py --current "${DATADIR}"/habitat/current_raw.tif \
-                                        --scenario "${DATADIR}"/habitat/restore_all.tif \
-                                        --area "${DATADIR}"/area-per-pixel.tif \
-                                        --scale 0.016666666666667 \
-                                        --output "${DATADIR}"/habitat/restore_all_diff_area.tif
+if [ ! -f "${DATADIR}"/habitat/restore_all_diff_area.tif ]; then
+    python3 ./prepare_layers/make_diff_map.py --current "${DATADIR}"/habitat/current_raw.tif \
+                                            --scenario "${DATADIR}"/habitat/restore_all.tif \
+                                            --area "${DATADIR}"/area-per-pixel.tif \
+                                            --scale "${PIXEL_SCALE}" \
+                                            --output "${DATADIR}"/habitat/restore_all_diff_area.tif
+fi
 
 # Generate urban all map
-python3 ./prepare_layers/make_constant_habitat.py --examplar "${DATADIR}"/habitat_maps/arable/lcc_1401.tif \
-                                                   --habitat_code 14.5 \
-                                                   --crosswalk "${DATADIR}"/crosswalk.csv \
-                                                   --output "${DATADIR}"/habitat_maps/urban
-
-python3 ./prepare_layers/make_constant_diff_map.py --current "${DATADIR}"/habitat/current_raw.tif \
-                                                   --habitat_code 14.5 \
-                                                   --crosswalk "${DATADIR}"/crosswalk.csv \
-                                                   --area "${DATADIR}"/area-per-pixel.tif \
-                                                   --scale 0.016666666666667 \
-                                                   --output "${DATADIR}"/habitat/urban_diff_area.tif
+if [ ! -d "${DATADIR}"/habitat_maps/urban ]; then
+    python3 ./prepare_layers/make_constant_habitat.py --examplar "${DATADIR}"/habitat_maps/arable/lcc_1401.tif \
+                                                    --habitat_code 14.5 \
+                                                    --crosswalk "${DATADIR}"/crosswalk.csv \
+                                                    --output "${DATADIR}"/habitat_maps/urban
+fi
+
+if [ ! -f "${DATADIR}"/habitat/urban_diff_area.tif ]; then
+    python3 ./prepare_layers/make_constant_diff_map.py --current "${DATADIR}"/habitat/current_raw.tif \
+                                                    --habitat_code 14.5 \
+                                                    --crosswalk "${DATADIR}"/crosswalk.csv \
+                                                    --area "${DATADIR}"/area-per-pixel.tif \
+                                                    --scale "${PIXEL_SCALE}" \
+                                                    --output "${DATADIR}"/habitat/urban_diff_area.tif
+fi
 
 # Fetch and prepare the elevation layers
-reclaimer zenodo --zenodo_id 5719984  --filename dem-100m-esri54017.tif --output "${DATADIR}"/elevation.tif
-gdalwarp -t_srs EPSG:4326 -tr 0.016666666666667 -0.016666666666667 -r max -co COMPRESS=LZW -wo NUM_THREADS=40 "${DATADIR}"/elevation.tif "${DATADIR}"/elevation-max.tif
-gdalwarp -t_srs EPSG:4326 -tr 0.016666666666667 -0.016666666666667 -r min -co COMPRESS=LZW -wo NUM_THREADS=40 "${DATADIR}"/elevation.tif "${DATADIR}"/elevation-min.tif
+if [ ! -f "${DATADIR}"/elevation.tif ]; then
+    reclaimer zenodo --zenodo_id 5719984  --filename dem-100m-esri54017.tif --output "${DATADIR}"/elevation.tif
+fi
+if [ ! -f "${DATADIR}"/elevation-max.tif ]; then
+    gdalwarp -t_srs EPSG:4326 -tr "${PIXEL_SCALE}" -"${PIXEL_SCALE}" -r max -co COMPRESS=LZW -wo NUM_THREADS=40 "${DATADIR}"/elevation.tif "${DATADIR}"/elevation-max.tif
+fi
+if [ ! -f "${DATADIR}"/elevation-min.tif ]; then
+    gdalwarp -t_srs EPSG:4326 -tr "${PIXEL_SCALE}" -"${PIXEL_SCALE}" -r min -co COMPRESS=LZW -wo NUM_THREADS=40 "${DATADIR}"/elevation.tif "${DATADIR}"/elevation-min.tif
+fi
 
 # Get species data per taxa from IUCN data
 for TAXA in "${TAXAS[@]}"
diff --git a/tests/test_food_layer.py b/tests/test_food_layer.py
new file mode 100644
index 0000000..96f5e0e
--- /dev/null
+++ b/tests/test_food_layer.py
@@ -0,0 +1,158 @@
+
+import numpy as np
+import pytest
+import yirgacheffe as yg
+from yirgacheffe.layers import RasterLayer
+from yirgacheffe.operators import DataType
+from yirgacheffe.window import Area, PixelScale
+
+from prepare_layers.make_food_current_map import TileInfo, process_tile
+
+@pytest.mark.parametrize("initial,crop_diff,pasture_diff,expected", [
+    (42, float("nan"), float("nan"), 42),
+
+    # Other habitat replacement
+    (42, 1.0, float("nan"), 1401),
+    (42, float("nan"), 1.0, 1402),
+    (42, 0.0, float("nan"), 42),
+    (42, float("nan"), 0.0, 42),
+    (42, -1.0, float("nan"), 42),
+    (42, float("nan"), -1.0, 42),
+
+    # Crop replacement
+    (1401, 1.0, float("nan"), 1401),
+    (1401, float("nan"), 1.0, 1401),
+    (1401, 0.0, float("nan"), 1401),
+    (1401, float("nan"), 0.0, 1401),
+    (1401, -1.0, float("nan"), 31),
+    (1401, -1.0, 1.0, 1402),
+    (1401, float("nan"), -1.0, 1401),
+
+    # Pasture replacement
+    (1402, 1.0, float("nan"), 1402),
+    (1402, 1.0, float("nan"), 1402),
+    (1402, float("nan"), 0.0, 1402),
+    (1402, 0.0, float("nan"), 1402),
+    (1402, float("nan"), -1.0, 31),
+    (1402, 1.0, -1.0, 1401),
+    (1402, -1.0, float("nan"), 1402),
+
+    # Don't touch urban
+    (1405, 1.0, float("nan"), 1405),
+    (1405, float("nan"), 1.0, 1405),
+    (1405, 0.0, float("nan"), 1405),
+    (1405, float("nan"), 0.0, 1405),
+    (1405, -1.0, float("nan"), 1405),
+    (1405, -1.0, 1.0, 1405),
+    (1405, float("nan"), -1.0, 1405),
+
+])
+def test_process_tile_all(initial, crop_diff, pasture_diff, expected) -> None:
+    with RasterLayer.empty_raster_layer(Area(-180, 90, 180, -90), PixelScale(1.0, -1.0), datatype=DataType.Int16) as current:
+        yg.constant(initial).save(current)
+        with RasterLayer.empty_raster_layer(Area(-180, 90, 180, -90), PixelScale(1.0, -1.0), datatype=DataType.Int16) as pnv:
+            yg.constant(31).save(pnv)
+
+            test_tile = TileInfo(
+                x_position=0,
+                y_position=0,
+                width=10,
+                height=12,
+                crop_diff=crop_diff,
+                pasture_diff=pasture_diff,
+            )
+
+            result = process_tile(current, pnv, test_tile)
+            expected = np.full((12, 10), expected, dtype=np.int16)
+            assert (result == expected).all()
+
+@pytest.mark.parametrize("initial,crop_diff,pasture_diff,expected_crop_count,expected_pasture_count,expected_pnv_count", [
+    (42, float("nan"), float("nan"), 0, 0, 0),
+    (42, 0.5, float("nan"), 50, 0, 0),
+    (42, float("nan"), 0.5, 0, 50, 0),
+    (42, -0.5, float("nan"), 0, 0, 0),
+    (42, float("nan"), -0.5, 0, 0, 0),
+
+    (1401, float("nan"), float("nan"), 100, 0, 0),
+    (1401, -0.5, float("nan"), 50, 0, 50),
+    (1401, float("nan"), -0.5, 100, 0, 0),
+
+    (1402, float("nan"), float("nan"), 0, 100, 0),
+    (1402, float("nan"), -0.5, 0, 50, 50),
+    (1402, -0.5, float("nan"), 0, 100, 0),
+])
+def test_partial_replacement(initial, crop_diff, pasture_diff, expected_crop_count, expected_pasture_count, expected_pnv_count) -> None:
+    with RasterLayer.empty_raster_layer(Area(-180, 90, 180, -90), PixelScale(1.0, -1.0), datatype=DataType.Int16) as current:
+        yg.constant(initial).save(current)
+        with RasterLayer.empty_raster_layer(Area(-180, 90, 180, -90), PixelScale(1.0, -1.0), datatype=DataType.Int16) as pnv:
+            yg.constant(31).save(pnv)
+
+            test_tile = TileInfo(
+                x_position=0,
+                y_position=0,
+                width=10,
+                height=10,
+                crop_diff=crop_diff,
+                pasture_diff=pasture_diff,
+            )
+
+            result = process_tile(current, pnv, test_tile)
+            crop_count = (result == 1401).sum()
+            assert crop_count == expected_crop_count
+            pasture_count = (result == 1402).sum()
+            assert pasture_count == expected_pasture_count
+            pnv_count = (result == 31).sum()
+            assert pnv_count == expected_pnv_count
+
+@pytest.mark.parametrize("initial,crop_diff,pasture_diff,expected_crop_count,expected_pasture_count,expected_pnv_count", [
+    (42, float("nan"), float("nan"), 0, 0, 0),
+    (42, 0.5, float("nan"), 50, 0, 0),
+    (42, float("nan"), 0.5, 0, 50, 0),
+    (42, -0.5, float("nan"), 0, 0, 0),
+    (42, float("nan"), -0.5, 0, 0, 0),
+
+    (1401, float("nan"), float("nan"), 50, 0, 0),
+    (1401, 1.0, float("nan"), 100, 0, 0),
+    (1401, 0.5, float("nan"), 100, 0, 0),
+    (1401, 0.1, float("nan"), 60, 0, 0),
+    (1401, -0.1, float("nan"), 40, 0, 10),
+    (1401, -0.5, float("nan"), 0, 0, 50),
+    (1401, -1.0, float("nan"), 0, 0, 50),
+    (1401, float("nan"), 1.0, 50, 50, 0),
+
+    (1405, float("nan"), float("nan"), 0, 0, 0),
+    (1405, 1.0, float("nan"), 50, 0, 0),
+    (1405, 0.5, float("nan"), 50, 0, 0),
+    (1405, 0.1, float("nan"), 10, 0, 0),
+    (1405, -0.1, float("nan"), 0, 0, 0),
+    (1405, -0.5, float("nan"), 0, 0, 0),
+    (1405, -1.0, float("nan"), 0, 0, 0),
+])
+def test_partial_initial_tile(initial, crop_diff, pasture_diff, expected_crop_count, expected_pasture_count, expected_pnv_count) -> None:
+    with RasterLayer.empty_raster_layer(Area(-180, 90, 180, -90), PixelScale(1.0, -1.0), datatype=DataType.Int16) as current:
+
+        # Cheating as Yirgacheffe doesn't have a neat way to provide numpy data straight to a layer
+        band = current._dataset.GetRasterBand(1)
+        vals = np.array([[initial, 22], [22, initial]])
+        all_vals = np.tile(vals, (90, 180))
+        band.WriteArray(all_vals, 0, 0)
+
+        with RasterLayer.empty_raster_layer(Area(-180, 90, 180, -90), PixelScale(1.0, -1.0), datatype=DataType.Int16) as pnv:
+            yg.constant(31).save(pnv)
+
+            test_tile = TileInfo(
+                x_position=0,
+                y_position=0,
+                width=10,
+                height=10,
+                crop_diff=crop_diff,
+                pasture_diff=pasture_diff,
+            )
+
+            result = process_tile(current, pnv, test_tile)
+            crop_count = (result == 1401).sum()
+            assert crop_count == expected_crop_count
+            pasture_count = (result == 1402).sum()
+            assert pasture_count == expected_pasture_count
+            pnv_count = (result == 31).sum()
+            assert pnv_count == expected_pnv_count
diff --git a/utils/footprint_of_humanity.py b/utils/footprint_of_humanity.py
index 4330f51..7399833 100644
--- a/utils/footprint_of_humanity.py
+++ b/utils/footprint_of_humanity.py
@@ -1,5 +1,6 @@
 # pylint: disable=C0301
 import argparse
+import os
 from pathlib import Path
 
 import numpy as np
@@ -49,6 +50,7 @@ def absolute(
     slimmed = pd.concat([slimmed_resident, slimmed_migratory])
     slimmed["extinction"] = slimmed.capped_current_persistence - slimmed.capped_scenario_persistence
 
+    os.makedirs(output_filename.parent, exist_ok=True)
     slimmed.to_csv(output_filename, index=False)
 
     for taxa in TAXA:
diff --git a/utils/raster_diff.py b/utils/raster_diff.py
new file mode 100644
index 0000000..d61055a
--- /dev/null
+++ b/utils/raster_diff.py
@@ -0,0 +1,51 @@
+import argparse
+import os
+from pathlib import Path
+
+import yirgacheffe as yg
+
+def raster_diff(
+    raster_a_path: Path,
+    raster_b_path: Path,
+    output_path: Path,
+) -> None:
+    os.makedirs(output_path.parent, exist_ok=True)
+
+    with yg.read_raster(raster_a_path) as raster_a:
+        with yg.read_raster(raster_b_path) as raster_b:
+            result = raster_a - raster_b
+            result.to_geotiff(output_path, parallelism=True)
+
+def main() -> None:
+    parser = argparse.ArgumentParser(description="Calculates the difference between two rasters")
+    parser.add_argument(
+        "--raster_a",
+        type=Path,
+        required=True,
+        dest="raster_a_path",
+        help="Left hand side of the difference."
+    )
+    parser.add_argument(
+        "--raster_b",
+        type=Path,
+        required=True,
+        dest="raster_b_path",
+        help="Right hands side of the difference"
+    )
+    parser.add_argument(
+        "--output",
+        type=Path,
+        required=True,
+        dest="output_filename",
+        help="Destination geotiff file for results."
+    )
+    args = parser.parse_args()
+
+    raster_diff(
+        args.raster_a_path,
+        args.raster_b_path,
+        args.output_filename,
+    )
+
+if __name__ == "__main__":
+    main()