接下來我們就試著用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，不然算出來的資料會剛好相反。
儲存圖層
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
參考網頁
沒有留言:
張貼留言