蒸氣壓差 (VPD)

蒸氣壓差的英文為Vapor pressure deficit，縮寫為VPD。

大氣當中所能保有的水蒸氣量為蒸氣壓，在這裡是把氣體濃度使用分壓方式表示，在特定的環境下大氣所能含有最多的蒸氣壓量就是飽和蒸汽壓 (saturated vapor pressure, SVP)，當氣溫越高的時候，大氣所能含有的水分就增加，溫度降低時，飽和蒸汽壓就降低 (我們可以從微觀的想法，溫度越高時水分子的動能越大，因此大氣當中更多的飽和蒸汽壓)。

實際蒸氣壓 (actual vapor pressure, AVP)是指大氣當中現有的水分含量，因此，蒸氣壓差(VPD)是指實際蒸氣壓和飽和蒸汽壓之間的差值，我們可以想為是大氣當中實際水蒸氣濃度和飽和濃度的差值，因此蒸氣壓差越大，就代表空氣越乾燥。

相對溼度(relative humidity, RH) 是實際蒸氣壓和飽和蒸汽壓的比值，相對溼度越低，空氣越乾燥。蒸氣壓差是濃度差的概念，而相對溼度是比例的概念，因此在科學上，蒸氣壓差比相對濕度更具有代表性。

關鍵概念

RH(%) = AVP/SVP*100
VPD = SVP - 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. New equations for computing vapor pressure and enhancement factor. American Meteorological Society 20:1527-1532.

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.
    
    Math will be calculated (b1-b2) / (b1 + b2).
    
    Parameters
    ----------
    b1, b2 : numpy arrays
        Two numpy arrays that will be used to calculate the normalized
        difference. Math will be calculated (b1-b2) / (b1+b2).
    
    Returns
    ----------
    n_diff : numpy array
        The element-wise result of (b1-b2) / (b1+b2) calculation. Inf values
        are set to nan. Array returned as masked if result includes nan values.

看起來他是利用numpy 資料進行向量運算，這個函示已經指定好第1個變數一定要是NIR，第2個變數是 Red band，不然算出來的資料會剛好相反。

儲存圖層

Python 需要建立出檔案，再把資料寫入到檔案當中。

先設定meta data黨，載入原本圖層的meta data

with rio.open(pine_path) as msi: 
    pine_data = msi.read() 
    pine_meta = msi.profile 
 

接著檢視 meta data

pine_meta

{'driver': 'GTiff', 'dtype': 'uint16', 'nodata': None, 'width': 9634, 'height': 4604, 'count': 6, 'crs': CRS.from_dict(init='epsg:3824'), 'transform': Affine(2.2059400000054369e-07, 0.0, 120.470607193754, 0.0, -2.0339900000014106e-07, 23.5296130451426), 'blockxsize': 256, 'blockysize': 256, 'tiled': True, 'compress': 'lzw', 'interleave': 'pixel'}

這邊可以看到，count 6 表示圖層有6個波段，資料型態是 16位元無符號整數 (uint 16)。我們必須把圖層改成只有1個波段，資料型態改為64位浮點數 (float64)。

pine_meta['count'] = 1 
pine_meta['dtype'] = "float64" 

接著，我們就可以建立圖層，利用rasterio把資料寫入檔案中。

pine_ndvi_outpath = os.path.join("pineapple","pineapple0808_ndvi.tif") 
with rio.open(pine_ndvi_outpath, 'w', **pine_meta) as dst: 
    dst.write(pine_ndvi, 1)

最後，我們來觀察一下earhpy的原始碼，https://earthpy.readthedocs.io/en/latest/_modules/earthpy/spatial.html

def normalized_diff(b1, b2):
    """Take two n-dimensional numpy arrays and calculate the normalized
    difference.

    Math will be calculated (b1-b2) / (b1 + b2).

    Parameters
    ----------
    b1, b2 : numpy arrays
        Two numpy arrays that will be used to calculate the normalized
        difference. Math will be calculated (b1-b2) / (b1+b2).

    Returns
    ----------
    n_diff : numpy array
        The element-wise result of (b1-b2) / (b1+b2) calculation. Inf values
        are set to nan. Array returned as masked if result includes nan values.

    Examples
    --------
    >>> import numpy as np
    >>> import earthpy.spatial as es
    >>> # Calculate normalized difference vegetation index
    >>> nir_band = np.array([[6, 7, 8, 9, 10], [16, 17, 18, 19, 20]])
    >>> red_band = np.array([[1, 2, 3, 4, 5], [11, 12, 13, 14, 15]])
    >>> ndvi = es.normalized_diff(b1=nir_band, b2=red_band)
    >>> ndvi
    array([[0.71428571, 0.55555556, 0.45454545, 0.38461538, 0.33333333],
           [0.18518519, 0.17241379, 0.16129032, 0.15151515, 0.14285714]])

    >>> # Calculate normalized burn ratio
    >>> nir_band = np.array([[8, 10, 13, 17, 15], [18, 20, 22, 23, 25]])
    >>> swir_band = np.array([[6, 7, 8, 9, 10], [16, 17, 18, 19, 20]])
    >>> nbr = es.normalized_diff(b1=nir_band, b2=swir_band)
    >>> nbr
    array([[0.14285714, 0.17647059, 0.23809524, 0.30769231, 0.2       ],
           [0.05882353, 0.08108108, 0.1       , 0.0952381 , 0.11111111]])
    """
    if not (b1.shape == b2.shape):
        raise ValueError("Both arrays should have the same dimensions")

    # Ignore warning for division by zero
    with np.errstate(divide="ignore"):
        n_diff = (b1 - b2) / (b1 + b2)

    # Set inf values to nan and provide custom warning
    if np.isinf(n_diff).any():
        warnings.warn(
            "Divide by zero produced infinity values that will be replaced "
            "with nan values",
            Warning,
        )
        n_diff[np.isinf(n_diff)] = np.nan

    # Mask invalid values
    if np.isnan(n_diff).any():
        n_diff = np.ma.masked_invalid(n_diff)

    return n_diff

