python+GDAL绘制遥感影像直方图

参考《Python地理空间分析指南》第六章
使用osgeo的gdal_array模块（处理遥感影像）和turtle模块（绘制直方图）、numpy模块（对数组进行操作，其实gdal_array模块中添加了numpy的引用，可以不导入numpy，但是这样做的话使用numpy还需要用gdal_array作为命名空间：gdal_array.numpy，使用起来会比较繁琐）

用到的遥感影像：

下面是代码实现，解释见注释
内联代码片

from osgeo import gdal_array
import turtle as t
import numpy
img='swap.tif'
img_arr=gdal_array.LoadFile(img)
#计算频率
gray=range(0,256)
hist=[]
#迭代图像中的不同波段
for item in img_arr:
    #迭代器迭代出每个像素
    im_iter=item.flat
    #对像素的灰度值进行排序
    rank=numpy.sort(im_iter)
    #找出灰度值的顺序号，不包括最大的，searchsorted（a，v）在数组a中插入数组v
    #返回插入位置下标组成的数组
    index=numpy.searchsorted(rank,gray)
    #附加最大值到最后
    n=numpy.append(index, [len(im_iter)])
    #错位相减得到每个灰度值的像素个数
    item_hist=n[1:]-n[:-1]
    #作为列表元素附加到hist中，这个元素即代表某个波段的频率数据
    hist.append(item_hist)


#隐藏画笔
t.pen(shown=False)
#绘制坐标轴
t.color('black')
#设置窗体的大小，用screensize设置的是画布大小及背景色，窗体和画布不是一个概念。
# 如果画布大于窗体，窗体会出现滚动条，如果画布小于窗体，画布会填充整个窗体。
t.Screen().setup(1000,700)
#绘制坐标轴
t.up()
axes=[(-400,250),(-400,-250),(400,-250)]
for point in range(0,3):
    t.goto(axes[point])
    t.down()
#绘制坐标轴的标题
t.up()
t.goto(0,-300)
t.write('gray value',align='center',font=('Arial',12,'bold'))
t.up()
t.goto(-400,270)
t.write('frequency',align='center',font=('Arial',12,'bold'))
#x轴刻度线
t.up()
x=-400;y=-250
t.goto(x,y)
for i in range(1,11):
    x+=80
    t.goto(x,y)
    t.down()
    t.goto(x,y+10)
    t.up()
    t.goto(x,y-20)
    t.write(i*25,align='center')
#按照频率最大值伸缩y轴
max=0
#计算三通道影像中的最大频率值
for item in hist:
    hmax=item.max()
    if hmax>max:
        max=hmax
#y轴刻度间隔
y_interval= max // 10
#y轴刻度线
t.up()
x=-400;y=-250
t.goto(x,y)
for i in range(1,11):
    y+=50
    t.goto(x, y)
    t.down()
    t.goto(x+10, y)
    t.up()
    t.goto(x-50, y-6)
    t.write(i * y_interval, align='left')
#画直方图
#将横纵坐标轴的值与屏幕像素相对应，即一个刻度值代表屏幕上的几个像素
x_screen=800/255;y_screen=500/max
#用不同颜色绘制三个波段的直方图
color=['red','green','blue'];a=0
for item in hist:
    x = -400;y = -250
    t.up()
    t.goto(x,y)
    t.color(color[a])
    t.down()
    for i in range(0,256):
        y=item[i]*y_screen-250
        x=x_screen*i-400
        t.goto(x,y)
    a+=1
#绘制结束后不关闭窗口
t.done()

运行结果：
后续可以将部分代码封装到函数中，再用起来更方便