python使用gdal将栅格文件转为shp_栅格转shp

python使用gdal将栅格文件转为shp

将栅格数据转为shp面

from osgeo import gdalconst, gdal, ogr, osr
import os


def raster2poly(raster, outshp):
    inraster = gdal.Open(raster)  # 读取路径中的栅格数据
    inband = inraster.GetRasterBand(1)  # 这个波段就是最后想要转为矢量的波段，如果是单波段数据的话那就都是1
    prj = osr.SpatialReference()
    prj.ImportFromWkt(inraster.GetProjection())  # 读取栅格数据的投影信息，用来为后面生成的矢量做准备

    drv = ogr.GetDriverByName("ESRI Shapefile")
    if os.path.exists(outshp):  # 若文件已经存在，则删除它继续重新做一遍
        drv.DeleteDataSource(outshp)
    Polygon = drv.CreateDataSource(outshp)  # 创建一个目标文件
    Poly_layer = Polygon.CreateLayer(raster[:-4], srs=prj, geom_type=ogr.wkbMultiPolygon)  # 对shp文件创建一个图层，定义为多个面类
    newField = ogr.FieldDefn('value', ogr.OFTReal)  # 给目标shp文件添加一个字段，用来存储原始栅格的pixel value
    Poly_layer.CreateField(newField)

    gdal.FPolygonize(inband, None, Poly_layer, 0)  # 核心函数，执行的就是栅格转矢量操作
    Polygon.SyncToDisk()
    Polygon = None


raster2poly("input_tif.tif", "output_shp.shp")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

将栅格数据转为shp点

# 002  转为shp点
from osgeo import ogr, osr, gdal
import numpy as np
import pandas as pd

inputTIF = "input_tif.tif"  # 待输入的TIF路径
NoData_value = -999  # 掩膜遮挡部分的值 NoData
result_path = 'output_shp.shp'  # shp结果文件的存放路径


def createPoint_strNosrs(reader, result_path):
    """
    创建无空间参考且字段类型全部为string的点shp文件
    :param reader: 待生成shp文件的表属性信息 为dict格式 其中最后两个key为点的位置信息  key='point'
    :param result_path: 待生成的shp 文件的存储地址
    :return:
    """
    # 设置shapefile驱动
    driver = ogr.GetDriverByName("ESRI Shapefile")

    # 创建数据源
    data_source = driver.CreateDataSource(result_path)

    # 创建图层
    layer = data_source.CreateLayer("point", geom_type=ogr.wkbPoint)  # 无坐标系文件

    # 初始化字段
    l_list = []
    l_col = list(reader[0].keys())
    for i in range(len(l_col) - 1):
        l_list.append({list(reader[0].keys())[i]: ogr.FieldDefn(str(l_col[i][:10]), ogr.OFTString)})
        ogr.FieldDefn(str(l_col[i][:10]), ogr.OFTString).SetWidth(120)  # 设置长度
        layer.CreateField(ogr.FieldDefn(str(l_col[i][:10]), ogr.OFTString))  # 创建字段

    for row in reader:
        # 创建特征
        feature = ogr.Feature(layer.GetLayerDefn())
        # 和设置字段内容进行关联  ,从数据源中写入数据

        # 设置创建字段
        for use_i in l_list:
            feature.SetField(str(list(use_i.keys())[0][:10]), row[list(use_i.keys())[0]])

        # 创建WKT 文本点
        wkt = row['point']

        # 生成实体点
        point = ogr.CreateGeometryFromWkt(wkt)

        # 使用点
        feature.SetGeometry(point)
        # 添加点
        layer.CreateFeature(feature)
        # 关闭 特征
        feature = None
    # 关闭数据
    data_source = None

# 读取栅格文件
data = gdal.Open(inputTIF, gdal.GA_ReadOnly)
data1= data.GetRasterBand(1).ReadAsArray()

#存储着栅格数据集的地理坐标信息
transform = data.GetGeoTransform()

# NoData_value = data1[0][0]  # 找到tif中代表NoData的值 以免报错
res = np.where(data1 != NoData_value )  # 获取栅格数据中 除去掩膜部分的点的坐标索引

lon = transform[3] + res[1] * transform[4] + res[0] * transform[5]
lat = transform[0] + res[1] * transform[1] + res[0] * transform[2]
values = data1[data1 != NoData_value ]  # 获取栅格数据中  除去掩膜部分的点的值

point = ['POINT (' + str(lat[k]) + ' ' + str(lon[k])+')' for k in range(len(lon))]


dataUse = pd.DataFrame({"value": values, "point": point})

reader = dataUse.to_dict(orient='records')
createPoint_strNosrs(reader, result_path)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79