conservation-python 教學

1這是什麼

Python 的空氣物理層其實不缺：psychrolib、CoolProp、 psychrochart、MetPy 都很成熟。缺的是保存科學那一層 —— 保存指數（preservation index）、壽命倍數（lifetime multiplier）、Zeng 黴菌等值線、 VTT 黴菌指數、木材平衡含水率，以及 Getty 那套「這筆讀數該加熱、加濕還是除濕」的 區間分類。這個專案補的就是這一塊。

它由兩個套件組成：

兩者都是 R 套件 ConSciR v0.3.0 的 Python 移植，GPL-3。函式名稱刻意跟 R 原版一致（calcDP、calcPI…）， 所以 R 那邊的文件和腳本可以直接對照過來。

2安裝

2.1 本機（Jupyter / 一般 Python）

兩個套件都還沒上 PyPI，所以直接從 GitHub 裝。先開一個虛擬環境 —— macOS 的 Homebrew Python 和新版 Linux 發行版會直接以 externally-managed-environment 拒絕你（PEP 668），不是壞掉，是保護機制。

python -m venv .venv
# Windows:  .venv\Scripts\activate
# macOS / Linux:
source .venv/bin/activate

pip install "git+https://github.com/Tai-ShengYeh/conservation-python.git#subdirectory=conscipy"
pip install matplotlib jupyterlab

macOS 上要用 python3，系統沒有 python 這個指令。

2.2 Google Colab

新開一個 notebook，第一格貼這個，然後 Shift+Enter：

!pip install -q "git+https://github.com/Tai-ShengYeh/conservation-python.git#subdirectory=conscipy"

matplotlib 和 pandas 是 Colab 內建的，不用裝。裝完直接 import conscipy as cp 就能用。

這份教學有一份對應的 notebook，可以直接開：

2.3 Orange Data Mining

Orange 本身從 orangedatamining.com/download 下載安裝。接著關鍵的一步是：附加元件必須裝進「正在跑 Orange 的那個 Python」， 不是你 $PATH 上那個。

最不會出錯的做法是用 Orange 自己的介面： Options ▸ Add-ons ▸ Add more…， 它一定裝在對的地方。但因為目前還沒上 PyPI，這條路要等發佈之後才有；現在請用終端機。

Windows（用 Orange 安裝目錄下的 python）：

pip install "git+https://github.com/Tai-ShengYeh/conservation-python.git#subdirectory=conscipy"
pip install "git+https://github.com/Tai-ShengYeh/conservation-python.git#subdirectory=orange3-conservation"

macOS，如果 Orange 是從 .dmg 裝的，它的 Python 在 app bundle 裡面：

/Applications/Orange3.app/Contents/MacOS/pip install "git+https://github.com/Tai-ShengYeh/conservation-python.git#subdirectory=conscipy"
/Applications/Orange3.app/Contents/MacOS/pip install "git+https://github.com/Tai-ShengYeh/conservation-python.git#subdirectory=orange3-conservation"

裝完重開 Orange。工具箱側邊欄會出現 Conservation Science 分類。

3第一個計算

先確認裝對了。這三個數字是從 ConSciR 的 README 抄過來的基準值， 移植的正確性就是拿它們比對的（第 21 節）。

import conscipy as cp

print("conscipy", cp.__version__)
print("calcDP(21.8, 36.8)     =", cp.calcDP(21.8, 36.8))
print("calcPI(21.8, 36.8)     =", cp.calcPI(21.8, 36.8))
print("calcRH_AH(23.8, 7.0524)=", cp.calcRH_AH(23.8, 7.0524))

實際輸出

conscipy 0.1.0
calcDP(21.8, 36.8)     = 6.383969538523943
calcPI(21.8, 36.8)     = 45.25849056382369
calcRH_AH(23.8, 7.0524)= 32.819637454815386

三個函式分別是：21.8 °C / 36.8 %RH 的露點是 6.38 °C； 同樣條件下的保存指數是 45.3 年；如果把這團空氣加熱到 23.8 °C （絕對濕度不變，7.0524 g/m³），相對濕度會掉到 32.8 %。

最後這個是保存工作最常問的問題之一 —— 「冬天開暖氣，濕度會掉到多少」。

4資料從哪來

先講結論：上游 R 套件那份 mydata 可以用， 這份教學從這裡開始就是用它。它是 ConSciR 的示範資料集，跟那個套件的其他部分一樣 是 GPL-3（DESCRIPTION 寫 License: GPL (>= 3)，資料沒有另外的授權聲明）， 附上出處即可使用。

但沒有必要把 3 MB 別人的資料複製進自己的 repo —— 直接從上游網址讀就好， Colab 上也一樣能跑：

URL = ("https://raw.githubusercontent.com/BhavShah01/ConSciR/"
       "master/inst/extdata/mydata.xlsx")

real = cp.tidy_TRHdata(URL, sheet="mydata")
print(real.head())

實際輸出

     Site    Sensor                    Date  Temp    RH
