Python+OGR库学习（八）：关于面矢量文件的一些小操作

对于OGR库对shp文件的基本读写没有问题，整理一些对geometry简单的小操作

代码功能

1、面矢量转线：主要理解面由ring构成，提取ring直接写入线矢量文件就OK（或者构建多线几何体，将所有ring都添加进去，完成后可以在arcmap中用拆分多部件要素拆分）
2、读取面矢量文件，输出四角多边形：获取图层四角范围，直接写入ring
3、输入面矢量，输出它的凸多边形：关键在于创建多面几何体，把每一个feature都写入，geomcolle.ConvexHull()构造凸多边形，写入图层
4、读取面矢量文件中，每个面的质心并保存为点矢量文件

代码及结果

1、面转线

##面转线方法2：一个面对应一个线要素
import ogr,osr,os
def pol2line(polyfn,linefn):
    driver = ogr.GetDriverByName('ESRI Shapefile')
    polyds = ogr.Open(polyfn,0)
    polyLayer = polyds.GetLayer()
    #创建输出文件
    if os.path.exists(linefn):
        driver.DeleteDataSource(linefn)
    lineds =driver.CreateDataSource(linefn)
    linelayer = lineds.CreateLayer(linefn,geom_type = ogr.wkbLineString)
    featuredefn = linelayer.GetLayerDefn()
    #获取ring到几何体
    #geomline = ogr.Geometry(ogr.wkbGeometryCollection)
    for feat in polyLayer:
        geom = feat.GetGeometryRef()
        ring = geom.GetGeometryRef(0)
        #geomcoll.AddGeometry(ring)
        outfeature = ogr.Feature(featuredefn)
        outfeature.SetGeometry(ring)
        linelayer.CreateFeature(outfeature)
        outfeature = None

if __name__ == '__main__':
    os.chdir(r'F:\Python+gdal\CookBook\data')
    pol2line('cache_towns.shp','poly2line_2.shp')

2、面转单多边形

##获取已有shp的四角范围，生成polygon写入新的shp
try:
    from osgeo import ogr
except ImportError:
    import ogr
import os,glob
os.chdir(r'F:\Python+gdal\CookBook\data')
# shps = glob.glob('*.shp')
# print(shps)#['ut_counties.shp']
driver = ogr.GetDriverByName('ESRI Shapefile')
inds = ogr.Open('cache_towns.shp',0)
inlayer = inds.GetLayer()

#空间参考
insrs = inlayer.GetSpatialRef()
#print(insrs)
extent = inlayer.GetExtent(True)
print(extent)#,(-114.047272999, -109.043206409, 36.991746304, 42.0023005849)
#创建环
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(extent[0],extent[2])
ring.AddPoint(extent[0],extent[3])
ring.AddPoint(extent[1],extent[3])
ring.AddPoint(extent[1],extent[2])
ring.CloseRings()
#创建多边形
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
#创建输出文件,如果已存在，删除
outfile = 'out_extent.shp'
if os.path.exists(outfile):
    driver.DeleteDataSource(outfile)
outds = driver.CreateDataSource(outfile)

outlayer = outds.CreateLayer('out_extent',srs = insrs,geom_type = ogr.wkbPolygon)
#输出图层添加字段属性
outfield = ogr.FieldDefn('id',ogr.OFTInteger)
outlayer.CreateField(outfield)
#获取属性表信息
outfielddefn = outlayer.GetLayerDefn()
#创建要素，添加几何体，写入字段，要素写入图层
feat = ogr.Feature(outfielddefn)
feat.SetGeometry(poly)
feat.SetField('id',1)

outlayer.CreateFeature(feat)
feat.Destroy()
#
inds = None
outds = None

3、面转凸多边形

import ogr,os
os.chdir(r'F:\Python+gdal\CookBook\data')

driver = ogr.GetDriverByName('ESRI Shapefile')
inds = ogr.Open('cache_towns.shp',0)
inlayer = inds.GetLayer()
#获取空间参考
insrs = inlayer.GetSpatialRef()
#创建几何体集合,写入多边形
geomcoll = ogr.Geometry(ogr.wkbGeometryCollection)
for feat in inlayer:
    geomcoll.AddGeometry(feat.GetGeometryRef())
#获取凸多边形
geomconvex = geomcoll.ConvexHull()
#创建输出文件
outfile = 'convexHull.shp'
if os.path.exists(outfile):
    driver.DeleteDataSource(outfile)
outds = driver.CreateDataSource(outfile)
outlayer = outds.CreateLayer('convexHull',srs = insrs,geom_type = ogr.wkbPolygon)
#创建表属性
idfield = ogr.FieldDefn('id',ogr.OFTInteger)
outlayer.CreateField(idfield)
#获取属性表信息
outfeatureddefn = outlayer.GetLayerDefn()
#创建要素并写入
outfeat = ogr.Feature(outfeatureddefn)
outfeat.SetGeometry(geomconvex)
outfeat.SetField('id',1)
outlayer.CreateFeature(outfeat)
feat = None

inds = None
outds = None

4、面质心转点矢量

import ogr,os
os.chdir(r'F:\Python+gdal\CookBook\data')

infile = 'cache_towns.shp'
driver = ogr.GetDriverByName('ESRI Shapefile')
inds = ogr.Open(infile,0)
inlayer = inds.GetLayer()

insrs =inlayer.GetSpatialRef()
#print(insrs)

#创建输出文件
outfile = 'Centroids_countries.shp'
if os.path.exists(outfile):
    driver.DeleteDataSource(outfile)
outds = driver.CreateDataSource(outfile)

outlayer = outds.CreateLayer('Centroids_countries',srs = insrs,geom_type = ogr.wkbPoint)
#定义属性表
infeaturedefn = inlayer.GetLayerDefn()
fieldcount = infeaturedefn.GetFieldCount()
for i in range(fieldcount):
    outfielddefn = infeaturedefn.GetFieldDefn(i)
    outlayer.CreateField(outfielddefn)
#获取输出图层属性定义
outfeaturedefn = outlayer.GetLayerDefn()
#遍历要素，获取中心点，创建输出要素，写入几何点，写入字段值，写入图层
for i in range((inlayer.GetFeatureCount())):
    feat = inlayer.GetFeature(i)
    outfeat = ogr.Feature(outfeaturedefn)
    #写入字段值
    for j in range(outfeaturedefn.GetFieldCount()):
        outfeat.SetField(outfeaturedefn.GetFieldDefn(j).GetNameRef(), feat.GetField(j))
    #获取几何，提取中心点
    geom = feat.GetGeometryRef()
    centroids = geom.Centroid()
    outfeat.SetGeometry(centroids)
    #添加到图层
    outlayer.CreateFeature(outfeat)
    # feat = None
    # outfeat = None
inds = None
outds = None

代码注重几何要素的变化，对字段值没有过多说明，后期对Geometry类多加学习

转载自：https://blog.csdn.net/weixin_42924891/article/details/85868278