Python GDAL在GIS数据处理中的应用.pdf

GI S的数据可以有多存在方式，可以将数据以某种方式存成TXT,EXCEL等,GI S数据的表现形式往往是衣服很漂亮的地图，而数据也是GI S的核心，GI S的数据往往是很多的，我们在做数据处理的时候往往希望能够批量处理而数据处理从某种意义上来说就是按照一定的规则和方法对其进行加工，通 过这样的操作使得这些数据成为符合我们项目需求的成果。数据处理中经常遇到 的就是数据转换，GI S的数据格式太多了Python作为一门很强大的开源的动态语言，以其灵活型和间接性征服了很 多人，GD AL将空间数据抽象，用一个抽象的模型来表达空间数据的实体，同时 它提供了很多操作数据的函数，0GR作为GD AL的分支提供了很多对矢量数据的 操作。这些心得是我在学习Python和GD AL的时候所做的一些记录，希望这些能给 大家学习上或工作中带来帮助，这些心得分对矢量数据和栅格数据的处理做了简 单的介绍，也通过两个完整的例子，说明了 GD AL的强大。语言是一中工具，用 于将人的大脑中的想法，告诉给计算机的工具。一门语言要学精通是很难的，我 常常在和朋友交流的时候说，语言，了解了它的特点，语法，能简单的写出几个 像样的函数，至于更深入的，用的时候借助GOOGLE等搜索引擎。GD AL作为空 间数据描述处理的工具，对它提供的函数，我只是用到了其中的冰山一角，我也 和很多人一样，是初学者，也需要不断的努力，这篇心得中的两个例子分别对矢 量数据的处理和栅格数据的处理做了介绍，这两个例子并不是多么的深奥，只是 为了体验，就像我前面说的，发现这种数据处理的规则，然后借助这样的工具，帮助我们取得学习上或者工作上的胜利在接触一个新的东西的时候，我喜欢从这新东西的一些名词和结构入手。记 得以前一位老师说过，一个新的东西，肯定会有自己的专业术语，然后就是体系 结构等，我们就从这个开始吧GD AL(Geospatial D ata Abstraction Library)是一个在 X/MI T 许可协议下 的开源栅格空间数据转换库。它利用抽象数据模型来表达所支持的各种文件格式。它还有一系列命令行工具来进行数据转换和处理。0GR是GD AL项目的一个分支，功能与GD AL类似，只不过它提供对矢量数据 的支持。有很多著名的GI S类产品都使用了 GD AL/0GR库，包括ESRI的ArcGI S,Google Earth和跨平台的GRASS GI S系统。利用GD AL/0GR库，可以使基于Linux的地理空间数据管理系统提供对矢量 和栅格文件数据的支持。1.GD ALgdal(Geospatial Data Abstraction Library)提供对多种栅格 数据的支持，包括 Arc/I nfo ASCH Grid(asc),GeoTiff(tiff),Erdas I magine I mages(img),ASCI I D EM(dem)等格式。GD AL使用抽象数据模型(Abstract D ata Model)来解析它所支持的数据格式,抽象数据模型包括数据集(D ataset),坐标系统，仿射地理坐标转换(Affine Ge。Transform),大地控制点(GCPs),元数据(Metadata),栅格波段(Raster Band),颜色表(Color Table),子数据集域(Subdatasets D omain),图像结构域(I mage_Structure D omain),XML 域(XML:D omains)。详细的结构描述请访问http:/www.gdal.org/gdal_datamodel.htmlGD AL的核心类结构设计如图所示:其中的类说明如下：GD ALMajorObject类：带有元数据的对象。GD ALD dataset类：通常是从一个栅格文件中提取的相关联的栅格波段集合 和这些波段的元数据；GD ALD dataset也负责所有栅格波段的地理坐标转换(georeferencing transform)和坐标系定义。GD ALD river类：文件格式驱动类，GD AL会为每一个所支持的文件格式创建 一个该类的实体，来管理该文件格式。GD ALD riverManager类：文件格式驱动管理类，用来管理GD ALD river类。2.OGROGR提供对矢量数据格式的读写支持，它所支持的文件格式包括：ESRI Shapefiles,S-57,SD TS,PostGI S,Oracle Spatial,Mapinfo mid/mif,Mapinfo TAB。Format NameCodeCreation Georeferencing:Arc/Info Binary Co verag eAVCBinjNo)Yes:Arc/Info E00(ASCII)Co verag eAVCEOONo|YesAtlas BNABNA|Yes岫:Auto CAD DXFDXFYes|NoCo iro na S巳parat巳d Vallie(csv)CSVYesNo|56DS/OPeNDAPDODSjNo|YesEERI P巳iso nal G巳o Databas巳PGeoNo)YesiESRI ArcDfSDENoiYesiESRI ShapefileESRI ShapefileYes|YesiFMEUbj 巳cts Gat巳wayFMEObjects GatewayiNo|YesGeo JSON|Yes)YesG&o c口nc已pt Ex pert-Geo co ncept(Yes|YesGoRSS-Geo RSSYes|YesfeMLGML|Yes|Yes丽GMT|Yes|YesiGPSBabelGPSBabel|Yes|YesIgpxGPXYes|Yes庇ASSGRASSNo|YesGPSTrackMaker(g tm,g tz)GPSTrackMaker|YesY esMydro g raphic Transfer Fo rmatHTFNoiYesnfo rniix DataBladeIDBiYesY esIlNTERLISZ/Interlis 1 and Interlis 2|Yes|Yes更详细的请访问 http:/gdal.org/ogr/ogr_formats.htmlOGR包括如下几部分：Geometry：类 Geometry(包括 OGRGeometry 等类)封装了 OpenGI S 的矢量数 据模型，并提供了一些几何操作，WKB(Well Knows Binary)和 WKT(Well Known Text)格式之间的相互转换，以及空间参考系统(投影)。Spatial Reference：类OGRSpatialReference封装了投影和基准面的定义。Feature：类OGRFeature封装了一*个完整Feature的定义，一*个完整的 Feature 包括Geometry 和 Geometry 的一系列属性。Feature D efinition：类 OGRFeatureD efn 里面封装了 feature 的属性，类 型、名称及其默认的空间参考系统等。一个OGRFeatureD efn对象通常与一个层(layer)对应。Layer：类OGRLayer是一个抽象基类，表示数据源类OGRD ataSource里面的 一层要素(Feature)0D ata Source：类OGRD ataSource是一个抽象基类，表示含有OGRLayer对象 的一个文件或一个数据库。D rivers：类OGRSFD river对应于每一个所支持的矢量文件格式。类 OGRSFD river 由类 OGRSFD riverRegistrar 来注册和管理。2)OGR 的 Geometry 模型OGR的Geometry模型是建立在OpenGI S的简单要素数据模型之上的。如下图OGRGeometryI；f；,OGRCurve OGRGeornetryCollection OGRPoint ijGRSurfaceOGRLineString OGRMultiLineString OGRMultiPoint OGRMultiPolygon OGRPolygonOGRLinearRing所示:图-OGR的Geometry模型关系图图-OpenGI S的简单要素数据模型由上面两图的对比，可以清楚的看到，OGR的Geomey模型是严格遵循 OpenGI S的简单要素数据规范的。OGR的Geometry模型不仅在继承体系上与 OpenGI S的简单要素数据模型一致，在函数接口上也向其靠拢，从基本的获取 Geometry 对象信息的方法如 D imension()、GeometryType()、SRI D()、Envelope()AsText()Boundary()等到判定空间未知关系的方法如 Equals(anotherGeometry:Geometry)、D isjoint(anotherGeometry:Geometry)、I ntersects(anotherGeometry:Geometry)、Touches(anotherGeometry:Geometry)等都是符合其标准的。SUBCLASSESSUPERCLASSt tST_MultiLineString ST MultiPolygonESRI的ST_Geometry数据存储模型3.安装ArcGI S 10 D esktop装了 Python2.6,那么我们就不用装这个了，然后去下面的 地址找到 gdal 的 Python 包 GD AL-1.6.1.win32-py2.6.exe,自动安装到 python 的安装目录 E:Python26 ArcGI S10.OLibsite-packagesc CM.GDAL下载下面的文件(GDAL Windows binaries)tmW of-lArrwoh gvetw4Index of/gdal/win32/l.6X 4nib I%!扈P，T，lp4iuxJ PatE DyetJnn缶3或一一 _.，二 U-；an 2009!051 S1K15 Fai；：的 a 30-Dec-JOOR MJfrom osgeo import gdalfrom osgeo.gdalconst import*dataset=gdal.Open(Z，C:dataliuyu.jpg,GA ReadOnly)dataset.GetD river().ShortN ame4.矢量篇对矢量数据的操作时通过OGR,对每种数据有一个驱动，因此在操作数据的时候,我们先要获取相应数据格式的驱动：如driver=ogr.GetD riverByN ame(ESRI Shapefile)，但是对于栅格数据略有不同；GD AL数据驱动，与0GR数据驱动 类似，需要先创建某一类型的数据驱动，再创建相应的栅格数据集。一次性注册 所有的数据驱动，但是只能读不能写：gdal.AllRegister().单独注册某一类型 的数据驱动，这样的话可以读也可以写，可以新建数据集：driver=gdal.GetD riverByN ame(HFA)driver.Register()要读取某种类型的数据，必须要先载入数据驱动，也就是初始化一个对象，让它“知道”某种数据结构。例如：driver=ogr.GetD riverByN ame(，ESRI Shapefile，)一般ESRI的shapefile都是填0的，如果不填的话默认也是0.4.1创建文件 datasource=driver.CreateD ataSource(C:liuyu.shp)shplayer=datasource.Cr eat eLayer(，liuyu，)geom_type=ogr.wkbPoint)_ liuyu.shp _ liuyu.shx,liuyu.dbf shplayer=datasource.CreateLayer(liuyu”,geom_type=ogr.wkbPoint)shplayer相当于我们的FeatureClass,因为我们将要素写入文件是通过 shplayer的，还记得C#和Engine是如何操作数据的要添加一个新字段，只能在layer里面加，而且还不能有数据添加的字段如果是字符串，还要设定宽度fieldD efn=ogr.FieldD efn(，id5,ogr.OFTString)fieldD efn.SetWidth(4)shplayer.CreateField(fieldD efn)添加一个新的feature,首先得完成上一步，把字段field都添加齐了然后从shpfile中读取相应的feature类型,并创建featurefieldD efn=ogr.FieldD efn(，id5,ogr.OFTString)f ieldD efn.SetWidth(4)shplayer.CreateField(fieldD efn)featureD efn=shplayer.feature=ogr.Feature(featureD efn)将 feature 写入 shplayer4.2建立点对象 feature.SetGeometry(point)shpfile.CreateFeature(feature)n=shplayer.GetFeatureCount()field=feature.SetField(id,liu)刚才说过Engine中操作数据是通过FeatureClass的，而通过上面的Python可 以看出Python中操作数据是通过这个shplayer的。/IFeatureClass CreateFeature Examplepublic void IFeatureClass_CreateFeature_Example(IFeatureClass featureclass)/Function is designed to work with polyline dataif(featureclass.ShapeType!=ESRI,ArcGIS.Geometry.esriGeometryType.esriGeometryPolyline)return;/create a geometry for the features shapeESRI.ArcGIS.Geometry.IPolyline polyline=newESRI.ArcGIS.Geometry.PolylineClass();ESRI.ArcGIS.Geometry.IPoint point=new ESRI.ArcGIS.Geometry.Pointclass();point.X=0;point.Y=0;polyline.FromPoint=point;point=new ESRI.ArcGIS.Geometry.Pointclass();point.X=10;point.Y=10;polyline.ToPoint=point;IFeature feature=featureClass.CreateFeature();/Apply the constructed shape to the new features shape feature.Shape=polyline;ISubtypes subtypes=(ISubtypes)featureclass;IRowSubtypes rowSubtypes=(IRowSubtypes)feature;if(subtypes.HasSubtype)/does the feature class have subtypes?(rowSubtypes.SubtypeCode=1;/in this example 1 represents the Primary Pipeline subtype)/initalize any default values that the feature has rowSubtypes.InitDefaultvalues();/Commit the default values in the feature to the databasefeature.Store();/update the value on a string field that indicates who installed the feature.feature.set_Value(feature.Fields.FindField(nInstalledByn),”K Johnston”);/Commit the updated values in the feature to the databasefeature.Store();如果用我们ArcGI S装了之后的那个Python环境敲写代码，一两句的话，还能忍受，但是如果大量的话，那可是很痛苦的，所以我们用Eclipse作为环境通过选择help（帮助）-software updates（软件更新）Find and I nsatllSearch for new Feature to I nstall 下一中添力口个 N ew Remote Site 为http: 直接从网站上下载安装的。这个和Eclipse安装其他插件一样，或者将这个插件 单独下载，然后将插件的文件复制到Eclipse的相关目录中.我们通过下面这个例子将上面提到的，做个小结。下面这个数据是很简单的，前面六位（1-6）表示X坐标，而紧接着X的六位（7-1 2）表示丫坐标，我们就将这个文件利用OGR弄成shp数据，下面我们看下这个数据：E points.匿|叵|区|File Edit Format View Help28138523045607D 人 22135572665697E 25780592415434A 31533212690568A 30643552333409C31765882304171C34180122337232E 33568792068170E41351762712735A 12238273635520A 26855272569021E40464382558364B18016822735241E vimport sysfrom osgeo import ogrogr.RegisterAll()driver=ogr.GetDriverByName(ESRI Shapefile)#C:txt2shp.shp this is wrongds=driver.CreateDataSource(C:liuyu.shp”)layer=ds.CreateLayer(rrliuyu”,geom_type=ogr.wkbPoint)#Create fieldsfieldl=ogr.FieldDefn(Num”,ogr.OFTInteger)field2=ogr.FieldDefn(nCharrogr.OFTString)field2.Setwidth(4)layer.CreateField(fieldl)layer.CreateField(field2)featuredefn=layer.GetLayerDefn()featu re=ogr.Featu re(featuredefn)#open txttxtfile=open(,rC:points.cspnr r)line=txtfile.readline()while line:cleanline=line.replace(r n r r r r)point=ogr.Geometry(ogr.wkbPoint)point.AddPoint(float(line0:6),float(line6:12)feature.SetGeometry(point)feature.SetField(0,int(line12:14)feature.SetField(1,line14:-1)layer.CreateFeature(feature)line=txtfile.readline()卜面的是图形建立新的儿何形状建立空的 geometry 对象：ogr.Geometry定义各种不同的geometry使用的方法是不一样的(point,line,polygon,etc)新建点，如我们上面的例子point=ogr.Geometry(ogr.wkbPoint)point.AddPoint(float(line0:6),float(line6:12)4.3构建面对于面的构建，如果用过AE的话，我们都知道，面是有一组ring构成册，而 ring又是由segment构成，segment又是由point构成，我大学的时候老师说过 这样一句话矢量数据的本质是point,而栅格数据的本质是pixel。我们看一下 AE中的这个图：那么OGR中的polygon呢？它的结构是什么样子？我们从下面的例子是否可以从例子中看到？0GR中Polygon是由ring构成的。#coding=gbk#面是由环构成的polygon.AddGeometry(outring)polygon.AddGeometry(inring)featuredefn=shpfile.GetLayerDefn()feature=ogr.Feature(featuredefn)faatmr.SatGomatry(polygon)print polygon.GetArea()F面这句话可以帮我们数数数polygon能有几个ring print polygon.GetGeometryCount()从polygon中读取ring,index的顺序和创建polygon时添加ring的顺序相同outring=polygon.GetGeometryRef(0)inring=polygon.GetGeometryRef(1)4.4 关于投影空间数据有一个东西非常重要,那就是投影,而0GR支持多种多样的Projections,GD AL 支持 WKT,PR0J.4,ESPG,USGS,ESRI.prj,使用 SpatialReference 对象可以从layer和Geometry中读取Projections,例如:spatialRef=layer.GetSpatialRef()spatialRef=geom.GetSpatialReference()投影信息一般存储在.pij文件中，如果没有这个文件，上述函数返回N one首先导入 osr 库，之后使用 osr.SpatialReference()创建 SpatialReference 对象，之后用下列语句向SpatialReference对象导入投影信息3 ImportFromEPSGQ ImportFromEPSGO ImportFromESRIO ImportFromMICoordSysO ImportFromPCIO心 ImportFromProj40 ImportFromUrlO ImportFromUSGSO。ImportFromWktO ImportFromXMLQ导出Projection,使用下面的语句可以导出为字符串0 ExportToMICoordSysQ ExportToPCIO ExportToPrettyWktO心 ExportToPrqj40 ExportToUSGSO ExportToWktO%ExportToXMLQ怎么给这刚才生成的文件添加投影信息？我们知道可以给一个Geometry指定投 影信息？想一下ESRI是怎么搞的？没错PRJ文件sr=osr.SpatialReference()sr,ImportFromEPSG(32612)#如果将geometry的投影者B改变了，那要一个个的循环，数据量太多的话，很麻烦，可以像ESR工 一样提供一个文件，写进去pr j f ile=open(rrC:liuyu.prj,rwr)sr.MorphToESRI()prj file.write(sr.ExportToWkt()prj file.close()print上面的代码执行后，我们就可以看到下面的文件名称liuyu.dbf囱 liuyu.prj liuyu.shp；,liuyu.shx厘 二 I I 回 w(v)PROJCSfWGS 1984 IfTU Zone I 5N*,GEOGCS-CCS WCS I 9H4DATUM L-D WGS I9H4,SPHEROID WGS I 9H4”,6378137,29R.滤7223563J.PR【MEM-Grccnwlch=.0,UNITLDcard.0.01 74B32926 I 9943295 .PROJ KCT ION I Trnniiverne MnrcntorJ.PAKABKTKKnt i tudr:of_oriRin*.Oj.PAKAMKTKK t centr nl mcr i di mi*.11 I J.VARA1KTER Kacnlc.factor*.0.9996 I.PAKAUKTEK(nl oe?_cnn t inn.500000|.rAKAHETEK I*ffal tie_northinR4*,0 I.UNIT I Metor”.1 I I怎么做投影变换？用这个方法TransformTo-Python的函数function,异常exception和模块module可以从任何一本python 教材上找到，在此不多赘述。5.栅格篇GD AL原生支持超过100种栅格数据类型，涵盖所有主流GI S与RS数据格式，如 下图：完整的支持列表可以参考http的支ww.g表1.or参formats list,html我们知道对于栅格数据，从数学的角度来考虑的话就是一个矩阵，从我们在学习 栅格数据和矢量数据的时候，我们知道栅格数据便于建立模型，这是因为他的结 构的原因。5.1栅格数据操作GD AL数据驱动，与OGR数据驱动类似，需要先创建某一类型的数据驱动，再创 建响应的栅格数据集。一次性注册所有的数据驱动，但是只能读不能写：gdal.AURegisterO单独注册某一类型的数据驱动，这样的话可以读也可以写，可以新建数据集：driver=gdal.GetD riverByN ame C GTiff)driver.Register()from osgeo import gdal from osgeo.gdalconst import*driver=gdal.GetD riverByN ame(JPEG)driver.Register()ds二gdal.Open(Z，C:dataliuyu.jpg,GA ReadOnly)band=ds.GetRasterBand(1)获取第一个波段 rown=band.XSize 没注意看这个到底是行还是列 cown=band.YSize a=band.ReadAsArray()将波段的数据用数组表示 aarray(5,4,4,.,53,55,57,5,5,5,48,49,51,6,6,7,45,45,46 ,81,81,82,6 1,6 2,6 4,81,82,85,6 7,6 6,6 6 ,87,87,87,72,6 8,6 5 ,dtype=uint8)通过这样，我们就可以对这个数组操作，然后将数组在写回去。Engine中是怎么读取栅格数据的？I RasterProps rasterProps=(I RasterProps)clipRaster;int dHeight 二 rasterProps.Height;当前栅格数据集的行数int dWidth二rasterProps.Width;当前栅格数据集的列数double dX=rasterProps.MeanCellSize().X;栅格的宽度double dY=rasterProps.MeanCellSize().Y;栅格的高度 lEnvelope extent=rasterProps.Extent;当前栅格数据集的范围 rstPixelType pixelType=rasterProps.PixelType;当前栅格像素类型 I Pnt pntSize=new PntClass();pntSize.SetCoords(dX,dY);I PixelBlock pixelBlock=clipRaster.CreatePixelBlock(pntSize);会在后 面说这个块，ESRI将栅格数据是按照块来存储的，有什么好处呢？I Pnt pnt=new PntClass();for(int i=0;i dHeight;i+)for(int j=0;j Arr=0*10*10 Arr0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 Arr 0 0 =1 Arr1,o,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,o,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1，0,0,0,0,0,0,0,0,0发现这样构造的数组用不成，他只不过是复制了 10行，借助强大的 google我找到了另外一种构造数组的方法：Arr=0 for column in range(ncols)for row in range(nrows)#calculate the valuesj=0;while j nrows:i=0;yl=ytop-(j+1)*(string.atof(cellsize)+0.5*(string.atof(cellsize)while incols:al=0 a2=0 xl=xleft+(i+1)*(string.atof(cellsize)-0.5*(string.atof(cellsize)n=0 while n PointCount:x=Xlistn y=Ylistn z=Zlistn n+=l;#利用平面距离的算法 d=idw_pdistance(xl,yl,x,yzstring.atof(Power)if dRadius:a2=a2+l/d al=al+z/dif a20:Arrji=al/a2i+=l;j+=l；当这个txt生成后，然后利用ArcGI S提供的工具rasterType=FLOAT”Execute ASCIIToRastergp.outputCoordinateSystem=SRinputPtsgp.ASCIIToRaster_conversion(outTxt,outRasterz rasterType)结果如下:I D W这个是第一步，还有其他的功能有待完善，稍后等我将这个完善了，继 续和大家分享