本节将介绍如何利用python完成对shp的基本操作

1.读取shp四至

import shapefile
sf = shapefile.Reader(r"E:\shp\1.shp")
#读取shp四至
min_x, min_y, max_x, max_y = sf.bbox

#读取每个图斑四至
shapes = sf.shapes()
arr = []
for i in range(0, len(shapes)):
arr.append(shapes[i].bbox)

利用GDAL ogr读取shp图版四至，并添加到属性表中。

import os
from osgeo import ogr

shp中添加四至

def shapes_boundary(shp_path):
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shp_path, 1)
layer = dataSource.GetLayer()
new_field1 = ogr.FieldDefn("minX", ogr.OFTReal)
layer.CreateField(new_field1)
new_field2= ogr.FieldDefn("minY", ogr.OFTReal)
layer.CreateField(new_field2)
new_field2 = ogr.FieldDefn("maxX", ogr.OFTReal)
layer.CreateField(new_field2)
new_field2 = ogr.FieldDefn("maxY", ogr.OFTReal)
layer.CreateField(new_field2)
t = int(layer.GetFeatureCount())
for i in range(0, t):
feat = layer.GetFeature(i)
geom = feat.GetGeometryRef()
minX = geom.GetEnvelope()[0]
minY = geom.GetEnvelope()[2]
maxX = geom.GetEnvelope()[1]
maxY = geom.GetEnvelope()[3]
feat.SetField("minX", minX)
feat.SetField("minY", minY)
feat.SetField("maxX", maxX)
feat.SetField("maxY", maxY)
layer.SetFeature(feat)
if __name__ == '__main__':
shp=r"H:\test\1.shp"
shapes_boundary(shp)

python利用ogr写入shp文件，定义shp文件属性字段（field）的数据格式为：

OFTInteger　　　　　　　# 整型
OFTIntegerList 　　　　# 整型list
OFTReal 　　　　　　　　# 双精度
OFTRealList 　　　　　 # 双精度list
OFTString 　　　　　　 # 字符
OFTStringList 　　　　 # 字符list
OFTWideString 　　　　 # 长字符
OFTWideStringList 　　# 长字符list
OFTBinary 　　　　　　　
OFTDate
OFTTime
OFTDateTime
OFTInteger64
OFTInteger64List
OFSTNone
OFSTBoolean
OFSTInt16
OFSTFloat32
OJUndefined

2.判断某个点是否在shp中

import shapefile
import shapely.geometry as geometry
import numpy as np
lons, lats = [], []

lons = np.linspace(123.67, 123.32, 50)

lats = np.linspace(23.48, 25.73, 25)

grid_lon, grid_lat = np.meshgrid(lons, lats)
flat_lon = grid_lon.flatten()
flat_lat = grid_lat.flatten()
flat_points = np.column_stack((flat_lon, flat_lat))
in_shape_points = []
sf = shapefile.Reader("E:/test/1.shp", encoding='gbk')
shapes = sf.shapes()
for pt in flat_points:
if geometry.Point(pt).within(geometry.shape(shapes[0])):
in_shape_points.append(pt)
print("The point is in shp")
else:
print("The point is not in shp")
print(in_shape_points)

3.gdal生成shp

import osgeo.ogr as ogr
import osgeo.osr as osr

def run():
driver = ogr.GetDriverByName("ESRI Shapefile")
data_source = driver.CreateDataSource("Gooise.shp")
srs = osr.SpatialReference()
srs.ImportFromEPSG(32631)
layer = data_source.CreateLayer("Gooise", srs, ogr.wkbPolygon)
feature = ogr.Feature(layer.GetLayerDefn())
#创建wkt文本
wkt = 'polygon((646080 5797460,648640 5797460,648640 5794900,646080 5794900,646080 5797460))'
polygon = ogr.CreateGeometryFromWkt(wkt)
feature.SetGeometry(polygon)
layer.CreateFeature(feature)
feature = None
data_source = None

4.shp拆分成多个shp

import osgeo.ogr as ogr
import osgeo.osr as osr

def create_shp(shp_name,wkt):
driver = ogr.GetDriverByName("ESRI Shapefile")
data_source = driver.CreateDataSource(f'E:/cq/tif/shp/{shp_name}.shp')
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
layer = data_source.CreateLayer(f'{shp_name}', srs, ogr.wkbMultiPolygon)
feature = ogr.Feature(layer.GetLayerDefn())
polygon = ogr.CreateGeometryFromWkt(wkt)
feature.SetGeometry(polygon)
layer.CreateFeature(feature)
feature = None
data_source = None

