算一个NDVI

GDAL算NDVI，其实本质就是转换成数组，按NDVI算法计算，在绘图的时候出了点小问题，一般如果要在python上绘制一下，最好把NDVI的值拉伸到8bit，这样可以直观的看到，但如果你要继续处理，最好还是按照原本的计划将它保留：

from osgeo import gdal_array as ga

import gdal, ogr, os, osr

import numpy as np

b3=r'D:\宁波四明山\aod-ningbo\imagery\2016-\1\LC08_L1TP_118039_20170824_20170912_01_T1_B4.TIF'

b4=r'D:\宁波四明山\aod-ningbo\imagery\2016-\1\LC08_L1TP_118039_20170824_20170912_01_T1_B5.TIF'

arr=ga.LoadFile(b3)

arr1=ga.LoadFile(b4)

ga.numpy.seterr(all="ignore")

ndvi=((arr1-arr)*1.0)/((arr1+arr)*1.0)

ndvi1=ga.numpy.nan_to_num(ndvi)

target=r'D:\Temp\ndvi1.tif'

out=ga.SaveArray(ndvi1,target,format = "GTiff",prototype = b4)

out=None

如果是arcpy，那就更方便了：

from arcpy import env

from arcpy.sa import*

arcpy.CheckOutExtension("spatial")

env.workspace = r'C:\Your\workspace'

input = r'C:\Your\raster.tif'

result = "outputName.tif"

# You may need to change the band combinations.

# This is for 4-band NAIP imagery or Landsat TM.

NIR = input + "\Band_4"

Red = input + "\Band_3"

NIR_out = "NIR.tif"

Red_out = "Red.tif"

arcpy.CopyRaster_management(NIR,NIR_out)

arcpy.CopyRaster_management(Red, Red_out)

Num = arcpy.sa.Float(Raster(NIR_out) - Raster(Red_out))

Denom = arcpy.sa.Float(Raster(NIR_out) + Raster(Red_out))

NIR_eq = arcpy.sa.Divide(Num, Denom)

NIR_eq.save(result)