0  London  External 2024-01-01 00:00:00.000   7.4  79.1
1  London  External 2024-01-01 00:15:00.000   7.3  79.5
2  London  External 2024-01-01 00:29:59.990   7.2  79.2
3  London  External 2024-01-01 00:44:59.985   6.8  80.2
4  London  External 2024-01-01 00:59:59.980   6.8  80.8

shape:       (104826, 5)
site:        ['London']
sensors:     ['External', 'Room 1', 'Room 1 Case']
date range:  2024-01-01 00:00:00 -> 2024-12-31 23:44:59.610000
interval:    0 days 00:15:00

the file holds 105408 rows; 582 have no
reading at all, and tidy_TRHdata drops them, leaving 104826.

三個感測器，這才是重點

它有 室外、室內、展示櫃三個量測點，一整年 15 分鐘一筆。 這組三聯資料一張表就說完了保存工作最核心的一件事 —— 緩衝：

實際輸出

             Temp                RH            
              min  mean   max   min  mean   max
Sensor                                         
External      0.7  13.4  32.6  20.4  68.5  96.8
Room 1       16.0  20.9  26.2  20.5  46.1  76.0
Room 1 Case  16.3  21.7  25.9  27.0  42.2  57.9

annual range (max - min):
             Temp    RH
Sensor                 
External     31.9  76.4
Room 1       10.2  55.5
Room 1 Case   9.6  30.9

室外全年溫度跨 31.9 °C、濕度跨 76.4 個百分點；進到室內收斂成 10.2 和 55.5； 再進到展示櫃只剩 9.6 和 30.9。每一層外殼都把外界的變動吃掉一截 —— 這就是為什麼要用展示櫃，而數字自己會講。

那 example_data() 還留著做什麼

套件裡另外有一個合成資料產生器。它不是因為上游資料不能用才存在 —— 早期版本的說明這樣寫，那是錯的，已經更正。它留著是因為產生器能給固定資料集 給不了的東西：

df = cp.example_data(n_days=90)     # 參數：n_days, freq, seed

實際輸出

     Site  Sensor                Date  Temp    RH
0  London  Room 1 2024-01-01 00:00:00  12.1  57.5
1  London  Room 1 2024-01-01 00:15:00  12.0  52.3
2  London  Room 1 2024-01-01 00:30:00  12.2  56.3

shape: (17280, 5)  sensors: ['Room 1', 'Room 2']

a = cp.example_data(n_days=2, seed=0)
b = cp.example_data(n_days=2, seed=0)
c = cp.example_data(n_days=2, seed=1)
print("seed 0 twice identical :", a.equals(b))
print("seed 0 vs seed 1 equal :", a.equals(c))

實際輸出

seed 0 twice identical : True
seed 0 vs seed 1 equal : False

5讀自己的紀錄器檔案

真的要做事的時候，資料來自 Meaco、Hanwell 或其他紀錄器的匯出檔。 tidy_TRHdata() 會把常見的欄位命名變體對齊成標準欄位 （Site / Sensor / Date / Temp / RH），也接受檔案路徑（CSV 或 Excel）。

raw = pd.DataFrame({
    "Date/Time":       ["2024-03-01 09:00", "2024-03-01 09:07", "2024-03-01 09:15"],
    "Temperature (C)": [18.4, 18.5, 18.6],
    "RH (%)":          [52.1, 51.8, 51.5],
    "Room":            ["Gallery 1"] * 3,
})
tidy = cp.tidy_TRHdata(raw, round_to="15min")   # 也可以 cp.tidy_TRHdata("logger.csv")

實際輸出

raw columns: ['Date/Time', 'Temperature (C)', 'RH (%)', 'Room']

   Site     Sensor                Date  Temp    RH
0  <NA>  Gallery 1 2024-03-01 09:00:00  18.4  52.1
1  <NA>  Gallery 1 2024-03-01 09:00:00  18.5  51.8
2  <NA>  Gallery 1 2024-03-01 09:15:00  18.6  51.5

dtypes:
Site              object
Sensor            object
Date      datetime64[ns]
Temp             float64
RH               float64
dtype: object

注意 Room 被認成 Sensor、Date/Time 被認成 Date 並轉成 datetime，而 round_to="15min" 把 09:07 這筆 對齊到 09:15。三個專用讀取器：

真的廠商匯出檔長什麼樣

上游那個 Excel 檔除了 mydata 之外，還附了 Meaco、 Hanwell、VAISALA 三張工作表 —— 正是這些讀取器存在的理由。 拿真檔案來跑：

cp.tidy_Meaco(URL, sheet="Meaco")
cp.tidy_Hanwell(URL, sheet="Hanwell")

實際輸出

--- tidy_Meaco: 10 readings ---
   Site Sensor                Date   Temp     RH
0  York  Store 2024-01-01 00:03:00  18.17  54.61
1  York  Store 2024-01-01 00:13:00  18.16  54.76
  site ['York']  sensor ['Store']
  2024-01-01 00:03:00 -> 2024-01-01 01:38:00

--- tidy_Hanwell: 7075 readings ---
   Site    Sensor                Date  Temp    RH
