14 How it works GDAL2Tiles源码阅读 _生成

引入
(以下简称g2t),这个历史悠久的切图脚本依然能发挥其功用,因为它稳定的做好了它应做的东西.相比前面说过的(以下简称g2m),我倒是更喜欢它,单文件脚本,运行只安装一个GDAL库足矣.同样因为有了g2m,我也是带着对比的心态提出几个问题:
与g2m对比,其算法有何差异精简
原始的g2t脚本近3000行,包含了详细的注释和一些其他的功能,分析起来会产生干扰,因此我精简掉全部注释和我所用不上的功能,保留了核心功能,方便分析和使用.
相比原来精简掉的内容:
相比原来修改的内容:
以下是修改后的脚本,不到500行,下面就基于这个修改后的核心脚本进行分析
from xml.etree import ElementTreeimport jsonfrom osgeo import gdal, osrfrom uuid import uuid4import sysimport shutilimport tempfileimport osimport mathfrom multiprocessing import Pool, Process, ManagerMAXZOOMLEVEL = 24TILESIZE = 256TILEDRIVER = 'PNG'TILEEXT = 'png'class GlobalMercator(object):def __init__(self):self.initialResolution = 2 * math.pi * 6378137 / TILESIZEself.originShift = 2 * math.pi * 6378137 / 2.0def MetersToLatLon(self, mx, my):lon = (mx / self.originShift) * 180.0lat = (my / self.originShift) * 180.0lat = 180 / math.pi * (2 * math.atan(math.exp(lat * math.pi / 180.0)) - math.pi / 2.0)return lat, londef PixelsToMeters(self, px, py, zoom):res = self.Resolution(zoom)mx = px * res - self.originShiftmy = py * res - self.originShiftreturn mx, mydef MetersToPixels(self, mx, my, zoom):res = self.Resolution(zoom)px = (mx + self.originShift) / respy = (my + self.originShift) / resreturn px, pydef PixelsToTile(self, px, py):tx = int(math.ceil(px / float(TILESIZE)) - 1)ty = int(math.ceil(py / float(TILESIZE)) - 1)return tx, tydef MetersToTile(self, mx, my, zoom):px, py = self.MetersToPixels(mx, my, zoom)return self.PixelsToTile(px, py)def TileBounds(self, tx, ty, zoom):minx, miny = self.PixelsToMeters(tx * TILESIZE, ty * TILESIZE, zoom)maxx, maxy = self.PixelsToMeters((tx + 1) * TILESIZE, (ty + 1) * TILESIZE, zoom)return (minx, miny, maxx, maxy)def Resolution(self, zoom):return self.initialResolution / (2**zoom)def check(status, message):if status:sys.stderr.write("运行出错: %s\n" % message)sys.exit(3)class TileDetail(object):tx = 0ty = 0tz = 0rx = 0ry = 0rxsize = 0rysize = 0wx = 0wy = 0wxsize = 0wysize = 0def __init__(self, **kwargs):for key in kwargs:if hasattr(self, key):setattr(self, key, kwargs[key])class TileJobInfo(object):srcFile = ""nbDataBands = 0outputFilePath = ""tminmax = []tminz = 0tmaxz = 0outGeoTrans = []def __init__(self, **kwargs):for key in kwargs:if hasattr(self, key):setattr(self, key, kwargs[key])class TileJobsMaker(object):def __init__(self, inputFile, outputFolder, options):self.dataBandsCount = 4self.vrtFilename = os.path.join(tempfile.mkdtemp(), str(uuid4()) + '.vrt')self.inputFile = inputFileself.outputFolder = outputFolderself.options = optionsminmax = self.options.zoom.split('-', 1)minmax.extend([''])zoom_min, zoom_max = minmax[:2]self.tminz = int(zoom_min)if zoom_max:self.tmaxz = int(zoom_max)else:self.tmaxz = int(zoom_min)def updateNoDataValue(self):def gdalVrtWarp(options, key, value):tb = ElementTree.TreeBuilder()tb.start("Option", {"name": key})tb.data(value)tb.end("Option")elem = tb.close()options.insert(0, elem)tempFile = tempfile.mktemp('-TileJobsMaker.vrt')self.warpedDataset.GetDriver().CreateCopy(tempFile, self.warpedDataset)with open(tempFile, 'r', encoding='utf-8') as f:vrtString = f.read()vrtRoot = ElementTree.fromstring(vrtString)options = vrtRoot.find("GDALWarpOptions")gdalVrtWarp(options, "INIT_DEST", "NO_DATA")gdalVrtWarp(options, "UNIFIED_SRC_NODATA", "YES")vrtString = ElementTree.tostring(vrtRoot).decode()with open(tempFile, 'w') as f:f.write(vrtString)correctedDataset = gdal.Open(tempFile)os.unlink(tempFile)correctedDataset.SetMetadataItem('NODATA_VALUES', '0 0 0 0')self.warpedDataset = correctedDatasetdef openData(self):gdal.AllRegister()inputDataset = gdal.Open(self.inputFile, gdal.GA_ReadOnly)check(not inputDataset, "数据无法打开")check(inputDataset.RasterCount == 0, "数据无波段")GetGeoTransform = inputDataset.GetGeoTransform()gcpCount = inputDataset.GetGCPCount()check(GetGeoTransform == (0.0, 1.0, 0.0, 0.0, 0.0, 1.0) and gcpCount == 0, "数据缺少空间信息")inputSrs = osr.SpatialReference()inputSrs.ImportFromWkt(inputDataset.GetProjection())outputSrs = osr.SpatialReference()outputSrs.ImportFromEPSG(3857)self.warpedDataset = gdal.AutoCreateWarpedVRT(inputDataset,inputSrs.ExportToWkt(), outputSrs.ExportToWkt())self.updateNoDataValue()self.warpedDataset.GetDriver().CreateCopy(self.vrtFilename, self.warpedDataset)outGeotrans = self.warpedDataset.GetGeoTransform()check((outGeotrans[2], outGeotrans[4]) != (0, 0), "不支持变形后的数据")self.ominx = outGeotrans[0]self.omaxx = outGeotrans[0] + self.warpedDataset.RasterXSize * outGeotrans[1]self.omaxy = outGeotrans[3]self.ominy = outGeotrans[3] - self.warpedDataset.RasterYSize * outGeotrans[1]self.mercator = GlobalMercator()self.tminmax = list(range(0, 32))for tz in range(0, 32):tminx, tminy = self.mercator.MetersToTile(self.ominx, self.ominy, tz)tmaxx, tmaxy = self.mercator.MetersToTile(self.omaxx, self.omaxy, tz)tminx, tminy = max(0, tminx), max(0, tminy)tmaxx, tmaxy = min(2**tz - 1, tmaxx), min(2**tz - 1, tmaxy)self.tminmax[tz] = (tminx, tminy, tmaxx, tmaxy)def makeMetadata(self):south, west = self.mercator.MetersToLatLon(self.ominx, self.ominy)north, east = self.mercator.MetersToLatLon(self.omaxx, self.omaxy)south, west = max(-85.05112878, south), max(-180.0, west)north, east = min(85.05112878, north), min(180.0, east)metadata = http://www.kingceram.com/post/{"south": south, "north": north, "west": west, "east": east}with open(os.path.join(self.outputFolder, 'metadata.json'), 'w') as f:json.dump(metadata, f)def makeBaseTiles(self):tminx, tminy, tmaxx, tmaxy = self.tminmax[self.tmaxz]tileDetails = []tz = self.tmaxzfor ty in range(tmaxy, tminy - 1, -1):for tx in range(tminx, tmaxx + 1):tilefilename = os.path.join(self.outputFolder, str(tz), str(tx), "%s.%s" % (ty, TILEEXT))if not os.path.exists(os.path.dirname(tilefilename)):os.makedirs(os.path.dirname(tilefilename))b = self.mercator.TileBounds(tx, ty, tz)rb, wb = self.geoQuery(b[0], b[3], b[2], b[1])rx, ry, rxsize, rysize = rbwx, wy, wxsize, wysize = wbtileDetails.append(TileDetail(tx=tx,ty=ty,tz=tz,rx=rx,ry=ry,rxsize=rxsize,rysize=rysize,wx=wx,wy=wy,wxsize=wxsize,wysize=wysize,))conf = TileJobInfo(srcFile=self.vrtFilename,nbDataBands=self.dataBandsCount,outputFilePath=self.outputFolder,tminmax=self.tminmax,tminz=self.tminz,tmaxz=self.tmaxz,)return conf, tileDetailsdef geoQuery(self, ulx, uly, lrx, lry):ds = self.warpedDatasetgeotran = ds.GetGeoTransform()rx = int((ulx - geotran[0]) / geotran[1] + 0.001)ry = int((uly - geotran[3]) / geotran[5] + 0.001)rxsize = int((lrx - ulx) / geotran[1] + 0.5)rysize = int((lry - uly) / geotran[5] + 0.5)wxsize, wysize = 4 * TILESIZE, 4 * TILESIZEwx = 0if rx < 0:rxshift = abs(rx)wx = int(wxsize * (float(rxshift) / rxsize))wxsize = wxsize - wxrxsize = rxsize - int(rxsize * (float(rxshift) / rxsize))rx = 0if rx + rxsize > ds.RasterXSize:wxsize = int(wxsize * (float(ds.RasterXSize - rx) / rxsize))rxsize = ds.RasterXSize - rxwy = 0if ry < 0:ryshift = abs(ry)wy = int(wysize * (float(ryshift) / rysize))wysize = wysize - wyrysize = rysize - int(rysize * (float(ryshift) / rysize))ry = 0if ry + rysize > ds.RasterYSize:wysize = int(wysize * (float(ds.RasterYSize - ry) / rysize))rysize = ds.RasterYSize - ryreturn (rx, ry, rxsize, rysize), (wx, wy, wxsize, wysize)class ProgressBar(object):def __init__(self, total_items, title):sys.stdout.write("%s 共%d张 \n" % (title, total_items))self.total_items = total_itemsself.nb_items_done = 0self.current_progress = 0self.STEP = 2.5def start(self):sys.stdout.write("0")def updateProgress(self, nb_items=1):self.nb_items_done += nb_itemsprogress = float(self.nb_items_done) / self.total_items * 100if progress >= self.current_progress + self.STEP:done = Falsewhile not done:if self.current_progress + self.STEP <= progress:self.current_progress += self.STEPif self.current_progress % 10 == 0:sys.stdout.write(str(int(self.current_progress)))if self.current_progress == 100:sys.stdout.write("\n")else:sys.stdout.write(".")else:done = Truesys.stdout.flush()class SingleProcessTiling(object):def __init__(self, inputFile, outputFolder, options):self.inputFile = inputFileself.outputFolder = outputFolderself.options = optionsself.total = 0self.tiling()self.createOverviewTiles()shutil.rmtree(os.path.dirname(self.tileJobInfo.srcFile))def tiling(self):tileDetails = self.workerTileDetails()tilecount = len(tileDetails)self.total += tilecountprogressBar = ProgressBar(tilecount, '切割顶层瓦片')progressBar.start()for tileDetail in tileDetails:self.createBaseTile(tileDetail)progressBar.updateProgress()def createBaseTile(self, tileDetail, queue=None):gdal.AllRegister()tileJobInfo = self.tileJobInfooutput = tileJobInfo.outputFilePathtilebands = tileJobInfo.nbDataBandsds = gdal.Open(tileJobInfo.srcFile, gdal.GA_ReadOnly)memDrv = gdal.GetDriverByName('MEM')outDrv = gdal.GetDriverByName(TILEDRIVER)alphaband = ds.GetRasterBand(1).GetMaskBand()tx = tileDetail.txty = tileDetail.tytz = tileDetail.tzrx = tileDetail.rxry = tileDetail.ryrxsize = tileDetail.rxsizerysize = tileDetail.rysizewx = tileDetail.wxwy = tileDetail.wywxsize = tileDetail.wxsizewysize = tileDetail.wysizequerysize = 4 * TILESIZEtilefilename = os.path.join(output, str(tz), str(tx), "%s.%s" % (ty, TILEEXT))dstile = memDrv.Create('', TILESIZE, TILESIZE, tilebands)data = alpha = Noneif rxsize != 0 and rysize != 0 and wxsize != 0 and wysize != 0:data = ds.ReadRaster(rx, ry, rxsize, rysize, wxsize, wysize, band_list=list(range(1, tilebands)))alpha = alphaband.ReadRaster(rx, ry, rxsize, rysize, wxsize, wysize)if data:dsquery = memDrv.Create('', querysize, querysize, tilebands)dsquery.WriteRaster(wx, wy, wxsize, wysize, data, band_list=list(range(1, tilebands)))dsquery.WriteRaster(wx, wy, wxsize, wysize, alpha, band_list=[tilebands])self.scaleQueryToTile(dsquery, dstile, tilefilename)del dsquerydel dsdel dataoutDrv.CreateCopy(tilefilename, dstile, strict=0)del dstileif queue:queue.put("tile %s %s %s" % (tx, ty, tz))def workerTileDetails(self):tileJobsMaker = TileJobsMaker(self.inputFile, self.outputFolder, self.options)tileJobsMaker.openData()tileJobsMaker.makeMetadata()conf, tileDetails = tileJobsMaker.makeBaseTiles()self.tileJobInfo = confreturn tileDetailsdef scaleQueryToTile(self, dsquery, dstile, tilefilename=''):tilebands = dstile.RasterCountfor i in range(1, tilebands + 1):res = gdal.RegenerateOverview(dsquery.GetRasterBand(i), dstile.GetRasterBand(i), 'average')check(res != 0, "概览生成失败 %s,%d" % (tilefilename, res))def createOverviewTiles(self):tileJobInfo = self.tileJobInfomemDriver = gdal.GetDriverByName('MEM')outDriver = gdal.GetDriverByName(TILEDRIVER)tilebands = tileJobInfo.nbDataBandstcount = 0for tz in range(tileJobInfo.tmaxz - 1, tileJobInfo.tminz - 1, -1):tminx, tminy, tmaxx, tmaxy = tileJobInfo.tminmax[tz]tcount += (1 + abs(tmaxx - tminx)) * (1 + abs(tmaxy - tminy))if tcount == 0:returnself.total += tcountprogressBar = ProgressBar(tcount, '切割下层瓦片')progressBar.start()for tz in range(tileJobInfo.tmaxz - 1, tileJobInfo.tminz - 1, -1):tminx, tminy, tmaxx, tmaxy = tileJobInfo.tminmax[tz]for ty in range(tmaxy, tminy - 1, -1):for tx in range(tminx, tmaxx + 1):tilefilename = os.path.join(self.outputFolder, str(tz), str(tx), "%s.%s" % (ty, TILEEXT))if not os.path.exists(os.path.dirname(tilefilename)):os.makedirs(os.path.dirname(tilefilename))dsquery = memDriver.Create('', 2 * TILESIZE, 2 * TILESIZE, tilebands)dstile = memDriver.Create('', TILESIZE, TILESIZE, tilebands)for y in range(2 * ty, 2 * ty + 2):for x in range(2 * tx, 2 * tx + 2):minx, miny, maxx, maxy = tileJobInfo.tminmax[tz + 1]if x >= minx and x <= maxx and y >= miny and y <= maxy:path = os.path.join(self.outputFolder, str(tz + 1), str(x), "%s.%s" % (y, TILEEXT))dsquerytile = gdal.Open(path, gdal.GA_ReadOnly)if (ty == 0 and y == 1) or (ty != 0 and(y % (2 * ty)) != 0):tileposy = 0else:tileposy = TILESIZEif tx:tileposx = x % (2 * tx) * TILESIZEelif tx == 0 and x == 1:tileposx = TILESIZEelse:tileposx = 0tempRaseter = dsquerytile.ReadRaster(0, 0, TILESIZE, TILESIZE)dsquery.WriteRaster(tileposx, tileposy, TILESIZE, TILESIZE, tempRaseter, band_list=list(range(1, tilebands + 1)))self.scaleQueryToTile(dsquery, dstile, tilefilename=tilefilename)outDriver.CreateCopy(tilefilename, dstile, strict=0)progressBar.updateProgress()class MultiProcessTiling(SingleProcessTiling):def __init__(self, inputFile, outputFolder, options):super().__init__(inputFile, outputFolder, options)def tiling(self):processes = self.options.processes or 1tileDetails = self.workerTileDetails()tilecount = len(tileDetails)self.total += tilecountmanager = Manager()queue = manager.Queue()pool = Pool(processes=processes)for tileDetail in tileDetails:pool.apply_async(self.createBaseTile, (tileDetail, ), {"queue": queue})p = Process(target=self.progressPrinter, args=[queue, tilecount, '切割顶层瓦片'])p.start()pool.close()pool.join()p.join()def progressPrinter(self, queue, nb_jobs, title):pb = ProgressBar(nb_jobs, title)pb.start()for _ in range(nb_jobs):queue.get()pb.updateProgress()queue.task_done()def process_args(argv):import argparseparser = argparse.ArgumentParser()parser.add_argument("-z", dest='zoom', metavar='切割级别', required=True, help="如 '12-20'")parser.add_argument("-p", dest='processes', metavar='进程数', type=int, default=1, help='默认单进程')parser.add_argument("-i", dest='input', metavar='tif文件', required=True)parser.add_argument("-o", dest='output', metavar="输出文件夹", help="可选,默认在输入文件夹下的同名文件夹")args = parser.parse_args(argv)inputFile = args.inputcheck(not os.path.isfile(inputFile), "%s不存在或非文件" % inputFile)outputFolder = args.outputif not outputFolder:tifname = os.path.basename(inputFile).split('.')[0]outputFolder = os.path.join(os.path.dirname(os.path.abspath(inputFile)), tifname)if not os.path.exists(outputFolder):os.makedirs(outputFolder)return inputFile, outputFolder, argsdef main():argv = gdal.GeneralCmdLineProcessor(sys.argv)inputFile, outputFolder, options = process_args(argv[1:])if options.processes == 1:SingleProcessTiling(inputFile, outputFolder, options)else:MultiProcessTiling(inputFile, outputFolder, options)if __name__ == '__main__':main()
运行流程
g2t的流程没有分支,非常清晰:

