python使用gdal对shp读取,新建和更新的实例
昨天要处理一个shp文件，读取里面的信息，做个计算然后写到后面新建的field里面。先写个外面网上都能找到的新建和读取吧。
1.读取shp文件
#-*-coding:cp936-*-
try:
fromosgeoimportgdal
fromosgeoimportogr
exceptImportError:
importgdal
importogr
defReadVectorFile():
#为了支持中文路径，请添加下面这句代码
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8","NO")
#为了使属性表字段支持中文，请添加下面这句
gdal.SetConfigOption("SHAPE_ENCODING","")
strVectorFile="E:\\Datum\\GDALCsTest\\Debug\\beijing.shp"
#注册所有的驱动
ogr.RegisterAll()
#打开数据
ds=ogr.Open(strVectorFile,0)
ifds==None:
print("打开文件【%s】失败！",strVectorFile)
return
print("打开文件【%s】成功！",strVectorFile)
#获取该数据源中的图层个数，一般shp数据图层只有一个，如果是mdb、dxf等图层就会有多个
iLayerCount=ds.GetLayerCount()
#获取第一个图层
oLayer=ds.GetLayerByIndex(0)
ifoLayer==None:
print("获取第%d个图层失败！\n",0)
return
#对图层进行初始化，如果对图层进行了过滤操作，执行这句后，之前的过滤全部清空
oLayer.ResetReading()
#通过属性表的SQL语句对图层中的要素进行筛选，这部分详细参考SQL查询章节内容
oLayer.SetAttributeFilter("\"NAME99\"LIKE\"北京市市辖区\"")
#通过指定的几何对象对图层中的要素进行筛选
#oLayer.SetSpatialFilter()
#通过指定的四至范围对图层中的要素进行筛选
#oLayer.SetSpatialFilterRect()
#获取图层中的属性表表头并输出
print("属性表结构信息：")
oDefn=oLayer.GetLayerDefn()
iFieldCount=oDefn.GetFieldCount()
foriAttrinrange(iFieldCount):
oField=oDefn.GetFieldDefn(iAttr)
print("%s:%s(%d.%d)"%(\
oField.GetNameRef(),\
oField.GetFieldTypeName(oField.GetType()),\
oField.GetWidth(),\
oField.GetPrecision()))
#输出图层中的要素个数
print("要素个数=%d",oLayer.GetFeatureCount(0))
oFeature=oLayer.GetNextFeature()
#下面开始遍历图层中的要素
whileoFeatureisnotNone:
print("当前处理第%d个:\n属性值：",oFeature.GetFID())
#获取要素中的属性表内容
foriFieldinrange(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("数据集关闭！")
2.新建shp文件
#-*-coding:cp936-*-
try:
fromosgeoimportgdal
fromosgeoimportogr
exceptImportError:
importgdal
importogr
defWriteVectorFile():
#为了支持中文路径，请添加下面这句代码
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8","NO")
#为了使属性表字段支持中文，请添加下面这句
gdal.SetConfigOption("SHAPE_ENCODING","")
strVectorFile="E:\\TestPolygon.shp"
#注册所有的驱动
ogr.RegisterAll()
#创建数据，这里以创建ESRI的shp文件为例
strDriverName="ESRIShapefile"
oDriver=ogr.GetDriverByName(strDriverName)
ifoDriver==None:
print("%s驱动不可用！\n",strDriverName)
return
#创建数据源
oDS=oDriver.CreateDataSource(strVectorFile)
ifoDS==None:
print("创建文件【%s】失败！",strVectorFile)
return
#创建图层，创建一个多边形图层，这里没有指定空间参考，如果需要的话，需要在这里进行指定
papszLCO=[]
oLayer=oDS.CreateLayer("TestPolygon",None,ogr.wkbPolygon,papszLCO)
ifoLayer==None:
print("图层创建失败！\n")
return
#下面创建属性表
#先创建一个叫FieldID的整型属性
oFieldID=ogr.FieldDefn("FieldID",ogr.OFTInteger)
oLayer.CreateField(oFieldID,1)
#再创建一个叫FeatureName的字符型属性，字符长度为50
oFieldName=ogr.FieldDefn("FieldName",ogr.OFTString)
oFieldName.SetWidth(100)
oLayer.CreateField(oFieldName,1)
oDefn=oLayer.GetLayerDefn()
#创建三角形要素
oFeatureTriangle=ogr.Feature(oDefn)
oFeatureTriangle.SetField(0,0)
oFeatureTriangle.SetField(1,"三角形")
geomTriangle=ogr.CreateGeometryFromWkt("POLYGON((00,200,1015,00))")
oFeatureTriangle.SetGeometry(geomTriangle)
oLayer.CreateFeature(oFeatureTriangle)
#创建矩形要素
oFeatureRectangle=ogr.Feature(oDefn)
oFeatureRectangle.SetField(0,1)
oFeatureRectangle.SetField(1,"矩形")
geomRectangle=ogr.CreateGeometryFromWkt("POLYGON((300,600,6030,3030,300))")
oFeatureRectangle.SetGeometry(geomRectangle)
oLayer.CreateFeature(oFeatureRectangle)
#创建五角形要素
oFeaturePentagon=ogr.Feature(oDefn)
oFeaturePentagon.SetField(0,2)
oFeaturePentagon.SetField(1,"五角形")
geomPentagon=ogr.CreateGeometryFromWkt("POLYGON((700,850,9015,8030,6515,700))")
oFeaturePentagon.SetGeometry(geomPentagon)
oLayer.CreateFeature(oFeaturePentagon)
oDS.Destroy()
print("数据集创建完成！\n")
3.更新
其实更新无非就是获取到field然后设置新值就可以了
其实用SetField()方法就行
importos,sys
fromosgeoimportgdal
fromosgeoimportogr
fromosgeoimportosr
importnumpy
importtransformer
#为了支持中文路径，请添加下面这句代码
pathname=sys.argv[1]
choose=sys.argv[2]
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8","NO")
#为了使属性表字段支持中文，请添加下面这句
gdal.SetConfigOption("SHAPE_ENCODING","")
#注册所有的驱动
ogr.RegisterAll()
#数据格式的驱动
driver=ogr.GetDriverByName('ESRIShapefile')
ds=driver.Open(pathname,update=1)
ifdsisNone:
print'Couldnotopen%s'%pathname
sys.exit(1)
#获取第0个图层
layer0=ds.GetLayerByIndex(0);
#投影
spatialRef=layer0.GetSpatialRef();
#输出图层中的要素个数
print'要素个数=%d'%(layer0.GetFeatureCount(0))
print'属性表结构信息'
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()
#下面开始遍历图层中的要素
whilefeatureisnotNone:
#获取要素中的属性表内容
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，已经设置了，所以不需要set。
网上的教程都是新建和读取，都没有提到这个，结果自己蠢到试了好久都没有发现问题在哪，以为是什么数据类型与设置字段属性不匹配，一头雾水哈哈哈。
补充知识：python使用GDAL生成shp文件
GDAL是一个开源的地理工具包，其支持基本所有的地理操作，其有python、java、c等语言包，是地理信息C端开发不可越过的工具，鉴于python语言的简单性，这里使用python中GDAL包来进行shp文件的生成，这里本质是利用ogc地理标准的坐标字符串来生成shp。
第一步：安装GDAL环境，建议下载后，本地安装，注意与python版本号要对应，可参考网上教程。
第二部：代码分析
引入GDAL工具包
importosgeo.ograsogr
importosgeo.osrasosr
注册驱动，这里是ESRIShapefile类型，并设置shp文件名称
driver=ogr.GetDriverByName("ESRIShapefile")
data_source=driver.CreateDataSource("ceshi.shp")
注入投影信息，这里使用的4326，表示经纬度坐标，根据情况可以自行更改
srs=osr.SpatialReference()
srs.ImportFromEPSG(4326)
这里定义的是，生成的要素类型，包括点、线、面
#ogr.wkbPoint点 #ogr.wkbLineString线 #ogr.wkbMultiPolygon面
这里的图层名称要与上面注册驱动的shp名称一致
layer=data_source.CreateLayer("ceshi",srs,ogr.wkbLineString)
这里设置要素的属性字段，其中设置了两个字段，分别是Name、data，其中ogr.OFTString表示字符串类型，其长度都是14字节，可自行设置宽度
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数据
其中各要素格式如下：
POINT(610) LINESTRING(34,1050,2025) POLYGON((11,51,55,15,11),(22,23,33,32,22)) MULTIPOINT(3.55.6,4.810.5) MULTILINESTRING((34,1050,2025),(-5-8,-10-8,-15-4)) MULTIPOLYGON(((11,51,55,15,11),(22,23,33,32,22)),((63,92,94,63))) GEOMETRYCOLLECTION(POINT(46),LINESTRING(46,710)) POINTZM(11560) POINTM(1180)
需要注意的是，这里应该与上面定义的生成要素的类型保持一致，最后是清空缓存，这里多说一句，字符串语法与postgis等开源gis一致，都遵循ogc国际标准
wkt='LINESTRING(34,1050,2025)' line=ogr.CreateGeometryFromWkt(wkt) feature.SetGeometry(line) layer.CreateFeature(feature) feature=None data_source=None
结果如下：
用arcgis打开
可以使用该方法，下载在线shp数据，只需要知道所需要素的geojson格式数据中坐标串即可。或者图像识别中获取的矢量边界赋予经纬度。
以上这篇python使用gdal对shp读取,新建和更新的实例就是小编分享给大家的全部内容了，希望能给大家一个参考，也希望大家多多支持毛票票。
声明：本文内容来源于网络，版权归原作者所有，内容由互联网用户自发贡献自行上传，本网站不拥有所有权，未作人工编辑处理，也不承担相关法律责任。如果您发现有涉嫌版权的内容，欢迎发送邮件至：czq8825#qq.com（发邮件时，请将#更换为@）进行举报，并提供相关证据，一经查实，本站将立刻删除涉嫌侵权内容。