遥感

python实现基于GDAL的哨兵2影像NDVI值计算

2019-12-25  本文已影响0人  GIS无情老博士

NDVI是什么

NDVI(归一化植被指数)是近红外波段的反射值与红光波段的反射值之差比上两者之和。即(NIR-R)/(NIR+R),NIR为近红外波段的反射值,R为红光波段的反射值。归一化植被指数是反映农作物长势和营养信息的重要参数之一。根据该参数,可以知道不同季节的农作物对氮的需求量, 对合理施用氮肥具有重要的指导作用。

计算代码

注意使用的是哨兵2影像,第8波段为红外波段,第4波段为红光波段

from osgeo import gdal
import os
import time
import numpy as np

img_root = "E:\\商丘yx"
img_type = (".img", ".dat", "tiff")
driver = gdal.GetDriverByName('GTiff')
result_name_temp = "temp2.tiff"
start = time.clock()

result_path = os.path.join(img_root, result_name_temp)
# 文件存在则删除文件
if os.path.exists(result_path):
    os.remove(result_path)


rater_file = "E:\\商丘yx\\t50smc_20190604t025551_b01.img"


def get_ndvi(path):  # 计算某一影像的ndvi值,返回二维数组
    dataset = gdal.Open(path)
    cols = dataset.RasterXSize  # 列数
    rows = dataset.RasterYSize  # 行数

    band8 = dataset.GetRasterBand(8).ReadAsArray(0, 0, cols, rows)
    band4 = dataset.GetRasterBand(4).ReadAsArray(0, 0, cols, rows)
    molecule = band8 - band4
    denominator = band8 + band4
    del dataset
    band = molecule / denominator
    band[band > 1] = 9999  # 过滤异常值
    return band


def compute_band(file):
    dataset = gdal.Open(file)
    cols = dataset.RasterXSize  # 列数
    rows = dataset.RasterYSize  # 行数
    # 生成影像
    target_ds = gdal.GetDriverByName('GTiff').Create(result_path, xsize=cols, ysize=rows, bands=1,
                                                     eType=gdal.GDT_Float32)
    target_ds.SetGeoTransform(dataset.GetGeoTransform())
    target_ds.SetProjection(dataset.GetProjection())
    del dataset
    band = get_ndvi(file)
    target_ds.GetRasterBand(1).SetNoDataValue(9999)
    target_ds.GetRasterBand(1).WriteArray(band)
    target_ds.FlushCache()


compute_band(rater_file)
elapsed = (time.clock() - start)
print("计算ndvi耗时:", elapsed)

计算前

计算前

计算结果

计算结果渲染
上一篇 下一篇

猜你喜欢

热点阅读