本文主要参考《python地理空间分析指南》。目的在于利用归一化植被指数NDVI来实现植被的分类制图。分类的方法比较简单，就是阈值切片。我对书中的内容进行了改进，并在分析中会指出书中存在的问题，欢迎进一步探讨。
二、数据进行实验的话数据很重要，本实验的数据是tiff格式的影像和一个shp矢量文件。原始数据的下载地址为：
http://git.io/v3fS9
下载好的数据放到项目文件夹下，直接解压即可。下面我们来看一下下载到的数据是什么样的：
下载到的数据虽然看着很多，但其实我们能够用到的只有两个，即原始影像数据文件farm.tif和农田的矢量文件field.shp，用ENVI软件打开这两个文件：
可以看到：原始数据中的farm.tif是3波段数据（打开波段查看器可以看的更清楚），而且skdls的田地和水体都是真彩色影像，所以可以推测出farm.tif的三个波段对应的是RGB三个波段；而field.shp对应的则是图中间的一小块区域；另外查看像元值可以知道，图像为8bit图像，各个波段的灰度范围都是0-255。同时该图像还带了空间投影及地理信息，是一幅比较典型的遥感影像。
三、实验这里需要重点强调的一点是，实验采用的植被指数NDVI的计算公式为：
即利用红光和近红外波段的反射率进行运算。
但是，我个人认为原实验中有两处小错误：①NDVI的计算公式利用的是反射率而非图像DN值，用DN值算下来的结果其实没有任何物理意义，实际上我们在进行计算的时候，根据需求，往往采用入瞳处的反射率或者地表真实反射率，计算这两种反射率都需要用到头文件；②原始图像应该是真彩色图像，所以是根本没有近红外波段的。
不过这只是一个简单的小实验，不影响我们的处理思路。
总体思路只有两步，计算NDVI和影像分类：
1.计算NDVI这一步主要的目的是利用gdal_array将农田影像载入numpy数组（gdal_array只能载入图像，而gdal可以获取图像的地理参照），然后分别获取近红外和红光波段；之后是将shp文件栅格化，填充并转化为numpy数组；最后是裁剪相应的波段，计算NDVI并保存。
先直接给出主要代码(ndvi.py):
from osgeo import gdal_arrayfrom PIL import Image, ImageDrawdef imageToArray(i): """ function: 将PIL图像库的数组转换为一个gdal_array图片 """ a = gdal_array.numpy.fromstring(i.tobytes(), 'b') a.shape = i.im.size[1], i.im.size[0] return adef world2Pixel(geoMatrix, x, y): """ function: 该函数使用gdal的geoMatrix方法，将目标图片的地理坐标转换成像素坐标 """ ulX = geoMatrix[0] ulY = geoMatrix[3] xDist = geoMatrix[1] yDist = geoMatrix[5] pixel = int((x - ulX) / xDist) line = int((ulY - y) / abs(yDist)) return (pixel, line)import gdalnumeric# 获取数据部分-----------------------------------------------------# 需要读取的影像source = "./NDVI-update/NDVI-update/farm.tif"# 最终生成的影像target = "ndvi.tif"# 将载入的数据转换成gdal_array数组srcArray = gdal_array.LoadFile(source)# 同时将源数据作为gdal图片载入，获取地理信息参照srcImage = gdal.Open(source)geoTrans = srcImage.GetGeoTransform()# 获取红光和近红外波段r = srcArray[1]ir = srcArray[2] # 这个应该是蓝光波段# 裁剪影像部分-----------------------------------------------------# 根据shapefile文件创建一个OGR图层field = ogr.Open("./NDVI-update/NDVI-update/field.shp")# 定义一个图层确保OGR格式完好lyr = field.GetLayer("field")# shapefile文件中只有一个多边形poly = lyr.GetNextFeature()# 将图层坐标转化为图片的像素坐标minX, maxX, minY, maxY = lyr.GetExtent()ulX, ulY = world2Pixel(geoTrans, minX, maxY) # 左上角坐标lrX, lrY = world2Pixel(geoTrans, maxX, minY) # 右下角坐标# 计算新图片的像素尺寸pxWidth = int(lrX - ulX)pxHeight = int(lrY - ulY)# 为遮罩层创建一个合适尺寸的空白图片clipped = gdal_array.numpy.zeros((3, pxHeight, pxWidth), gdal_array.numpy.uint8)# 裁剪波段rclip = r[ulY:lrY, ulX:lrX]irclip = ir[ulY:lrY, ulX:lrX]# 为图片创建一个新的geomatrix对象geoTrans = list(geoTrans)geoTrans[0] = minXgeoTrans[3] = maxY# 在8位黑白遮罩图片上把点映射为像素用于绘制区域边界points = []pixels = []# 获取多边形几何图形geom = poly.GetGeometryRef()pts = geom.GetGeometryRef(0)# 遍历几何图像，将点转换为易于管理的list对象for p in range(pts.GetPointCount()): points.append((pts.GetX(p), pts.GetY(p)))# 遍历点集并将之映射为像素，然后把这些像素添加到list中for p in points: pixels.append(world2Pixel(geoTrans, p[0], p[1]))# 用“L”模式创建一个栅格多边形的黑白图片，白色值=1rasterPoly = Image.new("L", (pxWidth, pxHeight), 1)# 创建一个PIL绘制对象rasterize = ImageDraw.Draw(rasterPoly)# 将像素点存储为图片上的多边形， 黑色=0rasterize.polygon(pixels, 0)# 将图片转换为gdal_array以便将其当作一个遮罩数组来用mask = imageToArray(rasterPoly)# 裁剪波段---------------------------------------------------------------# 使用遮罩分别裁剪红光和近红外波段rclip = gdal_array.numpy.choose(mask, (rclip, 0)).astype(gdal_array.numpy.uint8)irclip = gdal_array.numpy.choose(mask, (irclip, 0)).astype(gdal_array.numpy.uint8)# 计算NDVI---------------------------------------------------------------# 忽略裁剪过程中numpy遇到的任何NaN值gdal_array.numpy.seterr(all='ignore')# 按照公式计算ndvi = 1.0 * ((irclip - rclip) / (irclip + rclip + 1.0))# 将结果中的所有NaN值转换为0ndvi = gdal_array.numpy.nan_to_num(ndvi)# 将NDVI图像进行保存，并调整地理参照信息gdalnumeric.SaveArray(ndvi, target, format="GTiff", prototype=srcImage)# 并调整地理参照信息update = gdal.Open(target, 1)update.SetGeoTransform(list(geoTrans))update.GetRasterBand(1).SetNoDataValue(0.0)update = None直接执行上述程序，会生成计算好的NDVI影像，下面直接展示结果：
对于结果进行了分析，NDVI大部分的值还是介于0-1.2之间的，但是仍有像元的值为150。。。。。。考虑了一下，可能是程序中的计算公式，为了防止分母为0，在分母中又多加1，这样就有可能出现异常值。
2.影像分类这一步并没有用分类方法对其进行分类，而是使用棕色到深绿色7个类别对其进行分类。首先读取处理好的ndvi影像，然后进行直方图均衡化处理；之后是创建7个类别，对影像数据进行分类。
下面直接给出分类的代码（classification.py）：
from osgeo import gdal_array as gdimport operatorfrom functools import reducedef histogram(a, bins=list(range(256))): """ 多维数组的直方图函数 """ # 直方图的数组均衡化、排序和分割 fa = a.flat n = gd.numpy.searchsorted(gd.numpy.sort(fa), bins) n = gd.numpy.concatenate([n, [len(fa)]]) hist = n[1:] - n[:-1] return histdef stretch(a): """ 在gdal_array数组图像上执行直方图平衡化操作 """ hist = histogram(a) lut = [] for b in range(0, len(hist), 256): # 创建等间隔区间 step = reduce(operator.add, hist[b:b+256]) / 255 # 创建均衡化查找表 n = 0 for i in range(256): lut.append(n / step) n += hist[i+b] gd.numpy.take(lut, a, out=a) return a# 载入数据-----------------------------------------------------------# 原图片为NDVI脚本输出source = "./ndvi.tif"# 分类图片的目标文件名target = "ndvi_color.tif"# 将图片载入数组ndvi = gd.LoadFile(source).astype(gd.numpy.uint8)# 预解析NDVI# 执行直方图均衡化操作以便能够使用所有类别ndvi = stretch(ndvi)# 创建一个尺寸大小和NDVI相同的3波段图像rgb = gd.numpy.zeros((3, len(ndvi), len(ndvi[0])), gd.numpy.uint8)# 创建类-------------------------------------------------------------# 类别列表中NDVI的区间上限值，需要注意上限和下限值在列表两端classes = [58, 73, 110, 147, 184, 220, 255]# 颜色查找表。lut必须从深棕色到深绿色代表的RGB元组匹配lut = [[120, 69, 25], [255, 178, 74], [255, 237, 166], [173, 232, 94], [135, 181, 64], [3, 156, 0], [1, 100, 0]]# 第一种类别的起始值start = 1# 获取每个类别的区间值，获取区间值，然后使用遮罩过滤for i in range(len(classes)): mask = gd.numpy.logical_and(start <= ndvi, ndvi <= classes[i]) for j in range(len(lut[i])): rgb[j] = gd.numpy.choose(mask, (rgb[j], lut[i][j])) start = classes[i] + 1# 保存彩色的NDVI数据为geotiff图片output = gd.SaveArray(rgb, target, format="GTiff", prototype=source)output = None直接执行上述代码，最终的实验结果为：
个人感觉实验的结果并不好，一眼看过去，主要都是翠绿色，和一点墨绿色，其他颜色没有。。。虽然分了7类，不过这个分类的实验瑕疵挺大的。
四、分析
1.这个实验的主要问题表现在两个方面：①NDVI的计算公式利用的是反射率而非图像DN值，而原始实验中采用的是像元DN值②原始图像应该是真彩色图像，所以是根本没有近红外波段的。
2.关于分类，一般比较简单的是采用阈值分割，不过这个阈值切片的分类效果是真的很差。
3.实验的目的主要在于提供图像处理的思路，不过效果是真的不好。