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/