0  None  External 2025-08-12 00:06:03  22.8  60.2
1  None  External 2025-08-12 00:11:01  22.6  60.2
  site []  sensor ['External']
  2025-08-12 00:06:03 -> 2025-09-10 10:30:19

兩種格式差很多。Meaco 是長格式，每筆讀數一列，欄位叫 RECEIVER / TRANSMITTER / DATE / TEMPERATURE / HUMIDITY，讀取器只是改名。 Hanwell 則有一段檔頭中繼資料、感測器名稱單獨一行、資料區塊， 最後還有一段統計摘要和 ---- END OF DATA ---- 標記。

6單點與向量計算

所有 calc* 函式都是 NumPy 向量化的：給純量回純量， 給陣列或 Series 回同樣形狀的結果。單位一律沿用 R 原版 —— 溫度 °C、相對濕度 %（0–100）、氣壓 hPa。

import numpy as np

temps = np.array([15.0, 20.0, 25.0, 30.0])
print("calcPws    :", np.round(cp.calcPws(temps), 4))
print("calcDP @55%:", np.round(cp.calcDP(temps, 55.0), 4))

實際輸出

input Temp : [15. 20. 25. 30.]
calcPws    : [17.0517 23.3834 31.6853 42.4513]
calcDP @55%: [ 6.0302 10.6855 15.3345 19.9773]

7四種飽和水氣壓公式

calcPws 提供四種算法。要選哪個？在保存工作的溫度範圍內， 差異小到無關緊要 —— 但知道差多少，比不知道好：

for m in ("Buck", "IAPWS", "Magnus", "VAISALA"):
    print(f"{m:<8} {cp.calcPws(20.0, method=m):.4f} hPa")

實際輸出

Buck     23.3834 hPa
IAPWS    23.3919 hPa
Magnus   23.3344 hPa
VAISALA  23.3925 hPa
spread   0.0581 hPa (0.248 % of the mean)

四者相差約 0.25 %，遠小於任何紀錄器的量測不確定度 （一般是 ±2–3 %RH）。預設的 Buck 在 0 °C 以下會自動切換到冰面係數， 一般情況照用即可。

8反函數

保存工作真正要問的通常是反過來的問題：要把濕度降到 50 %， 空氣要加熱到幾度？ 三個反函式負責這件事，而且它們能繞回去 —— 這本身就是一種驗證。

T, RH = 21.8, 36.8
dp = cp.calcDP(T, RH)
ah = cp.calcAH(T, RH)
print(f"  RH back from DP  {cp.calcRH_DP(T, dp):.6f} %")
print(f"  RH back from AH  {cp.calcRH_AH(T, ah):.6f} %")
print(f"  T needed for {RH} % at that DP: {cp.calcTemp(RH, dp):.6f} degC")
print(f"  same, Buck bisection:           {cp.calcTemp(RH, dp, method='Buck'):.6f} degC")

實際輸出

T=21.8 degC, RH=36.8 %
  dew point           6.383970 degC
  RH back from DP     36.800000 %   (should be 36.8)
  absolute humidity   7.052415 g/m3
  RH back from AH     36.800000 %   (should be 36.8)
  T needed for 36.8 % at that DP: 21.800000 degC
  same, Buck bisection:           21.794722 degC

calcTemp 有兩種解法。Magnus 有解析解，直接算。 Buck 沒有，R 版用純量的 uniroot() 逐點解； 這裡改成向量化二分法，所以可以直接吃整個陣列。 兩種方法在上面的例子差在小數第四位，是二分法容差造成的。

9濕度欄位

實務上不會一個一個算，而是給一張表、長出一堆欄位。 add_* 系列做的就是這件事。

out = cp.add_humidity_calcs(df.head(3))
print(out.round(3).T)

實際輸出（轉置後比較好讀）

                            0                    1                           2
Site                   London               London                      London
Sensor                 Room 1               Room 1                      Room 1
Date      2024-01-01 00:00:00  2024-01-01 00:15:00  2024-01-01 00:29:59.990000
Temp                     21.8                 21.8                        21.8
RH                       36.8                 36.7                        36.6
Pws                    26.121               26.121                      26.121
Pw                      9.613                9.586                        9.56
DP                      6.384                6.344                       6.305
AH                      7.052                7.033                       7.014
AD                      1.192                1.192                       1.192
MR                      5.957                5.941                       5.925
SH                      0.856                0.856                       0.856
Enthalpy               37.157               37.115                      37.074

一次補上 Pws、Pw、DP、AH、AD、MR、SH、Enthalpy 八欄。 可以用 method= 換飽和水氣壓公式、P_atm= 換氣壓 （高海拔場館會需要）。

10保存風險欄位

year = cp.example_data(n_days=365)
out = cp.add_conservation_calcs(year)
print(out[["Temp","RH","Mould_LIM","Mould_risk","Mould_rate",
           "Mould_index","PreservationIndex","Lifetime","EMC_wood"]].head(3).round(3).T)

實際輸出

                          0         1         2
