GDAL python教程（2）——几何形状geometry与投影projection
==================================================================    


建立新的几何形状
----------------------------------

建立空的geometry对象：ogr.Geometry

定义各种不同的geometry使用的方法是不一样的(point, line, polygon, etc)

 
新建点point，使用方法AddPoint( <x>, <y>, [<z>])。其中的z坐标一般是省略的，默认值是0

例如：

::

    point = ogr.Geometry(ogr.wkbPoint)
    point.AddPoint(10,20)

 

新建line
~~~~~~~~~~~~~~~~~~~~~~~~~

使用AddPoint(<x>, <y>, [<z>])添加点

使用SetPoint(<index>, <x>, <y>, [<z>])更改点的坐标

例如下面这段代码，更改了0号点的坐标：

::

    line = ogr.Geometry(ogr.wkbLineString)
    line.AddPoint(10,10)
    line.AddPoint(20,20)
    line.SetPoint(0,30,30)   #(10,10) -> (30,30)

 

统计所有点的数目

::

    print line.GetPointCount()

读取0号点的x坐标和y坐标

::

    print line.GetX(0)
    print line.GetY(0)

 

新建多边形，首先要新建环(ring)，然后把环添加到多边形对象中。

 

如何创建一个ring？先新建一个ring对象，然后向里面逐个添加点。

::

    ring = ogr.Geometry(ogr.wkbLinearRing)
    ring.AddPoint(0,0)
    ring.AddPoint(100,0)
    ring.AddPoint(100,100)
    ring.AddPoint(0,100)

结束的时候，用CloseRings关闭ring，或者将最后一个点的坐标设定为与第一个点相同。

::

    ring.CloseRings()
    ring.AddPoint(0,0)

 

下面举一个例子，创建一个方框。这是个polygon对象，又例外两层ring构成。

::

    outring = ogr.Geometry(ogr.wkbLinearRing)
    outring.AddPoint(0,0)
    outring.AddPoint(100,0)
    outring.AddPoint(100,100)
    outring.AddPoint(0,100)
    outring.AddPoint(0,0)

    inring = ogr.Geometry(ogr.wkbLinearRing)inring = ogr.Geometry(ogr.wkbLinearRing)
    inring.AddPoint(25,25)
    inring.AddPoint(75,25)
    inring.AddPoint(75,75)
    inring.AddPoint(25,75)
    inring.CloseRings()

    polygon = ogr.Geometry(ogr.wkbPolygon)
    polygon.AddGeometry(outring)
    polygon.AddGeometry(inring)

最后三句话比较重要，就是先建立一个polygon对象，然后添加外层ring和内层ring

下面这句话可以帮你数数你的polygon能有几个ring

::

    print polygon.GetGeometryCount()

从polygon中读取ring，index的顺序和创建polygon时添加ring的顺序相同

::

    outring = polygon.GetGeometryRef(0)
    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。

 

不要删除一个已存在的Feature的Geometry，会把python搞崩溃的

只能删除脚本运行期间创建的Geometry，比方说手工创建出来的，或者调用其他函数自动创建的。就算这个Geometry已经用来创建别的Feature，你还是可以删除它。

例如：Polygon.Destroy()

 

关于投影Projections，使用SpatialReference对象

多种多样的Projections，GDAL支持WKT, PROJ.4, ESPG, USGS, ESRI.prj

可以从layer和Geometry中读取Projections，例如：

::

    spatialRef = layer.GetSpatialRef()
    spatialRef = geom.GetSpatialReference()

投影信息一般存储在.prj文件中，如果没有这个文件，上述函数返回None

 

建立一个新的Projection：

首先导入osr库，之后使用osr.SpatialReference()创建SpatialReference对象

之后用下列语句向SpatialReference对象导入投影信息

* ImportFromWkt(<wkt>)
* ImportFromEPSG(<epsg>)
* ImportFromProj4(<proj4>)
* ImportFromESRI(<proj_lines>)
* ImportFromPCI(<proj>, <units>, <parms>)
* ImportFromUSGS(<proj_code>, <zone>)
* ImportFromXML(<xml>)

导出Projection，使用下面的语句可以导出为字符串

* ExportToWkt()
* ExportToPrettyWkt()
* ExportToProj4()
* ExportToPCI()
* ExportToUSGS()
* ExportToXML()


对一个几何形状Geometry进行投影变换，要先初始化两个Projection，然后创建一个CoordinateTransformation对象，用它来做变换

::

    sourceSR = osr.SpatialReference()
    sourceSR.ImportFromEPSG(32612) #UTM 12N WGS84
    targetSR = osr.SpatialReference()
    targetSR.ImportFromEPSG(4326) #Geo WGS84
    coordTrans = osr.CoordinateTransformation(sourceSR, targetSR)
    geom.Transform(coordTrans)

但是这段代码很让人蛋疼！在windows里面跑不通。老外的论坛里面有讨论，说在linux里面没问题，windows死活不行，哎。。。

http://n2.nabble.com/PROJ-4-EPSG-28992-td2033665.html

 

另外还有几个要注意的地方：

要在适当的时候编辑Geometry，投影变换之后最好就不要再动了吧。

对一个数据源DataSource里面的所有Geometry做投影变换，你得一个一个来。用个循环吧。

 

将你的投影写入.prj文件，其实很简单。首先MorphToESRI()，转成字符串，然后开个文本文件往里面写就行了。例如：

::

    targetSR.MorphToESRI()
    file = open('test.prj', 'w')
    file.write(targetSR.ExportToWkt())
    ffile.close()