def run():
driver = ogr.GetDriverByName('ESRI Shapefile')
fileName = "E:/cq/中小河流.shp"
dataSource = driver.Open(fileName, 0)
layer = dataSource.GetLayer(0)
print("空间参考 :{0}".format(layer.GetSpatialRef()))
for i in range(0, layer.GetFeatureCount()):
feat = layer.GetFeature(i)
wkt = feat.geometry()
print(wkt)
create_shp(i, str(wkt))

if __name__ == '__main__':
run()

5.基于gdal的面矢量面积计算

import ogr

def area(shpPath):
'''计算面积'''
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shpPath, 1)
layer = dataSource.GetLayer()
new_field = ogr.FieldDefn("Area", ogr.OFTReal)
new_field.SetWidth(32)
new_field.SetPrecision(16) # 设置面积精度,小数点后16位
layer.CreateField(new_field)
for feature in layer:
geom = feature.GetGeometryRef()
area = geom.GetArea() # 计算面积
# m_area = (area/(0.0089**2))*1e+6 # 单位由十进制度转为米
# print(m_area)
feature.SetField("Area", area) # 将面积添加到属性表中
layer.SetFeature(feature)
dataSource = None

6.使用面积和Value值过滤矢量图层

import sys
from osgeo import ogr, osr, gdal

过滤矢量图层

def guolv(shp,mask,out,name):
'''
:param shp:shp路径
:param mask:栅格路径，主要用来取投影信息
:param out:输出路径
:param name:文件名
:return:None
'''
driver = ogr.GetDriverByName('ESRI Shapefile')
dataSource = driver.Open(shp, 0) # 0是只读，1可写
if dataSource is None:
print('could not open')
sys.exit(1)
# 获取图层
layer = dataSource.GetLayer(0)
t = int(layer.GetFeatureCount())
drv = ogr.GetDriverByName('ESRI Shapefile')
Polygon = drv.CreateDataSource(out)
data = gdal.Open(mask, gdal.GA_ReadOnly)
prj = osr.SpatialReference()
prj.ImportFromWkt(data.GetProjection())
# oLayer = Polygon.CreateLayer(name, srs=prj, geom_type=ogr.wkbMultiPolygon)
oLayer = Polygon.CreateLayer(name, srs=prj, geom_type=ogr.wkbPolygon)
oDefn = oLayer.GetLayerDefn() # 定义要素
gardens = ogr.Geometry(ogr.wkbPolygon)
oFieldID = ogr.FieldDefn("ID", ogr.OFTInteger) # 创建一个叫ID的整型属性
oLayer.CreateField(oFieldID, 1)
new_field1 = ogr.FieldDefn("minX", ogr.OFTReal)
oLayer.CreateField(new_field1)
new_field2= ogr.FieldDefn("minY", ogr.OFTReal)
oLayer.CreateField(new_field2)
new_field3 = ogr.FieldDefn("maxX", ogr.OFTReal)
oLayer.CreateField(new_field3)
new_field4 = ogr.FieldDefn("maxY", ogr.OFTReal)
oLayer.CreateField(new_field4)
feature = ogr.Feature(oLayer.GetLayerDefn())
ID=0
for i in range(0, t):
feat = layer.GetFeature(i)
wkt = (feat.geometry())
geom = feat.GetGeometryRef()
area = geom.GetArea()
m_area = (area / (0.0089 ** 2)) * 1e+6
v = feat.GetField('Value')
# 面积过滤
if m_area > 10000 and v==1:
ID = ID+1
polygon = ogr.CreateGeometryFromWkt(str(wkt)) ## 生成面
minX = geom.GetEnvelope()[0]
minY = geom.GetEnvelope()[2]
maxX = geom.GetEnvelope()[1]
maxY = geom.GetEnvelope()[3]

        feature.SetField("minX", float(minX))  
        feature.SetField("minY", float(minY))  
        feature.SetField("maxX", float(maxX))  
        feature.SetField("maxY", float(maxY))

        feature.SetGeometry(polygon)  ## 设置面  
        feature.SetField(0, ID)  
        oLayer.CreateFeature(feature)