Sensor             External  External  External
Temp                    7.4       7.3       7.2
RH                     79.1      79.5      79.2
Mould_LIM             84.52    84.649    84.779
Mould_risk          No risk   No risk   No risk
Mould_rate              0.0       0.0       0.0
Mould_index             0.0       0.0       0.0
PreservationIndex   139.388   140.614    143.11
Lifetime              0.857     0.855     0.856
EMC_wood             16.128    16.269    16.167

readings above the Zeng isopleth, per sensor:
  External       8497 / 34825  (24.4 %)
  Room 1            4 / 35127  (0.0 %)
  Room 1 Case       0 / 34874  (0.0 %)

peak VTT mould index, per sensor:
Sensor
External       1.4668
Room 1         0.0052
Room 1 Case    0.0015

三個感測器的結果差得很開，而且方向合理：室外有將近四分之一的讀數超過 Zeng 等值線、VTT 指數爬到 1.47；室內只有 4 筆越線；展示櫃一筆都沒有。 這是一份表現正常的建築該有的樣子。

11目標區間與調整量

這是 Getty《Tools for the Analysis of Collection Environments》那套框架的核心： 給一個目標溫濕度區間，把每一筆讀數分類，並算出要花多少力氣把它拉回區間內。

out = cp.add_humidity_adjustments(year, LowT=16, HighT=25, LowRH=40, HighRH=60)
print(out["zone"].value_counts())
print((out["zone"].value_counts() / 4).round(1))   # 讀數間隔 15 分鐘 → 換算成小時

實際輸出

Room 1, hours per zone over the year:
zone
Within               5502
Hum or cooling       2190
Dehum or heating      552
Hum only              392
Dehum only             70
Cooling and dehum      54
Cooling only           21

first three readings needing action:
   Temp    RH            zone  dTemp_TRHadj  dRH_TRHadj
0  21.8  36.8  Hum or cooling           0.0         3.2
1  21.8  36.7  Hum or cooling           0.0         3.3
2  21.8  36.6  Hum or cooling           0.0         3.4

zone 有十一種可能值，說的是「該做什麼」：Within、 Heating only、Dehum or heating、Dehum only、 Hum only、Cooling and dehum、Heating and hum、 Heating and dehum、Hum or cooling、Cooling only、 Cooling and hum。

「or」那幾種（Dehum or heating、Hum or cooling）是有意義的 —— 它們代表兩種手段都可以達成目標，你可以挑比較省能源的那個。 這也是為什麼函式會同時算出三組調整量：

把 dTemp_TRHadj 和 dAH_TRHadj 接到 第 17 節的 HVAC 函式，就能把「分區」變成「度電」。

12時間變數

out = cp.add_time_vars(df.head(2))
print(out[["Date","day","hour","weekday","month","year","Summer"]])

實際輸出

                 Date        day  hour weekday month  year  Summer
0 2024-01-01 00:00:00 2024-01-01     0     Mon   Jan  2024  Winter
1 2024-01-01 00:15:00 2024-01-01     0     Mon   Jan  2024  Winter

13繪圖

三個繪圖函式，都回傳 matplotlib 的 Figure， 所以你可以繼續改它、或存成檔案。需要 pip install matplotlib。

fig = cp.graph_TRH(year)                    # 溫濕度時間序列，依 Sensor 分面
fig = cp.graph_psychrometric(year)          # 心理量圖
fig = cp.graph_TRHbivariate(year)           # 溫濕度二維密度
fig.savefig("chart.png", dpi=150)

心理量圖的縱軸可以換成任何一個計算量，用 y_func=：

fig = cp.graph_psychrometric(year, y_func="calcPI")   # 縱軸改成保存指數

14兩個黴菌模型

套件裡有兩套彼此獨立的黴菌模型，回答的問題不一樣。

Zeng 等值線 —— 「現在這一刻危險嗎」

calcMould_Zeng 是無記憶的：只看當下的溫濕度， 回答「這個溫度下，相對濕度超過多少會開始長黴」。

for lim in (0, 0.1, 0.5, 1, 2, 3, 4):
    print(f"{lim} mm/day  20 degC: {cp.calcMould_Zeng(20.0, 50.0, LIM=lim):.2f}")

# label=True 反過來：給實際讀數，回推生長速率
cp.calcMould_Zeng(25.0, 90.0, label=True)

實際輸出

critical RH (%) at which each growth rate starts:
     0 mm/day  20 degC:  75.62   25 degC:  74.50
   0.1 mm/day  20 degC:  77.09   25 degC:  75.97
   0.5 mm/day  20 degC:  80.84   25 degC:  79.72
     1 mm/day  20 degC:  83.59   25 degC:  82.47
     2 mm/day  20 degC:  86.59   25 degC:  85.47
     3 mm/day  20 degC:  89.11   25 degC:  87.99
     4 mm/day  20 degC:  91.10   25 degC:  89.98

