使用ogr裁剪向量資料

使用ogr裁剪向量資料

由來：​ 近期有個需求，內容是這樣的：我們有兩個向量資料，現在要求以乙個向量檔案為底板，按字段對另乙個向量檔案進行分割，生成若干小的shpfile檔案

分析：​ 經過分析之後，將步驟拆解如下：

首先確保兩個shpfile投影座標系統一

​ 如果出現不統一的情況，那麼用arcgis的工具project進行投影轉換。（data management tools--projections and transformations--project），鏈結如下

其次編寫按屬性字段分割shpfile，生成若干小shpfile的**

接著編寫根據layer建立shpfile的**（>使用ogr中拷貝方法建立新的shpfile

最後進行迴圈，對每次生成的shpfile和原先的shpfile進行空間查詢，得到layer後生成需要的shpfile

核心**：

from osgeo import ogr, osr

import os
def createshpbylayer(shp, layer, filetype):
'''根據layer建立shpfile
'''driver = ogr.getdriverbyname("esri shapefile")
ds = driver.createdatasource(shp)
pt_layer = ds.copylayer(layer, 'layername')
ds.destroy()
def splitshp(shpfile, outpath, splitfield):
'''按屬性字段分割shpfile
'''ds = ogr.open(shpfile)
lyr = ds.getlayer(0)
for feat in lyr:
cityname = feat.getfield(splitfield)  # 以欄位名為檔名稱
outshp = os.path.join(outpath, str(cityname)+'.shp')
geom = feat.getgeometryref()
driver = ogr.getdriverbyname("esri shapefile")
outds = driver.createdatasource(outshp)
outlyr = outds.createlayer("layername", lyr.getspatialref(), ogr.wkbmultipolygon)
outlyr.createfields(lyr.schema)  # 建立字段屬性
outfeat = ogr.feature(lyr.getlayerdefn())
for i in range(feat.getfieldcount()):
val = feat.getfield(i)
outfeat.setfield(i, val)
outfeat.setgeometry(geom)
outlyr.createfeature(outfeat)
outds = none
def splitshpbyshp(covershp, splitshp, outpath):
'''@desc: 將乙個shp按屬性分割後，再用來分割另乙個shp（空間查詢）
'''ds = ogr.open(covershp)
lyr = ds.getlayer(0)
splitds = ogr.open(splitshp)
splitlyr = splitds.getlayer(0)
for feat in lyr:
cityname = feat.getfield("name")  # 以欄位名為檔名稱
outshp = os.path.join(outpath, str(cityname)+'.shp')
geom = feat.getgeometryref()
driver = ogr.getdriverbyname("memory")
outds = driver.createdatasource("temp")
outlyr = outds.createlayer("layername", lyr.getspatialref(), ogr.wkbmultipolygon)
outlyr.createfields(lyr.schema)  # 建立字段屬性
outfeat = ogr.feature(lyr.getlayerdefn())
for i in range(feat.getfieldcount()):
val = feat.getfield(i)
# val = val.encode('utf8') # 將unicode編碼轉化為中文，還有點問題
outfeat.setfield(i, val)
outfeat.setgeometry(geom)
outlyr.createfeature(outfeat)
outds = none
splitlyr.setspatialfilter(geom) # 按空間查詢
createshpbylayer(outshp, splitlyr, "esri shapefile")
print("按空間查詢的要素數量："+ str(splitlyr.getfeaturecount()))

轉出來的shpfile屬性欄位中，中文部分亂碼，目前未能解決。

向量裁剪向量

也不知道為啥，向量裁剪向量這麼普通的東西這麼難找，趕緊放出來讓大家用用 import os import numpy as np import geopandas as gpd import warnings warnings.filterwarnings ignore geoseries.notn...

使用向量面裁剪柵格資料的對齊問題

最近湊巧有幾個比較多的柵格裁剪問題，整理如下 我們只有由於柵格與向量資料的儲存模型不相同，這就導致柵格資料的像元無法與向量資料的點等同，從而導致裁切後的對齊問題，放大資料我們就能發現，如下圖可以說明 其中黑白色為柵格資料，每個正方形代表乙個像元，紅色區域為向量面資料。我們按照預設設定執行 raste...

狸貓的面試 專案描述 向量裁剪

專案描述 要求 在4s內，用單執行緒，完成對百萬級的向量圖形的裁剪 輸入 1.乙個多邊形 儲存在xml檔案中，順序儲存多邊形所有的端點 2.乙個向量圖形 line或者circle line儲存其兩個端點ab，circle儲存其圓心o和半徑r。實際有一百萬條直線和一百萬個圓 輸出 輸出裁剪後的圖形。例...