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(发邮件时,请将#更换为@)进行举报,并提供相关证据,一经查实,本站将立刻删除涉嫌侵权内容。