python 使用GDAL实现栅格tif转矢量shp的方式小结

前言

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

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

核心函数

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，非常大，打开非常卡

  • 使用Polygonize

合并之后的矢量数据有48MB，相对第一种方法数据量大大减少

到此这篇关于python 使用GDAL实现栅格tif转矢量shp的文章就介绍到这了,更多相关python栅格tif转矢量shp内容请搜索我们以前的文章或继续浏览下面的相关文章希望大家以后多多支持我们！

赞 (0)

相关推荐

  • 在Python中用GDAL实现矢量对栅格的切割实例

    概述: 本文讲述如何在Python中用GDAL实现根据输入矢量边界对栅格数据的裁剪. 效果: 裁剪前 矢量边界 裁剪后 实现代码: # -*- coding: utf-8 -*- """ @author lzugis @date 2017-06-02 @brief 利用shp裁剪影像 """ from osgeo import gdal, gdalnumeric, ogr from PIL import Image, ImageDraw impo

  • python 矢量数据转栅格数据代码实例

    这篇文章主要介绍了python 矢量数据转栅格数据代码实例,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友可以参考下 投影包osr与proj4的使用 osr投影转换示例 from osgeo import osr,ogr #定义投影 #wgs84 source=osr.SpatialReference() source.ImportFromEPSG(4326) #google target=osr.SpatialReference() target.Imp

  • python安装gdal的两种方法

    1.不用手动下载文件,直接执行以下命令即可 conda install gdal 2.首先,下载gdal的whl文件  链接, 官网下载比较慢,GDAL-2.2.4-cp27-cp27m-win_amd64.whl 链接: https://pan.baidu.com/s/1prPHLJKwoKK505i5qTVZ7g 提取码: egj6 有百度云可以下载,然后放入本机目录. 这里目录有两种,一是放入anaconda安装目录的Scripts目录,我的是D:\anaconda\Scripts目录:二

  • 利用python GDAL库读写geotiff格式的遥感影像方法

    如下所示: from osgeo import gdal import numpy as np def read_tiff(inpath): ds=gdal.Open(inpath) row=ds.RasterXSize col=ds.RasterYSize band=ds.RasterCount geoTransform=ds.GetTransform() proj=ds.GetTransform() data=np.zeros([row,col,band]) for i in range(b

  • python使用gdal对shp读取,新建和更新的实例

    昨天要处理一个shp文件,读取里面的信息,做个计算然后写到后面新建的field里面.先写个外面网上都能找到的新建和读取吧. 1.读取shp文件 #-*- coding: cp936 -*- try: from osgeo import gdal from osgeo import ogr exceptImportError: import gdal import ogr defReadVectorFile(): # 为了支持中文路径,请添加下面这句代码 gdal.SetConfigOption(

  • python 使用GDAL实现栅格tif转矢量shp的方式小结

    python 使用GDAL实现栅格tif转矢量shp的方式小结

    前言 目前有一张tif格式的栅格影像,需要在web地图上进行展示,使用动态切片WMS的方式,渲染速度比较慢,而且大的时候会出现模糊的问题.并且后面需要做多期影像的切换,渲染与加载效率也值得关注. 计划是使用栅格转矢量的方式,将栅格数据转为矢量shp文件,然后进行矢量切片,使用Mapbox进行前端动态渲染.在网上查询了很多资料,有人说使用d3-contour在node.js中生成或者使用rasterio在python中进行转换,整体过程都比较麻烦,很不易实现.最终选定了使用GDAL进行栅格转矢量的

  • python与C、C++混编的四种方式(小结)

    混编的含义有两种, 一种是在python里面写C 一种是C里面写python 本文主要是进行简化,方便使用. ##################################################################################################### 第一种.Python调用C动态链接库(利用ctypes) pycall.c /***gcc -o libpycall.so -shared -fPIC pycall.c*/ #inclu

  • GDAL 矢量属性数据修改方式(python)

    Case:需要给一个现有的shp数据创建一个字段,并将属性表中原有的一个文本类型的属性转换为整型后填入新创建的字段. Problem:新字段创建成功,但是赋值操作无效,即无法成功给字段写入值. solution:对字段进行赋值后需要,重新写入Feature,否则赋值无效,即layer0.SetFeature(feature). 特别注意:在对数据进行读写操作,一定要以读写的方式打开,即Open(filePath,1),该方法的原型为Open(pszName,int bUpdate = false

  • Python叠加两幅栅格图像的实现方法

    目的 现有两幅栅格图像,一个是某地区道路栅格图,一个是某地区土地利用类型图,需要将道路叠加到土地利用类型图中,即叠加后,重合的像元值以道路图为准,其余的像元值仍是土地利用类型图原有的像元值. 图1 道路信息图 图2 土地利用类型图 图3 结果图 具体实现 from gdalconst import * from osgeo import gdal import osr import sys import copy #叠加两个栅格图像(一个道路栅格图,一个土地利用类型图),两幅图像重叠的像元值都是

  • 使用Python和GDAL给图片加坐标系的实现思路(坐标投影转换)

    使用Python和GDAL给图片加坐标系的实现思路(坐标投影转换)

    ** 使用Python和GDAL给图片加坐标系 ** 假设你已经知道arcgis地理配准(如下图内容),懂一点python. ** -目的和背景 1.从地图网站获得一张PNG格式的截图,已知坐标系为WGS84和左上角坐标.arcgis地理配准再定义投影即可给它加上原图的坐标系. 2.假设有上千张图片,可用Python和GDAL给图片加坐标系. -实现思路 1.使用GDAL需要知道待投影图片的地理坐标信息.仿射矩阵参数. 仿射矩阵参数是干什么的?见:https://zhuanlan.zhihu.c

  • 基于Python实现nc批量转tif格式

    基于Python实现nc批量转tif格式

    由于做项目需要运用到netCDF格式的气象数据,而ArcGIS中需要用栅格影像进行处理,对于较多的文件,ArcGIS一个个手动转换过于繁琐,因此我们采用Python进行转换,当然也可以采用matlab进行转换. 首先需要安装下面几个库: import os import netCDF4 as nc import numpy as np from osgeo import gdal, osr, ogr import glob 我们可以在下面网址中寻找对应python安装版本的安装包,下载后,在py

  • Anaconda下Python中GDAL模块的下载与安装过程

    本文介绍在Anaconda环境下,安装Python中栅格.矢量等地理数据处理库GDAL的方法. 需要注意的是,本文介绍基于conda install命令直接联网安装GDAL库的方法:这一方法有时不太稳定,且速度较慢.因此,如果有需要,大家可以参考Anaconda环境GDAL库基于whl文件的配置方法(https://www.jb51.net/article/280638.htm)这篇文章中的方法,可以更快速地配置GDAL库. 首先,我们打开“Anaconda Prompt (Anaconda)”

随机推荐