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

You may also like...

退出移动版