Lesson 8.2.1: Converting model predictions into disease probability maps

flowchart TD
    A[Input Image] --> B[Deep Learning Model]
    B --> C[Raw logits]
    C --> D[Softmax → Probability map]
    D --> E[Optional Post‑processing]
    E --> F[Thresholding → Binary disease map]
    F --> G[Output Disease Probability Map]

Concept: From Scores to Probabilities

In computer vision for medical imaging, a deep learning model often outputs a score for each voxel or pixel in a scan. That score is usually a real‑valued number that reflects how strongly the model believes that the location belongs to a particular class, such as tumour tissue versus normal brain. However, for clinical reporting and downstream decision‑support, we rarely want a raw score; we need a probability that a voxel truly belongs to the disease class. Converting a score to a probability is a small but crucial mathematical step that bridges the gap between the model’s internal representation and a clinically interpretable metric.

Why Probabilities Matter

Probabilities have several desirable properties:

• Interpretability: A 0.8 probability is immediately understood as an 80 % chance that the voxel is diseased.
• Calibration: Probabilities can be calibrated against ground truth data so that, for example, all voxels whose model outputs a probability of 0.7 actually contain disease in 70 % of cases.
• Decision‑making: Thresholds can be chosen to optimise sensitivity or specificity, or to balance them via ROC analysis.
• Aggregation: Probabilities can be summed or averaged over patches, slices, or the entire volume to produce region‑level or patient‑level metrics.

Raw scores, on the other hand, are often unbounded or symmetric around zero and do not directly convey certainty. For instance, the model might output a score of 2.5 for a lesion and 0.5 for normal tissue; without a mapping to probabilities, it is unclear what those numbers mean in absolute terms.

Common Score‑to‑Probability Transformations

Two transformations dominate the literature: sigmoid and softmax. The choice depends on the number of classes and the structure of the network’s final layer.

Sigmoid for Binary or One‑vs‑Rest Problems

A sigmoid function maps any real‑valued input to the open interval (0, 1). It is defined as:

def sigmoid(x):
    return 1 / (1 + math.exp(-x))

Applied voxel‑wise to a score s, the sigmoid yields a probability p = sigmoid(s). The function has an S‑shaped curve that asymptotically approaches 0 and 1, ensuring that extreme scores translate into probabilities close to the limits.

Softmax for Multi‑Class Problems

When the model outputs a vector of scores [s1, s2, …, sk] for \(k\) classes, the softmax function normalises them into a probability distribution:

import math
def softmax(scores):
    exp_scores = [math.exp(s) for s in scores]
    sum_exp = sum(exp_scores)
    return [e / sum_exp for e in exp_scores]

Each output probability sums to 1 across all classes, making it straightforward to interpret which class is most likely.

Temperature Scaling

Sometimes a model is overconfident; its probabilities are too close to 0 or 1. Temperature scaling introduces a scalar T > 0 to soften the distribution:

def softmax_temperature(scores, T=1.0):
    exp_scores = [math.exp(s / T) for s in scores]
    sum_exp = sum(exp_scores)
    return [e / sum_exp for e in exp_scores]

A higher temperature (>1) flattens the distribution, while a lower temperature (<1) sharpens it. Temperature is typically tuned on a validation set to improve calibration.

Implementing the Transformation in Practice

Below is a concise example that demonstrates how to convert a 3‑D numpy array of logits (raw scores) into a probability map using the sigmoid function. The code is written in a way that is directly usable in a typical PyTorch or TensorFlow inference pipeline.

import numpy as np
import math

def sigmoid(x):
    "Return the logistic sigmoid of x."
    return 1.0 / (1.0 + math.exp(-x))

def logits_to_probability_map(logits):
    "
    Convert a 3-D array of logits to a probability map.

    Parameters
    ----------
    logits : np.ndarray
        3‑D array of shape (Z, Y, X) containing raw model scores.

    Returns
    -------
    probs : np.ndarray
        3‑D array of shape (Z, Y, X) with values in [0, 1].
    "
    # Vectorise the sigmoid to operate on the whole array efficiently
    vectorised_sigmoid = np.vectorize(sigmoid)
    probs = vectorised_sigmoid(logits)
    return probs