參考網頁

https://www.earthdatascience.org/courses/use-data-open-source-python/multispectral-remote-sensing/vegetation-indices-in-python/calculate-NDVI-python/

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，但是不選擇就會有錯誤訊息，這根本就是：「沒關係，我沒有強迫你一定要選擇哦。」決定不選就發現爆掉了，有夠機歪的啦。

算完之後預設是灰階的圖片，這裡有一個小技巧，在圖層屬性的Symbology，選擇 single band pseudocolor，選擇紅黃綠的色階，把最小值設定為0，最大值設定為1。

將~將~將~，一張看起來舒服的NDVI照片就出現了。

單一波段匯出

先前測試的結果，轉成xyz之後檔案會變得非常大(大概有5G，根本打不開)，我們就選定中間都是鳳梨的區域來進行。

先數化圖層

再從Raster -> Extraction -> Clip Raster by Mask Layer，記得選擇多光譜的圖層和數化的polygon進行裁剪。

套疊RGB一下，我們感興趣的區域就是這麼大

重點就是 Rearrange Band

選擇Tool box 裡面的 Rearrange Bands，可以將波段分開來

我們就從藍光開始做起，選擇Band 1

副檔名選擇.xyz

經過來漫長的等待後，就順利的產製檔案

用notepad打開來看看

果然資料就是用XYZ的方式呈現，我有發現，QGIS會把圖片修正為長方形，雖然原本我們所切割的圖層面積較小，但它似乎會用圖層最邊界的XY座標進行計算，目前還想不到破解方法。

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()

Matplotlib

Matplotlib的imshow函數可以進行繪圖，我們使用matplotlib.pyplot(縮寫為plt) 進行繪圖，imshow這一個函數只能畫出單一圖層，我們就選擇藍光波段來畫，在圖後面加上[0]。

fig, ax = plt.subplots()
ax.imshow(pine0808[0])
ax.set_title("Pineapple Multispectral image \n Band 1 Blue")
plt.show()

Python 的計數是從0開始，我們就可以依此類推：
[0] - Blue
[1] - Green
[2] - Red
[3] - Red Edge
[4] - Near IR

Earthpy

earthpy 函數可以畫單一圖層，也可以進行RGB 的套疊

先用ep.plot_bands()來畫單一波段的圖層

ep.plot_bands(pine0808[0],
              title="Pineapple Multispectral image \n Band 1 Blue",
              cbar=True,vmin=100,vmax=2000)
plt.show()

上面的vmax和vmin可以調整顯示的波段值

ep.plot_bands()把所有的波段畫出來，話說不知道為什麼，正射影像會多一個波段的資訊，但是那實際是空的。

titles = ["Blue Band", "Green Band", "Red Band","Red Edge", "Near Infrared (NIR)","plot"]
ep.plot_bands(pine0808,
              figsize=(12, 5),
              cols=3,
              title=titles,
              cbar=False,
              vmax=15000)
plt.show()

ep.plot_rgb()，可以選擇波段，畫出具有紅、藍、綠三個顏色，我們先嘗試使用原本的顏色來畫圖

ep.plot_rgb(pine0808,
            rgb=[2, 1, 0],
            title="RGB Composite image - pineapple",
            figsize=(15, 8),
            stretch=True)
plt.show()

畫出來的圖比qgis自然許多，但是圖片偏暗，應該和各個波段的訊號有關。

另一個在植物光譜上面很有意思的方法，是用紅外光取代紅光的波段。將上面指令稍作修改rgb=[4, 1, 0]

參考網站

https://www.earthdatascience.org/courses/use-data-open-source-python/multispectral-remote-sensing/intro-multispectral-data/

https://github.com/earthlab/earthpy

https://earthpy.readthedocs.io/en/latest/gallery_vignettes/index.html