基于Python-GDAL的遥感影像镶嵌脚本

#!/usr/bin/env python ############################################################################### # $Id$ # # Project: InSAR Peppers # Purpose: Module to extract data from many rasters into one output. # Author: Frank Warmerdam, warmerdam@pobox.com # ############################################################################### # Copyright (c) 2000, Atlantis Scientific Inc. (www.atlsci.com) # Copyright (c) 2009-2011, Even Rouault <even dot rouault at mines-paris dot org> # # This library is free software; you can redistribute it and/or # modify it under the terms of the GNU Library General Public # License as published by the Free Software Foundation; either # version 2 of the License, or (at your option) any later version. # # This library is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # Library General Public License for more details. # # You should have received a copy of the GNU Library General Public # License along with this library; if not, write to the # Free Software Foundation, Inc., 59 Temple Place - Suite 330, # Boston, MA 02111-1307, USA. ############################################################################### # changes 29Apr2011 # If the input image is a multi-band one, use all the channels in # building the stack. # anssi.pekkarinen@fao.org import math import sys import time from osgeo import gdal try: progress = gdal.TermProgress_nocb except: progress = gdal.TermProgress __version__ = '$id$'[5:-1] verbose = 0 quiet = 0 # ============================================================================= def raster_copy( s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n, t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n, nodata=None ): if verbose != 0: print('Copy %d,%d,%d,%d to %d,%d,%d,%d.' % (s_xoff, s_yoff, s_xsize, s_ysize, t_xoff, t_yoff, t_xsize, t_ysize )) if nodata is not None: return raster_copy_with_nodata( s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n, t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n, nodata ) s_band = s_fh.GetRasterBand( s_band_n ) m_band = None # Works only in binary mode and doesn't take into account # intermediate transparency values for compositing. if s_band.GetMaskFlags() != gdal.GMF_ALL_VALID: m_band = s_band.GetMaskBand() elif s_band.GetColorInterpretation() == gdal.GCI_AlphaBand: m_band = s_band if m_band is not None: return raster_copy_with_mask( s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n, t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n, m_band ) s_band = s_fh.GetRasterBand( s_band_n ) t_band = t_fh.GetRasterBand( t_band_n ) data = s_band.ReadRaster( s_xoff, s_yoff, s_xsize, s_ysize, t_xsize, t_ysize, t_band.DataType ) t_band.WriteRaster( t_xoff, t_yoff, t_xsize, t_ysize, data, t_xsize, t_ysize, t_band.DataType ) return 0 # ============================================================================= def raster_copy_with_nodata( s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n, t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n, nodata ): try: import numpy as Numeric except ImportError: import Numeric s_band = s_fh.GetRasterBand( s_band_n ) t_band = t_fh.GetRasterBand( t_band_n ) data_src = s_band.ReadAsArray( s_xoff, s_yoff, s_xsize, s_ysize, t_xsize, t_ysize ) data_dst = t_band.ReadAsArray( t_xoff, t_yoff, t_xsize, t_ysize ) nodata_test = Numeric.equal(data_src,nodata) to_write = Numeric.choose( nodata_test, (data_src, data_dst) ) t_band.WriteArray( to_write, t_xoff, t_yoff ) return 0 # ============================================================================= def raster_copy_with_mask( s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n, t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n, m_band ): try: import numpy as Numeric except ImportError: import Numeric s_band = s_fh.GetRasterBand( s_band_n ) t_band = t_fh.GetRasterBand( t_band_n ) data_src = s_band.ReadAsArray( s_xoff, s_yoff, s_xsize, s_ysize, t_xsize, t_ysize ) data_mask = m_band.ReadAsArray( s_xoff, s_yoff, s_xsize, s_ysize, t_xsize, t_ysize ) data_dst = t_band.ReadAsArray( t_xoff, t_yoff, t_xsize, t_ysize ) mask_test = Numeric.equal(data_mask, 0) to_write = Numeric.choose( mask_test, (data_src, data_dst) ) t_band.WriteArray( to_write, t_xoff, t_yoff ) return 0 # ============================================================================= def names_to_fileinfos( names ): """ Translate a list of GDAL filenames, into file_info objects. names -- list of valid GDAL dataset names. Returns a list of file_info objects. There may be less file_info objects than names if some of the names could not be opened as GDAL files. """ file_infos = [] for name in names: fi = file_info() if fi.init_from_name( name ) == 1: file_infos.append( fi ) return file_infos # ***************************************************************************** class file_info: """A class holding information about a GDAL file.""" def init_from_name(self, filename): """ Initialize file_info from filename filename -- Name of file to read. Returns 1 on success or 0 if the file can't be opened. """ fh = gdal.Open( filename ) if fh is None: return 0 self.filename = filename self.bands = fh.RasterCount self.xsize = fh.RasterXSize self.ysize = fh.RasterYSize self.band_type = fh.GetRasterBand(1).DataType self.projection = fh.GetProjection() self.geotransform = fh.GetGeoTransform() self.ulx = self.geotransform[0] self.uly = self.geotransform[3] self.lrx = self.ulx + self.geotransform[1] * self.xsize self.lry = self.uly + self.geotransform[5] * self.ysize ct = fh.GetRasterBand(1).GetRasterColorTable() if ct is not None: self.ct = ct.Clone() else: self.ct = None return 1 def report( self ): print('Filename: '+ self.filename) print('File Size: %dx%dx%d' % (self.xsize, self.ysize, self.bands)) print('Pixel Size: %f x %f' % (self.geotransform[1],self.geotransform[5])) print('UL:(%f,%f) LR:(%f,%f)' % (self.ulx,self.uly,self.lrx,self.lry)) def copy_into( self, t_fh, s_band = 1, t_band = 1, nodata_arg=None ): """ Copy this files image into target file. This method will compute the overlap area of the file_info objects file, and the target gdal.Dataset object, and copy the image data for the common window area. It is assumed that the files are in a compatible projection ... no checking or warping is done. However, if the destination file is a different resolution, or different image pixel type, the appropriate resampling and conversions will be done (using normal GDAL promotion/demotion rules). t_fh -- gdal.Dataset object for the file into which some or all of this file may be copied. Returns 1 on success (or if nothin