label=True gives the rate implied by an actual reading:
  20.0 degC / 70.0 %RH -> 0.0 mm/day
  20.0 degC / 85.0 %RH -> 1.0 mm/day
  25.0 degC / 90.0 %RH -> 5.0 mm/day
  28.0 degC / 95.0 %RH -> 5.0 mm/day

VTT 指數 —— 「累積下來到什麼程度了」

VTT（Viitanen）模型是遞迴的：每一步的黴菌指數取決於前一步。 指數從 0 跑到 6（0 = 沒有生長，6 = 完全覆蓋），濕的時候往上爬， 乾的時候會回落。這是它跟 Zeng 最大的差別 —— 它有記憶。

合成資料在這裡還有用嗎？有，但角色不同

真實資料告訴你這棟建築做了什麼；它不會告訴你這棟建築 可能變成什麼樣。要問「如果除濕機壞掉三個月會怎樣」， 只能自己造出那個情境 —— 從同一組真實讀數衍生就行，不必另外發明產生器：

damp = room.copy()
damp["RH"] = np.clip(damp["RH"] + 25.0, 5, 99)     # 同一個房間，濕 25 個百分點

實際輸出

Real data covers what the building did, not what it might do.
A generator reaches the rest -- here the same room, 25 %RH damper:

  Room 1 as measured     RH max  76.0 %   at risk      4   peak index 0.0052
  Room 1 + 25 %RH        RH max  99.0 %   at risk  11415   peak index 1.3293

同一個房間、同一年的溫度，只把濕度抬高 25 個百分點，越線讀數就從 4 筆變成 11,415 筆，VTT 指數從 0.005 升到 1.33。這是真實資料做不到的事 —— 它只有一種結果。

四個材料敏感度等級

dt = cp.dt_days_from_times(outside["Date"])
for s in cp.VTT_SENSITIVITY:                 # very / sensitive / medium / resistant
    m = cp.mould_index_series(outside["Temp"], outside["RH"],
                              dt_days=dt, sensitivity=s)
    print(f"{s:<10} -> {m[-1]:.4f}")

實際輸出

final VTT index after a year outdoors:
  very       (A=1.0, B=7.0, C=2.0)  ->  1.4668
  sensitive  (A=0.3, B=6.0, C=1.0)  ->  1.4646
  medium     (A=0.0, B=5.0, C=1.5)  ->  1.4578
  resistant  (A=0.0, B=3.0, C=1.0)  ->  1.3733

the same four classes inside Room 1:
  very       ->  0.0052
  sensitive  ->  0.0052
  medium     ->  0.0052
  resistant  ->  0.0052

15時間、遞迴與順序

VTT 是遞迴，這帶來三個要注意的地方，其中第一個最容易致命。

15.1 每筆讀數代表多少時間

模型的生長速率是每天的 —— dM_dt 裡那個 7 * 就是把 Viitanen 的反應時間從「週」換算成「天」。 所以遞迴每走一步，必須知道那一步代表多久。15 分鐘一筆的資料， 每步是 1/96 天，不是一天。

dt = cp.dt_days_from_times(outside["Date"])       # 每筆讀數代表幾天
m  = cp.mould_index_series(outside["Temp"], outside["RH"], dt_days=dt)

實際輸出

External, 34825 readings over 365 days
  median step   15 minutes
  largest gap   49.8 hours

final mould index, real elapsed time   1.4668
final mould index, one day per reading 6.0000   <- saturated

resampled hourly, same year            1.4701

差別是決定性的：正確計時得到 1.47，把每筆讀數都當成一天則直接撞上限 6.0。而把同一年重新取樣成每小時一筆，答案是 1.4701 —— 幾乎不變。 這正是重點：結果應該取決於條件持續了多久，而不是紀錄器多久記一次。

15.2 上游的 element-wise 問題

R 版的 calcMould_VTT() 用 mapply 逐列套用， 而 M_prev 預設是 0 —— 也就是每一列都假設前一刻的黴菌指數是零。 那只有在「只算一個時間步」的時候才正確。整條序列這樣算，結果會嚴重低估。

15.3 順序

實際輸出

readings                   34825
final index, time order    1.466842
final index, shuffled      1.466282
difference                 0.038 %

Small, but not zero, and only small while growth stays well below
M_max. Once k2 starts biting, when the index rises, order matters:
  as measured  peak 1.4668   readings above 1.0: 5387
  shuffled     peak 1.4663   readings above 1.0: 10952

upstream element-wise (M_prev = 0 every row):
  max 0.034676
cumulative recursion:
  max 1.466842

把資料打亂再跑，答案幾乎一樣但不完全一樣。 早期版本的這份教學說它「完全順序無關」，那是因為當時用的示範資料指數一直很低， k2 幾乎恆等於 1，每一步的增量就與 M 無關。 一旦指數爬起來、k2 開始起作用，順序就會有影響。

實務結論不變，而且更強：排序是呼叫端的責任。 Orange 元件會自己依時間欄排序，函式庫不會。

16分組陷阱

whole = cp.add_conservation_calcs(real)                    # 三個感測器一起
for name, g in real.groupby("Sensor"):
    print(name, cp.add_conservation_calcs(g)["Mould_index"].max())

