How it works(14) GDAL2Tiles源码阅读

本文主要是介绍How it works(14) GDAL2Tiles源码阅读，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

引入

gdal2tiles(以下简称g2t),这个历史悠久的切图脚本依然能发挥其功用,因为它稳定的做好了它应做的东西.相比前面说过的gdal2mbtiles(以下简称g2m),我倒是更喜欢它,单文件脚本,运行只安装一个GDAL库足矣.同样因为有了g2m,我也是带着对比的心态提出几个问题:

从表现来看,g2t更慢
- 慢的原因是什么
- 可以采用g2m加速吗
与g2m对比,其算法有何差异

精简

原始的g2t脚本近3000行,包含了详细的注释和一些其他的功能,分析起来会产生干扰,因此我精简掉全部注释和我所用不上的功能,保留了核心功能,方便分析和使用.

相比原来精简掉的内容:

与KML输出相关的功能和变量
与生成HTML脚本相关的功能和变量
只保留了生成web墨卡托投影瓦片的功能
只保留了"average"一种重采样算法
限定只能生成png格式的瓦片

相比原来修改的内容:

将相关的功能整合进类中
修改命名规范为小驼峰

以下是修改后的脚本,不到500行,下面就基于这个修改后的核心脚本进行分析

from xml.etree import ElementTree
import json
from osgeo import gdal, osr
from uuid import uuid4
import sys
import shutil
import tempfile
import os
import math
from multiprocessing import Pool, Process, ManagerMAXZOOMLEVEL = 24
TILESIZE = 256
TILEDRIVER = '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 = {"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的调用可以发现其实际运行流程如下

# 打开数据
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 = {"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)

这里最重要的部分就是坐偏移量的计算,生成的几组变量在后面切图的时候都会用到:

控制如何在源图上找到瓦片所需位置的数据
- rx/ry:该张瓦片的左上角相对于原始tif的像素偏移
- rxsize/rysize:该张瓦片的长/宽对应于原始tif多少像素.不过由于瓦片都是正方形且是固定值,所以这两个值在一般情况下都相等且等于同一值
控制如何把对应的数据绘制到瓦片上.有个绘图
- wx/wy:窗口偏移,大多数情况为0
- wxsize/wysize:实际绘制时的大小,大多数情况为4倍瓦片尺寸

特殊情况:在所有tif的边缘处,都会出现瓦片冗余,这些地方需要单独处理:

x最小的瓦片,rx会小于0,即取到不在tif上的位置,需要将rx置零,缩小rxsize到有数据的范围,保证只在有数据的地方取数据.同时调整wx,wxsize,将数据绘制在对应的地方
x最大的瓦片,rx+rxsize会超过tif本身宽度,需要缩小rxsize,同样保证只在有数据的地方取数据.调整wxsize,只在有数据的地方绘制,其他无数据的地方已经在前面设为透明了,无需绘制.

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 1
tileDetails = 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()

切瓦片的时候充分利用核心资源无异是非常合算的.在没有阅读之前,我构思的多线程是固定若干个处理进程,然后给每个进程单独分配任务.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)