用python的方法进行gdal编程.pdf

分享心得，相互交流，共同进步Email：sishui198醉意人间在接触一个新的东西的时候，我喜欢从这新东西的一些名词和结构入手。记得以前一位老师说过，一个新的东西，肯定会有自己的专业术语，然后就是体系结构等，我们就从这个开始吧GDAL(Geospatial Data Abstraction Library)是一个在 X/MIT 许可协议下的开源栅格空间数据转换库。它利用抽象数据模型来表达所支持的各种文件格式。它还有一系列命令行工具来进行数据转换和处理。OGR 是 GDAL 项目的一个分支，功能与 GDAL 类似，只不过它提供对矢量数据的支持。有很多著名的GIS类产品都使用了GDAL/OGR库，包括ESRI的ArcGIS，Google Earth 和跨平台的 GRASS GIS 系统。利用 GDAL/OGR 库，可以使基于 Linux 的地理空间数据管理系统提供对矢量和栅格文件数据的支持。1.GDAL1.GDALGDAL（Geospatial Data Abstraction Library）提供对多种栅格数据的支持，包括 Arc/Info ASCII Grid(asc)，GeoTiff(tiff)，Erdas Imagine Images(img)，ASCII DEM(dem)等格式。GDAL 使用抽象数据模型(Abstract Data Model)来解析它所支持的数据格式，抽象数据模型包括数据集(Dataset)，坐标系统，仿射地理坐标转换(Affine Geo Transform)，大地控制点(GCPs)，元数据(Metadata)，栅格波段(Raster Band)，颜色表(Color Table)，子数据集域(Subdatasets Domain)，图像结构域(Image_Structure Domain)，XML 域(XML:Domains)。详细的结构描述请访问www.gdal.org/gdal_datamodel.html分享心得，相互交流，共同进步Email：sishui198醉意人间GDAL 的核心类结构设计如图所示：其中的类说明如下：GDALMajorObject 类：带有元数据的对象。GDALDdataset 类：通常是从一个栅格文件中提取的相关联的栅格波段集合和这些波段的元数据；GDALDdataset 也负责所有栅格波段的地理坐标转换(georeferencing transform)和坐标系定义。GDALDriver 类：文件格式驱动类，GDAL 会为每一个所支持的文件格式创建一个该类的实体，来管理该文件格式。GDALDriver 类：文件格式驱动类，GDAL 会为每一个所支持的文件格式创建一个该类的实体，来管理该文件格式。GDALDriverManager 类：文件格式驱动管理类，用来管理 GDALDriver 类。2.OGR2.OGROGR 提供对矢量数据格式的读写支持，它所支持的文件格式包括：ESRI Shapefiles，S-57，SDTS，PostGIS，Oracle Spatial，Mapinfo mid/mif，Mapinfo TAB。分享心得，相互交流，共同进步Email：sishui198醉意人间更详细的请访问 gdal.org/ogr/ogr_formats.html1)OGR 体系结构OGR 包括如下几部分：Geometry：类 Geometry(包括 OGRGeometry 等类)封装了 OpenGIS 的矢量数据模型，并提供了一些几何操作，WKB(Well Knows Binary)和 WKT(Well KnownText)格式之间的相互转换，以及空间参考系统(投影)。Spatial Reference：类OGRSpatialReference封装了投影和基准面的定义。Feature：类 OGRFeature 封装了一个完整 Feature 的定义，一个完整的Feature 包括一个 Geometry 和 Geometry 的一系列属性。分享心得，相互交流，共同进步Email：sishui198醉意人间Feature Definition：类 OGRFeatureDefn 里面封装了 feature 的属性，类型、名称及其默认的空间参考系统等。一个 OGRFeatureDefn 对象通常与一个层(layer)对应。Layer：类 OGRLayer 是一个抽象基类，表示数据源类 OGRDataSource 里面的一层要素(Feature)。Data Source：类 OGRDataSource 是一个抽象基类，表示含有 OGRLayer 对象的一个文件或一个数据库。Drivers：类 OGRSFDriver 对应于每一个所支持的矢量文件格式。类OGRSFDriver 由类 OGRSFDriverRegistrar 来注册和管理。Drivers：类 OGRSFDriver 对应于每一个所支持的矢量文件格式。类OGRSFDriver 由类 OGRSFDriverRegistrar 来注册和管理。2)OGR 的 Geometry 模型OGR 的 Geometry 模型是建立在 OpenGIS 的简单要素数据模型之上的。如下图所示：图-OGR 的 Geometry 模型关系图分享心得，相互交流，共同进步Email：sishui198醉意人间图-OpenGIS 的简单要素数据模型由上面两图的对比，可以清楚的看到，OGR 的 Geometry 模型是严格遵循OpenGIS 的简单要素数据规范的。OGR 的 Geometry 模型不仅在继承体系上与OpenGIS 的简单要素数据模型一致，在函数接口上也向其靠拢，从基本的获取Geometry 对象信息的方法如 Dimension()、GeometryType()、SRID()、Envelope()、AsText()、Boundary()等到判定空间未知关系的方法如Equals(anotherGeometry:Geometry)、Disjoint(anotherGeometry:Geometry)、Intersects(anotherGeometry:Geometry)、Touches(anotherGeometry:Geometry)等都是符合其标准的。分享心得，相互交流，共同进步Email：sishui198醉意人间 ESRI 的 ST_Geometry 数据存储模型GDAL/OGR 是非常著名的开源 GIS 库，GDAL 是对栅格数据进行操作，而 OGR 是对矢量数据进行操作，它们相当于一个通用数据访问库。甚至 ESRI 的产品中都用了此库。安装ArcGIS 10 Desktop 装了 Python2.6，那么我们就不用装这个了，然后去下面的地址找到 gdal 的 Python 包 GDAL-1.6.1.win32-py2.6.exe，自动安装到 python的安装目录 E:Python26ArcGIS10.0Libsite-packages分享心得，相互交流，共同进步Email：sishui198醉意人间下载下面的文件(GDAL Windows binaries)分享心得，相互交流，共同进步Email：sishui198醉意人间解压后按照 README_EXE.TXT 的说明添加环境变量，我的路径为 C:gdalwin32-1.6添加数据路径的变量为了验证是否安装成功，可以执行下列命令:from osgeo import gdalfrom osgeo.gdalconst import*分享心得，相互交流，共同进步Email：sishui198醉意人间dataset=gdal.Open(C:dataliuyu.jpg,GA_ReadOnly)dataset.GetDriver().ShortName对矢量数据的操作时通过 OGR，对每种数据有一个驱动，因此在操作数据的时候，我们先要获取相应数据格式的驱动：如 driver=ogr.GetDriverByName(ESRI Shapefile)，但是对于栅格数据略有不同；GDAL 数据驱动，与 OGR 数据驱动类似，需要先创建某一类型的数据驱动，再创建相应的栅格数据集。一次性注册所有的数据驱动，但是只能读不能写：gdal.AllRegister().单独注册某一类型的数据驱动，这样的话可以读也可以写，可以新建数据集：driver=gdal.GetDriverByName(HFA)driver.Register()driver.Register()导入库：from osgeo import ogr要读取某种类型的数据，必须要先载入数据驱动，也就是初始化一个对象，让它“知道”某种数据结构。from osgeo import ogrdriver=ogr.GetDriverByName(ESRI Shapefile)数据驱动 driver 的 open()方法返回一个数据源对象open(,)其中 update 为 0 是只读，为 1 是可写例如：from osgeo import ogrdriver=ogr.GetDriverByName(ESRI Shapefile)inputs=C:Export_Output.shp分享心得，相互交流，共同进步Email：sishui198醉意人间datasource=driver.Open(inputs,0)if datasource is None:print could not open读取数据层layer=datasource.GetLayer()一般 ESRI 的 shapefile 都是填 0 的，如果不填的话默认也是 0.#读取数据layer=datasource.GetLayer()n=layer.GetFeatureCount()n255feature=layer.GetNextFeature()#注意里面没有 FeatureClass 这个说法FID=feature.GetFID()FID0读出上下左右边界extent=layer.GetExtent()print extent:,extentprint ul:,extent0,extent3分享心得，相互交流，共同进步Email：sishui198醉意人间print lr:,extent1,extent2读取字段的个数fieldscount=feature.GetFieldCount()fieldscount23释放内存feature.Destroy()关闭数据源，相当于文件系统操作中的关闭文件datasource.Destroy()#创建文件datasource=driver.CreateDataSource(C:liuyu.shp)shplayer=datasource.CreateLayer(liuyu,geom_type=ogr.wkbPoint)shplayer=datasource.CreateLayer(liuyu,geom_type=ogr.wkbPoint)shplayer 相当于我们的 FeatureClass，因为我们将要素写入文件是通过shplayer 的，还记得 C#和 Engine 是如何操作数据的要添加一个新字段，只能在 layer 里面加，而且还不能有数据添加的字段如果是字符串，还要设定宽度fieldDefn=ogr.FieldDefn(id,ogr.OFTString)分享心得，相互交流，共同进步Email：sishui198醉意人间fieldDefn.SetWidth(4)shplayer.CreateField(fieldDefn)添加一个新的 feature，首先得完成上一步，把字段 field 都添加齐了然后从 shpfile 中读取相应的 feature 类型，并创建 featurefieldDefn=ogr.FieldDefn(id,ogr.OFTString)fieldDefn.SetWidth(4)shplayer.CreateField(fieldDefn)0featureDefn=shplayer.FeatureDefn()feature=ogr.Feature(featureDefn)将 feature 写入 shplayer建立点对象point=ogr.Geometry(ogr.wkbPoint)point.AddPoint(10,10)feature.SetGeometry(point)0point.GetX()10.0创建要素分享心得，相互交流，共同进步Email：sishui198醉意人间shpfile.CreateFeature(feature)0 n=shplayer.GetFeatureCount()n1field=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=new ESRI.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();分享心得，相互交流，共同进步Email：sishui198醉意人间/Apply the constructed shape to the new features shapefeature.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 hasrowSubtypes.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(InstalledBy),K Johnston);/Commit the updated values in the feature to the databasefeature.Store();如果用我们 ArcGIS 装了之后的那个 Python 环境敲写代码，一两句的话，还能忍受，但是如果大量的话，那可是要死人的，所以我们用 Eclipse 作为环境通过选择help（帮助）-software updates（软件更新）Find and InsatllSearch for new Feature to Install 下一中添加一个 New Remote Site为 ,然后就可以点击下一步自动升级了,这个是直接从网站上下载安装的。这个和 Eclipse 安装其他插件一样，或者将这个插件单独下载，然后将插件的文件复制到 Eclipse 的相关目录中当我们用这个 Eclipse 调试的时候，调试速度太慢了，看着半天没反应，真难受，看到下面的红字pydev debugger:warning:psyco not available for speedups(the debugger 分享心得，相互交流，共同进步Email：sishui198醉意人间will still work correctly,but a bit slower)pydev debugger:starting于是在网上找到 psyco（我就这个为调试加速器），下载一个适合自己平台的，安装就可以了，我没有装，我的代码不是很多最多就 200 行左右，我相信我的机器，哈哈下面这个数据是很简单的，前面六位（1-6）表示 X 坐标，而紧接着 X 的六位（7-12）表示 Y 坐标，我们就将这个文件利用 OGR 弄成 shp 数据，下面我们看下这个数据：分享心得，相互交流，共同进步Email：sishui198醉意人间import sysfrom osgeo import ogrogr.RegisterAll()driver=ogr.GetDriverByName(ESRI Shapefile)#C:txt2shp.shp this is wrongds=driver.CreateDataSource(C:liuyu.shp)layer=ds.CreateLayer(liuyu,geom_type=ogr.wkbPoint)#Create fieldsfield1=ogr.FieldDefn(Num,ogr.OFTInteger)field2=ogr.FieldDefn(Char,ogr.OFTString)field2.SetWidth(4)layer.CreateField(field1)layer.CreateField(field2)featuredefn=layer.GetLayerDefn()feature=ogr.Feature(featuredefn)#open txttxtfile=open(C:points.csp,r)line=txtfile.readline()while line:cleanline=line.replace(n,)point=ogr.Geometry(ogr.wkbPoint)分享心得，相互交流，共同进步Email：sishui198醉意人间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()下面是生成后的的数据表下面的是图形分享心得，相互交流，共同进步Email：sishui198醉意人间建立新的几何形状建立空的 geometry 对象：ogr.Geometry定义各种不同的 geometry 使用的方法是不一样的(point,line,polygon,etc)新建点，如我们上面的例子point=ogr.Geometry(ogr.wkbPoint)point.AddPoint(float(line0:6),float(line6:12)新建 line分享心得，相互交流，共同进步Email：sishui198醉意人间使用 AddPoint(,)添加点使用 SetPoint(,)更改点的坐标读取 0 号点的 x 坐标和 y 坐标print line.GetX(0)print line.GetY(0)新建面对于面的构建，如果用过 AE 的话，我们都知道，面是有一组 ring 构成册，而ring 又是由 segment 构成，segment 又是由 point 构成，我大学的时候老师说过这样一句话矢量数据的本质是 point，而栅格数据的本质是 pixel。我们看一下AE 中的这个图：那么OGR中的polygon呢？它的结构是什么样子？我们从下面的例子是否可以从例子中看到？#coding=gbkCreated on 2010-11-9author:刘宇分享心得，相互交流，共同进步Email：sishui198醉意人间from osgeo import ogrogr.RegisterAll()driver=ogr.GetDriverByName(ESRI Shapefile)#创建一个文件ds=driver.CreateDataSource(C:test.shp)shpfile=ds.CreateLayer(test,geom_type=ogr.wkbPolygon)#创建内环out=ogr.Geometry(ogr.wkbLinearRing)out.AddPoint(0,0)out.AddPoint(1000,0)out.AddPoint(1000,1000)out.AddPoint(0,1000)out.CloseRings()#创建外环inring=ogr.Geometry(ogr.wkbLinearRing)inring.AddPoint(250,250)inring.AddPoint(750,250)inring.AddPoint(750,750)inring.AddPoint(250,750)inring.AddPoint(250,250)polygon=ogr.Geometry(ogr.wkbPolygon)#面是由环构成的polygon.AddGeometry(out)polygon.AddGeometry(inring)featuredefn=shpfile.GetLayerDefn()feature=ogr.Feature(featuredefn)feature.SetGeometry(polygon)print polygon.GetArea()下面这句话可以帮我们数数数 polygon 能有几个 ringprint polygon.GetGeometryCount()从 polygon 中读取 ring，index 的顺序和创建 polygon 时添加 ring 的顺序相同outring=polygon.GetGeometryRef(0)分享心得，相互交流，共同进步Email：sishui198醉意人间inring=polygon.GetGeometryRef(1)创建复合几何形状 multi geometry例如 MultiPoint,MultiLineString,MultiPolygon用 AddGeometry 把普通的几何形状加到复合几何形状中，例如：multipoint=ogr.Geometry(ogr.wkbMultiPoint)point=ogr.Geometry(ogr.wkbPoint)point=ogr.Geometry(ogr.wkbPoint)point.AddPoint(10,10)multipoint.AddGeometry(point)point.AddPoint(20,20)multipoint.AddGeometry(point)读取 MultiGeometry 中的 Geometry，方法和从 Polygon 中读取 ring 是一样的，可以说 Polygon 是一种内置的 MultiGeometry。关于投影 Projections关于投影 Projections空间数据有一个东西非常重要，那就是投影，而 OGR 支持多种多样的 Projections，GDAL 支持 WKT,PROJ.4,ESPG,USGS,ESRI.prj，使用 SpatialReference 对象可以从 layer 和 Geometry 中读取 Projections，例如：spatialRef=layer.GetSpatialRef()spatialRef=geom.GetSpatialReference()投影信息一般存储在.prj 文件中，如果没有这个文件，上述函数返回 None分享心得，相互交流，共同进步Email：sishui198醉意人间建立一个新的 Projection：首先导入 osr 库，之后使用 osr.SpatialReference()创建 SpatialReference 对象之后用下列语句向 SpatialReference 对象导入投影信息导出 Projection，使用下面的语句可以导出为字符串我们对刚才的的例子加深下，因为我们上面的没有投影，现在我们给它指定一个投影#coding=gbkCreated on 2010-11-9author:刘宇from osgeo import ogr分享心得，相互交流，共同进步Email：sishui198醉意人间from osgeo import osrogr.RegisterAll()driver=ogr.GetDriverByName(ESRI Shapefile)#创建一个文件ds=driver.CreateDataSource(C:test.shp)shpfile=ds.CreateLayer(test,geom_type=ogr.wkbPolygon)#创建内环out=ogr.Geometry(ogr.wkbLinearRing)out.AddPoint(0,0)out.AddPoint(1000,0)out.AddPoint(1000,1000)out.AddPoint(0,1000)out.CloseRings()#创建外环inring=ogr.Geometry(ogr.wkbLinearRing)inring.AddPoint(250,250)inring.AddPoint(750,250)inring.AddPoint(750,750)inring.AddPoint(250,750)inring.AddPoint(250,250)polygon=ogr.Geometry(ogr.wkbPolygon)#面是由环构成的polygon.AddGeometry(out)polygon.AddGeometry(inring)featuredefn=shpfile.GetLayerDefn()feature=ogr.Feature(featuredefn)feature.SetGeometry(polygon)print polygon.GetArea()print polygon.GetGeometryCount()if(polygon.GetSpatialReference()=None):printOKsr=osr.SpatialReference()sr.ImportFromEPSG(32612)polygon.AssignSpatialReference(sr)/指定一个SpatialReference分享心得，相互交流，共同进步Email：sishui198醉意人间有什么问题？如果我的 geometry 很多，那我是不是要遍历每一个？想一下 ESRI是怎么搞的？没错 PRJ 文件#coding=gbkCreated on 2010-11-9author:刘宇from osgeo import ogrfrom osgeo import osrogr.RegisterAll()driver=ogr.GetDriverByName(ESRI Shapefile)#创建一个文件ds=driver.CreateDataSource(C:liuyu.shp)shpfile=ds.CreateLayer(liuyu,geom_type=ogr.wkbPolygon)#创建内环out=ogr.Geometry(ogr.wkbLinearRing)out.AddPoint(0,0)out.AddPoint(1000,0)out.AddPoint(1000,1000)out.AddPoint(0,1000)out.CloseRings()#创建外环inring=ogr.Geometry(ogr.wkbLinearRing)inring.AddPoint(250,250)inring.AddPoint(750,250)inring.AddPoint(750,750)inring.AddPoint(250,750)inring.AddPoint(250,250)polygon=ogr.Geometry(ogr.wkbPolygon)#面是由环构成的polygon.AddGeometry(out)polygon.AddGeometry(inring)featuredefn=shpfile.GetLayerDefn()feature=ogr.Feature(featuredefn)feature.SetGeometry(polygon)print polygon.GetArea()分享心得，相互交流，共同进步Email：sishui198醉意人间print polygon.GetGeometryCount()if(polygon.GetSpatialReference()=None):printOKsr=osr.SpatialReference()sr.ImportFromEPSG(32612)#如果将geometry的投影都改变了，那要一个个的循环，数据量太多的话，很麻烦，可以像ESRI一样提供一个文件，写进去prjfile=open(C:liuyu.prj,w)#转成字符串sr.MorphToESRI()prjfile.write(sr.ExportToWkt()prjfile.close()printOK上面的代码执行后，我们就可以看到下面的文件怎么做投影变换？用这个方法怎么做投影变换？用这个方法TransformTo()简单的地理数据处理 geoprocessing 可以对比 ArcGIS 提供的简单的地理数据处理 geoprocessing 可以对比 ArcGIS 提供的多边形的：交：poly3.Intersection(poly2)并：poly3.Union(poly2)差：poly3.Difference(poly2)补：poly3.SymmetricDifference(poly2)geometry 的：分享心得，相互交流，共同进步Email：sishui198醉意人间.Buffer()给geometry 加 buffer，就是把点线变成多边形，变粗了.Equal()两个 geometry 相等吗？.Distance()返回两个 geometry 之间的最短距离.GetEnvelope()信封，有意思，其实就是用一个方框框住这个几何形状，返回四个角的坐标(minx,maxx,miny,maxy)python 的函数 function，异常 exception 和模块 module 可以从任何一本 python教材上找到，在此不多赘述。栅格篇栅格篇GDAL 原生支持超过 100 种栅格数据类型，涵盖所有主流 GIS 与 RS 数据格式，如下图：完整的支持列表可以参考 www.gdal.org/formats_list.html我们知道对于栅格数据，从数学的角度来考虑的话就是一个矩阵，从我们在学习栅格数据和矢量数据的时候，我们知道栅格数据便于建立模型，这是因为他的结构的原因，我们看一下几个相关的模块。分享心得，相互交流，共同进步Email：sishui198醉意人间相关模块1.Numeric：高速的数组处理，对栅格数据尤其重要2.NumPy：下一代的 Numeric3.更强大的 gis 库 www.gispython.org/导入 GDAL 支持库from osgeo import gdal,gdalconstGDAL 数据驱动，与 OGR 数据驱动类似，需要先创建某一类型的数据驱动，再创建响应的栅格数据集。一次性注册所有的数据驱动，但是只能读不能写：gdal.AllRegister()单独注册某一类型的数据驱动，这样的话可以读也可以写，可以新建数据集：driver=gdal.GetDriverByName(GTiff)driver.Register()from osgeo import gdalfrom osgeo.gdalconst import*driver=gdal.GetDriverByName(JPEG)driver.Register()15ds=gdal.Open(C:dataliuyu.jpg,GA_ReadOnly)band=ds.GetRasterBand(1)/获取第一个波段 rown=band.XSize/没注意看这个到底是行还是列rown3888cown=band.YSizecown分享心得，相互交流，共同进步Email：sishui198醉意人间2592 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,.,61,62,64,81,82,85,.,67,66,66,87,87,87,.,72,68,65,dtype=uint8)Engine 中是怎么读取栅格数据的？IRasterProps rasterProps=(IRasterProps)clipRaster;int dHeight=rasterProps.Height;/当前栅格数据集的行数int dWidth=rasterProps.Width;/当前栅格数据集的列数double dX=rasterProps.MeanCellSize().X;/栅格的宽度double dY=rasterProps.MeanCellSize().Y;/栅格的高度IEnvelope extent=rasterProps.Extent;/当前栅格数据集的范围rstPixelType pixelType=rasterProps.PixelType;/当前栅格像素类型IPnt pntSize=new PntClass();pntSize.SetCoords(dX,dY);IPixelBlock pixelBlock=clipRaster.CreatePixelBlock(pntSize);/会在后面说这个块，ESRI 将栅格数据是按照块来存储的，有什么好处呢？IPnt pnt=new PntClass();for(int i=0;i dHeight;i+)for(int j=0;j dWidth;j+)pnt.SetCoords(i,j);clipRaster.Read(pnt,pixelBlock);if(pixelBlock!=null)object obj=pixelBlock.GetVal(0,0,0);分享心得，相互交流，共同进步Email：sishui198醉意人间MessageBox.Show(Convert.