本節將介紹如何利用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)