文章插图
从原始数据获取所需的最高级别的瓦片,更低级的瓦片只需从这些最高级瓦片一层一层生成.
这样速度更快:因为最高级的瓦片只能利用gdal从原始tif中获取,其速度受tif尺寸影响很大,且从tif上取得级别越低,单次所取范围越大,速度也越慢.举个实际的例子,从原始tif上获取某位置17级的瓦片的时间将远远大于从原始tif获取4张对应位置的18级瓦片,并将其合成的时间.
整体架构 1. 生成切割顶级瓦片文物列表
生成任务列表的功能被完全封装在这个类中,观察的调用可以发现其实际运行流程如下

# 打开数据tileJobsMaker = TileJobsMaker(self.inputFile, self.outputFolder, self.options)tileJobsMaker.openData()# 生成元数据文件tileJobsMaker.makeMetadata()# 生成任务列表conf, tileDetails = tileJobsMaker.makeBaseTiles()

1.1 打开数据

def __init__(self, inputFile, outputFolder, options):# 默认只支持包含RGB波段的数据self.dataBandsCount = 4# 流程采用gdal的vrt模式,加快运行(类似g2m)self.vrtFilename = os.path.join(tempfile.mkdtemp(), str(uuid4()) + '.vrt')self.inputFile = inputFileself.outputFolder = outputFolderself.options = options# 格式化缩放范围minmax = self.options.zoom.split('-', 1)minmax.extend([''])zoom_min, zoom_max = minmax[:2]self.tminz = int(zoom_min)if zoom_max:self.tmaxz = int(zoom_max)else:self.tmaxz = int(zoom_min)def openData(self):gdal.AllRegister()inputDataset = gdal.Open(self.inputFile, gdal.GA_ReadOnly)# 强制将数据投影到web墨卡托投影,因为一般我们都是把瓦片发布为互联网服务的,3857无疑是最方便的inputSrs = osr.SpatialReference()inputSrs.ImportFromWkt(inputDataset.GetProjection())outputSrs = osr.SpatialReference()outputSrs.ImportFromEPSG(3857)# 投影.在这里,操作的数据其实已经变为vrt格式self.warpedDataset = gdal.AutoCreateWarpedVRT(inputDataset,inputSrs.ExportToWkt(), outputSrs.ExportToWkt())# 强制将nodata值设为透明self.updateNoDataValue()# 将vrt格式的数据集写入到指定的vrt文件中,供后期使用self.warpedDataset.GetDriver().CreateCopy(self.vrtFilename, self.warpedDataset)# 计算元数据的四至,因为已经投影为web墨卡托,这里的单位是米outGeotrans = self.warpedDataset.GetGeoTransform()self.ominx = outGeotrans[0]self.omaxx = outGeotrans[0] + self.warpedDataset.RasterXSize * outGeotrans[1]self.omaxy = outGeotrans[3]self.ominy = outGeotrans[3] - self.warpedDataset.RasterYSize * outGeotrans[1]# GlobalMercator是一个封装了各种投影方法的类self.mercator = GlobalMercator()self.tminmax = list(range(0, 32))# 计算每一缩放级,瓦片的行列号范围# 这里可以只计算用户指定的缩放范围,但其实影响不大for tz in range(0, 32):tminx, tminy = self.mercator.MetersToTile(self.ominx, self.ominy, tz)tmaxx, tmaxy = self.mercator.MetersToTile(self.omaxx, self.omaxy, tz)tminx, tminy = max(0, tminx), max(0, tminy)tmaxx, tmaxy = min(2**tz - 1, tmaxx), min(2**tz - 1, tmaxy)self.tminmax[tz] = (tminx, tminy, tmaxx, tmaxy)def updateNoDataValue(self):"""更新nodata的地方为透明"""def gdalVrtWarp(options, key, value):"""vrt文件修改"""tb = ElementTree.TreeBuilder()tb.start("Option", {"name": key})tb.data(value)tb.end("Option")elem = tb.close()options.insert(0, elem)tempFile = tempfile.mktemp('-TileJobsMaker.vrt')self.warpedDataset.GetDriver().CreateCopy(tempFile, self.warpedDataset)# 直接通过修改vrt文件的方式添加nodata值with open(tempFile, 'r', encoding='utf-8') as f:vrtString = f.read()vrtRoot = ElementTree.fromstring(vrtString)options = vrtRoot.find("GDALWarpOptions")# 设定数据集的每一个像素初始值都为no_datagdalVrtWarp(options, "INIT_DEST", "NO_DATA")# 当所有波段都符合no_data时,将整个波段都视为no_data,而不将每个波段独立对待gdalVrtWarp(options, "UNIFIED_SRC_NODATA", "YES")vrtString = ElementTree.tostring(vrtRoot).decode()with open(tempFile, 'w') as f:f.write(vrtString)# 加载修改后的vrt文件correctedDataset = gdal.Open(tempFile)os.unlink(tempFile)# 设置no_data值为透明(RGBA四个波段都是0)correctedDataset.SetMetadataItem('NODATA_VALUES', '0 0 0 0')self.warpedDataset = correctedDataset

