Phase 3 - Protocol Standards

From Raw Data to Actionable Intelligence

Phase 3 of the WIA standard addresses the critical question: how should raw satellite imagery, weather station data, and sensor measurements be processed into standardized drought indices? While Phases 1 and 2 define what drought data looks like and how to access it, Phase 3 ensures that different implementations calculate drought indices consistently using validated scientific methods.

Without processing protocol standards, two systems could claim WIA compliance while producing completely different drought assessments for the same location and time, defeating the purpose of standardization.

Satellite Data Acquisition Standards

The foundation of modern drought monitoring is satellite remote sensing. The WIA standard specifies protocols for acquiring and preprocessing data from major Earth observation satellites.

Supported Satellite Platforms

Satellite	Sensor	Products	WIA Standard Resolution
MODIS Terra/Aqua	MOD13Q1/MYD13Q1	NDVI, EVI	250m, 16-day composite
Landsat 8/9	OLI, TIRS	NDVI, LST	30m, 16-day revisit
Sentinel-2	MSI	NDVI, Moisture Index	10-20m, 5-day revisit
SMAP	L-band Radiometer	Soil Moisture	9km, 2-3 day revisit
SMOS	MIRAS	Soil Moisture	25km, 3-day revisit

Atmospheric Correction Protocol

Raw satellite reflectance values are affected by atmospheric conditions. The WIA standard requires atmospheric correction to obtain surface reflectance suitable for vegetation index calculation.

Standard Method: 6S Radiative Transfer Model

The Second Simulation of a Satellite Signal in the Solar Spectrum (6S) is the WIA standard method for atmospheric correction:

// Atmospheric Correction Pseudocode

function atmosphericCorrection(rawReflectance, sensorGeometry, atmosphericConditions) {
  // Input parameters
  const {
    solarZenithAngle,
    sensorZenithAngle,
    relativeAzimuth
  } = sensorGeometry;

  const {
    waterVapor_cm,
    ozone_atm_cm,
    aerosolOpticalDepth,
    surfacePressure_mb
  } = atmosphericConditions;

  // Run 6S model
  const correction = run6SModel({
    geometry: sensorGeometry,
    atmosphere: atmosphericConditions,
    wavelength: sensorBands
  });

  // Apply correction
  const surfaceReflectance = (rawReflectance - correction.pathReflectance) /
                              correction.transmittance;

  return surfaceReflectance;
}

// Quality requirements
// - RMSE vs ground targets: < 5%
// - Residual atmospheric effects: < 2% reflectance units

Alternative Methods (If 6S Unavailable)

Method	Applicability	Accuracy	Computational Cost
DOS (Dark Object Subtraction)	Quick processing	Not specified (ground-target validation required)	Very Low
QUAC (Quick Atmospheric Correction)	No atmospheric data needed	Not specified (ground-target validation required)	Low
FLAASH	Multispectral and hyperspectral imagery	Not specified (ground-target validation required)	High
Sen2Cor (Sentinel-2)	Sentinel-2 specific	Not specified (ground-target validation required)	Medium

Cloud Masking Standards

Clouds obscure surface conditions and must be identified and masked before drought index calculation. The WIA standard adopts the Fmask (Function of mask) algorithm for cloud detection.

Fmask Cloud Detection Algorithm

function cloudMask(surfaceReflectance, brightnessTemp, metadata) {
  // Stage 1: Potential cloud pixels
  const potentialClouds = identifyPotentialClouds({
    blue: surfaceReflectance.blue,
    green: surfaceReflectance.green,
    red: surfaceReflectance.red,
    nir: surfaceReflectance.nir,
    swir1: surfaceReflectance.swir1,
    swir2: surfaceReflectance.swir2,
    thermal: brightnessTemp
  });

  // Stage 2: Cloud probability
  const cloudProbability = calculateCloudProbability(potentialClouds, {
    temperature: brightnessTemp,
    whiteness: calculateWhiteness(surfaceReflectance),
    brightness: calculateBrightness(surfaceReflectance),
    variability: calculateSpectralVariability(surfaceReflectance)
  });

  // Stage 3: Cloud shadows
  const shadowMask = detectCloudShadows({
    cloudMask: cloudProbability > 0.5,
    solarAzimuth: metadata.solarAzimuth,
    solarZenith: metadata.solarZenith,
    cloudHeight: estimateCloudHeight(brightnessTemp)
  });

  // Final mask
  return {
    cloudMask: cloudProbability > 0.5,
    cloudShadowMask: shadowMask,
    cloudProbability: cloudProbability,
    quality: assessMaskQuality(cloudProbability, shadowMask)
  };
}