實際輸出

mydata holds three sensors, one after another:
  ['External', 'Room 1', 'Room 1 Case']

add_conservation_calcs(real)      peak index 1.4736
  ...on 'External' alone           peak index 1.4668
  ...on 'Room 1' alone           peak index 0.0052
  ...on 'Room 1 Case' alone           peak index 0.0015

Split by sensor, and sort, before the recursion runs:
  grouped                         peak index 1.4668

每個感測器都從前一個結束的地方接著長，而不是從零開始。 更糟的是接縫處時間會倒退一年，那一步被算成沒有經過時間 —— 所以錯誤的方式不只是數字偏高，而是根本沒有意義。正確做法是先排序、再分組：

fixed = pd.concat(
    [cp.add_conservation_calcs(g)
     for _, g in real.sort_values(["Sensor", "Date"]).groupby("Sensor")],
    ignore_index=True)

用 concat 而不是 groupby(...).apply(...) 是刻意的 —— 後者需要 include_groups= 參數，而那個參數在 pandas 2.2 已經標記棄用。

Orange 的 Conservation Risk 元件有處理這件事 —— 它有 Time column 和 Group by 兩個設定，會先排序再分組跑， 然後把結果寫回原本的列位置。這是目前函式庫層和元件層行為不一致的地方； 用函式庫時要自己補上。

17HVAC 負荷

把「需要多少調整」換算成「需要多少能量」。

cp.calcCoolingPower(28, 21, 70, 50, 2.0)      # 28°C/70% → 21°C/50%，2 m³/s
cp.calcSensibleHeating(21, 28, 2.0)           # 顯熱負荷
cp.calcTotalHeating(21, 28, 50, 70, 2.0)      # 全熱負荷
cp.calcSensibleHeatRatio(21, 28, 50, 70, 2.0) # 顯熱比 (%)
cp.calcCoolingCapacity(3000)                  # 移除 3 kW 電力負載所需冷氣容量
cp.calcOffCoilDewPoint(26, 6, 12)             # 出風溫度

實際輸出

Cool 28 degC/70 %RH outside air to 21 degC/50 %RH at 2 m3/s:
  cooling power     0.0697 kW
  sensible heating  16.8223 kW
  total heating     71.7508 kW
  sensible heat ratio 23.45 %

Cooling capacity for a 3 kW lighting load: 4.3714 kW
Off-coil temperature, 26 degC on-coil, 6/12 degC chilled water: 10.7000 degC

18五個元件

裝好之後，工具箱側邊欄會多一個 Conservation Science 分類。 以下截圖都是真的把元件建起來、送 30 天示範資料進去之後截下來的。

為什麼結果欄位分成兩種

元件有一個刻意的設計：連續型結果放進 X（特徵）， 類別型結果放進 metas。這樣新增的數值欄位會自動被 PCA、PLS、迴歸吃進去；而分區和風險標籤留在 metas，可以拿來著色或分群， 但不會扭曲距離與投影。

另外，每個 DiscreteVariable 都是用完整的可能值清單建立的， 不是只用當批出現過的值。所以兩個檔案分開處理之後仍然共用同一套編碼， 可以直接 concatenate。

19組成工作流

典型的一條線是這樣：

Logger File ──▶ Humidity Calculations ──▶ Conservation Risk ──▶ Envelope Zones
                                                    │
                                                    ▼
                                         Psychrometric Chart ──▶ Data Table

因為輸出就是一般的 Orange Table，可以直接接到畫布上任何其他元件 —— Data Table 看數字、Distributions 看分布、Scatter Plot、PCA、 Hierarchical Clustering 都行。這是用 Orange 而不是自己寫腳本的主要理由。

套件附了一份現成的工作流。裝好之後從 Help ▸ Example Workflows ▸ Conservation Science Examples ▸ collection-environment 開啟， 在 Logger File 裡勾 Use synthetic demo data，整條線就會跑起來， 不需要任何自己的檔案。

20模擬資料證明什麼、不證明什麼

回到你最初的問題。答案是可以，但要把兩件事分清楚：

要驗證公式，需要另外兩種手段：對照已知數值（第 21 節） 和內部一致性（第 22 節）。三者合起來才算完整。

21對照上游數值

移植的正確性靠的是把 Python 的結果跟 ConSciR README 印出來的數字對比， 相對容差 1e-6。這就是 conscipy/tests/test_parity.py 在做的事：

實際輸出

call                            conscipy  ConSciR README    rel. diff
calcDP(21.8, 36.8)               6.38397         6.38397     7.23e-08
calcAH(21.8, 36.8)               7.05242         7.05242     6.72e-07
calcPI(21.8, 36.8)              45.25849        45.25850     2.08e-07
calcRH_AH(23.8, 7.0524)         32.81964        32.81970     1.91e-06

四項的相對差都在 1e-6 到 1e-8 之間 —— 差距來自 README 只印了六位有效數字， 不是計算不一致。

python -m pytest conscipy/tests -q

