Please, someone help me to calculate Potential Evapotranspiration (PET) using raster temperature?

by

Summary

This postmortem analyzes a failure that occurred when attempting to compute Potential Evapotranspiration (PET) using the Thornthwaite method in R. The issue stemmed from incorrectly deriving per‑pixel latitude from a raster stack, leading to invalid PET outputs and runtime errors.

Root Cause

The root cause was the misuse of:

  • yFromCell(Temperature_stack, 1:nrow(Temperature_stack))

This call returns one latitude value per row, not one latitude value per pixel.
However, the Thornthwaite function expects:

  • A single latitude value, or
  • A vector of latitudes matching the number of time steps, not the number of pixels.

This mismatch caused the PET computation to fail.

Why This Happens in Real Systems

Real geospatial pipelines often break due to:

  • Incorrect assumptions about raster geometry
  • Confusing “rows” with “cells”
  • Passing spatial vectors where statistical functions expect temporal vectors
  • Using raster stacks without checking dimensionality

These mistakes are common when combining spatial data structures with climatological formulas that were never designed for per‑pixel computation.

Real-World Impact

Such errors can lead to:

  • Incorrect PET maps, causing flawed hydrological modeling
  • Silent propagation of wrong climate indicators
  • Wasted compute time on large raster stacks
  • Misleading scientific conclusions if unnoticed

Example or Code (if necessary and relevant)

Below is a corrected pattern showing how senior engineers compute PET per pixel by iterating or vectorizing properly:

library(raster)
library(SPEI)

tstack <- stack(list.files("C:/data", full.names = TRUE, pattern = "tif"))
lat_raster <- init(tstack, 'y')

pet <- calc(tstack, fun = function(ts) {
  thornthwaite(ts, lat = mean(ts, na.rm = TRUE))
})

How Senior Engineers Fix It

Experienced engineers avoid this class of bugs by:

  • Inspecting raster dimensions before applying climate formulas
  • Separating spatial and temporal vectors
  • Using calc() or overlay() to apply functions per pixel
  • Deriving latitude as a raster, not a vector
  • Validating assumptions about what each function expects

They also test with:

  • A single pixel
  • A single time slice
  • A known latitude

before scaling up.

Why Juniors Miss It

Junior engineers often miss this because:

  • They assume “one value per cell” when the function expects “one value per time step”
  • They rely on trial‑and‑error instead of reading function signatures
  • They misunderstand how raster stacks represent space vs. time
  • They expect PET libraries to be spatially aware, when most are not

The result is a subtle but common mismatch between spatial data structures and climatological algorithms, leading to confusing failures.

Leave a Comment