Example usage
logits = np.random.randn(64, 128, 128)  # random scores for a 64x128x128 volume
prob_map = logits_to_probability_map(logits)
print(prob_map.shape, prob_map.min(), prob_map.max())

Line‑by‑line explanation:

• import numpy as np and import math bring in the numerical libraries required for array operations and the exponential function.
• The sigmoid function implements the standard logistic function. It takes a single scalar x and returns a value between 0 and 1.
• logits_to_probability_map is a helper that accepts a 3‑D numpy array of logits. The docstring explains the purpose, parameters, and return type for clarity.
• Inside the function, np.vectorize(sigmoid) creates a vectorised version of the sigmoid that can be applied element‑wise across the entire array without explicit loops.
• The vectorised sigmoid is then applied to logits, producing probs, a probability map with the same shape.
• The example at the bottom generates a random 3‑D volume, converts it, and prints its shape and extremal values to confirm that all probabilities lie between 0 and 1.

Performance Considerations

In production, you would typically rely on the deep learning framework’s built‑in sigmoid (e.g. torch.sigmoid or tf.nn.sigmoid) to utilise GPU acceleration. The numpy example is primarily for educational purposes and can be easily swapped out:

import torch

def logits_to_probs_torch(logits_tensor):
    return torch.sigmoid(logits_tensor)

Using the framework’s native function ensures that the operation is differentiable and can be part of a larger computational graph if required for further optimisation.

Calibration: Aligning Predicted Probabilities with Reality

Even with a mathematically correct transformation, the model’s output probabilities may not be well calibrated. Calibration curves plot the empirical accuracy against predicted probabilities, revealing systematic over‑ or under‑confidence. Common calibration techniques include:

1. Platt Scaling: Fit a logistic regression model to map logits to probabilities on a validation set.
2. : A non‑parametric method that fits a piecewise‑constant function to the logits.
3. : Compute the weighted average difference between predicted probabilities and observed frequencies across probability bins.

After calibration, the probability map will better reflect the true risk of disease, which is critical when a clinician makes treatment decisions based on the map.

From Voxel‑Level Probabilities to Clinical Decision Support

Once a voxel‑wise probability map is produced, several downstream steps can be applied to translate the map into actionable information:

• Thresholding: Choose a probability threshold (e.g. 0.5) to generate a binary segmentation mask.
• Connected‑Component Analysis: Identify contiguous disease regions to compute lesion volumes.
• Statistical Summaries: Calculate the mean probability within a region of interest or the area under the ROC curve across patients.
• Visualisation: Overlay the probability map on the original scan with colour‑coding (e.g. red for high probability, blue for low) to aid visual inspection.

Each of these steps relies fundamentally on the reliability of the probability estimates. A well‑calibrated, sigmoid‑transformed probability map forms the foundation upon which robust, clinically useful tools can be built.

Common Pitfalls and How to Avoid Them

• Ignoring Numerical Stability: When logits are very large or very small, computing math.exp(-x) can overflow or underflow. Frameworks handle this internally, but in pure Python you should clip logits or use numerically stable formulations.
• Assuming Sigmoid Suffices for Multi‑Class: If you have more than two classes but still use a sigmoid per class, the probabilities will not sum to one. Use softmax instead.
• Overlooking Calibration: Raw sigmoid outputs are often overconfident. Always evaluate calibration on a held‑out dataset.
• Mis‑specifying Temperature: An inappropriate temperature can severely distort the probability distribution. Tune T on a validation set and monitor calibration metrics.

By rigorously applying the appropriate mathematical transformation, ensuring numerical stability, and validating calibration, clinicians and researchers can trust the probability maps produced by deep learning models. This confidence is essential for integrating such models into the diagnostic workflow and for ultimately improving patient outcomes.

Additional Example: Using NumPy’s Vectorised Sigmoid

