python：根据RGB波段计算相应的指数

import os
import numpy as np
from osgeo import gdal
from osgeo.gdalconst import GDT_Byte, GDT_Float32, GDT_UInt16
import os
from skimage import io as sio
from sklearn.preprocessing import normalize
img_path='I:/xxx/aa/img/'
lbl_path='I:/xxx/aa/relbl/'

img_paths = os.listdir(img_path)
os.chdir(r'I:/xxx/aaa/img_proj/')
for path in (img_paths):    
    RGBimg=gdal.Open(str(img_path+path))
    band=RGBimg.GetRasterBand(1)
    img_proj=sio.imread(str(img_path+path))
    B=img_proj[:,:,0]
    G=img_proj[:,:,1]
    R=img_proj[:,:,2]
    NIR=img_proj[:,:,3]

    B=B.astype(np.float16)
    G=G.astype(np.float16)
    R=R.astype(np.float16)
    NIR=NIR.astype(np.float16)

    ndvi=(NIR-R)/(NIR+R+1e-5)
    ndwi=(G-NIR)/(G+NIR+1e-5)

    max_ndvi=ndvi.max()
    max_ndwi=ndwi.max()
    
    min_ndvi=ndvi.min()
    min_ndwi=ndwi.min()
    ndvi=(ndvi-min_ndvi)/(max_ndvi-min_ndvi+1e-5)
    ndwi=(ndwi-min_ndwi)/(max_ndwi-min_ndwi+1e-5)    
    img_proj=img_proj[:,:,0:3]

#读取标签
    lbl_name=lbl_path+path.replace('aaa', 'bbb')
    lbl_name=lbl_name.replace('2017', '2016')
    lbl=sio.imread(lbl_name)


    ndvi_name=path.replace('naip', 'ndvi')
    ndwi_name=path.replace('naip', 'ndwi')
    img_proj_name=path
    lbl_proj_name=path.replace('naip', 'nlcd')

    #利用gdal创建结果栅格（带投影信息）
    gtiff_driver =gdal.GetDriverByName('GTiff')
    result = gtiff_driver.Create(ndvi_name,band.XSize, band.YSize, 1, GDT_Float32)
    result1 = gtiff_driver.Create(ndwi_name,band.XSize, band.YSize, 1, GDT_Float32)
    result2 = gtiff_driver.Create(img_proj_name,band.XSize, band.YSize, 3, GDT_Byte) 
    result3 = gtiff_driver.Create(lbl_proj_name,band.XSize, band.YSize, 1, GDT_Byte)       
    result.SetProjection(RGBimg.GetProjection())
    result1.SetProjection(RGBimg.GetProjection())  
    result2.SetProjection(RGBimg.GetProjection())   
    result3.SetProjection(RGBimg.GetProjection())  

    result.SetGeoTransform(RGBimg.GetGeoTransform()) 
    result1.SetGeoTransform(RGBimg.GetGeoTransform()) 
    result2.SetGeoTransform(RGBimg.GetGeoTransform())    
    result3.SetGeoTransform(RGBimg.GetGeoTransform())

    result.WriteArray(ndvi)
    result1.WriteArray(ndwi)
    result3.WriteArray(lbl)

    for i in range(3):
        result2.GetRasterBand(i+1).WriteArray(img_proj[:,:,i])


    result.FlushCache()
    result1.FlushCache()
    result2.FlushCache()
    result3.FlushCache()