在Python中使用GDAL為衛片制作縮略圖

衛片和航片文件一般都很大，

有時候需要制作一個縮略圖，放在導航窗體中做底圖

這個文件不用太清晰，只要反映整體面貌即可

在ArcMap，PCI和Erdas中，可以直接截屏，或者使用重投影等功能很方便地完成，

最近在用python和GDAL，於是乎，查了查資料，並不復雜，

首先要確定python安裝了numpy庫，代碼

1. import gdal
2. import numpy
3. from gdalconst import *
4. dataset = gdal.Open("F://數據//Raster//earth.img")
5. width = dataset.RasterXSize
6. height = dataset.RasterYSize
7. bw = 0.1
8. bh = 0.1
9. x = int(width*bw)
10. y = int(height*bh)
11. datas = []
12. for i in range(3):
13. band = dataset.GetRasterBand(i+1)
14. data = band.ReadAsArray(0,0,width,height,x,y)
15. datas.append(numpy.reshape(data,(1,-1)))
16. datas = numpy.concatenate(datas)
17. driver = gdal.GetDriverByName("GTiff")
18. tods = driver.Create("F://數據//Raster//tpix1.tif",x,y,3,options=["INTERLEAVE=PIXEL"])
19. tods.WriteRaster(0,0,x,y,datas.tostring(),x,y,band_list=[1,2,3])

以上的這個方法，用到了numpy的矩陣運算，不過使用band的ReadAsArray

方法取得的數據，采用的是最鄰近采樣，雖然縮略圖的要求不高，但是這個

采樣方式總讓人不爽，查了資料，使用gdal.ReprojectImage

可以設置采樣方式，代碼如下：

1. # -*- coding: cp936 -*-
2. import gdal
3. import numpy
4. from gdalconst import *
5. dataset = gdal.Open("F://數據//Raster//earth.img")
6. width = dataset.RasterXSize
7. height = dataset.RasterYSize
8. bw = 0.1
9. bh = 0.1
10. x = int(width*bw)
11. y = int(height*bh)
12. driver = gdal.GetDriverByName("GTiff")
13. tods = driver.Create("F://數據//Raster//tpix1.tif",x,y,3,options=["INTERLEAVE=PIXEL"])
14. tods.SetGeoTransform(dataset.GetGeoTransform())
15. tods.SetProjection(dataset.GetProjection())
16. gdal.ReprojectImage(dataset,tods,dataset.GetProjection(),tods.GetProjection(),GRA_Cubic)

但這個生成的不是縮略圖，而是從大圖中裁剪出來的一部分