Below is a compact NumPy implementation that leverages broadcasting to apply the sigmoid function to an entire array without explicit loops or np.vectorize. This approach is often faster for large volumes.

import numpy as np

def vectorised_sigmoid(logits):
    ""Return a probability map using NumPy broadcasting.""
    return 1.0 / (1.0 + np.exp(-logits))

Example usage
logits = np.random.randn(32, 64, 64)
probs = vectorised_sigmoid(logits)
print(probs.shape, probs.min(), probs.max())

Implementation/Code: Building a Probability Raster

In this section we step through the complete workflow that transforms raw model outputs—normally a list or array of class‑probabilities—into a tidy raster that can be visualised in any GIS platform. The process is split into three logical blocks:

• Preparing the reference grid and spatial metadata
• Reshaping the model output into a 2‑D array that matches the grid
• Writing the probability layer to disk with correct georeferencing information

All of the code below is written in Python 3.10+, using the most common geospatial libraries: rasterio for raster I/O, numpy for array manipulation, and geopandas to handle the reference vector that defines the study area. The example assumes a binary classifier that outputs a single probability value per pixel for the presence of a disease.

1. Setting up the reference grid

Before we can write a raster we must know its spatial extent, resolution, and coordinate reference system (CRS). Usually this information comes from an existing reference raster or from a vector layer that bounds the study area. In the example we read a vector shapefile that contains the bounding polygon of the target region.

import geopandas as gpd
import rasterio
from rasterio.transform import from_origin
import numpy as np

Load the vector with the study area boundary
boundary = gpd.read_file('data/study_area.shp')

Ensure the boundary is in a projected CRS that preserves area (e.g. EPSG:3857)
boundary = boundary.to_crs(epsg=3857)

Calculate the bounding box of the boundary
minx, miny, maxx, maxy = boundary.total_bounds

Define raster resolution (e.g. 30 metres per pixel)
pixel_size = 30

Compute the number of rows and columns
width = int(np.ceil((maxx - minx) / pixel_size))
height = int(np.ceil((maxy - miny) / pixel_size))

Build an affine transform that maps pixel coordinates to map units
transform = from_origin(minx, maxy, pixel_size, pixel_size)

Line‑by‑line explanation:

1. Import the required libraries.
2. Read the shapefile containing the study area into a GeoDataFrame.
3. Re‑project the geometry to EPSG:3857 (Web Mercator) which is a common projected CRS that preserves distances and is suitable for rasterisation.
4. Extract the minimum and maximum coordinates that bound the geometry.
5. Set the desired pixel size.
6. Calculate how many columns (width) and rows (height) are needed to cover the bounding box at the chosen resolution. np.ceil ensures we round up so no area is left out.
7. Create an affine transform that tells rasterio how to convert from pixel indices to map coordinates.

At this stage we have the skeleton of the raster: its shape, resolution, and CRS.

2. Reshaping the model output

The model typically returns a flat array of probabilities, one per pixel in the reference grid, ordered row‑major. We need to reshape this into a 2‑D numpy array that matches the raster dimensions. In practice the model might output a 3‑D array (batch, height, width) or a flat list; the code below handles the most common case.

# Mock model output: a flat list of probabilities
In practice this would come from your model inference pipeline
flat_probs = np.random.uniform(low=0.0, high=1.0, size=width * height)

Reshape into a 2‑D array (height rows, width columns)
prob_array = flat_probs.reshape((height, width))

Explanation:

1. Generate a random array to simulate model output. Replace this with the real probabilities from your classifier.
2. Use reshape to turn the flat vector into a 2‑D matrix. The first dimension is height (rows), the second is width (columns). numpy will fill the array in row‑major order, which aligns with how rasterio writes data.

In many cases you might need to mask out areas outside the study boundary. rasterio.features.rasterise can create a mask from the vector geometry. Below we show how to apply such a mask, setting pixels outside the polygon to np.nan so they are ignored in visualisations.

from rasterio.features import rasterize

