在Python中使用GDAL为卫片制作缩略图
在Python中使用GDAL为卫片制作缩略图
卫片和航片文件一般都很大,
有时候需要制作一个缩略图,放在导航窗体中做底图
这个文件不用太清晰,只要反映整体面貌即可
在ArcMap,PCI和Erdas中,可以直接截屏,或者使用重投影等功能很方便地完成,
最近在用python和GDAL,于是乎,查了查资料,并不复杂,
首先要确定python安装了numpy库,代码
- import gdal
- import numpy
- from gdalconst import *
- dataset = gdal.Open("F://数据//Raster//earth.img")
- width = dataset.RasterXSize
- height = dataset.RasterYSize
- bw = 0.1
- bh = 0.1
- x = 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 gdal
- import numpy
- from gdalconst import *
- dataset = gdal.Open("F://数据//Raster//earth.img")
- width = dataset.RasterXSize
- height = dataset.RasterYSize
- bw = 0.1
- bh = 0.1
- x = 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)
但这个生成的不是缩略图,而是从大图中裁剪出来的一部分
评论暂时关闭