在python中使用GDAL为卫片制作缩略图

iojsioi

　　卫片和航片文件一般都很大，
　　有时候需要制作一个缩略图，放在导航窗体中做底图
　　这个文件不用太清晰，只要反映整体面貌即可
　　在ArcMap，PCI和Erdas中，可以直接截屏，或者使用重投影等功能很方便地完成，
　　最近在用python和GDAL，于是乎，查了查资料，并不复杂，
　　首先要确定python安装了numpy库，代码
　　import gdalimport numpyfrom gdalconst import *dataset = gdal.Open("F://数据//Raster//earth.img")width = dataset.RasterXSizeheight = dataset.RasterYSizebw = 0.1bh = 0.1x = int(width*bw)y = int(height*bh)datas = []for i in range(3):band = dataset.GetRasterBand(i+1)data = band.ReadAsArray(0,0,width,height,x,y)datas.append(numpy.reshape(data,(1,-1)))datas = numpy.concatenate(datas)driver = gdal.GetDriverByName("GTiff")tods = driver.Create("F://数据//Raster//tpix1.tif",x,y,3,options=["INTERLEAVE=PIXEL"])tods.WriteRaster(0,0,x,y,datas.tostring(),x,y,band_list=[1,2,3])
　　以上的这个方法，用到了numpy的矩阵运算，不过使用band的ReadAsArray
　　方法取得的数据，采用的是最邻近采样，虽然缩略图的要求不高，但是这个
　　采样方式总让人不爽，查了资料，使用gdal.ReprojectImage
　　可以设置采样方式，代码如下：
　　# -*- coding: cp936 -*-import gdalimport numpyfrom gdalconst import *dataset = gdal.Open("F://数据//Raster//earth.img")width = dataset.RasterXSizeheight = dataset.RasterYSizebw = 0.1bh = 0.1x = int(width*bw)y = int(height*bh)driver = gdal.GetDriverByName("GTiff")tods = driver.Create("F://数据//Raster//tpix1.tif",x,y,3,options=["INTERLEAVE=PIXEL"])tods.SetGeoTransform(dataset.GetGeoTransform())tods.SetProjection(dataset.GetProjection())gdal.ReprojectImage(dataset,tods,dataset.GetProjection(),tods.GetProjection(),GRA_Cubic)
　　但这个生成的不是缩略图，而是从大图中裁剪出来的一部分
　　于是乎，还没找到方法，先探索到这里，等日后再议了