Create a mask array where inside the boundary is 1, outside is 0
mask = rasterize(
    [(geom, 1) for geom in boundary.geometry],
    out_shape=(height, width),
    transform=transform,
    fill=0,
    dtype=np.uint8
)

Apply the mask: set probabilities outside the boundary to NaN
prob_array = np.where(mask == 1, prob_array, np.nan)

Now prob_array contains the disease probability for every pixel inside the study area, with NaN values indicating no data.

3. Writing the probability raster to disk

With the data array and spatial metadata ready, we can write the raster. rasterio supports various data types; for probability values a 32‑bit floating point is appropriate. We also embed a descriptive name and a colour map reference to aid interpretation in GIS viewers.

# Define output file path
output_path = 'output/disease_probability.tif'

Build the metadata dictionary
meta = {
    'driver': 'GTiff',
    'height': height,
    'width': width,
    'count': 1,            # single band
    'dtype': 'float32',
    'crs': 'EPSG:3857',
    'transform': transform,
    'nodata': np.nan,
    'compress': 'lzw',
    'description': 'Probability of disease presence'
}

Write the raster
with rasterio.open(output_path, 'w', **meta) as dst:
    dst.write(prob_array.astype('float32'), 1)

Line‑by‑line:

1. Specify where the resulting GeoTIFF will be stored.
2. Construct a dictionary of metadata:
   • driver specifies the file format.
   • Dimension and band count.
   • Data type float32.
   • CRS matches the reference grid.
   • Affine transform from earlier.
   • nodata set to np.nan so that GIS software recognises missing values.
   • Compression to reduce file size.
   • A human‑readable description.
3. Open the file for writing, passing the metadata.
4. Write the probability array to the first (and only) band, converting to float32 to match the declared dtype.

After running the script, the file disease_probability.tif can be loaded into QGIS, ArcGIS, or any other GIS application. When displayed, the probability values can be colour‑graded using a suitable colour ramp (e.g. green for low risk, red for high risk). The NaN cells outside the boundary will appear transparent, ensuring a clean visualisation.

Common pitfalls and best practices

Extending the workflow

Once the basic probability raster is in place, several enhancements are common in production pipelines:

1. Thresholding – Convert probabilities to binary presence/absence maps by applying a threshold (e.g. 0.5). This can be done in numpy with prob_array > 0.5 and then writing a uint8 raster.
2. Uncertainty quantification – If the model produces a distribution of predictions (e.g. via Monte Carlo dropout), write separate raster layers for mean, variance, or confidence intervals.
3. Multi‑class probability rasters – For multi‑class problems, reshape the output to (n_classes, height, width) and write each class as a separate band. Include a colour table that maps each band to a distinct colour.
4. Geoprocessing automation – Wrap the entire process in a Python function or a Snakemake rule, enabling reproducible runs across multiple study areas.

By following the steps above you will reliably convert any model’s raw probability outputs into a GIS‑ready raster that can be analysed, visualised, and shared with stakeholders. The key lies in aligning spatial metadata, handling missing data correctly, and preserving numerical precision throughout the pipeline.

Alternative approach: clipping with rasterio.mask

Instead of building a mask array manually, you can clip the probability raster as soon as it is written, using rasterio.mask.mask. This is handy when the model output is already in raster form or when you want to keep the workflow tidy.

from rasterio.mask import mask

Assuming `prob_array` has been written to disk already
with rasterio.open(output_path, 'r+') as src:
    out_image, out_transform = mask(src, boundary.geometry, crop=True, nodata=np.nan)
    # out_image is a 3‑D array (bands, rows, cols); keep the first band
    src.update_tags(1, description='Clipped probability')
    src.write(out_image[0].astype('float32'), 1)

Clipping in this way ensures the final file only contains data within the study boundary, potentially reducing file size.

3.1 Alternative Implementation: Random Forest Classifier

Below is a compact Python script that replaces the logistic regression with a random forest classifier from scikit-learn. This approach captures non‑linear relationships and interactions between covariates, which can be important for complex disease systems.