1.2 生成元数据文件

def makeMetadata(self):# 以经纬度的方式计算瓦片的四至,供前端地图跳转使用south, west = self.mercator.MetersToLatLon(self.ominx, self.ominy)north, east = self.mercator.MetersToLatLon(self.omaxx, self.omaxy)south, west = max(-85.05112878, south), max(-180.0, west)north, east = min(85.05112878, north), min(180.0, east)metadata = http://www.kingceram.com/post/{"south": south, "north": north, "west": west, "east": east}with open(os.path.join(self.outputFolder, 'metadata.json'), 'w') as f:json.dump(metadata, f)

当然可以尽可能的将有用的信息写入元数据中.
1.3 生成任务列表

def makeBaseTiles(self):# 最高缩放级别的瓦片行列号范围tminx, tminy, tmaxx, tmaxy = self.tminmax[self.tmaxz]tileDetails = []tz = self.tmaxz# y倒序,从左上角开始切图for ty in range(tmaxy, tminy - 1, -1):for tx in range(tminx, tmaxx + 1):# tz,yx,ty对应着行列号中的x/y/ztilefilename = os.path.join(self.outputFolder, str(tz), str(tx), "%s.%s" % (ty, TILEEXT))if not os.path.exists(os.path.dirname(tilefilename)):os.makedirs(os.path.dirname(tilefilename))# 计算该瓦片的投影经纬度范围b = self.mercator.TileBounds(tx, ty, tz)# 获取该瓦片具体的各种偏移参数rb, wb = self.geoQuery(b[0], b[3], b[2], b[1])rx, ry, rxsize, rysize = rbwx, wy, wxsize, wysize = wbtileDetails.append(TileDetail(tx=tx,ty=ty,tz=tz,rx=rx,ry=ry,rxsize=rxsize,rysize=rysize,wx=wx,wy=wy,wxsize=wxsize,wysize=wysize,))conf = TileJobInfo(srcFile=self.vrtFilename,nbDataBands=self.dataBandsCount,outputFilePath=self.outputFolder,tminmax=self.tminmax,tminz=self.tminz,tmaxz=self.tmaxz,)return conf, tileDetailsdef geoQuery(self, ulx, uly, lrx, lry):ds = self.warpedDataset# geotran[0/3]是tif左上角点x/y# geotran[1/5]是像源宽/高geotran = ds.GetGeoTransform()# 计算该瓦片的左上角在源图上的x/y像素偏移量rx = int((ulx - geotran[0]) / geotran[1] + 0.001)ry = int((uly - geotran[3]) / geotran[5] + 0.001)# 计算该瓦片在源图上的像素宽度rxsize = int((lrx - ulx) / geotran[1] + 0.5)rysize = int((lry - uly) / geotran[5] + 0.5)# 窗口尺寸.4倍于瓦片尺寸,提高缩放重采样时的瓦片效果wxsize, wysize = 4 * TILESIZE, 4 * TILESIZE# 窗口偏移wx = 0# 特殊情况下的修正if rx < 0:rxshift = abs(rx)wx = int(wxsize * (float(rxshift) / rxsize))wxsize = wxsize - wxrxsize = rxsize - int(rxsize * (float(rxshift) / rxsize))rx = 0if rx + rxsize > ds.RasterXSize:wxsize = int(wxsize * (float(ds.RasterXSize - rx) / rxsize))rxsize = ds.RasterXSize - rxwy = 0if ry < 0:ryshift = abs(ry)wy = int(wysize * (float(ryshift) / rysize))wysize = wysize - wyrysize = rysize - int(rysize * (float(ryshift) / rysize))ry = 0if ry + rysize > ds.RasterYSize:wysize = int(wysize * (float(ds.RasterYSize - ry) / rysize))rysize = ds.RasterYSize - ryreturn (rx, ry, rxsize, rysize), (wx, wy, wxsize, wysize)