整套跑起來是 19 個測試。Orange 元件那邊另外有 39 個， 用的是 Orange 自己的 WidgetTest 框架，所以測到的是真的訊號傳遞， 不只是計算函式。

其中 conscipy/tests/test_upstream_quirks.py 這三個比較特別： 它們釘住的是下一節那兩個問題與分組陷阱的 目前行為。這樣做的用意是讓文件和程式一起壞掉 —— 哪天有人把 calcCoolingPower 改成物理上正確的版本， 測試會失敗並指向這份教學該更新的段落，而不是讓網頁默默描述一個已經不存在的行為。

22內部一致性

第三種手段：不需要任何外部基準，只需要物理上必然成立的關係。 在合成資料上跑一遍，任何一條不成立就代表有 bug。

df = cp.add_humidity_calcs(cp.example_data(n_days=30))
T, R = df["Temp"].to_numpy(), df["RH"].to_numpy()

assert (df["DP"] <= T).all()                                    # 露點不可能高於乾球
assert (df["Pw"] <= df["Pws"]).all()                            # 分壓不可能高於飽和壓
assert np.allclose(cp.calcRH_DP(T, df["DP"]), R, atol=1e-6)     # 反函數繞得回來
assert np.allclose(cp.calcRH_AH(T, df["AH"]), R, atol=1e-3)
assert (df["SH"] < df["MR"]).all()                              # 比濕必小於混合比
assert df["AD"].between(1.0, 1.3).all()                         # 空氣密度的物理範圍

實際輸出

rows checked: 35127
dew point never above dry bulb      : True
Pw never above Pws                  : True
RH recovered from DP within 1e-6    : True
RH recovered from AH within 1e-3    : True
specific humidity < mixing ratio    : True
air density in 1.0-1.3 kg/m3        : True  (range 1.1701-1.2173)

這種檢查特別值得寫成測試，因為它不需要你去找基準值 —— 而且會抓到「數值看起來很合理但其實錯了」的那種 bug。 下一節那五個發現，多半就是這樣冒出來的。

但也要看清楚它的極限。上面 SH < MR 那一條是通過的， 而 calcSH 的值其實錯了 8 倍 —— 只檢查大小關係的斷言抓不到尺度錯誤。 細節在第 23.4 節。

23五個發現

照著上面那套方法做這份教學的時候，跑出五個彼此矛盾、或與文獻不符的地方。 四個追到上游 R 套件，conscipy 是忠實移植（NOTICE 寫明公式沿用原版）， 所以它們不是移植錯誤；一個是這份移植自己的，已經修掉。

23.1 calcCoolingPower 的單位

第 17 節那兩個數字差了 1000 倍。手算一次就看得出來為什麼：

實際輸出

air density rho          1.1605 kg/m3
enthalpy in  h1          70.8741 kJ/kg
enthalpy out h2          40.8384 kJ/kg
rho * V * (h1 - h2)      69.7144 kW   <- kg/s * kJ/kg = kW

calcCoolingPower(...)    0.0697 kW
calcTotalHeating(...)    71.7508 kW

ratio between the two    1029.2

calcEnthalpy 回傳的是 kJ/kg。質量流率是 kg/s。 兩者相乘已經是 kW，不需要再除以 1000。ConSciR 的 R 原始碼是這樣寫的：

coolingPowerWatts <- massFlowRate * (h1 - h2)
coolingPowerKW    <- coolingPowerWatts / 1000

變數名 coolingPowerWatts 洩漏了假設 —— 以為 h 的單位是 J/kg。 同一個檔案裡的 calcTotalHeating 就沒有那個除法。

實務建議：要冷卻負荷的時候用 abs(calcTotalHeating(...))。

23.2 VTT 臨界相對濕度的正負號

ConSciR 與 conscipy 寫的是 -0.00267·T³ + 0.160·T² + 3.13·T + 100， 而 Hukka & Viitanen (1999) 原始論文是 − 3.13·T：

實際輸出

 T (degC)   ConSciR / conscipy   Hukka & Viitanen
      5.0               119.32              88.02
     10.0               144.63              82.03
     15.0               173.94              80.04
     19.9               204.61              80.03
     20.0               205.24              80.04
     20.1                80.00              80.00
     25.0                80.00              80.00

A critical RH above 100 % can never be exceeded by a real reading,
so frac = (RHcrit - RH)/(RHcrit - 100) turns positive and M_max grows.
That is why a 12 degC / 53 %RH room accumulates a VTT index at all,
while the Zeng isopleth on the same readings reports no risk.

看 T = 20.0 和 20.1 那兩列。模型在 20 °C 以上固定用 80 %， 以下用多項式。正確的係數會讓多項式在 20 °C 剛好收斂到 80.04，跟另一段接得上；用 +3.13T 則是 205.24 掉到 80，中間裂了 125 個百分點。

後果很具體：RH_crit 大於 100 % 時， frac = (RH_crit − RH) / (RH_crit − 100) 變成正的大數， M_max 跟著變大，於是模型在又冷又乾的空氣裡預測出黴菌生長。 20 °C 以上不受影響，Zeng 等值線是獨立實作也沒有這個問題。