import geopandas as gpd
import rasterio
import numpy as np
import pandas as pd
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import roc_auc_score
from rasterio.transform import from_origin

1. Load covariate rasters
rasters = {
    'temp': rasterio.open('temp.tif'),
    'rain': rasterio.open('rain.tif'),
    'elev': rasterio.open('elev.tif'),
    'land': rasterio.open('land.tif'),
    'popd': rasterio.open('popd.tif')
}

2. Read hotspot points
hotspots = gpd.read_file('hotspots.shp')
hotspots = hotspots.to_crs(rasters['temp'].crs)

3. Extract raster values at hotspot locations
def extract_values(gdf, raster, name):
    coords = [(x, y) for x, y in zip(gdf.geometry.x, gdf.geometry.y)]
    values = [val[0] for val in raster.sample(coords)]
    gdf[name] = values

for key in rasters:
    extract_values(hotspots, rasters[key], key)

4. Prepare design matrix and response
X = hotspots[['temp', 'rain', 'elev', 'land', 'popd']].values
y = hotspots['label'].values

5. Fit Random Forest
rf = RandomForestClassifier(
    n_estimators=200,
    max_depth=10,
    random_state=42,
    n_jobs=-1
)
rf.fit(X, y)

6. Model evaluation
y_pred_proba = rf.predict_proba(X)[:, 1]
print(f'Cross‑validated AUC: {roc_auc_score(y, y_pred_proba):.3f}')

7. Predict across the study area
meta = rasters['temp'].meta
rows, cols = meta['height'], meta['width']
transform = meta['transform']

Stack rasters into a 2‑D array (rows*cols, covariates)
stack = np.vstack([
    rasters['temp'].read(1).ravel(),
    rasters['rain'].read(1).ravel(),
    rasters['elev'].read(1).ravel(),
    rasters['land'].read(1).ravel(),
    rasters['popd'].read(1).ravel()
]).T

mask = np.any(np.isnan(stack), axis=1)
stack_clean = stack[~mask]

Predict probabilities
probs = rf.predict_proba(stack_clean)[:, 1]

Re‑inject masked values
full_probs = np.full((rows * cols), np.nan)
full_probs[~mask] = probs
prob_raster = full_probs.reshape((rows, cols))

8. Write probability raster
with rasterio.open(
    'malaria_prob_rf.tif',
    'w',
    driver='GTiff',
    height=rows,
    width=cols,
    count=1,
    dtype=rasterio.float32,
    crs=meta['crs'],
    transform=transform,
) as dst:
    dst.write(prob_raster.astype(rasterio.float32), 1)

This script demonstrates how to swap the modelling backbone while keeping the spatial workflow unchanged. The resulting GeoTIFF, malaria_prob_rf.tif, can be compared directly with the logistic‑regression output to assess model performance and spatial bias.

Pitfalls: Common Errors and How to Avoid Them

When converting the raw outputs of a predictive model into spatial probability maps, a range of subtle yet consequential errors can creep in. These mistakes often arise from a mismatch between the statistical assumptions of the model, the preprocessing of the data, and the visualisation pipeline. Below we dissect the most frequent pitfalls, analyse why they occur, and present concrete strategies to eliminate them.

1. Ignoring Class Imbalance During Training

In many disease mapping studies, the proportion of infected individuals is far lower than that of healthy ones. If the model is trained on such imbalanced data without adjusting for it, the output probabilities will be systematically biased toward the majority class. As a result, the resulting map will under‑estimate risk in genuinely high‑prevalence areas.

• Solution: Apply class‑weighting or resampling techniques during model fitting. Most scikit‑learn estimators accept a class_weight argument; in neural networks, use a weighted loss function.
• Verification: Examine the calibration curve of the model on a held‑out set; a perfectly calibrated model will produce a diagonal line in the plot.

2. Mis‑scaled Input Features

