延續上一講的雙域視角(巨孔隙/微孔隙)，這一講的參數會分成兩組來看:定義「桶子形狀」的靜態參數(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...
接下來我們就試著用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/
留言
張貼留言