gdal.py

gdal.py封装了所需的与gdal相关的操作.主要用来进行切割前的处理.
主要操作有投影变换与抽取波段:

def preprocess(inputfile, outputfile, band=None, spatial_ref=None,
               resampling=None, compress=None, **kwargs):
    # 所有要进行的预处理动作
    functions = []
    dataset = Dataset(inputfile)
  # 指定了抽取的波段或栅格数大于1就要进行波段抽取
    if band is not None and dataset.RasterCount > 1:
    # 制作偏函数
    # 因为extract_color_band这个函数需要band和inputfile两个参数
    # 制作成将band参数打包的偏函数,最终只需infutfile一个参数即可.
    functions.append(
        ('Extracting band {0}'.format(band),
         partial(extract_color_band, band=band))
    )

    # 如果需要进行投影转换
    if spatial_ref is not None and \
            dataset.GetSpatialReference() != spatial_ref:
        # 制作只需要infutfile为参数的偏函数
        functions.append(
            ('Reprojecting to EPSG:{0}'.format(spatial_ref.GetEPSGCode()),
             partial(warp,
                     spatial_ref=spatial_ref, resampling=resampling))
        )

    if not functions:
        # 不需要转换的话只需把源文件和临时文件连接起来就行了,不作处理
        rmfile(outputfile, ignore_missing=True)
        srcfile = os.path.relpath(inputfile, os.path.dirname(outputfile))
        os.symlink(srcfile, outputfile)
        return inputfile
    # 需要转换的话就进行转换
    return pipeline(inputfile=inputfile, outputfile=outputfile,
                    functions=functions, compress=compress, **kwargs)


def pipeline(inputfile, outputfile, functions, **kwargs):
    tmpfiles = []
    try:
        previous = inputfile
        # 遍历所有设定好的预处理动作
        for name, f in functions:
          # 每个动作都会生成vrt文件
            vrt = f(previous)
            # 获取临时产生的文件的路径
            current = vrt.get_tempfile(suffix='.vrt', prefix='gdal')
            tmpfiles.append(current)
            previous = current.name
        # 最终执行所有的vrt文件
        return vrt.render(outputfile=outputfile, **kwargs)
    finally:
      # 执行完毕清理临时文件
        for f in tmpfiles:
            f.close()

这其中有两个值得拿出来说一下:

• 偏函数封装
• GDAL的VRT模式

偏函数

偏函数,相当于把函数和该函数的部分或全部调用参数打包到一个新函数中去,这样调用新函数时就能用更少的参数,甚至不需要参数直接调用.
在上面的代码中,制作了两个偏函数:

• extract_color_band(inputfile, band)
• warp(inputfile, spatial_ref=None, cmd=GDALWARP, resampling=None, maximum_resolution=None)

可以发现,他们拥有相同且作用独立于其他参数的参数"inputfile".
制作成偏函数的他们,将其他调用函数打包进偏函数内部后,都形成了一个只接受inputfile的新函数,这样就可以在同一个流程内统一调度执行.

GDAL的VRT

为什么一定要统一执行呢?因为这里使用了GDAL的VRT文件模式.

VRT模式是GDAL的延迟处理模式.对于一系列的GDAL操作,它延迟了所有计算密集型处理,只会生成一个vrt描述文件,记录要进行什么操作,只有等需要处理时,一并处理.这种工作流式的调度不大会减少处理时间,但简化了操作流程,也大大节省了硬盘空间.

以extract_color_band为例:

def extract_color_band(inputfile, band):
    dataset = Dataset(inputfile)
  # 待执行的命令
    command = [
        'gdal_translate',
        '-q', # 安静模式
        '-of', 'VRT', # 输出到VRT
        '-b', band, # 指定波段
        inputfile,
        '/vsistdout'
    ]
    # 返回执行后的VRT文件
    return VRT(check_output_gdal([str(e) for e in command]))


def check_output_gdal(*popenargs, **kwargs):
  """
  执行所有命令
  """
    p = Popen(stderr=PIPE, stdout=PIPE, *popenargs, **kwargs)
    stdoutdata, stderrdata = p.communicate()
    return stdoutdata