When raster layers are fed into the model, each band often carries a different numerical range (e.g. NDVI in [−1, 1] versus temperature in [−50, 50]). Without normalising these inputs, the optimiser will give undue weight to features with larger absolute values, distorting the probability estimates.

• Solution: Standardise all covariates to zero mean and unit variance, or scale them to a common range such as [0, 1]. Keep the scaler fitted on the training set and apply the same transformation to the prediction tiles.
• Tip: Store the scaler parameters in a JSON file so they can be re‑used during inference.

3. Overlooking Spatial Autocorrelation

Most machine‑learning algorithms assume independence between samples. In spatial epidemiology, this assumption is rarely true; nearby locations often exhibit similar disease risk due to unmeasured environmental drivers. Ignoring this can lead to overconfident probability estimates.

• Solution: Incorporate a spatial random effect (e.g. via a Gaussian Process or Conditional Autoregressive model) or use a mixed‑effects framework. Alternatively, add a spatial lag of the predictors to the feature set.
• Check: Compute the Moran’s I statistic on the residuals; a significant value indicates remaining spatial structure.

4. Inaccurate Raster Alignment and Projection

Probability maps are only meaningful if they align perfectly with the reference geographic coordinate system. A small mis‑registration between the model output and the target grid can produce misleading hotspots.

• Solution: Use geospatial libraries (e.g. rasterio) to reproject all rasters to the same CRS before feeding them into the model. Validate the alignment by overlaying a vector boundary.
• Validation step: Run a quick overlay test: rasterio.plot.show((raster1, raster2)) and confirm that the visual alignment matches.

5. Over‑Thresholding Probabilities

Converting continuous probabilities to binary maps by applying a hard threshold (e.g. 0.5) without considering the prevalence or cost of false positives/negatives can produce artefacts. A single pixel may be misclassified simply because the threshold was arbitrary.

• Solution: Use a cost‑sensitive threshold derived from the ROC curve or from a decision‑analytic framework such as Youden’s J statistic. Alternatively, retain the continuous probability field and let subsequent analyses (e.g. risk zoning) decide thresholds.
• Example threshold optimisation: optimal_threshold = roc_curve(y_true, y_pred)[1][np.argmax(roc_curve(y_true, y_pred)[2] - roc_curve(y_true, y_pred)[1])]

6. Neglecting Uncertainty Propagation

Most practitioners present a single deterministic probability map, masking the inherent uncertainty of the model. This can be misleading, especially in policy contexts where risk thresholds are critical.

• Solution: Quantify predictive uncertainty using techniques such as Monte‑Carlo dropout, ensemble methods, or Bayesian neural networks. Visualise the standard deviation or credible intervals as an overlay or separate thematic layer.
• Practical tip: Store the full distribution of predictions for each pixel in a NetCDF file for downstream analysis.

7. Failure to Re‑Calibrate Post‑Training

Even a well‑trained model may produce poorly calibrated probabilities due to changes in data distribution between training and deployment. This is particularly common when the model is applied to a new region.

• Solution: Employ post‑hoc calibration methods such as Platt scaling or isotonic regression on a calibration set drawn from the target area. Libraries like scikit‑learn provide CalibratedClassifierCV for this purpose.
• Check calibration: Plot a reliability diagram; the closer the points lie to the diagonal, the better the calibration.

8. Over‑Fitting to a Limited Dataset

High‑dimensional models with too many predictors relative to the number of observations tend to over‑fit, producing probability maps that look convincing but are not generalisable.

• Solution: Apply dimensionality reduction (e.g. PCA) or feature selection based on mutual information. Use cross‑validation to monitor generalisation performance.
• Rule of thumb: For logistic regression, keep the number of predictors below one‑third of the number of positive cases.

9. Poor Handling of Missing Data

Spatial datasets often contain gaps due to cloud cover or sensor failure. If missing values are not properly imputed, the model may produce NaNs or biased probabilities in those areas.

• Solution: Use spatial interpolation (e.g. kriging) or machine‑learning imputers (e.g. KNNImputer) that respect spatial structure. Always flag imputed areas in the final map.
• Implementation hint: In rasterio, use rasterio.fillnodata with a suitable method before feeding the data into the model.

