当前位置：首页 > article >正文

【GIS教程】使用GDAL实现栅格转矢量（GeoJSON、Shapefile）- 附完整代码

article 2024/12/25 13:50:11

文章目录

一、应用场景
- 1、GeoJSON
- 2、ESRI Shapefile
- 3、GDAL
二、基本思路
- 1、数据准备
- 2、重投影（可选）
- 3、创建空的矢量图层
- 4、栅格转矢量
三、完整代码
四、总结
五、拓展（使用ArcGIS工具进行栅格转矢量）

一、应用场景

TIFF格式的影像数据包含了丰富的地理信息，但是直接使用栅格数据进行分析和规划存在一定的局限性，因为栅格数据不利于进行空间查询和属性分析。

因此我们可以将这些栅格TIFF数据转换为矢量格式，以便更好地进行空间分析和数据集成。

1、GeoJSON

GeoJSON是一种基于JSON（JavaScript Object Notation）的格式，用于在Web应用程序中编码各种地理数据结构。它能够描述几何、拓扑关系、属性和地理空间数据的其他特征，非常适合于矢量化处理和空间分析。

GeoJSON 的关键特性：

1、基于JSON： GeoJSON 是一种轻量级的数据交换格式，易于人阅读和编写，同时也易于机器解析和生成。
2、几何类型： GeoJSON支持多种几何类型，包括点（Point）、线（LineString）、线环（LinearRing）、多边形（Polygon）、多点（MultiPoint）、多线（MultiLineString）、多多边形（MultiPolygon）和几何集合（GeometryCollection）。
3、特征（Feature）： GeoJSON中的一个“特征”（Feature）是一个包含几何对象和属性的对象。特征可以代表地图上的一个实体，如一个城市或一条河流。
4、特征集合（FeatureCollection）： 特征集合是一组特征的集合，可以用来表示多个地理实体。
5、属性： 每个特征都可以有一个与之关联的属性对象，该对象存储了关于特征的信息，例如名称、类型或其他描述性数据。
6、 坐标参考系统： GeoJSON 对象可以包含一个“坐标参考系统”（CRS）属性，用来指定坐标系统，最常见的是WGS84（EPSG:4326）。

2、ESRI Shapefile

SHP（Shapefile）是一种广泛使用的矢量数据存储格式，由美国ESRI开发，它主要用于存储地理信息系统（GIS）中的点、线或多边形数据，一个完整的Shapefile由以下文件组成：

.shp文件：这是Shapefile的主文件，包含几何图形的几何信息，即空间坐标数据

.shx文件：这是Shapefile的索引文件，用于存储.shp文件中要素的位置，加快数据访问速度

.dbf文件：这是Shapefile的属性数据文件，存储每个几何图形的属性信息，如名称、类型等，采用dBase IV格式

除了这三个必须的文件外，Shapefile还可以包括以下可选文件：

.prj文件：包含地图投影的信息，定义了数据的坐标系和投影信息

.sbn和.sbx文件：这些是Shapefile的空间索引文件，提供了额外的空间查询性能，用于支持空间索引的快速查找

.xml文件：包含Shapefile的元数据信息

.cpg文件：存储的是属性表编码，例如ANSI或UTF-8，如果SHP打开属性表乱码很可能就是CPG文件出问题了

3、GDAL

GDAL（Geospatial Data Abstraction Library）是一个开源的地理空间数据处理库，它提供了强大的数据转换功能。通过编写一段GDAL命令或脚本，可以将栅格TIFF文件转换为矢量格式。

二、基本思路

1、数据准备

输入数据： 栅格数据可以是遥感图像，DEM地形，也可以是经过栅格计算后提取的坡度坡向等数据。格式为tif。栅格数据的像元值通常为整数，每个整数代表栅格的一种颜色或类别。

输出数据： 需要定义输出的矢量数据的格式，一般而言，栅格转面要素的情况最常见。在转换过程中，可以选择一个字段，该字段的值将被用来指定输出面要素的属性。默认情况下，这个字段是“Value”，它包含了每个栅格像元中的值。

栅格转矢量的输出格式可以是GeoJSON，也可以是shapefile。

这里使用DEM数据作为输入数据。
在这里插入图片描述

2、重投影（可选）

数据重投影是一个很常见的处理方法，将图像或数据集从一个地理坐标系统转换到另一个坐标系统。这项技术对于确保数据的一致性和可用性至关重要，尤其是在处理来自不同来源和不同坐标系统的数据时。

之前博客【GIS教程】使用Arcpy+GDAL批量投影并批量生成COG文件-附完整代码中介绍了使用Arcpy进行批量投影的方法，本篇补充使用GDAL函数修改栅格数据的坐标系。

使用函数gdal.Warp：

# 打开源TIF文件
input_tif = r'C:\yourpath\Desktop\input.tif' #源栅格数据路径
prj_tif = r'C:\yourpath\Desktop\prj_output.tif' #投影后的栅格数据路径
raster_dataset = gdal.Open(input_tif, gdal.GA_ReadOnly) # 读取源数据

# 设置目标坐标系为WGS84（EPSG:4326）
target_srs = osr.SpatialReference()
target_srs.ImportFromEPSG(4326)

