Python计算方程出现RuntimeWarning(double_scalars)错误的解决咨询

运行以下计算ED的代码时出现警告：

T = 15
H = range(0, 100)
for hum in H:
  ED = 0.942 * (hum ** 0.679) + 11 * exp((hum - 100) / 10) + 0.18 * (21.1 - T) ** (1 - exp(-0.115 * hum))

报错信息：

RuntimeWarning: invalid value encountered in double_scalars
  ED = 0.942 * (hum ** 0.679) + 11 * exp((hum - 100) / 10) + 0.18 * (21.1 - T) ** (1 - exp(-0.115 * hum))

该方程在计算器或Excel中可正常运行，推测错误源于复杂的双精度标量运算部分。以下是完整代码片段（该方程为大型Python代码的一部分）：

import numpy as np
from netCDF4 import Dataset
import os
from numpy.ma.core import MaskedArray
import gdal
from math import exp
from math import log

nlat = 600
nlon = 465

years = range(1996,2005)
for year in years:
    starting_FFMC_input = np.zeros(shape = (nlat, nlon))

    #Temperature at 2m above ground
    temperature_at_2m_raw: MaskedArray()
    temperature_at_2m_input: MaskedArray()
    with Dataset(dir + "WRFoutputTemperature.nc") as file_temperature_at_2m:
        temperature_at_2m_raw = file_temperature_at_2m.variables['T2MEAN'][number_of_days_in_timeperiod * 8 : (number_of_days_in_timeperiod + days_in_year) * 8 - 1, :, :] - 273.15
        dimsizes_temperature_at_2m = temperature_at_2m_raw.shape
        temperature_at_2m_input = temperature_at_2m_raw[4 : (dimsizes_temperature_at_2m[0]-4) : 8, :, :] #3hourly, source data starts from 01-12-1995
    del temperature_at_2m_raw


    #Relative humidity at 2m above ground
    relative_humidity_at_2m_raw: MaskedArray()
    relative_humidity_at_2m_input: MaskedArray()
    with Dataset(dir + 'WRFoutputHumidity') as file_relative_humidity_at_2m:
        relative_humidity_at_2m_raw = file_relative_humidity_at_2m.variables['HURS'][number_of_days_in_timeperiod * 8 : (number_of_days_in_timeperiod + days_in_year) * 8 - 1, :, :] * 100
        dimsizes_relative_humidity_at_2m = relative_humidity_at_2m_raw.shape
        relative_humidity_at_2m_input = relative_humidity_at_2m_raw[4 : (dimsizes_relative_humidity_at_2m[0]-4) : 8, :, :] #3hourly, source data starts from 01-12-1995
    del relative_humidity_at_2m_raw

    index_i = range(0, 13132)
    for i in index_i:
        index_x = coords[i,0]
        index_y = coords[i,1]

        days = range(0, days_in_year-1)
        for day in days:
            temperature = temperature_at_2m_input[day, index_y, index_x]
            humidities = relative_humidity_at_2m_input[day, index_y, index_x]

            if humidities > 100:
                humidities = 100
            elif humidities < 0:
                humidities = 0
            else:
                humidities = humidities

            ED = 0.942 * (humidities ** 0.679) + 11 * exp((humidities - 100) / 10) + 0.18 * ((21.1 - temperature) ** (1 - exp(-0.115 * humidities))

温度和湿度数据示例：

Temperature is: 12.101776
Humidity is: 75.87772
Temperature is: 11.955383
Humidity is: 90.84902
Temperature is: 11.471436
Humidity is: 99.00755
Temperature is: 13.040009
Humidity is: 95.4102

错误原因分析 #

1. 负数的非整数次幂无意义：当21.1 - temperature为负数（即温度高于21.1℃），且指数1 - exp(-0.115 * humidities)为非整数时，实数范围内无法计算负数的非整数次幂，Python会返回无效值并触发警告。Excel可能通过自动取绝对值等方式规避了这个问题。
2. 语法错误：原代码中ED的计算表达式末尾缺失一个闭合括号，会导致语法报错。
3. 函数适配问题：使用math.exp处理numpy的MaskedArray数组，可能存在标量与数组运算的适配问题，建议改用numpy的矢量化函数。

改写后的代码 #

针对上述问题，修改ED计算部分的代码如下：

# 替换math.exp为np.exp，适配数组运算
exponent = 1 - np.exp(-0.115 * humidities)
base = 21.1 - temperature

# 处理负数底数的情况，确保计算结果为有效实数（参考Excel逻辑取绝对值）
term_3 = 0.18 * np.where(base >= 0, base ** exponent, np.abs(base) ** exponent)

# 计算ED，修复语法错误
ED = 0.942 * (humidities ** 0.679) + 11 * np.exp((humidities - 100) / 10) + term_3

额外说明 #

• 如果ED的物理定义中，当温度高于21.1℃时该部分有特殊计算逻辑，请根据实际公式调整np.where中的处理方式，而非直接取绝对值。
• 可以提前检查温度数据中是否存在高于21.1℃的记录，这是触发警告的直接场景。

内容的提问来源于stack exchange，提问作者Elsri