Python——栅格转矢量

Python栅格转矢量

在利用python进行栅格转矢量操作时，我们可以选择gdal、shapely等多个库包进行

现提出两种栅格转矢量方法，其一利用shapely+geopandas+rasterio，其二利用gdal进行

总体来说，前者更为简单，且矢量相同属性合并时，建议使用geopandas，因为gdal真的很麻烦（就目前所使用的的方法来说。。。

如果小伙伴有针对gdal相同属性矢量合并功能较为简单的代码的话，欢迎留言，一起学习啦~~

利用shapely+rasterio+geopandas进行栅格转矢量

#!usr/bin/env python
# -*- coding: utf-8 -*-
"""
date：2022/11/17
author：甲戌_Tr
email: liu_xxxi@163.com
"""

import rasterio
from rasterio import features
from shapely.geometry import shape
import geopandas as gpd

def raster2polygon(inraster,outshp,nodata=0):
    out_shp = gpd.GeoDataFrame(columns=['DN','geometry'])

    with rasterio.open(inraster) as f:
        image = f.read(1)
        img_crs = f.crs
        image[image == f.nodata] = nodata

        i = 0
        for coords, value in features.shapes(image, transform=f.transform):
            if value != nodata:
                geom = shape(coords)
                out_shp.loc[i] = [value,geom]
                i += 1

    out_shp.set_geometry('geometry',inplace=True)
    out_shp = out_shp.dissolve(by='DN')

    out_shp.to_file(outshp,crs=img_crs)

    return out_shp

利用gdal进行栅格转矢量

#!usr/bin/env python
# -*- coding: utf-8 -*-
"""
date：2022/11/17
author：甲戌_Tr
email: liu_xxxi@163.com
"""

from osgeo import gdal,ogr,osr

def gdal_raster2polygon(inraster,outshp):
    raster_ds = gdal.Open(inraster)
    band = raster_ds.GetRasterBand(1)

    raster_proj = osr.SpatialReference(wkt=raster_ds.GetProjection())
    epsg_num = raster_proj.GetAttrValue('AUTHORITY',1)
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(int(epsg_num))

    # srs.ImportFromWkt(raster_ds.GetProjection())

    driver = ogr.GetDriverByName("memory")

    out_ds = driver.CreateDataSource('out')

    out_layer = out_ds.CreateLayer('out',srs=srs,geom_type=ogr.wkbPolygon)

    field = ogr.FieldDefn('DN',ogr.OFTInteger)
    out_layer.CreateField(field)

    gdal.Polygonize(band,band,out_layer,0,[],callback=None)

    # 对相同DN值的矢量进行合并
    multipoly = True

    driver = ogr.GetDriverByName('ESRI Shapefile')

    union_ds = driver.CreateDataSource(outshp)
    union_layer = union_ds.CreateLayer(outshp,srs=srs,geom_type=ogr.wkbPolygon)


    ###分级字段创建
    fd = ogr.FieldDefn('DN', ogr.OFTInteger)
    union_layer.CreateField(fd)  # 创建字段

    unique_DN = []
    for f in out_layer:
        f_value = f.GetField("DN")
        if f_value not in unique_DN:

            unique_DN.append(int(f_value))

    out_layer.ResetReading()

    defn = out_layer.GetLayerDefn()

    for i in unique_DN:
        field_sql = f"DN={str(i)}"

        out_layer.SetAttributeFilter(field_sql)
        multi = ogr.Geometry(ogr.wkbMultiPolygon)

        for feat in out_layer:
            if feat.geometry():
                wkt = feat.geometry().ExportToWkt()
                multi.AddGeometryDirectly(ogr.CreateGeometryFromWkt(wkt))

        union = multi.UnionCascaded()
        geo_name = union.GetGeometryName()

        if geo_name != 'MULTIPOLYGON':
            multipoly = False

        if multipoly is False:
            for geom in union:
                poly = ogr.CreateGeometryFromWkb(geom.ExportToWkb())
                feat = ogr.Feature(defn)
                feat.SetGeometry(poly)
                union_layer.CreateFeature(feat)
                union_layer.ResetReading()
        else:
            out_feat = ogr.Feature(defn)
            out_feat.SetGeometry(union)
            out_feat.SetField('DN', i)
            union_layer.CreateFeature(out_feat)
            union_layer.ResetReading()

    union_layer.SetAttributeFilter(None)
    for ft in union_layer:
        if ft.GetField ('DN') ==None:
            union_layer.DeleteFeature(ft.GetFID())

    union_layer = None
    union_ds = None

欢迎指正~~