helper.py
把相关的模块梳理的差不多了,可以进入整个g2m的调度中心helper.py,而helper.py的核心,就是warp_mbtiles:

def warp_mbtiles(inputfile, outputfile, metadata, colors=None, band=None,
                 spatial_ref=None, resampling=None,
                 min_resolution=None, max_resolution=None, fill_borders=None,
                 zoom_offset=None, renderer=None, pngdata=None):
    # 生成临时文件,存储进行过预处理后的源文件
    with NamedTemporaryFile(suffix='.tif') as tempfile:
        # gdal加载源文件,读取相关信息
        dataset = Dataset(inputfile)
        # 预处理(投影变换/波段提取)
        warped = preprocess(inputfile=inputfile, outputfile=tempfile.name,
                            band=band, spatial_ref=spatial_ref,
                            resampling=resampling, compress='LZW')
        # 制作偏函数
        preprocessor = partial(resample_after_warp,
                               whole_world=dataset.IsWholeWorld())
        # 切割瓦片到mbtiles中
        return image_mbtiles(inputfile=warped, outputfile=outputfile,
                             metadata=metadata,
                             min_resolution=min_resolution,
                             max_resolution=max_resolution,
                             colors=colors, renderer=renderer,
                             preprocessor=preprocessor,
                             fill_borders=fill_borders,
                             zoom_offset=zoom_offset,
                             pngdata=pngdata)


def resample_after_warp(pyramid, colors, whole_world, **kwargs):
    resolution = pyramid.dataset.GetNativeResolution()
    if whole_world:
        # 如果当前的地图范围非常大,覆盖了整个地球,则需要重采样,因为gdal有时会将其投影到一个局部投影下
        pyramid.dataset.resample_to_world()
    else:
        pyramid.dataset.resample(resolution=resolution)
    colorize(pyramid=pyramid, colors=colors)
    pyramid.dataset.align_to_grid(resolution=resolution)
    return pyramid

def image_mbtiles(inputfile, outputfile, metadata,
                  min_resolution=None, max_resolution=None, fill_borders=None,
                  zoom_offset=None, colors=None, renderer=None,
                  preprocessor=None, pngdata=None):

    if pngdata is None:
        pngdata = dict()
    if renderer is None:
        renderer = PngRenderer(**pngdata)
  # 创建mbtiles
    with MbtilesStorage.create(filename=outputfile,
                               metadata=metadata,
                               zoom_offset=zoom_offset,
                               renderer=renderer) as storage:
        # 创建切割器
        pyramid = TmsPyramid(inputfile=inputfile,
                             storage=storage,
                             min_resolution=min_resolution,
                             max_resolution=max_resolution)
        # 切割前进行颜色处理
        if preprocessor is None:
            preprocessor = colorize
        pyramid = preprocessor(**locals())
    # 正式切割瓦片
        pyramid.slice(fill_borders=fill_borders)
        # 切割完成后补充mbtiles的元数据
        if zoom_offset is None:
            zoom_offset = 0
        if min_resolution is None:
            min_resolution = pyramid.resolution
        if max_resolution is None:
            max_resolution = pyramid.resolution
        metadata = storage.mbtiles.metadata
        metadata['x-minzoom'] = min_resolution + zoom_offset
        metadata['x-maxzoom'] = max_resolution + zoom_offset

可以重新总结一下全部的流程了:

核心流程.png

总结

看完源码回头再审视g2m,依然能感觉到深深地无力.主要有3点:

• gdal与vips的运用
• python中比较高级的方法
• 数字图形学的知识

这3点的熟练运用可以看出g2m的作者非常的专业.
当然,也有些问题:

• 没有进度的反馈
• 单图像体积过大/切得瓦片数量极多的时候,使用时间可能呈指数级别增加.

当然,在小图片/少数量瓦片时,g2m的表现相当令人满意.实测1张500Mb的rgb波段tiff地图,在i3-3220-8G-500G机械硬盘的测试环境下,切割20级以下的时候,基本能1分钟切100Mb的成果(因为没有实时反馈,所以是按照mbtiles文件的体积增加量来计算的).