模式學習日誌

發表文章

GLYCIM模式待討論事項

1. 株高：印象中這邊是原始碼有點問題 2. 模擬的莢數 >> 實際莢數： 3. 模擬的分支數 >> 實際分支數： rosetta 程式需要有visual studio才能開啟(還要有C++套件)，有可能是當初building的時候有問題 1. 株高 GLYCIM程式中，與株高相關的參數為PARM(21)與PARM(22)，模式假設株高是Vstage的指數關係，也就是株高 = Parm(21)*VSTAGE^PARM(22) 為了表示增加速率，我們把上式進行微分可獲得 d株高/dvstage = Parm(21)*Parm(22)*Vstage^(Parm(22)-1) 程式碼在2750行 !YA PDMH=(PARM(21)+(PARM(22)*(VSTAGE+PDV/2.0)**1.37))*PDV*SLOW*HTFATR IF (VSTAGE.GT.0.0) THEN F1=PARM(21)*PARM(22)*(VSTAGE)**(PARM(22)-1.) ELSE F1=0.0 END IF IF (VSTAGE+PDV.GT.0.0) THEN F2=PARM(21)*PARM(22)*(VSTAGE+PDV)**(PARM(22)-1.) ELSE F2=0.0 END IF PDMH=0.5*(F1+F2)*PDV 這裡的計算方法很值得學習，它分別計算 Vstage (這一個時間點)和 Vstage + dv (下一個時間點所對應的株高增長速度，再計算平均植以取得較準確的生長速率值。值得注意的是，在第4270行的時候開始累計株高，這邊使用Maturity group和 R stage作為條件 ! CALCULATE CHANGES IN MAINSTEM AND B...

日射量資料剖析

Diurnal Variation 從 Charles-Edward et al., (1986) 建議可以使用half-wave sine response (De Vries, 1955)，或Full-wave sine function (Charles-Edwards and Acock, 1977)，2Dsoil 是使用half wave sine function (P42, Equation 5.11)，GLYCIM 模式使用的公式應該也是 Half wave sine function。 Half wave sine function $I(t) = \tfrac{\pi S}{2h} \times \sin(\tfrac{\pi t}{h}) $ , 0 < t < h Full wave sine function $I(t) = S \times [1 + \sin (\tfrac{2 \pi t}{h} + \tfrac{3 \pi}{2})/h]$ 其中，S是整天的日射量累積值，h為日照時數，t為日照開始的小時數。輻射基本概念參考 Campbell and Norman (1998) An introduction to environmental biophysics. P 148. 回憶一下能量和波長轉換的方程式 $e = \tfrac{hc}{\lambda}$ h是普朗克常數($6.63 \times 10^{-34} $ Js)，$\lambda$是波長(m)。光合有效輻射(photosynthetically active radiation, PAR) 是指 400 - 700 nm 的輻射量，對於太陽在海平面、太陽天頂角為60度的情況下，在400-700 nm 波段內的太陽輻射，中位波長約為550 nm，如果有一個光子的波長是550 nm，那麼它的能量為 $e = \tfrac{6.63 \times 10^{-34} (Js) \times 3 \times 10^8 (m/s)}{550 \times 10^{-9} (m)} = 3.6 \times 10^{-19} $ (J)...

氣溫日資料轉換為小時資料

氣溫日資料轉換為小時資料 Modeling temporal variation in air temperature From: Campbell, G.S. and J.M. Norman (1997) An introduction to environmental biophysics. Equation (2.2) and (2.3) page23 dimentionless diurnal temperature function Gamma(t) = 0.44 – 0.46 sin(w * t + 0.9) + 0.11 sin(2wt + 0.9) Where w = pi/12, t is time of day in hours The hourly temperature calculated from following equations T(t) = Tmax,i-1 Gamma(t) + Tmin,i [1-Gamma(t) ] t = 0- 5 T(t) = Tmax,i Gamma(t) + Tmin,i [1-Gamma(t) ] t = 5 - 14 T(t) = Tmax,i Gamma(t) + Tmin,i+1 [1-Gamma(t) ] t = 14-0 Here, Tmax is the daily maximum temperature and Tmin is the minimum temperature. The subscript I represents the present day; i - 1 is the previous day, and i + 1 is the next day. 檔案下載 https://drive.google.com/file/d/1AqnJnoA3M0SyjoCvRWWSDPR2OSH0U6Iq/view?usp=sharing

日照時長 Day length

太陽赤緯 (solar declination- δ delta) $\sin\delta = 0.39785 \sin[278.97 + 0.9856J + 1.9165 \sin(356.6+0.9856J)]$ 其中，J 是積日(Julian day, DOY) 日長 (day length - hd) $h_d = \cos^{-1}(\frac{\cos\psi - \sin\phi \sin\delta}{\cos\phi\cos\delta})$ 其中， ϕ phi, 是緯度 δ delta, 是太陽赤緯 $\psi$ psi, 是日出或日落時的天頂角，如果選擇民用曙暮光 (civil twilight)，其定義為太陽在地平線下6度，因此 ψ 應該使用96度，計算上要輸入 $\psi = 96 \pi /180$ 如果使用太陽在地平線就算起，選擇 $\psi$ 為90度，因為 $\cos\psi = 0$，因此公式又可以改成 $h_d = \cos^{-1}(-\frac{\sin\phi \sin\delta}{\cos\phi\cos\delta}) = \cos^{-1}(-\tan\phi \tan\delta)$ Excel檔下載連結 https://drive.google.com/file/d/1Oj5n220PiloLYAqG-HdYLg8TZAjmXpkO/view?usp=sharing 參考文獻：出自於 Campbell and Norman. 1998. An introduction to environmental biophysics. 第170頁

git

打開git bash，他預設的應該是在users 底下的資料夾我們就要移動到程式碼所在的資料夾建立倉庫(repository)，這邊要留意的是如果要移動到d槽，要先輸入 cd /d 才會到D槽下面，如果按d: 的話是不會有反應的。接著就如同cmd一樣，cd 到我們要的資料夾下面。 cd Python cd 'maizsim phenol' 如果資料夾內還沒有git的專案資料，就必須建立倉庫，輸入 git init 檢查資料夾，就會有一個.git的資料夾出現了這裡面的main.py, out.csv, phenology.py, __pycache__, READMD.md 是我所建立的檔案，我們可以透過 git add . 將所有檔案提交到倉庫裡面裡面，使用 git status 檢查裡面的檔案狀態接下來我們要把檔案加到repository裡面，使用再檢視一下狀態(git status) 接下來就使用commit，可以在這裡加入註解 git commit -m "modify variable" -m 後面的""可以加上我們任何需要的註解接著我們就要把專案上傳到github裡面去，使用 git remote add origin https://github.com/Chuchung0604/MAIZSIM_phenolstage 如果專案已經有連線的話，就可以跳過這個步驟，add 和 commit 完後就可以push 到github上去，push的指令如下 git push -u origin master https://www.itread01.com/content/1548500237.html https://ithelp.ithome.com.tw/users/20004901/ironman/525

蒸氣壓差 (VPD)

蒸氣壓差的英文為Vapor pressure deficit，縮寫為VPD。大氣當中所能保有的水蒸氣量為蒸氣壓，在這裡是把氣體濃度使用分壓方式表示，在特定的環境下大氣所能含有最多的蒸氣壓量就是飽和蒸汽壓 (saturated vapor pressure, SVP)，當氣溫越高的時候，大氣所能含有的水分就增加，溫度降低時，飽和蒸汽壓就降低 (我們可以從微觀的想法，溫度越高時水分子的動能越大，因此大氣當中更多的飽和蒸汽壓)。實際蒸氣壓 (actual vapor pressure, AVP)是指大氣當中現有的水分含量，因此，蒸氣壓差(VPD)是指實際蒸氣壓和飽和蒸汽壓之間的差值，我們可以想為是大氣當中實際水蒸氣濃度和飽和濃度的差值，因此蒸氣壓差越大，就代表空氣越乾燥。相對溼度(relative humidity, RH) 是實際蒸氣壓和飽和蒸汽壓的比值，相對溼度越低，空氣越乾燥。蒸氣壓差是濃度差的概念，而相對溼度是比例的概念，因此在科學上，蒸氣壓差比相對濕度更具有代表性。關鍵概念 $$RH = \frac{AVP}{SVP} \times 100$$ $$VPD = VSP - AVP$$ 如何計算大氣蒸氣壓差 AVPD? 1. 計算飽和蒸氣壓(SVP)：使用 Buck (1981) 的公式8 $$SVP = [1.0007+3.46\times 10^{-6}\times P ]\times 6.1121\times e^{17.502 T / (T+240.97) }$$ 其中SVP的單位為kPa，P為大氣壓力，1 atm為101.325 kPa，T為溫度( °C )。 2. 計算實際蒸氣壓(AVP) $$AVP= RH \times SVP $$ 3. 計算蒸氣壓差(VPD) $$ VPD = SVP-AVP$$ 或 $$ VPD = SVP \times (1-RH)$$ 如何計算葉片蒸氣壓差 LVPD? 葉片蒸氣壓差的計算方法和大氣蒸氣壓差的算法類似，區別在於先利用葉片溫度計算葉片飽和蒸汽壓差(LSVP)，再跟大氣實際蒸氣壓差(AVP = RH * SVP ) 計算葉片蒸氣壓差 LVDP = LSVP - AVP 參考資料 https://pulsegrow.com/blogs/learn/vpd Buck A.L. 1981. ...

Python 解析多光譜 (2)

接下來我們就試著用Python畫出NDVI的作品一樣的，先載入必要的套件，並設定working directory import os import matplotlib.pyplot as plt import numpy as np import rasterio as rio import geopandas as gpd import earthpy as et from earthpy import spatial as es from earthpy import plot as ep path = "D:\Python\multispectral" os.chdir(path) 接著，利用rasterio載入影像，用shape函式查看影像基本特性 pine_path = os.path.join("pineapple","pineapple0808.tif") with rio.open(pine_path) as msi: pine0808 = msi.read() pine0808.shape (6, 4604, 9634) 再來就利用earthpy 的 normalized_diff 函式來算 NDVI，再使用ep.plot_bands來畫圖 pine_ndvi = es.normalized_diff(pine0808[4], pine0808[2]) ep.plot_bands(pine_ndvi, cmap='PiYG', scale=True, vmin=-1,vmax=1, title="NDVI of Pineapple 0808") plt.show() 我們來觀察一下 es.normalized_diff 這一個函式吧。 help(es.normalized_diff) normalized_diff(b1, b2) Take two n-dimensional numpy arrays and calculate the normalized difference. Mat...