// Quality Requirements:
// - Cloud detection accuracy: > 95%
// - Commission error (false positives): < 3%
// - Omission error (missed clouds): < 5%

NDVI Calculation Protocol

The Normalized Difference Vegetation Index is fundamental to vegetation-based drought monitoring. The WIA standard specifies exact calculation methods and quality requirements.

Standard NDVI Calculation

function calculateNDVI(surfaceReflectance, cloudMask) {
  // Extract spectral bands (atmospherically corrected)
  const red = surfaceReflectance.red;     // 0.63-0.69 μm
  const nir = surfaceReflectance.nir;     // 0.76-0.90 μm

  // Apply cloud mask
  const validPixels = cloudMask.cloudMask.map(
    (isCloud, i) => !isCloud && !cloudMask.cloudShadowMask[i]
  );

  // Calculate NDVI where valid
  const ndvi = validPixels.map((valid, index) => {
    if (!valid) return NaN;

    const nirVal = nir[index];
    const redVal = red[index];

    // Avoid division by zero
    if (Math.abs(nirVal + redVal) < 0.001) return NaN;

    // Standard NDVI formula
    return (nirVal - redVal) / (nirVal + redVal);
  });

  // Quality control
  const qualityNDVI = applyQualityControl(ndvi, {
    validRange: [-1.0, 1.0],
    expectedVegetationRange: [0.2, 0.9],
    adjacencyEffect: checkAdjacentPixels(ndvi)
  });

  return {
    ndvi: qualityNDVI,
    validPixelPercent: calculateValidPercent(qualityNDVI),
    meanNDVI: calculateMean(qualityNDVI),
    stddevNDVI: calculateStdDev(qualityNDVI)
  };
}

// Precision requirement: 0.01 NDVI units
// Cross-sensor consistency: ±0.02 NDVI units

NDVI Anomaly Calculation

Drought detection requires comparing current NDVI to historical norms:

function calculateNDVIAnomaly(currentNDVI, historicalData, location, date) {
  // Extract same period from historical record
  const samePeriod = historicalData.filter(record =>
    isSamePeriod(record.date, date, 16)  // window: ±16 days
  );

  // Calculate statistics
  const historicalMean = calculateMean(samePeriod.map(r => r.ndvi));
  const historicalStdDev = calculateStdDev(samePeriod.map(r => r.ndvi));

  // Anomaly calculation
  const absoluteAnomaly = currentNDVI - historicalMean;
  const standardizedAnomaly = absoluteAnomaly / historicalStdDev;
  const percentile = calculatePercentile(currentNDVI, samePeriod.map(r => r.ndvi));

  // Drought classification based on anomaly
  let droughtStatus;
  if (standardizedAnomaly < -2.0) droughtStatus = "severe_drought";
  else if (standardizedAnomaly < -1.5) droughtStatus = "moderate_drought";
  else if (standardizedAnomaly < -1.0) droughtStatus = "mild_drought";
  else droughtStatus = "normal";

  return {
    ndvi: currentNDVI,
    anomaly: absoluteAnomaly,
    zscore: standardizedAnomaly,
    percentile: percentile,
    historicalMean: historicalMean,
    droughtStatus: droughtStatus
  };
}

// Historical baseline requirement: Minimum 10 years of data
// Window matching: ±16 days from current date

Evapotranspiration Calculation Standards

Evapotranspiration (ET) quantifies water loss, essential for water balance drought indices like PDSI. The WIA standard adopts the FAO-56 Penman-Monteith equation.

FAO-56 Penman-Monteith Method