10. Assuming Homogeneous Model Performance Across Space

It is tempting to apply a single probability map across the entire study area, assuming the model performs equally well everywhere. In reality, model error can vary spatially due to unmeasured covariates or differing data quality.

• Solution: Partition the study area into sub‑regions based on environmental strata or administrative boundaries, train local models, or at least evaluate spatially varying performance metrics.
• Evaluation: Compute the confusion matrix for each sub‑region and visualise the spatial distribution of sensitivity and specificity.

Practical Implementation Example

Below is a concise Python snippet demonstrating how to convert raw model logits into calibrated probabilities, apply a cost‑sensitive threshold, and save the result as a GeoTIFF. The code uses scikit‑learn for calibration and rasterio for geospatial handling.

import numpy as np
import rasterio
from rasterio.transform import from_origin
from sklearn.calibration import CalibratedClassifierCV
from sklearn.linear_model import LogisticRegression

Load pre‑processed raster patches (features) and corresponding labels
X = np.load('features.npy')          # shape: (n_samples, n_features)
y = np.load('labels.npy')            # shape: (n_samples,)

Train a base classifier
base_clf = LogisticRegression(max_iter=1000)
base_clf.fit(X, y)

Calibrate the probabilities using Platt scaling
calibrated_clf = CalibratedClassifierCV(base_clf, method='sigmoid', cv=5)
calibrated_clf.fit(X, y)

Predict probabilities on a new raster block
X_new = np.load('new_features.npy')   # same feature engineering as training
probabilities = calibrated_clf.predict_proba(X_new)[:, 1]

Determine optimal threshold via Youden's J
from sklearn.metrics import roc_curve
fpr, tpr, thresholds = roc_curve(y, calibrated_clf.predict_proba(X)[:, 1])
youden_index = tpr - fpr
optimal_thr = thresholds[np.argmax(youden_index)]

Apply threshold to obtain binary map
binary_map = (probabilities >= optimal_thr).astype(np.uint8)

Save as GeoTIFF
transform = from_origin(-180, 90, 0.01, 0.01)  # example geotransform
with rasterio.open(
    'probability_map.tif',
    'w',
    driver='GTiff',
    height=binary_map.shape[0],
    width=binary_map.shape[1],
    count=1,
    dtype=binary_map.dtype,
    crs='+proj=longlat',
    transform=transform,
) as dst:
    dst.write(binary_map, 1)

Line‑by‑line explanation

• import numpy as np – Load NumPy for numerical operations.
• import rasterio – Rasterio handles reading and writing GeoTIFFs.
• from rasterio.transform import from_origin – Create an affine transform for georeferencing.
• from sklearn.calibration import CalibratedClassifierCV – Wrapper for probability calibration.
• from sklearn.linear_model import LogisticRegression – Base classifier.
• Load training features and labels from pre‑saved NumPy arrays.
• Instantiate and fit the logistic regression model.
• Wrap the base classifier in CalibratedClassifierCV with Platt scaling (method='sigmoid') and 5‑fold cross‑validation.
• Fit the calibrated model on the same data.
• Predict raw probabilities on a new raster block.
• Compute the ROC curve on training data to obtain false positive/true positive rates.
• Calculate Youden's J statistic (tpr - fpr) to identify the threshold that maximises sensitivity + specificity.
• Determine the optimal threshold from the ROC array.
• Generate a binary map by thresholding the probabilities.
• Define the geotransform (here a placeholder of 0.01 degree resolution).
• Write the binary map to a GeoTIFF with appropriate metadata (CRS, dimensions, datatype).

By following the guidelines above, you minimise the risk of producing misleading disease probability maps. The combination of rigorous preprocessing, model calibration, spatially aware validation, and transparent uncertainty quantification ensures that the final visualisations are both scientifically robust and policy‑relevant.

Additional Example: Bayesian Neural Network for Uncertainty Estimation

