视频1 视频21 视频41 视频61 视频文章1 视频文章21 视频文章41 视频文章61 推荐1 推荐3 推荐5 推荐7 推荐9 推荐11 推荐13 推荐15 推荐17 推荐19 推荐21 推荐23 推荐25 推荐27 推荐29 推荐31 推荐33 推荐35 推荐37 推荐39 推荐41 推荐43 推荐45 推荐47 推荐49 关键词1 关键词101 关键词201 关键词301 关键词401 关键词501 关键词601 关键词701 关键词801 关键词901 关键词1001 关键词1101 关键词1201 关键词1301 关键词1401 关键词1501 关键词1601 关键词1701 关键词1801 关键词1901 视频扩展1 视频扩展6 视频扩展11 视频扩展16 文章1 文章201 文章401 文章601 文章801 文章1001 资讯1 资讯501 资讯1001 资讯1501 标签1 标签501 标签1001 关键词1 关键词501 关键词1001 关键词1501 专题2001
pythongdal教程之:几何形状geometry与投影projection
2020-11-27 14:27:10 责编:小采
文档


建立空的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死活不行,哎。。

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

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

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

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

targetSR.MorphToESRI()

file = open('test.prj', 'w')

file.write(targetSR.ExportToWkt())

ffile.close()

下载本文
显示全文
专题