function calculateET0(weatherData) {
  // Weather inputs
  const {
    temperatureMin_C,
    temperatureMax_C,
    relativeHumidityMin_percent,
    relativeHumidityMax_percent,
    windSpeed2m_ms,
    solarRadiation_MJm2day,
    latitude,
    elevation_m,
    date
  } = weatherData;

  // Derived parameters
  const tempMean = (temperatureMin_C + temperatureMax_C) / 2;
  const delta = saturationVaporPressureSlope(tempMean);  // kPa/°C
  const P = atmosphericPressure(elevation_m);  // kPa
  const gamma = psychrometricConstant(P);  // kPa/°C

  // Saturation vapor pressure
  const es = (saturationVaporPressure(temperatureMax_C) +
              saturationVaporPressure(temperatureMin_C)) / 2;

  // Actual vapor pressure
  const ea = (saturationVaporPressure(temperatureMin_C) * relativeHumidityMax_percent / 100 +
              saturationVaporPressure(temperatureMax_C) * relativeHumidityMin_percent / 100) / 2;

  // Net radiation
  const Rn = netRadiation(solarRadiation_MJm2day, latitude, date, tempMean, ea);

  // Soil heat flux (negligible for daily calculations)
  const G = 0;

  // Wind speed adjustment to 2m height (if needed)
  const u2 = windSpeed2m_ms;

  // FAO-56 Penman-Monteith equation
  const numerator = 0.408 * delta * (Rn - G) +
                     gamma * (900 / (tempMean + 273)) * u2 * (es - ea);
  const denominator = delta + gamma * (1 + 0.34 * u2);

  const ET0 = numerator / denominator;  // mm/day

  return {
    et0_mm_day: Math.max(0, ET0),
    temperature_c: tempMean,
    vaporPressureDeficit_kPa: es - ea,
    netRadiation_MJm2day: Rn,
    quality: assessETQuality(weatherData)
  };
}

// Supporting functions
function saturationVaporPressure(temp_C) {
  return 0.6108 * Math.exp((17.27 * temp_C) / (temp_C + 237.3));
}

function saturationVaporPressureSlope(temp_C) {
  const es = saturationVaporPressure(temp_C);
  return (4098 * es) / Math.pow(temp_C + 237.3, 2);
}

function atmosphericPressure(elevation_m) {
  return 101.3 * Math.pow((293 - 0.0065 * elevation_m) / 293, 5.26);
}

function psychrometricConstant(P_kPa) {
  return 0.000665 * P_kPa;
}

// Accuracy requirement: RMSE < 15% vs eddy covariance measurements

Crop-Specific ET Calculation

Converting reference ET to crop-specific ET using crop coefficients:

function calculateETc(ET0, cropType, growthStage, stressConditions) {
  // Crop coefficient lookup
  const Kc = getCropCoefficient(cropType, growthStage);

  // Stress coefficient (reduces ET under water stress)
  const Ks = calculateStressCoefficient(stressConditions);

  // Crop ET
  const ETc = ET0 * Kc * Ks;

  return {
    etc_mm_day: ETc,
    et0_mm_day: ET0,
    cropCoefficient: Kc,
    stressCoefficient: Ks,
    cropType: cropType,
    growthStage: growthStage
  };
}

// Crop coefficient database (example values)
const cropCoefficients = {
  wheat: {
    initial: 0.30,
    development: 0.75,
    mid_season: 1.15,
    late_season: 0.40
  },
  corn: {
    initial: 0.30,
    development: 0.80,
    mid_season: 1.20,
    late_season: 0.60
  },
  alfalfa: {
    cutting_period: 1.20,
    regrowth: 0.40
  }
};

// Real implementations build the crop coefficient database from validated sources such as the FAO-56 crop coefficient tables

Soil Moisture Processing Protocol

For satellite-derived soil moisture (SMAP, SMOS), the WIA standard specifies retrieval algorithms and quality control:

L-Band Radiometry Soil Moisture Retrieval

