遥感指数生成（python）

文章目录

• ​​相关代码​​
• ​​2. 生成NDVI​​
• ​​3. 生成NDWI​​

相关代码

import os, sys, time
import numpy as np
from osgeo import ogr
from osgeo import gdal
from osgeo import gdal_array as ga

def read_img(filename):
    dataset=gdal.Open(filename)

    im_width = dataset.RasterXSize
    im_height = dataset.RasterYSize

    im_geotrans = dataset.GetGeoTransform()
    im_proj = dataset.GetProjection()
    im_data = dataset.ReadAsArray(0,0,im_width,im_height)

    del dataset 
    return im_proj,im_geotrans,im_width, im_height,im_data

def ndvi(im_geotrans,im_proj, im_redBand, im_nirBand, im_width, im_height, outfile):
    driver = gdal.GetDriverByName("GTiff")
    dataset = driver.Create(outfile, im_width, im_height, 1, gdal.GDT_Int16)
    dataset.SetGeoTransform(im_geotrans)
    dataset.SetProjection(im_proj)
    ga.numpy.seterr(all='ignore')
    ndvi = ((im_nirBand - im_redBand) * 1.0)/((im_nirBand + im_redBand) * 1.0)
    ndvi = ndvi.astype(np.float32)
    img_min = 0
    img_max = 255
    ndvi = stretch_n(ndvi, img_min, img_max)
    ndvi = ndvi.astype(np.uint8)
    dataset.GetRasterBand(1).WriteArray(ndvi)
    dataset.GetRasterBand(1).SetNoDataValue(0)
    return ndvi

2. 生成NDVI

if __name__ == "__main__":
    import glob
    imgList = ['new_area.tif']

    for input in imgList:
        imgName = os.path.split(input)[-1].split('.')[0]
        im_proj, im_geotrans, im_width, im_height, im_data = read_img(input)
        data = gdal.Open(input)

        r = data.GetRasterBand(1).ReadAsArray()
        g = data.GetRasterBand(2).ReadAsArray()
        ir = data.GetRasterBand(4).ReadAsArray()

        ndvi_ = ndvi(im_geotrans, im_proj, r, ir, im_width, im_height, outfile=f'{imgName}_ndvi.tif')

3. 生成NDWI

if __name__ == "__main__":
    import glob
    imgList = ['new_area.tif']

    for input in imgList:
        imgName = os.path.split(input)[-1].split('.')[0]
        im_proj, im_geotrans, im_width, im_height, im_data = read_img(input)
        data = gdal.Open(input)

        r = data.GetRasterBand(1).ReadAsArray()
        g = data.GetRasterBand(2).ReadAsArray()
        ir = data.GetRasterBand(4).ReadAsArray()

        ndwi_ = ndvi(im_geotrans, im_proj, g, ir, im_width, im_height, outfile=f'{imgName}_ndwi.tif')