python使用gdal對shp讀取,新建和更新的實例
更新時間：2020年03月10日 14:54:30 作者：lionhenryzxxy
這篇文章主要介紹了python使用gdal對shp讀取,新建和更新的實例，具有很好的參考價值，希望對大家有所幫助。一起跟隨小編過來看看吧
昨天要處理一個shp文件，讀取里面的信息，做個計算然后寫到后面新建的field里面。先寫個外面網(wǎng)上都能找到的新建和讀取吧。
1.讀取shp文件
#-*- coding: cp936 -*-
try:
from osgeo import gdal
from osgeo import ogr
exceptImportError:
import gdal
import ogr
defReadVectorFile():
# 為了支持中文路徑，請?zhí)砑酉旅孢@句代碼
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8","NO")
# 為了使屬性表字段支持中文，請?zhí)砑酉旅孢@句
gdal.SetConfigOption("SHAPE_ENCODING","")
strVectorFile ="E:\\Datum\\GDALCsTest\\Debug\\beijing.shp"
# 注冊所有的驅(qū)動
ogr.RegisterAll()
#打開數(shù)據(jù)
ds = ogr.Open(strVectorFile, 0)
if ds == None:
print("打開文件【%s】失??！", strVectorFile)
return
print("打開文件【%s】成功！", strVectorFile)
# 獲取該數(shù)據(jù)源中的圖層個數(shù)，一般shp數(shù)據(jù)圖層只有一個，如果是mdb、dxf等圖層就會有多個
iLayerCount = ds.GetLayerCount()
# 獲取第一個圖層
oLayer = ds.GetLayerByIndex(0)
if oLayer == None:
print("獲取第%d個圖層失??！\n", 0)
return
# 對圖層進行初始化，如果對圖層進行了過濾操作，執(zhí)行這句后，之前的過濾全部清空
oLayer.ResetReading()
# 通過屬性表的SQL語句對圖層中的要素進行篩選，這部分詳細參考SQL查詢章節(jié)內(nèi)容
oLayer.SetAttributeFilter("\"NAME99\"LIKE \"北京市市轄區(qū)\"")
# 通過指定的幾何對象對圖層中的要素進行篩選
#oLayer.SetSpatialFilter()
# 通過指定的四至范圍對圖層中的要素進行篩選
#oLayer.SetSpatialFilterRect()
# 獲取圖層中的屬性表表頭并輸出
print("屬性表結(jié)構(gòu)信息：")
oDefn = oLayer.GetLayerDefn()
iFieldCount = oDefn.GetFieldCount()
for iAttr in range(iFieldCount):
oField =oDefn.GetFieldDefn(iAttr)
print( "%s: %s(%d.%d)" % ( \
oField.GetNameRef(),\
oField.GetFieldTypeName(oField.GetType() ), \
oField.GetWidth(),\
oField.GetPrecision()))
# 輸出圖層中的要素個數(shù)
print("要素個數(shù) = %d", oLayer.GetFeatureCount(0))
oFeature = oLayer.GetNextFeature()
# 下面開始遍歷圖層中的要素
while oFeature is not None:
print("當(dāng)前處理第%d個: \n屬性值：", oFeature.GetFID())
# 獲取要素中的屬性表內(nèi)容
for iField inrange(iFieldCount):
oFieldDefn =oDefn.GetFieldDefn(iField)
line = " %s (%s) = " % ( \
oFieldDefn.GetNameRef(),\
ogr.GetFieldTypeName(oFieldDefn.GetType()))
ifoFeature.IsFieldSet( iField ):
line = line+ "%s" % (oFeature.GetFieldAsString( iField ) )
else:
line = line+ "(null)"
print(line)
# 獲取要素中的幾何體
oGeometry =oFeature.GetGeometryRef()
# 為了演示，只輸出一個要素信息
break
print("數(shù)據(jù)集關(guān)閉！")
2.新建shp文件
#-*- coding: cp936 -*-
try:
from osgeo import gdal
from osgeo import ogr
exceptImportError:
import gdal
import ogr
defWriteVectorFile():
# 為了支持中文路徑，請?zhí)砑酉旅孢@句代碼
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8","NO")
# 為了使屬性表字段支持中文，請?zhí)砑酉旅孢@句
gdal.SetConfigOption("SHAPE_ENCODING","")
strVectorFile ="E:\\TestPolygon.shp"
# 注冊所有的驅(qū)動
ogr.RegisterAll()
# 創(chuàng)建數(shù)據(jù)，這里以創(chuàng)建ESRI的shp文件為例
strDriverName = "ESRIShapefile"
oDriver =ogr.GetDriverByName(strDriverName)
if oDriver == None:
print("%s 驅(qū)動不可用！\n", strDriverName)
return
# 創(chuàng)建數(shù)據(jù)源
oDS =oDriver.CreateDataSource(strVectorFile)
if oDS == None:
print("創(chuàng)建文件【%s】失??！", strVectorFile)
return
# 創(chuàng)建圖層，創(chuàng)建一個多邊形圖層，這里沒有指定空間參考，如果需要的話，需要在這里進行指定
papszLCO = []
oLayer =oDS.CreateLayer("TestPolygon", None, ogr.wkbPolygon, papszLCO)
if oLayer == None:
print("圖層創(chuàng)建失??！\n")
return
# 下面創(chuàng)建屬性表
# 先創(chuàng)建一個叫FieldID的整型屬性
oFieldID =ogr.FieldDefn("FieldID", ogr.OFTInteger)
oLayer.CreateField(oFieldID, 1)
# 再創(chuàng)建一個叫FeatureName的字符型屬性，字符長度為50
oFieldName =ogr.FieldDefn("FieldName", ogr.OFTString)
oFieldName.SetWidth(100)
oLayer.CreateField(oFieldName, 1)
oDefn = oLayer.GetLayerDefn()
# 創(chuàng)建三角形要素
oFeatureTriangle = ogr.Feature(oDefn)
oFeatureTriangle.SetField(0, 0)
oFeatureTriangle.SetField(1, "三角形")
geomTriangle =ogr.CreateGeometryFromWkt("POLYGON ((0 0,20 0,10 15,0 0))")
oFeatureTriangle.SetGeometry(geomTriangle)
oLayer.CreateFeature(oFeatureTriangle)
# 創(chuàng)建矩形要素
oFeatureRectangle = ogr.Feature(oDefn)
oFeatureRectangle.SetField(0, 1)
oFeatureRectangle.SetField(1, "矩形")
geomRectangle =ogr.CreateGeometryFromWkt("POLYGON ((30 0,60 0,60 30,30 30,30 0))")
oFeatureRectangle.SetGeometry(geomRectangle)
oLayer.CreateFeature(oFeatureRectangle)
# 創(chuàng)建五角形要素
oFeaturePentagon = ogr.Feature(oDefn)
oFeaturePentagon.SetField(0, 2)
oFeaturePentagon.SetField(1, "五角形")
geomPentagon =ogr.CreateGeometryFromWkt("POLYGON ((70 0,85 0,90 15,80 30,65 15,700))")
oFeaturePentagon.SetGeometry(geomPentagon)
oLayer.CreateFeature(oFeaturePentagon)
oDS.Destroy()
print("數(shù)據(jù)集創(chuàng)建完成！\n")
3.更新
其實更新無非就是獲取到field然后設(shè)置新值就可以了
其實用SetField()方法就行
import os,sys
from osgeo import gdal
from osgeo import ogr
from osgeo import osr
import numpy
import transformer
# 為了支持中文路徑，請?zhí)砑酉旅孢@句代碼
pathname = sys.argv[1]
choose = sys.argv[2]
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "NO")
# 為了使屬性表字段支持中文，請?zhí)砑酉旅孢@句
gdal.SetConfigOption("SHAPE_ENCODING", "")
# 注冊所有的驅(qū)動
ogr.RegisterAll()
# 數(shù)據(jù)格式的驅(qū)動
driver = ogr.GetDriverByName('ESRI Shapefile')
ds = driver.Open(pathname, update=1)
if ds is None:
print 'Could not open %s'%pathname
sys.exit(1)
# 獲取第0個圖層
layer0 = ds.GetLayerByIndex(0);
# 投影
spatialRef = layer0.GetSpatialRef();
# 輸出圖層中的要素個數(shù)
print '要素個數(shù)=%d'%(layer0.GetFeatureCount(0))
print '屬性表結(jié)構(gòu)信息'
defn = layer0.GetLayerDefn()
fieldindex = defn.GetFieldIndex('x')
xfield = defn.GetFieldDefn(fieldindex)
#新建field
fieldDefn = ogr.FieldDefn('newx', xfield.GetType())
fieldDefn.SetWidth(32)
fieldDefn.SetPrecision(6)
layer0.CreateField(fieldDefn,1)
fieldDefn = ogr.FieldDefn('newy', xfield.GetType())
fieldDefn.SetWidth(32)
fieldDefn.SetPrecision(6)
layer0.CreateField(fieldDefn,1)
feature = layer0.GetNextFeature()
# 下面開始遍歷圖層中的要素
while feature is not None:
# 獲取要素中的屬性表內(nèi)容
x = feature.GetFieldAsDouble('x')
y = feature.GetFieldAsDouble('y')
newx, newy = transformer.begintrans(choose, x, y)
feature.SetField('newx', newx)
feature.SetField('newy', newy)
layer0.SetFeature(feature)
feature = layer0.GetNextFeature()
feature.Destroy()
ds.Destroy()
這里我其實想說最重要的是這個SetFeature()，就是你更新好了field的feature一定要重新set一下，不然是根本起不到任何改變的。新建的時候有createfeature，已經(jīng)設(shè)置了，所以不需要set。
網(wǎng)上的教程都是新建和讀取，都沒有提到這個，結(jié)果自己蠢到試了好久都沒有發(fā)現(xiàn)問題在哪，以為是什么數(shù)據(jù)類型與設(shè)置字段屬性不匹配，一頭霧水哈哈哈。
補充知識：python使用GDAL生成shp文件
GDAL是一個開源的地理工具包，其支持基本所有的地理操作，其有python、java、c等語言包，是地理信息C端開發(fā)不可越過的工具，鑒于python語言的簡單性，這里使用python中GDAL包來進行shp文件的生成，這里本質(zhì)是利用ogc地理標(biāo)準(zhǔn)的坐標(biāo)字符串來生成shp。
第一步：安裝GDAL環(huán)境，建議下載后，本地安裝，注意與python版本號要對應(yīng)，可參考網(wǎng)上教程。
第二部：代碼分析
引入GDAL工具包
import osgeo.ogr as ogr
import osgeo.osr as osr
注冊驅(qū)動，這里是ESRI Shapefile類型，并設(shè)置shp文件名稱
driver = ogr.GetDriverByName("ESRI Shapefile")
data_source = driver.CreateDataSource("ceshi.shp")
注入投影信息，這里使用的4326，表示經(jīng)緯度坐標(biāo)，根據(jù)情況可以自行更改
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
這里定義的是，生成的要素類型，包括點、線、面
#ogr.wkbPoint 點 #ogr.wkbLineString 線 #ogr.wkbMultiPolygon 面
這里的圖層名稱要與上面注冊驅(qū)動的shp名稱一致
layer = data_source.CreateLayer("ceshi", srs, ogr.wkbLineString)
這里設(shè)置要素的屬性字段，其中設(shè)置了兩個字段，分別是Name、data，其中ogr.OFTString表示字符串類型，其長度都是14字節(jié)，可自行設(shè)置寬度
field_name = ogr.FieldDefn("Name", ogr.OFTString)
field_name.SetWidth(14)
layer.CreateField(field_name)
field_name = ogr.FieldDefn("data", ogr.OFTString)
field_name.SetWidth(14)
layer.CreateField(field_name)
在生成的字段名中插入要素值，即屬性表中每行的值
feature = ogr.Feature(layer.GetLayerDefn())
feature.SetField("Name", "ceshi")
feature.SetField("data", "1.2")
核心部分，生成line數(shù)據(jù)
其中各要素格式如下：
POINT(6 10) LINESTRING(3 4,10 50,20 25) POLYGON((1 1,5 1,5 5,1 5,1 1),(2 2,2 3,3 3,3 2,2 2)) MULTIPOINT(3.5 5.6, 4.8 10.5) MULTILINESTRING((3 4,10 50,20 25),(-5 -8,-10 -8,-15 -4)) MULTIPOLYGON(((1 1,5 1,5 5,1 5,1 1),(2 2,2 3,3 3,3 2,2 2)),((6 3,9 2,9 4,6 3))) GEOMETRYCOLLECTION(POINT(4 6),LINESTRING(4 6,7 10)) POINT ZM (1 1 5 60) POINT M (1 1 80)
需要注意的是，這里應(yīng)該與上面定義的生成要素的類型保持一致，最后是清空緩存，這里多說一句，字符串語法與postgis等開源gis一致，都遵循ogc國際標(biāo)準(zhǔn)
wkt = 'LINESTRING(3 4,10 50,20 25)' line = ogr.CreateGeometryFromWkt(wkt) feature.SetGeometry(line) layer.CreateFeature(feature) feature = None data_source = None
結(jié)果如下：
用arcgis打開
可以使用該方法，下載在線shp數(shù)據(jù)，只需要知道所需要素的geojson格式數(shù)據(jù)中坐標(biāo)串即可?；蛘邎D像識別中獲取的矢量邊界賦予經(jīng)緯度。
以上這篇python使用gdal對shp讀取,新建和更新的實例就是小編分享給大家的全部內(nèi)容了，希望能給大家一個參考，也希望大家多多支持腳本之家。
相關(guān)文章
這篇文章主要為大家詳細介紹了Python中函數(shù)的基本使用教程，文中示例代碼介紹的非常詳細，具有一定的參考價值，感興趣的小伙伴們可以參考一下，希望能夠給你帶來幫助2022-04-04
python數(shù)據(jù)可視化使用pyfinance分析證券收益示例詳解
這篇文章主要為大家介紹了python數(shù)據(jù)可視化使用pyfinance分析證券收益的示例詳解及pyfinance中returns模塊的應(yīng)用，有需要的朋友可以借鑒參考下，希望能夠有所幫助2021-11-11
隨著互聯(lián)網(wǎng)技術(shù)的發(fā)展,越來越多的人選擇通過網(wǎng)絡(luò)平臺購票,特別是在中國,12306作為官方火車票預(yù)訂平臺,承擔(dān)了巨大的訪問量,對于熱門線路或者節(jié)假日出行,往往會出現(xiàn)一票難求的情況,因此,一些技術(shù)愛好者嘗試?yán)镁幊陶Z言如Python來開發(fā)搶票腳本2025-01-01
Python3.5編程實現(xiàn)修改IIS WEB.CONFIG的方法示例
這篇文章主要介紹了Python3.5編程實現(xiàn)修改IIS WEB.CONFIG的方法,涉及Python針對xml格式文件的讀寫以及節(jié)點操作相關(guān)技巧,需要的朋友可以參考下2017-08-08
這篇文章主要介紹了python實現(xiàn)給數(shù)組按片賦值的方法,實例分析了Python在指定位置進行賦值的相關(guān)技巧,需要的朋友可以參考下2015-07-07
python中zip和unzip數(shù)據(jù)的方法
這篇文章主要介紹了python中zip和unzip數(shù)據(jù)的方法,實例分析了Python中zlib模塊的相關(guān)使用技巧,需要的朋友可以參考下2015-05-05
這篇文章主要和大家介紹一下異常的最后一個知識點——斷言 ，斷言是判斷一個表達式，在表達式為 False 的時候觸發(fā)異常。本文將通過示例詳細介紹一下斷言，需要的可以參考一下2022-03-03
Python Pillow Image.save 保存為jpg圖片壓縮問題
Pillow 庫支持多種圖片格式，Pillow 能夠很輕松地實現(xiàn)圖片格式之間的轉(zhuǎn)換。本文就來詳細的介紹一下Image.save的具體使用，感興趣的可以了解一下2021-11-11
PyCharm可以說是當(dāng)今最流行的一款Python?IDE了,下面這篇文章主要給大家介紹了關(guān)于pycharm安裝中文插件的2種方法,文中通過圖文介紹的非常詳細,需要的朋友可以參考下2023-06-06
python編程使用PyQt制作預(yù)覽窗口游戲中的小地圖
這篇文章主要為大家介紹了python使用PyQt制作預(yù)覽窗口游戲中的小地圖實現(xiàn)示例過程，有需要的朋友可以借鑒參考下，希望能夠有所幫助，祝大家多多進步2021-10-10