目前有一张tif格式的栅格影像,需要在web地图上进行展示,使用动态切片WMS的方式,渲染速度比较慢,而且大的时候会出现模糊的问题。并且后面需要做多期影像的切换,渲染与加载效率也值得关注。

计划是使用栅格转矢量的方式,将栅格数据转为矢量shp文件,然后进行矢量切片,使用Mapbox进行前端动态渲染。在网上查询了很多资料,有人说使用 d3-contour 在node.js中生成或者使用 rasterio 在python中进行转换,整体过程都比较麻烦,很不易实现。最终选定了使用GDAL进行栅格转矢量的方法,代码比较简单。
原始tif影像(12.8MB)如下:
原始tif

GDAL中栅格转矢量的函数主要是以下两个,二者的参数没有任何区别,只是功能有区别:

  • FPolygonize (*args, **kwargs)

    FPolygonize(Band srcBand, Band maskBand, Layer outLayer, int iPixValField, char options=None, GDALProgressFunc callback=0, void * callback_data=None) -> int

    将每个像元转成一个矩形。

  • Polygonize (*args, **kwargs) **

    Polygonize(Band srcBand, Band maskBand, Layer outLayer, int iPixValField, char ** options=None, GDALProgressFunc callback=0, void * callback_data=None) -> int

    将每个像元转成一个矩形,然后将相似的像元进行合并。

    from osgeo import gdal, ogr, osr
    import os
    import datetime
    import numpy as np
    path = "Z_NAFP20210727.tif"
    if __name__ == '__main__':
        start_time = datetime.datetime.now()
        inraster = gdal.Open(path)  # 读取路径中的栅格数据
        inband = inraster.GetRasterBand(1)  # 这个波段就是最后想要转为矢量的波段,如果是单波段数据的话那就都是1
        prj = osr.SpatialReference()
        prj.ImportFromWkt(inraster.GetProjection())  # 读取栅格数据的投影信息,用来为后面生成的矢量做准备
        outshp = path[:-4] + ".shp"  # 给后面生成的矢量准备一个输出文件名,这里就是把原栅格的文件名后缀名改成shp了
        drv = ogr.GetDriverByName("ESRI Shapefile")
        if os.path.exists(outshp):  # 若文件已经存在,则删除它继续重新做一遍
            drv.DeleteDataSource(outshp)
        Polygon = drv.CreateDataSource(outshp)  # 创建一个目标文件
        Poly_layer = Polygon.CreateLayer(path[:-4], srs=prj, geom_type=ogr.wkbMultiPolygon)  # 对shp文件创建一个图层,定义为多个面类
        newField = ogr.FieldDefn('value', ogr.OFTReal)  # 给目标shp文件添加一个字段,用来存储原始栅格的pixel value,浮点型,
        Poly_layer.CreateField(newField)
        gdal.Polygonize(inband, None, Poly_layer, 0)  # 核心函数,执行的就是栅格转矢量操作
        # gdal.FPolygonize(inband, None, Poly_layer, 0)  # 只转矩形,不合并
        Polygon.SyncToDisk()
        Polygon = None
        end_time = datetime.datetime.now()
        print("Succeeded at", end_time)
        print("Elapsed Time:", end_time - start_time)  # 输出程序运行所需时间
    
    • 使用FPolygonize
      转换之后的矢量数据有270MB,非常大,打开非常卡
      FPolygonize转换效果

    • 使用Polygonize
      合并之后的矢量数据有48MB,相对第一种方法数据量大大减少
      Polygonize转换效果

    最近不少小伙伴问我怎么搞前端插值分析,我在github上查找了一些资料,目前最常用的方式是webgis框架+idw(反距离权重算法)+d3-contour的方式实现,这种方式是比较简单同时基本能满足一般的气象分析,比如局部的平均降水发布图等,下来我将详细的讲述具体的实现方式 GDAL是一个开源的地理工具包,其支持基本所有的地理操作,其有python、java、c等语言包,是地理信息C端开发不可越过的工具,鉴于python语言的简单性,这里使用pythonGDAL包来进行shp文件的生成,这里本质是利用ogc地理标准的坐标字符串来生成shp。 第一步:安装GDAL环境,建议下载后,本地安装,注意与python版本号要对应,可参考网上教程。 第二部:代码分析 引入GDA... 直接到下列链接下载即可,按照说明,将lib目录添加的系统环境变量中即可windows下GDAL322的库-深度学习文档类资源-CSDN下载包含目录中设置include目录 库目录中设置lib的路径 附加依赖项中设置gdal_i.lib GDAL读取数钱需要注册一下驱动(用于编码解码图像的驱动),同时可以设置一下支持中文路径。加载数据时需要注意,GA_Update和GA_ReadOnly两种模式。 2.1 获取图像的尺寸 这里所有的波段的size都是一样的 2.2 获取图像的通道数