跳到主要內容

發表文章

目前顯示的是 9月, 2020的文章

Tipping Bucket Model (3): Parameters

延續上一講的雙域視角(巨孔隙/微孔隙),這一講的參數會分成兩組來看:定義「桶子形狀」的靜態參數(LL、DUL、SAT、bulk density),跟決定「桶子怎麼漏水」的動態參數(SWCON、𝐾𝑠𝑎𝑡)。 3.1 靜態參數:桶子的三條水位線 (1) LL(Lower Limit,凋萎點下限) 這是植物完全無法吸取的水分下限,操作型定義通常對應到 −1.5 MPa(15 bar)的基質勢。在微孔隙/巨孔隙框架下,LL 描述的完全是微孔隙域——這個含水量下,水被吸附在極細孔隙的表面,基質勢極負,只存在於微孔隙的毛細跟吸附作用力範圍內,巨孔隙在這個含水量早就空了。 DSSAT 通常用 LL(而不是永久凋萎點 wilting point)這個名稱,是因為它其實是「特定作物根系可達到的下限」,理論上依作物種類、根系穿透力略有差異,但實務上常直接借用土壤的物理凋萎點。 (2) DUL(Drained Upper Limit,排水上限) 這是整套 bucket 理論裡最關鍵的一條線——它同時也是上一講講的「巨孔隙排水」的觸發點:𝑆𝑊> 𝐷𝑈𝐿 SW>DUL 才會啟動𝐷𝑅𝐴𝐼𝑁𝐿 方程。 物理意義上,DUL 是「重力排水已經(近似)結束、只剩微孔隙毛細力撐著水」的那個含水量,通常對應基質勢在 −10 到 −33 kPa 之間(視質地而定,砂土接近 −10 kPa,黏土接近 −33 kPa)。用雙域語言講:DUL 就是巨孔隙域清空、微孔隙域仍飽滿的那條分界線。這也是為什麼 DUL 不是一個絕對的物理常數,而是跟排水時間、量測方式高度相關的操作型定義——現場量測通常是「飽和後排水1–2天」的含水量,這個「1–2天」本身就是承認巨孔隙排水需要一點時間,但排水速率遠快於微孔隙。 (3) SAT(Saturation,飽和含水量) 這是巨孔隙+微孔隙全部孔隙都灌滿水時的含水量,約等於總孔隙度 𝜙 。注意上一講 Emerman 論文提到的關鍵假設——「巨孔隙域跟微孔隙域必須有相同的飽和含水量才能讓雙域模型的擬合站得住腳」。這句話反過來提醒你:SAT 這個看似簡單的參數,其實隱含了「兩個域共享同一個總孔隙空間」的假設,在真實土壤中不一定成立(比如強烈結構化的黏土,巨孔隙可能只佔總孔隙度很小一部分)。 (4) Bulk...

蒸氣壓差 (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...

QGIS 操作多光譜

每次發現多光譜的影像無法直接用相片瀏覽器打開,一整個就是眼神死😪😪😪,在俊毅的幫忙之下開始使用QGIS來看影像,自己可以疊合出NDVI的照片,總算知道自己沒有那麼笨,可以做一點基礎的工作。今天的目的是完成打開影像、套疊NDVI、轉成XYZ格式的練習,之後來試試使用Python進行。 這裡面最主要是使用Tool box 裡面的兩個工具,Raster calculator 可以進行波段的計算,長出具有特色的波段圖層(例如NDVI, NBI....),Rearrange Bands 把我們感興趣的波段匯出,這裡我選擇xyz產出習慣處理的資料格式。 有關影像 這邊使用的是正射過的影像為案例,鏡頭是MicaSense的RedEdge-MX,屬於5波段的多光譜相機,總共有藍光(475 nm, 20 nm width)、綠光(560 nm, 20 nm width)、紅光(668 nm, 10 nm width)、紅邊(717nm, 10 nm width)、近紅外光(840 nm, 40 nm width)。 操作的時候會有5個band,band 1 是波長最短的藍光,依序到波長最長的進紅外光(band 5)。 打開影像 用QGIS打開影像,一般的RGB影像(左圖)看起來就是舒舒服服,可是多光譜影像(右圖)長得麻麻喳喳的,越看越懷疑自己眼睛業障重。 我發現在QGIS可以設定波段的呈現,但似乎怎麼調都怎麼怪。  Raster calculator - 波段疊加 在Tool box 裡面有一個 raster calculator,可以用來進行波段的套疊,打開後就可以選擇波段進行計算,我們發現這裡面有一個NDVI的功能,點了add之後就會跑出一個對話視窗,選擇NIR是第5個波段,Red是第3個波段。 點選確認之後,就會發現Expression 的地方會變成 (0808多光譜正射@5 - 0808多光譜正射@3) / (0808多光譜正射@5 + 0808多光譜正射@3) 之後也可以直接輸入公式,來呈現我們要的波長疊合方式。 科普一下,NDVI = (NIR - R) / (NRI + R) ,越健康的植物越會反射紅外光,因此越接近1。 這邊有一個bug,Reference layer(s)的地方一定要輸入參考圖層,作為CRS的選擇,雖然它叫做optional,但是不選擇...

Python 解析多光譜 (1)

跟著Earth lab的內容逐步練習 基本概念 通常在表達波段的時候,都是選用波段的中間值表示,例如820-830 nm,我們就會說是825 nm。 高光譜和多光譜的差異在於波段的數量,從GIS geography可以找到以下的解釋 Multispectral: 3-10 wider bands. Hyperspectral: Hundreds of narrow bands 資料出處: https://gisgeography.com/multispectral-vs-hyperspectral-imagery-explained/ 使用Python 操作多光譜 rasterio.open()  函數用來打開多光譜影像 stack() 函數可用來引入多光譜資料 earthpy的plot_rgb() 函數用來指定波段 開始吧 首先載入所有必要的元件 import os import matplotlib.pyplot as plt import numpy as np import rasterio as rio import geopandas as gpd import earthpy as et import earthpy.plot as ep 需要額外安裝的套件包括rasterio, geopandas, earthypy,可以從anaconda power shell prompt 安裝。 接著設定路徑 path = "D:\Python\multispectral" os.chdir(path) plt.rcParams['figure.figsize'] = (10, 10) plt.rcParams['axes.titlesize'] = 20 接著我們載入影像,使用rasterio(縮寫為rio)的open function,把它取名為msi,這裡是multispectral image 的意思,把圖檔命名為pine0808 pine_path = os.path.join("pineapple","pineapple0808.tif") with rio.open(pine_path) as msi: pine0808 = msi.read...