这里最重要的部分就是坐偏移量的计算,生成的几组变量在后面切图的时候都会用到:
控制如何把对应的数据绘制到瓦片上.有个绘图
特殊情况:在所有tif的边缘处,都会出现瓦片冗余,这些地方需要单独处理:
g2t本可以不生成任务列表,边算边切瓦片的.强制生成任务列表的目的只有一个:可以多进程切割.
2. 切割顶级瓦片

def createBaseTile(self, tileDetail, queue=None):gdal.AllRegister()# 获取任务参数tileJobInfo = self.tileJobInfooutput = tileJobInfo.outputFilePathtilebands = tileJobInfo.nbDataBands# 打开vrt文件作为数据集,读取模式,在多进程模式下不会出现锁?ds = gdal.Open(tileJobInfo.srcFile, gdal.GA_ReadOnly)memDrv = gdal.GetDriverByName('MEM')outDrv = gdal.GetDriverByName(TILEDRIVER)alphaband = ds.GetRasterBand(1).GetMaskBand()tx = tileDetail.txty = tileDetail.tytz = tileDetail.tzrx = tileDetail.rxry = tileDetail.ryrxsize = tileDetail.rxsizerysize = tileDetail.rysizewx = tileDetail.wxwy = tileDetail.wywxsize = tileDetail.wxsizewysize = tileDetail.wysize# '窗口'数据集尺寸就是4倍瓦片大小querysize = 4 * TILESIZEtilefilename = os.path.join(output, str(tz), str(tx), "%s.%s" % (ty, TILEEXT))# 最终要写入的数据集dstile = memDrv.Create('', TILESIZE, TILESIZE, tilebands)data = alpha = Noneif rxsize != 0 and rysize != 0 and wxsize != 0 and wysize != 0:# 根据上文获取到的参数,读取每一张瓦片对应的数据data = ds.ReadRaster(rx, ry, rxsize, rysize, wxsize, wysize, band_list=list(range(1, tilebands)))# alpha波段直接创建alpha = alphaband.ReadRaster(rx, ry, rxsize, rysize, wxsize, wysize)if data:# 所谓的'窗口'数据集dsquery = memDrv.Create('', querysize, querysize, tilebands)# 先将数据读取到窗口数据集中dsquery.WriteRaster(wx, wy, wxsize, wysize, data, band_list=list(range(1, tilebands)))dsquery.WriteRaster(wx, wy, wxsize, wysize, alpha, band_list=[tilebands])# 重采样到目标数据集self.scaleQueryToTile(dsquery, dstile, tilefilename)del dsquerydel dsdel data# 目标数据集导出png图片到指定位置outDrv.CreateCopy(tilefilename, dstile, strict=0)del dstile# 如果是多进程模式则向主进程传递进度if queue:queue.put("tile %s %s %s" % (tx, ty, tz))def scaleQueryToTile(self, dsquery, dstile, tilefilename=''):"""从'窗口'数据集使用'average'算法重采样到目标数据集"""tilebands = dstile.RasterCountfor i in range(1, tilebands + 1):res = gdal.RegenerateOverview(dsquery.GetRasterBand(i), dstile.GetRasterBand(i), 'average')

