[Python] GDAL/OGR操作矢量數(shù)據(jù)（shp隘谣、GeoJSON）
GeoDoer 2019-12-30 12:17:51? 2196? 收藏 15
分類專欄： # GIS|Python
版權(quán)
GDAL項(xiàng)目旨于地理數(shù)據(jù)抽象模型對(duì)地理數(shù)據(jù)文件進(jìn)行讀寫管理蝉稳；而其項(xiàng)目下有兩大類模塊：GDAL和OGR
OGR提供操作矢量數(shù)據(jù)的API淋纲，GDAL模塊提供柵格數(shù)據(jù)的API
【相關(guān)鏈接】
1.GDAL/OGR例子合集
2.GDAL官網(wǎng)API
3.網(wǎng)上的學(xué)習(xí)筆記抓艳，教程來自于猶他州立大學(xué)
4.中國OSGeo教程
文章目錄
OGR中的數(shù)據(jù)結(jié)構(gòu)
driver解析數(shù)據(jù)的驅(qū)動(dòng)
datasource管理數(shù)據(jù)源
layer管理圖層
feature管理要素類
從Layer中獲得feature
更新要素
刪除要素
關(guān)閉要素
創(chuàng)建要素
geometry管理要素的幾何屬性
Geometry的幾何屬性
從Feature獲取Geometry
常用方法
創(chuàng)建Geometry
field管理要素的屬性字段
創(chuàng)建字段
更改屬性定義
其他
經(jīng)驗(yàn)
wkt
查詢
屬性查詢
空間查詢
SQL查詢
其他
設(shè)置編碼
示例
打開shp
刪除數(shù)據(jù)
遍歷所有屬性值
創(chuàng)建
創(chuàng)建點(diǎn)狀數(shù)據(jù)集
創(chuàng)建線狀數(shù)據(jù)集
創(chuàng)建多邊形數(shù)據(jù)集
創(chuàng)建屬性
創(chuàng)建點(diǎn)Geometry
創(chuàng)建線Geometry
創(chuàng)建多邊形Geometry
拷貝
datasource層次的拷貝
layer層次的拷貝
feature層次的拷貝
OGR中的數(shù)據(jù)結(jié)構(gòu)
driver解析數(shù)據(jù)的驅(qū)動(dòng)
# 【方法一】從已有數(shù)據(jù)源中獲取驅(qū)動(dòng)變量
ds = ogr.open(r'D:\....\...shp')
driver = ds.GetDriver()
# 【方法二】通過名稱獲取
json_driver = ogr.GetDrverByName('GeoJSON')
? ? # 獲得驅(qū)動(dòng)名的兩種方法
? ? # 1. OGR網(wǎng)站上有介紹，通過GDAL/ORG自帶的ogrinfo
? ? # 2. 代碼中提供的print_drivers函數(shù)來獲取驅(qū)動(dòng)程序的名字
1
2
3
4
5
6
7
8
datasource管理數(shù)據(jù)源
（1）創(chuàng)建數(shù)據(jù)源
有了驅(qū)動(dòng)程序之后，提供數(shù)據(jù)源名稱來使用它創(chuàng)建一個(gè)空的數(shù)據(jù)源玷或。新創(chuàng)建的數(shù)據(jù)源會(huì)自動(dòng)打開等待寫入
# 創(chuàng)建一個(gè)功能齊全的SpatialLite數(shù)據(jù)源儡首，而不是使用SQLite
fn = "filename"
driver = ogr.GetDriverByName("SQLite")
# 創(chuàng)建新的數(shù)據(jù)源時(shí)，不能覆蓋現(xiàn)有的數(shù)據(jù)源偏友。如果你的代碼可能會(huì)覆蓋現(xiàn)有數(shù)據(jù)蔬胯，那么在創(chuàng)建新數(shù)據(jù)之前需要?jiǎng)h除舊數(shù)據(jù)
if os.path.exists(fn):
driver.DeleteDataSource(fn)
ds= driver.CreateDataSource(fn,
options=[ #創(chuàng)建選項(xiàng)：傳遞一個(gè)字符串列表（參數(shù)在OGR網(wǎng)站上有文檔介紹）
'SPATIALITE=yes'
]
)
if ds is None:
#【注意】如果數(shù)據(jù)源沒有創(chuàng)建成功，那么CreateDataSource函數(shù)會(huì)返回為空位他。如果為空對(duì)象氛濒，之后使用會(huì)報(bào)AttributeError錯(cuò)誤
sys.exit('Could not create {0}'.formaat(json_fn ) )
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
（2）刪除數(shù)據(jù)源
datasource.Destory() #關(guān)閉
del datasource # 刪除datasource
""" 刪除ds變量強(qiáng)制關(guān)閉文件，并將所有的編輯成果寫入磁盤中
【注意】刪除圖層變量并不會(huì)觸發(fā)這個(gè)操作鹅髓，必須關(guān)閉數(shù)據(jù)源才行
【備注】如果想保持?jǐn)?shù)據(jù)源打開的話舞竿，可以通過圖層對(duì)象或者數(shù)據(jù)源對(duì)象調(diào)用ds.SyncToDisk()
【警告】為了使你的編輯寫入到磁盤中，必須關(guān)閉文件或者調(diào)用SyncToDisk函數(shù)窿冯。如果沒有這么做骗奖，并且在交互環(huán)境中還打開數(shù)據(jù)源，那么你會(huì)很失望地發(fā)現(xiàn)創(chuàng)建了一個(gè)空的數(shù)據(jù)集
"""
1
2
3
4
5
6
7
layer管理圖層
layer = datasource.GetLayer(0) #根據(jù)下標(biāo)來獲取圖層（對(duì)于shp而言只有一個(gè)圖層）
layer.GetGeomType() #圖層幾何類型
n = layer.GetFeatureCount() #要素?cái)?shù)量
extent = layer.GetExtent() #上下左右邊界
readedNum = GetFeaturesRead() #已經(jīng)讀取多少條Feature
layer.ResetReading() #重置內(nèi)部feature指針醒串，指向FID=0的要素
layer.GetSpatialRef() #使用wkt文本表示空間參考系統(tǒng)的實(shí)例
layer.schema #獲取FieldDefn對(duì)象列表（屬性名）
1
2
3
4
5
6
7
8
feature管理要素類
從Layer中獲得feature
# 1. 【獲得圖層中的要素】
feature = layer.GetFeature(0)
# 2. 【獲取圖層中所有要素】
feat = layer.GetNextFeature()
while feat:
feat = layer.GetNextFeature()
layer.ResetReading()
1
2
3
4
5
6
7
更新要素
要素更新操作也大同小異执桌，只不過你操作的是圖層中的現(xiàn)有要素，而不是一個(gè)空白的對(duì)象
# 【例子】為圖層的每一個(gè)要素添加一個(gè)唯一的ID
layer.CreateField( # 為圖層添加一個(gè)ID字段
ogr.FieldDefn('ID', ogr.OFTInteger)
)
n = 1
for feat in layer:
feat.SetField('ID', n) #設(shè)置字段的值
layer.SetFeature(feat) #將數(shù)值傳遞給SetFeature函數(shù)芜赌，來更新圖層中的要素信息（這里不是CreateFeature函數(shù)）
n += 1
1
2
3
4
5
6
7
8
9
刪除要素
刪除要素更容易仰挣，需要知道的就是待刪除要素的FID
#【例子】刪除字段City_Name為"Seattle"的要素
# 找到目標(biāo)的FID進(jìn)行刪除
for feat in layer:
if feat.GetField('City_Name')=='Seattle':
layer.DeleteFeature(feat.GetFID() ) #刪除操作
"""
【但是】
1. 某些特定的數(shù)據(jù)格式使用此方式并不能完全刪除要素。
2. 你可能看不出缠沈，有時(shí)要素只是被標(biāo)記為刪除了膘壶，而不是被完全拋出，它仍然還潛伏著洲愤。
3. 正因如此香椎，你不會(huì)看到其他要素被分配到剛才刪除的FID
4. 這意味著，如果已經(jīng)刪除了很多要素禽篱，在文件中可能存在大量不必要的已用空間
【所以】刪除這些要素可以回收對(duì)應(yīng)的空間
【如何回收對(duì)應(yīng)空間】如果你有一些關(guān)系數(shù)據(jù)庫的經(jīng)驗(yàn)畜伐，應(yīng)該熟悉這一點(diǎn)。它類似于在微軟的Access數(shù)據(jù)庫上運(yùn)行壓縮和修復(fù)或在PostgreSQL數(shù)據(jù)庫中使用重組功能（VACUUM）
【具體做法】打開數(shù)據(jù)源躺率，然后執(zhí)行一條SQL語句來壓縮數(shù)據(jù)庫（不同的數(shù)據(jù)格式回收的方法不同）
"""
ds.ExecuteSQL('REPACK' + layer.GetName() ) #shapefile格式回收方式
ds.ExecuteSQL('RECOMPUTE EXTENT ON' + layer.GetName() ) #確甭杲纾空間范圍更新
# 當(dāng)現(xiàn)有要素發(fā)生改變或被刪除后，它不會(huì)更新元數(shù)據(jù)中的空間范圍
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
關(guān)閉要素
feature.Destory() #關(guān)閉
1
創(chuàng)建要素
往現(xiàn)有圖層中添加新要素和往全新的圖層中添加要素的操作一樣悼吱，具體請(qǐng)看示例
【步驟】
1.基于圖層字段創(chuàng)建一個(gè)空要素慎框，填充它
2.然后把它插入到該圖層
geometry管理要素的幾何屬性
Geometry的幾何屬性
【org.Geometry中的幾何類型】ogr.wkbPoint、org.wkbLineString后添、org.wkbPolygon等笨枯，請(qǐng)看wkt
# 獲得Geometry的幾何屬性
layer.GetGeomType() == ogr.wkbPolygon
layer.GetGeometryName()
1
2
3
從Feature獲取Geometry
# 兩種方法
geometry = feature.geometry()
geometry = feat.GetGeometryRef()
1
2
3
常用方法
geom.GetGeometryName() #要素類型
geom.GetPointCount() #要素點(diǎn)個(gè)數(shù)
geom.GetX(0);geom.GetY(0);geom.GetZ(0) #獲得第一個(gè)坐標(biāo)點(diǎn)的X、Y、Z
print(geom) #打印出所有點(diǎn)
geom.ExportToWkt() #導(dǎo)出WKT
1
2
3
4
5
創(chuàng)建Geometry
請(qǐng)看示例
field管理要素的屬性字段
【字段類型常量】
有一些不同的屬性字段類型馅精，但并不是每一種字段類型都支持所有的數(shù)據(jù)格式严嗜。這時(shí)候用于描述各種數(shù)據(jù)格式的在線文檔就派上用場(chǎng)了。
字段數(shù)據(jù)類型 OGR常量
Integer OFTInteger
List of integers OFTIntegerList
Floating point number OFTReal
List of floating point numbers OFTRealList
String OFTString
List of strings OFTStringList
Date OFTDate
Time of day OFTTime
Date and time OFTDateTime
創(chuàng)建字段
在Layer中新建屬性字段
# 創(chuàng)建并添加第一個(gè)屬性字段
coord_fld = ogr.FieldDefn('X', ogr.OFTReal)
# 設(shè)置限制
coord_fld.SetWidth(8)
coord_fld.SetPrecision(3)
out_lyr.CreateField(coord_fld)
# 重用FieldDefn對(duì)象來創(chuàng)建第二個(gè)字段
coor_fld.SetName('Y')
out_lyr.CreateField(coord_fld)
1
2
3
4
5
6
7
8
9
更改屬性定義
【函數(shù)】AlterFieldDefn(iField, field_def, nFlags)
【作用】用新字段的規(guī)則替換現(xiàn)有的字段
【參數(shù)說明】
iField是想改變的屬性字段對(duì)應(yīng)的索引值
field_def是新屬性字段的定義對(duì)象
nFlags是一個(gè)整數(shù)洲敢，是下表中的一個(gè)或多個(gè)常數(shù)的總和
【用于指明字段定義的哪個(gè)屬性可以更改的標(biāo)記漫玄，如果想使用多個(gè)，一起添加即可】
需要更改的字段屬性 OGR常量
Field name only ALTER_NAME_FLAG
Field type only ALTER_TYPE_FLAG
Field width and/or precision only ALTER_WIDTH_PRECISION_FLAG
All of the above ALTER_ALL_FLAG
# 【例子一】 更改原字段的名字
index = layer.GetLayerDefn().GetFieldIndex('NAME') #獲取字段的索引值
fld_defn = ogr.FieldDefn('City_Name', ogr.OFTString) #創(chuàng)建新屬性的字段定義
layer.AlterFieldDefn(index, fld_defn, ogr.ALTER_NAME_FLAG)
# 【例子二】 更改字段的多個(gè)屬性压彭，例如名稱和精度
lyr_defn = layer.GetLayerDefn()
index = lyr_defn.GetFieldIndex('X')
width = lyr_defn.GetFieldDefn(index).GetWidth()
fld_defn = ogr.FieldDefn('X_coord', ogr.OFTReal)
fld_defn.SetWidth(width)
fld_defn.SetPrecision(4)
flag = ogr.ALTER_NAME_FLAG + ogr.ALTER_WIDTH_PRECISION_FLAG
layer.AlterFieldDefn(index, fld_defn, flag)
"""【注意】在此例子中創(chuàng)建新的字段定義時(shí)使用的是原始字段寬度
如果沒有為原始數(shù)據(jù)設(shè)置足夠大的字段寬度睦优，結(jié)果將不正確
為了解決這個(gè)問題，需要使用與原始數(shù)據(jù)一樣的字段寬度
為了使更改的精度生效壮不，所有記錄必須重寫
為數(shù)據(jù)設(shè)置一個(gè)比它自身更高的精度并不會(huì)提高精度汗盘，因?yàn)閿?shù)據(jù)不能憑空產(chǎn)生，但如果設(shè)置的精度不夠高询一，精度就會(huì)降低
"""
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
其他
feature.GetFieldCount() # 屬性總數(shù)
feature.keys() # 屬性名
feat.GetField('AREA') # 獲取字段
feature.GetFieldAsString("FID") # 獲取字段成string
1
2
3
4
經(jīng)驗(yàn)
1.沒能為GeoJSON文件設(shè)置一個(gè)精度
2.想在shapefile中設(shè)置字段精度隐孽，必須設(shè)置字段寬度
3.設(shè)置字段的默認(rèn)值：FieldDefn.SetDefault("我是默認(rèn)值")
wkt
WKB（二進(jìn)制Well-KnownBinary），用于不同軟件程序間進(jìn)行幾何要素類型轉(zhuǎn)換的一種二進(jìn)制表示標(biāo)準(zhǔn)家凯。
因?yàn)樗嵌M(jìn)制格式，所以人們無法直接閱讀獲取其表示的內(nèi)容如失，但是熟知文本格式（Well-Known Text绊诲，WKT）可以閱讀。
幾何要素類型 對(duì)應(yīng)OGR常量
Point wkbPoint
Multipoint wkbMultiPoint
Line wkbLineString
Multiline wkbMultiLineString
Polygon wkbPolygon
Multipolygon wkbMultiPolygon
Unknown geometry type wkbUnknown（圖層如果有多種幾何類型就返回wkbUnknown）
No geometry wkbNone
查詢
屬性查詢
from osgeo import ogr
import os
shpfile = r'C:\tmp\data.shp'
ds = ogr.Open(shpfile)
layer = ds.GetLayer(0) #得到圖層
lyr_count = layer.GetFeatureCount()
print(lyr_count) #原先的要素總數(shù)
layer.SetAttributeFilter("AREA > 800000") #屬性查詢條件
lyr_count = layer.GetFeatureCount() #查詢過后褪贵，layer被更新
print(lyr_count) #滿足條件的layer
# 將該layer保存成shp
driver = ogr.GetDriverByName("ESRI Shapefile")
extfile = r'C:\tmp\data.shp'
if os.access( extfile, os.F_OK ):
driver.DeleteDataSource( extfile )
newds = driver.CreateDataSource(extfile)
layernew = newds.CreateLayer('rect',None,ogr.wkbPolygon)
# 遍歷掂之，復(fù)制
feat = layer.GetNextFeature()
while feat is not None:
layernew.CreateFeature(feat)
feat = layer.GetNextFeature()
newds.Destroy()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
空間查詢
SQL查詢
from osgeo import ogr
driver = ogr.GetDriverByName("ESRI Shapefile")
world_shp = r'C:\tmp\test.shp'
world_ds = ogr.Open(world_shp)
world_layer = world_ds.GetLayer()
world_layer_name = world_layer.GetName()
result = world_ds.ExecuteSQL("select * from %s where prov_id = '22' order by BNDRY_ID desc" % (world_layer_name)) # ) # ExecuteSQL是基于數(shù)據(jù)集進(jìn)行的，而不是圖層
resultFeat = result.GetNextFeature ()
out_shp = r'C:\tmp\test\test_sql_result.shp'
create_shp_by_layer(out_shp, result) #保存結(jié)果
# 對(duì)查詢結(jié)果進(jìn)行遍歷
while resultFeat :
? ? print resultFeat.GetField('BNDRY_ID')
? ? resultFeat = result.GetNextFeature ()
#執(zhí)行下一條SQL語句之前一定要先釋放
world_ds.ReleaseResultSet(result)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
其他
設(shè)置編碼
# 全局設(shè)置
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES") #支持文件名稱及路徑內(nèi)的中文
gdal.SetConfigOption("SHAPE_ENCODING", "GBK")? ? ? ? #支持屬性字段中的中文
# 為圖層指定編碼
out_layer = ds.CreateLayer( "CommercialHousing",
? ? srs = in_layer.GetSpatialRef(),
? ? geom_type= ogr.wkbPoint,
? ? options? = [
? ? ? ? "ENCODING=UTF-8"
? ? ]
)
1
2
3
4
5
6
7
8
9
10
11
示例
打開shp
from osgeo import ogr
inshp_path = 'D:\mycode\GISandPython\0data\park_point_shp\xiamen_20181128_park.shp'
driver = ogr.GetDriverByName('ESRI Shapefile') #查找一個(gè)特定的驅(qū)動(dòng)程序
datasource = driver.Open(inshp_path, 0) #0只讀脆丁，1可寫
dir( datasource) #使用Python的內(nèi)省函數(shù)dir()查看所有方法
1
2
3
4
5
刪除數(shù)據(jù)
# 刪除矢量數(shù)據(jù)文件
if os.path.exists(out_shp)
driver.DeleteDataSource(out_shp)
ds2 = driver.CreateDataSource(out_shp)
1
2
3
4
遍歷所有屬性值
#【遍歷所有屬性值】
for i in range(feature.GetFieldCount() ):
print( feature.GetField(i) )
### 【查看表的結(jié)構(gòu)世舰，各個(gè)字段的名稱等信息】在layer附加信息中看
layerdef = layer.GetLayerDefn()
for i in range(layerdef.GetFieldCount() ):
defn = layerdef.GetFieldDefn(i)
print(defn.GetName(), defn.GetWidth(), defn.GetType(), defn.GetPrecision() )
1
2
3
4
5
6
7
8
創(chuàng)建
創(chuàng)建點(diǎn)狀數(shù)據(jù)集
from osgeo import ogr
import os,math
driver = ogr.GetDriverByName("ESRI Shapefile")
extfile = 'point_demo.shp'
point_coors = [[300,450], [750, 700], [1200, 450], [750, 200], [750, 450]]
print point_coors
driver = ogr.GetDriverByName("ESRI Shapefile")
if os.access( extfile, os.F_OK ):
? ? driver.DeleteDataSource( extfile )
newds? = driver.CreateDataSource(extfile)
layernew = newds.CreateLayer('point',None,ogr.wkbPoint)
for aa in point_coors:
? ? wkt = 'POINT (' + str(aa[0]) + ' ' + str(aa[1]) + ')'
? ? geom = ogr.CreateGeometryFromWkt(wkt)
? ? feat = ogr.Feature(layernew.GetLayerDefn())
? ? feat.SetGeometry(geom)
? ? layernew.CreateFeature(feat)
newds.Destroy()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
創(chuàng)建線狀數(shù)據(jù)集
from osgeo import ogr
import os,math
driver = ogr.GetDriverByName("ESRI Shapefile")
extfile = 'line_demo.shp'
point_coors = [300,450, 750, 700, 1200, 450, 750, 200]
print point_coors
driver = ogr.GetDriverByName("ESRI Shapefile")
if os.access( extfile, os.F_OK ):
? ? driver.DeleteDataSource( extfile )
newds? = driver.CreateDataSource(extfile)
layernew = newds.CreateLayer('point',None,ogr.wkbLineString)
wkt = 'LINESTRING (%f %f,%f %f,%f %f,%f %f,%f %f)' % (point_coors[0],point_coors[1],
? ? point_coors[2],point_coors[3], point_coors[4],point_coors[5],
? ? point_coors[6],point_coors[7], point_coors[0],point_coors[1])
print(wkt)
geom = ogr.CreateGeometryFromWkt(wkt)
feat = ogr.Feature(layernew.GetLayerDefn())
feat.SetGeometry(geom)
layernew.CreateFeature(feat)
newds.Destroy()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
創(chuàng)建多邊形數(shù)據(jù)集
from osgeo import ogr
import os,math
driver = ogr.GetDriverByName("ESRI Shapefile")
extfile = 'rect_demo.shp'
if os.access( extfile, os.F_OK ):
driver.DeleteDataSource( extfile )
extent = [400, 1100, 300, 600]
newds = driver.CreateDataSource(extfile)
layernew = newds.CreateLayer('rect',None,ogr.wkbPolygon)
width = math.fabs(extent[1]-extent[0])
height = math.fabs(extent[3]-extent[2])
tw = width/2
th = width/2
extnew = extent[0]+tw
wkt = 'POLYGON ((%f %f,%f %f,%f %f,%f %f,%f %f))' % (extent[0],extent[3],
? ? extent[1],extent[3], extent[1],extent[2],
? ? extent[0],extent[2], extent[0],extent[3])
geom = ogr.CreateGeometryFromWkt(wkt)
feat = ogr.Feature(layernew.GetLayerDefn())
feat.SetGeometry(geom)
layernew.CreateFeature(feat)
newds.Destroy()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
創(chuàng)建屬性
from osgeo import ogr
import os,math
driver = ogr.GetDriverByName("ESRI Shapefile") #shp驅(qū)動(dòng)器
extfile = 'rect_field_demo.shp' #輸出的shp路徑
if os.access( extfile, os.F_OK ): #文件是否已存在
? ? driver.DeleteDataSource( extfile )
extent = [400, 1100, 300, 600] #范圍
newds = driver.CreateDataSource(extfile) #按照驅(qū)動(dòng)器創(chuàng)建數(shù)據(jù)源
layernew = newds.CreateLayer('rect',None,ogr.wkbPolygon) #創(chuàng)建圖層
# 創(chuàng)建字段
fieldcnstr = ogr.FieldDefn("fd",ogr.OFTString) #創(chuàng)建字段(字段名，類型)
fieldcnstr.SetWidth(32) #設(shè)置寬度
layernew.CreateField(fieldcnstr) #將字段設(shè)置到layer
fieldf = ogr.FieldDefn("f",ogr.OFTReal)
layernew.CreateField(fieldf)
# wkt格式的polygon
wkt = 'POLYGON ((%f %f,%f %f,%f %f,%f %f,%f %f))' % (extent[0],extent[3],
? ? extent[1],extent[3], extent[1],extent[2],
? ? extent[0],extent[2], extent[0],extent[3])
geom = ogr.CreateGeometryFromWkt(wkt) #根據(jù)wkt創(chuàng)建geometry
feat = ogr.Feature( layernew.GetLayerDefn() ) #從layer中得到feature的要求槽卫，創(chuàng)建一個(gè)Feature對(duì)象
feat.SetField('fd',"這里是字段的值") #設(shè)置字段值
feat.SetGeometry(geom) #設(shè)置要素的geometry屬性
layernew.CreateFeature(feat) #在圖層中創(chuàng)建該feature
newds.Destroy() #刪除（包括保存）
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
創(chuàng)建點(diǎn)Geometry
from osgeo import ogr
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(10,20) #向Point添加多個(gè)點(diǎn)跟压，不會(huì)報(bào)錯(cuò)，但最終只會(huì)用最后添加的點(diǎn)
1
2
3
創(chuàng)建線Geometry
from osgeo import ogr
line = ogr.Geometry(ogr.wkbLineString) # 創(chuàng)建
line.AddPoint(10, 20) # 添加點(diǎn)
line.SetPoint(1, 50, 50) # 修改點(diǎn)坐標(biāo)
line.GetPointCount() # 得到所有點(diǎn)數(shù)目
line.GetX(0); line.GetY(0) # 讀取0號(hào)點(diǎn)的坐標(biāo)
1
2
3
4
5
6
創(chuàng)建多邊形Geometry
from osgeo ipmort ogr
# 創(chuàng)建環(huán)
ring = ogr.Geometry(ogr.wkbLinearRing) #創(chuàng)建
ring.AddPoint(10,10) #添加點(diǎn)
ring.CloseRings() #閉合 或 添加的最后一個(gè)點(diǎn)和第一個(gè)相同
# 將環(huán)添加到多邊形中
polygon = ogr.Geometry(ogr.wkbPolygon)
polygon.AddGeometry(outring)
polygon.AddGeometry(inring)
ring = polygon.GetGeometryRef(0) #得到多邊形內(nèi)的幾何對(duì)象
1
2
3
4
5
6
7
8
9
10
拷貝
datasource層次的拷貝
from osgeo import ogr
import os,math
inshp = "C:\tmp\test.shp"
ds = ogr.Open(inshp)
driver = ogr.GetDriverByName("ESRI Shapefile")
outputfile = 'C:\tmp\test_copy.shp'
if os.access( outputfile, os.F_OK) : #已經(jīng)有該文件
driver.DeleteDataSource(outputfile)? #刪除
pt_cp = driver.CopyDataSource(ds, outputfile) #根據(jù)數(shù)據(jù)源ds創(chuàng)建數(shù)據(jù)歼培，返回指針
pt_cp.Release() #釋放指針震蒋，將數(shù)據(jù)寫入到磁盤
1
2
3
4
5
6
7
8
9
10
layer層次的拷貝
from osgeo import ogr
import os,math
inshp = "C:\tmp\test.shp"
ds = ogr.Open(inshp)
driver = ogr.GetDriverByName("ESRI Shapefile")
outputfile = "C:\tmp\test_copy2.shp"
if os.access( outputfile, os.F_OK):
driver.DeleteDataSourse( outputfile )
layer = ds.GetLayer()
newds = driver.CreateDataSource(outputfile) #首先得有數(shù)據(jù)源
pt_layer? = newds.CopyLayer(layer, 'abcd') #復(fù)制圖層，返回指針
newds.Destroy() #對(duì)newds進(jìn)行Destroy()操作躲庄，才能將數(shù)據(jù)寫入磁盤
1
2
3
4
5
6
7
8
9
10
11
12
feature層次的拷貝
from osgeo import ogr
import os,math
inshp = 'C:\tmp\test.shp'
ds = ogr.Open(inshp)
driver = ogr.GetDriverByName("ESRI Shapefile")
outputfile = 'c:\tmp\test.shp'
if os.access( outputfile, os.F_OK ):
? ? driver.DeleteDataSource( outputfile )
newds = driver.CreateDataSource(outputfile) #按outputfile的數(shù)據(jù)源創(chuàng)建一個(gè)newds
layernew = newds.CreateLayer('worldcopy',None,ogr.wkbLineString) #創(chuàng)建一個(gè)圖層
layer = ds.GetLayer() #得到原數(shù)據(jù)源的圖層
extent = layer.GetExtent() #圖層范圍
# 遍歷舊圖層的所有feature查剖，進(jìn)行復(fù)制
feature = layer.GetNextFeature()
while feature is not None:
? ? layernew.CreateFeature(feature)
? ? feature = layer.GetNextFeature()
newds.Destroy() # 包括了將數(shù)據(jù)flush到磁盤
————————————————
版權(quán)聲明：本文為CSDN博主「GeoDoer」的原創(chuàng)文章，遵循CC 4.0 BY-SA版權(quán)協(xié)議噪窘，轉(zhuǎn)載請(qǐng)附上原文出處鏈接及本聲明笋庄。
原文鏈接：https://blog.csdn.net/summer_dew/article/details/87930241