Below is a short snippet that demonstrates how to construct a simple Bayesian neural network in Pyro, predict probabilistic outcomes, and visualise the mean and standard deviation of the predictive distribution over a raster field.

import torch
import pyro
import pyro.distributions as dist
from pyro.infer import Predictive
import rasterio
import numpy as np

Define a Bayesian neural network with one hidden layer
class BNN(torch.nn.Module):
    def __init__(self, input_dim, hidden_dim=32):
        super().__init__()
        self.fc1 = torch.nn.Linear(input_dim, hidden_dim)
        self.fc2 = torch.nn.Linear(hidden_dim, 1)

    def forward(self, x):
        h = torch.relu(self.fc1(x))
        return torch.sigmoid(self.fc2(h))

Wrap the network in a Pyro model
def model(x, y=None):
    weight1 = pyro.sample('weight1', dist.Normal(0., 1.).expand([input_dim, hidden_dim]).to_event(2))
    bias1   = pyro.sample('bias1',   dist.Normal(0., 1.).expand([hidden_dim]).to_event(1))
    weight2 = pyro.sample('weight2', dist.Normal(0., 1.).expand([hidden_dim, 1]).to_event(2))
    bias2   = pyro.sample('bias2',   dist.Normal(0., 1.).expand([1]).to_event(1))

    def net(x):
        h = torch.relu(torch.matmul(x, weight1) + bias1)
        return torch.sigmoid(torch.matmul(h, weight2) + bias2)

    with pyro.plate('data', x.shape[0]):
        obs = pyro.sample('obs', dist.Bernoulli(net(x)), obs=y)
    return obs

Assume features have been pre‑processed and loaded
X = torch.tensor(np.load('features.npy'), dtype=torch.float32)
input_dim = X.shape[1]

Perform posterior predictive sampling
predictive = Predictive(model, num_samples=200)
samples = predictive(X)

Compute mean and std dev of predictions
mean_pred = samples['obs'].float().mean(0).numpy()
std_pred  = samples['obs'].float().std(0).numpy()

Save mean and std dev maps as separate GeoTIFFs
with rasterio.open('bnn_mean.tif', 'w', driver='GTiff', height=mean_pred.shape[0],
                   width=mean_pred.shape[1], count=1, dtype='float32',
                   crs='+proj=longlat', transform=from_origin(-180, 90, 0.01, 0.01)) as dst:
    dst.write(mean_pred, 1)

with rasterio.open('bnn_std.tif', 'w', driver='GTiff', height=std_pred.shape[0],
                   width=std_pred.shape[1], count=1, dtype='float32',
                   crs='+proj=longlat', transform=from_origin(-180, 90, 0.01, 0.01)) as dst:
    dst.write(std_pred, 1)

Approach Comparison Table

Further Reading

• World Health Organisation (2018) — Global Atlas of Disease Prevalence. A comprehensive guide to mapping disease burden using model outputs. This atlas demonstrates best practices for converting statistical predictions into spatial probability surfaces.
• Waller, L.E. & Gotway, C.A. (2004) — Applied Spatial Statistics for Public Health Data. The book offers detailed methodology for spatial interpolation and risk mapping, including case studies that transform model predictions into probability maps.
• Wang, H., et al. (2019) — Spatial Epidemiology: Methods and Applications. This peer‑reviewed article reviews modern techniques for converting predictive models into risk surfaces, emphasising validation and uncertainty quantification.
• Open Geospatial Consortium (OGC) (2021) — Web Coverage Service (WCS) Specification. Provides the technical standard for delivering continuous probability maps from model servers, essential for interoperable disease probability mapping.
• Jiang, Y., & Smith, R. (2022) — Bayesian Spatial Models in Epidemiology. The chapter discusses hierarchical modelling and the generation of probability maps, with practical examples in R and Python.
• Centres for Disease Control and Prevention (CDC) (2020) — Spatial Epidemiology and Risk Mapping. Official documentation outlining CDC’s approach to creating and interpreting disease probability maps from predictive models.