单进程调用:

# 产生任务tileDetails = self.workerTileDetails()# 每个任务串行执行for tileDetail in tileDetails:self.createBaseTile(tileDetail)

多进程调用:

processes = self.options.processes or 1tileDetails = self.workerTileDetails()manager = Manager()queue = manager.Queue()# 同时只允许指定个数的进程在运行pool = Pool(processes=processes)for tileDetail in tileDetails:pool.apply_async(self.createBaseTile, (tileDetail, ), {"queue": queue})pool.close()pool.join()

【14How it works GDAL2Tiles源码阅读】切瓦片的时候充分利用核心资源无异是非常合算的.在没有阅读之前,我构思的多线程是固定若干个处理进程,然后给每个进程单独分配任务.g2t的实现方式更加简单,实现也清晰.因为不断有进程创建/结束,我无法确定创建/结束进程的会耽误多长时间,不过实际来看应该影响不大.
3. 生成下层瓦片

 def createOverviewTiles(self):tileJobInfo = self.tileJobInfomemDriver = gdal.GetDriverByName('MEM')outDriver = gdal.GetDriverByName(TILEDRIVER)tilebands = tileJobInfo.nbDataBandsfor tz in range(tileJobInfo.tmaxz - 1, tileJobInfo.tminz - 1, -1):tminx, tminy, tmaxx, tmaxy = tileJobInfo.tminmax[tz]for ty in range(tmaxy, tminy - 1, -1):for tx in range(tminx, tmaxx + 1):# 遍历所有底层瓦片tilefilename = os.path.join(self.outputFolder, str(tz), str(tx), "%s.%s" % (ty, TILEEXT))if not os.path.exists(os.path.dirname(tilefilename)):os.makedirs(os.path.dirname(tilefilename))dsquery = memDriver.Create('', 2 * TILESIZE, 2 * TILESIZE, tilebands)dstile = memDriver.Create('', TILESIZE, TILESIZE, tilebands)# 每级的行列号都是它大一级的2分之1# 每张瓦片都与比它大一级的4张对应瓦片所示范围相同,所以根据这4张瓦片就能拼接出本级瓦片for y in range(2 * ty, 2 * ty + 2):for x in range(2 * tx, 2 * tx + 2):minx, miny, maxx, maxy = tileJobInfo.tminmax[tz + 1]# 只拼接有数据的瓦片if x >= minx and x <= maxx and y >= miny and y <= maxy:path = os.path.join(self.outputFolder, str(tz + 1), str(x), "%s.%s" % (y, TILEEXT))# 读取4张中的每一张瓦片dsquerytile = gdal.Open(path, gdal.GA_ReadOnly)# 把4张瓦片放到对应的位置if (ty == 0 and y == 1) or (ty != 0 and(y % (2 * ty)) != 0):tileposy = 0else:tileposy = TILESIZEif tx:tileposx = x % (2 * tx) * TILESIZEelif tx == 0 and x == 1:tileposx = TILESIZEelse:tileposx = 0# 读取瓦片再在'窗口'中绘制tempRaseter = dsquerytile.ReadRaster(0, 0, TILESIZE, TILESIZE)dsquery.WriteRaster(tileposx, tileposy, TILESIZE, TILESIZE, tempRaseter, band_list=list(range(1, tilebands + 1)))# 重采样后写入本地对应文件self.scaleQueryToTile(dsquery, dstile, tilefilename=tilefilename)outDriver.CreateCopy(tilefilename, dstile, strict=0)

切下层瓦片时并没有采用多进程模式,我猜可能是已经够快了,事实也确实如此.或许如果也采用多进程模式,整体运行还会更快?应在在瓦片数量比较多的情况下才有所体现吧.
总结
只用GDAL使得g2t部署简单,但也限制了其性能.如果把读取tif和生成png都替换成类似g2m的vips,或许能获得和g2m一样的效率.