23.3 calcLM：三個錯誤，其中兩個互相遮蔽

壽命倍數應該說的是「這個環境下東西能撐多久，相對於 20 °C / 50 %RH」。 Michalski 的經驗法則是降 5 °C 大約可以讓壽命加倍。 實際跑出來是這樣：

實際輸出

Lifetime multiplier, as ConSciR and conscipy compute it:
   5.0 degC / 50 %RH  ->  0.997790
  10.0 degC / 50 %RH  ->  0.998552
  15.0 degC / 50 %RH  ->  0.999288
  20.0 degC / 50 %RH  ->  1.000000
  25.0 degC / 50 %RH  ->  1.000688
  30.0 degC / 50 %RH  ->  1.001354

A 5 degC drop should roughly double the lifetime (Michalski);
here 20 -> 15 degC changes it by -0.07 %.

                                      25C/50%   15C/50%   20C/30%
as written      EA=100 J   -  1/3      1.0007    0.9993    1.1856
units fixed     EA=100 kJ  -  1/3      1.9899    0.4907    1.1856
units and sign  EA=100 kJ  +  1/3      0.5025    2.0380    1.1856
published form  EA=100 kJ  +  1.3      0.5025    2.0380    1.9427

降 5 °C 只讓壽命變了 −0.07 %，而且方向還是反的 —— 溫度愈高壽命愈長。原因是三個錯誤疊在一起：

• 單位：EA=100 被當成 J/mol，但文獻是 100 kJ/mol，差 1000 倍。指數項因此塌成 1.000x， 溫度幾乎沒有作用。
• 正負號：程式寫的是 exp(−EA/R · (1/T − 1/T₀))， 文獻是 +。
• 相對濕度指數：程式是 (50/RH)^(1/3)， 文獻是 ^1.3。

在上游修正之前，不要把 calcLM 的輸出當成可用的保存年限指標。 calcPI 是另一套獨立的公式（IPI 的 Arrhenius 式），不受影響。

23.4 calcSH 的單位換算

實際輸出

calcMR(20.0, 50.0)              7.260814 g/kg dry air
calcSH(20.0, 50.0)              0.878947   <- MR / (1 + MR)

q = w / (1 + w) holds for w in kg/kg. With w in g/kg:
  w  = 7.260814 g/kg = 0.007260814 kg/kg
  q  = 0.007208475 kg/kg = 7.208475 g/kg

So the label says g/kg but the value is 8.2x too small.
The ordering check SH < MR still passes: True  -- which is why it never caught this.

比濕與混合比的關係 q = w / (1 + w) 成立的前提是 w 的單位是 kg/kg。但 calcMR 回傳的是 g/kg， 直接代進去就少了 1000 這個尺度。結果標成 g/kg，數值卻小了 8.2 倍。

影響範圍比函式庫大：add_humidity_calcs() 每次都會寫出 SH 欄，Orange 的 Humidity Calculations 元件可以勾選它， Psychrometric Chart 還可以拿它當縱軸 —— 也就是會畫出一條標著 g/kg 但小了 8 倍的曲線。

23.5 缺溫度被當成「無風險」

實際輸出

A reading with no temperature, at 90 %RH:
  calcMould_Zeng(nan, 90)             nan
  calcMould_Zeng(nan, 90, label=True) 0.0   <- reads as 'no growth'
  calcMould_Zeng(20, nan, label=True) nan

mydata has 582 rows with no reading. tidy_TRHdata drops those,
but a logger that reports humidity while its thermistor is faulty
would not be dropped, and would be scored as safe.

calcMould_Zeng(..., label=True) 只檢查相對濕度是不是 NaN， 沒有檢查溫度。溫度缺失時所有比較都是 False，結果就停在 0.0 —— 而 0.0 的意思是「沒有生長」。接下來 DataFrame 和 Orange 元件會把它編成 No risk，跟一筆真正安全的讀數完全無法區分。

正確的行為應該是：溫度或濕度任一缺失，數值結果為 NaN、類別結果為 missing， 而不是 No risk。tidy_TRHdata() 會丟掉兩者都缺的列， 所以只有「濕度正常但溫度感測器壞掉」這種情況會踩到 —— 但那正是最需要被標記出來的情況。

24參考資料

• 原始碼與 issue：github.com/Tai-ShengYeh/conservation-python
• 上游 R 套件：ConSciR · 說明文件
• Orange Data Mining：orangedatamining.com
• Cosaert, A., Beltran, V., et al. (2022). Tools for the Analysis of Collection Environments. Getty Conservation Institute.
• Hukka, A. & Viitanen, H. (1999). A mathematical model of mould growth on wooden material. Wood Science and Technology, 33, 475–485. doi:10.1007/s002260050131

引用時請一併引用上游套件與它實作的框架：

Shah, B., Cosaert, A., Beltran, V., et al. ConSciR: Tools for Conservation
Science. R package. https://bhavshah01.github.io/ConSciR/