# 使用gdal.Warp进行坐标系转换，并生成重投影的TIF文件
gdal.Warp(prj_tif, raster_dataset, dstSRS=target_srs)

# 关闭数据集
raster_dataset = None

得到的prj_tif即为重投影后的栅格数据。

3、创建空的矢量图层

在GDAL中创建空的矢量图层是为了初始化一个具有特定几何类型和属性结构的空间数据容器，以便后续能够将几何要素和属性数据写入其中，从而构建完整的矢量数据集。

使用函数data_source.CreateLayer，需要定义矢量格式、坐标系、图层以及属性字段。

#定义坐标系
prj = osr.SpatialReference()
prj.ImportFromEPSG(4326)

# 创建一个空的矢量图层:GeoJOSN
output_file = r'C:\yourpath\Desktop\result.geojson'
driver = ogr.GetDriverByName('GeoJSON') 

data_source = driver.CreateDataSource(output_file) #这里添加输出文件路径
# polygons为geojson的"name"参数，geom_type为类型，prj为定义的坐标系。
layer = data_source.CreateLayer('polygons', geom_type=ogr.wkbPolygon, srs = prj)
# 定义矢量属性，不可省略
layer.CreateField(ogr.FieldDefn('value', ogr.OFTReal))


# 创建一个空的矢量图层:Shapefile
output_file = r'C:\yourpath\Desktop\result' #路径为文件夹
driver = ogr.GetDriverByName('ESRI Shapefile') 

data_source = driver.CreateDataSource(output_file) #这里添加输出文件路径
# polygons为geojson的"name"参数，geom_type为类型，prj为定义的坐标系。
layer = data_source.CreateLayer('polygons', geom_type=ogr.wkbPolygon, srs = prj)
# 定义矢量属性，不可省略
layer.CreateField(ogr.FieldDefn('value', ogr.OFTReal))

4、栅格转矢量

使用函数gdal.Polygonize：将每个像元转成一个矩形，然后将相似的像元进行合并。

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

#读取栅格波段
tmp_dataset = gdal.Open(input_tif, gdal.GA_ReadOnly) # 读取栅格数据
tmp_band = tmp_dataset.GetRasterBand(1) # 最后想要转为矢量的波段，单波段是1
mask_band = tmp_band.GetMaskBand() # 定义掩膜范围 

# 将栅格数据转换为矢量多边形
# outLayer 是定义好的矢量的图层Layer
gdal.Polygonize(srcBand=tmp_band, maskBand=mask_band, outLayer=layer, iPixValField=0)

输出结果：GeoJSON
在这里插入图片描述
输出结果：Shapefile

三、完整代码

# 栅格转geojson实验代码
import os
from osgeo import gdal, ogr, osr

def reproject_raster(input_tif, prj_tif, target_epsg):
    # 打开源TIF文件
    raster_dataset = gdal.Open(input_tif, gdal.GA_ReadOnly)

    # 设置目标坐标系
    target_srs = osr.SpatialReference()
    target_srs.ImportFromEPSG(target_epsg)

    # 使用gdal.Warp进行坐标系转换，并生成临时重投影的TIF文件
    gdal.Warp(prj_tif, raster_dataset, dstSRS=target_srs)

    # 关闭数据集
    raster_dataset = None

def create_empty_vector(output_file, target_epsg):
    prj = osr.SpatialReference()
    prj.ImportFromEPSG(target_epsg)
    # 创建一个空的矢量图层
    driver = ogr.GetDriverByName('GeoJSON') #修改 GeoJSON 或 ESRI Shapefile
    data_source = driver.CreateDataSource(output_file)
    layer = data_source.CreateLayer('polygons', geom_type=ogr.wkbPolygon,srs = prj)
    layer.CreateField(ogr.FieldDefn('value', ogr.OFTReal))
    return data_source, layer

def raster_to_vector(input_tif, output_file):
    # 对input_tif进行重投影
    tmp_tif = "temp.tif" # 中间的重投影文件，视为过程文件
    target_epsg = 4326
    reproject_raster(input_tif, tmp_tif, target_epsg)

    # 打开重投影后的TIF文件
    tmp_dataset = gdal.Open(tmp_tif, gdal.GA_ReadOnly)
    tmp_band = tmp_dataset.GetRasterBand(1)
    mask_band = tmp_band.GetMaskBand()

    # 创建空矢量图层
    data_source, layer = create_empty_vector(output_file, target_epsg)

    # 将栅格数据转换为矢量多边形
    gdal.Polygonize(srcBand=tmp_band, maskBand=mask_band, outLayer=layer, iPixValField=0)

    # 关闭数据源和数据集
    data_source.Destroy()
    tmp_dataset = None
    os.remove(tmp_tif)  # 移除中间文件


# 使用函数
input_tif = r'C:\yourpath\Desktop\input.tif'
output_file = r"C:\yourpath\Desktop\output.geojson"
raster_to_vector(input_tif, output_file)

四、总结

GDAL栅格转矢量是将栅格数据（像素表示的连续数据）转换为矢量数据（离散的几何对象）的过程，通常涉及重投影gdal.Warp、创建空矢量图层data_source.CreateLayer、使用gdal.Polygonize函数将栅格像素转换为矢量多边形，并最终编辑和保存转换结果，以便进行空间分析和GIS应用。