function retrieveSoilMoisture(brightnessTemperature, landCover, soilTexture, vegetation) {
  // Input: L-band brightness temperature (1.4 GHz)
  const Tb_h = brightnessTemperature.horizontal;
  const Tb_v = brightnessTemperature.vertical;

  // Remove vegetation attenuation
  const tau = vegetationOpticalDepth(vegetation.waterContent, vegetation.type);
  const Tsoil_h = removeVegetationEffect(Tb_h, tau, vegetation.temperature);
  const Tsoil_v = removeVegetationEffect(Tb_v, tau, vegetation.temperature);

  // Dielectric model inversion
  const soilDielectric = calculateDielectric(Tsoil_h, Tsoil_v, soilTexture);

  // Convert dielectric to volumetric soil moisture
  const soilMoisture = dielectricToMoisture(soilDielectric, soilTexture);

  // Quality control
  const quality = assessSoilMoistureQuality({
    landCover: landCover,
    vegetation: vegetation,
    retrieval: soilMoisture,
    RFI_check: radioFrequencyInterferenceCheck(Tb_h, Tb_v)
  });

  return {
    soil_moisture_percent: soilMoisture * 100,
    depth_cm: [0, 5],  // L-band penetration depth
    quality_flag: quality.flag,
    uncertainty_percent: quality.uncertainty,
    retrieval_algorithm: "tau-omega_model"
  };
}

// Accuracy target: 0.04 m³/m³ (4% volumetric)
// Validation: Against in-situ sensor networks

Quality Control and Validation

All WIA-compliant drought indices must undergo rigorous quality control:

Required Quality Checks

Quality Check	Method	Threshold	Action if Failed
Range validation	Check values within physical limits	Index-specific ranges	Flag as invalid
Spatial consistency	Compare with neighboring pixels	Z-score < 3.0	Flag for review
Temporal consistency	Compare with previous values	Change < 2 std dev	Flag rapid changes
Cross-validation	Compare multiple indices	Agreement > 80%	Lower confidence score
Ground truth comparison	Validate against field data	RMSE within targets	Adjust algorithms

Uncertainty Quantification

WIA-compliant systems must report uncertainty estimates:

function quantifyUncertainty(droughtIndex, sourceData, processingChain) {
  // Measurement uncertainty
  const sensorUncertainty = getSensorSpecifications(sourceData.satellite).uncertainty;

  // Processing uncertainty
  const atmosphericCorrectionError = 0.03;  // 3% reflectance
  const cloudMaskError = 0.02;  // 2% misclassification
  const algorithmError = getAlgorithmUncertainty(processingChain.algorithm);

  // Validation uncertainty
  const validationRMSE = getValidationStatistics(droughtIndex.type).rmse;

  // Combined uncertainty (root sum of squares)
  const totalUncertainty = Math.sqrt(
    Math.pow(sensorUncertainty, 2) +
    Math.pow(atmosphericCorrectionError, 2) +
    Math.pow(cloudMaskError, 2) +
    Math.pow(algorithmError, 2) +
    Math.pow(validationRMSE, 2)
  );

  return {
    uncertainty: totalUncertainty,
    components: {
      sensor: sensorUncertainty,
      processing: Math.sqrt(Math.pow(atmosphericCorrectionError, 2) +
                            Math.pow(cloudMaskError, 2)),
      algorithm: algorithmError,
      validation: validationRMSE
    },
    confidenceLevel: uncertaintyToConfidence(totalUncertainty)
  };
}

Data Fusion Protocols

Combining multiple data sources improves drought assessment accuracy:

Multi-Sensor Fusion

function fuseDroughtIndices(indices, weights, uncertainties) {
  // Weighted average based on inverse uncertainty
  const precisionWeights = uncertainties.map(u => 1 / (u * u));
  const totalWeight = precisionWeights.reduce((a, b) => a + b, 0);

  const fusedValue = indices.reduce((sum, index, i) => {
    return sum + (index * precisionWeights[i] / totalWeight);
  }, 0);

  // Fused uncertainty
  const fusedUncertainty = Math.sqrt(1 / totalWeight);

  return {
    fusedIndex: fusedValue,
    uncertainty: fusedUncertainty,
    sourceIndices: indices,
    weights: precisionWeights.map(w => w / totalWeight)
  };
}