>>> from env_helper import info; info()
页面更新时间： 2023-06-23 11:55:14
运行环境：
Linux发行版本: Debian GNU/Linux 12 (bookworm)
操作系统内核: Linux-6.1.0-9-amd64-x86_64-with-glibc2.36
Python版本: 3.11.2
4.3. 使用 OGR 创建矢量数据¶
我们以Shapefile为例，使用以下函数创建新的数据源（DataSource）：
driver.CreateDataSource(<filename>)
需注意的是，应先判断这个要被创建的文件（filename）不能存在，否则会出错。 通常需要先判断一下文件是否已存在，若存在的话，需使用其他文件名，或将已存在的矢量数据直接删除。
4.3.1. 删除矢量数据¶
大多数的GIS数据格式，如shapefile，mapinfo tab等，都不是单一文件。 像
Shapefile，除了最基本的shp文件外，还有保存属性的dbf文件。
因此，在删除GIS数据时，可使用 os
模块删除，不过需要查找相关的文件进行全部删除，较之OGR的删除数据函数
DeleteDataSource()要烦琐得多。
>>> import os
>>> from osgeo import ogr
>>> driver=ogr.GetDriverByName('ESRI Shapefile')
>>> out_shp = '/tmp/xx_world_borders.shp'
>>> if os.path.exists(out_shp):
>>> driver.DeleteDataSource(out_shp)
>>> dir(out_shp)[:6] + ['... ...'] + dir(out_shp)[-3:]
['__add__',
'__class__',
'__contains__',
'__delattr__',
'__dir__',
'__doc__',
'... ...',
'translate',
'upper',
'zfill']
然后，使用下面的语句创建数据：
>>> ds =driver.CreateDataSource(out_shp)
>>> layer=ds.CreateLayer('test',geom_type=ogr.wkbPoint)
若想添加一个新字段，只能在layer里加，而且还得保证layer里没有数据。
若添加的字段是字符串，还需设定宽度。
>>> fieldDefn = ogr.FieldDefn( 'id',ogr.OFTString)
>>> fieldDefn.SetWidth(4)
>>> layer.CreateField(fieldDefn)
0
若添加一个新的feature，首先得完成上一步，把字段field都添加完整。
然后从layer中读取相应的feature类型，并创建feature：
>>> featureDefn = layer.GetLayerDefn()
>>> feature = ogr.Feature(featureDefn)
设定几何形状：
>>> point = ogr.Geometry(ogr.wkbPoint)
>>> point.SetPoint(0,123,123)
>>> feature.SetGeometry(point)
0
设定某字段数值：
>>> feature.SetField('id', 23)
将feature写入layer：
>>> layer.CreateFeature(feature)
0
最后，从内存中清除 ds，将数据写入磁盘中。
>>> ds.Destroy()
4.3.2. 使用OGR创建要素几何形状（Geometry）¶
创建空的geometry对象：ogr.Geometry。定义不同的geometry方法也是不同的(point,
line, polygon,
etc)。ogr中提供了各种类型，常用的有：ogr.wkbPoint，
ogr.wkbLineString，ogr.wkbPolygon。在创建这三种不同要素时，需注意使用的格式。
OGR中定义的几何类型
OGR关键字 | 名称 | 编码值 |
|---|---|---|
wkbUnknown | 未知类型 | 0 |
wkbPoint | 点 | 1 |
wkbLineString | 线 | 2 |
wkbPolygon | 多边形 | 3 |
wkbMultiPoint | 多点 | 4 |
wkbMultiLineString | 多线 | 5 |
wkbMultiPolygon | 多部分多边形 | 6 |
wkbGeometryCollection | 几何形状集合 | 7 |
wkbNone | 无 | 100 |
wkbLinearRing | 线环 | 101 |
wkbPoint25D | ２.5维点 | 0x80000001 |
wkbLineString25D | 2.5维线 | 0x80000002 |
wkbPolygon25D | 2.5维多边形 | 0x80000003 |
wkbMultiPoint25D | 2.5维多点 | 0x80000004 |
wkbMultiLineString25D | 2.5维多线 | 0x80000005 |
wkbMultiPolygon25D | 2.5维多部分多边形 | 0x80000006 |
wkbGeometryCollection25D | 2.5维几何形状集合 | 0x80000007 |
4.3.3. 创建点要素¶
下面我们先来看一下如何创建一个点。注意，此时只是在内存中将对象创建出来， 跟文件没有关系。
新建点point，使用方法AddPoint( <x>, <y>, [<z>])。其中的z坐标一般是省略的，默认值是0
例如：
>>> from osgeo import ogr
>>> point = ogr.Geometry(ogr.wkbPoint)
>>> point.AddPoint(10,20)
>>> print(point)
POINT Z (10 20 0)
注意，若向point中添加多个点，并不会出错，
但结果只会是最后添加的点。
4.3.4. 创建线要素¶
新建Line方法与点基本一致。与点不同的是，线需要添加多个点。
使用AddPoint(<x>, <y>, [<z>])添加点
使用SetPoint(<index>, <x>, <y>, [<z>])更改点的坐标
例如下面这段代码，更改了0号点的坐标：
>>> from osgeo import ogr
>>> line = ogr.Geometry(ogr.wkbLineString)
>>> line.AddPoint(10,10)
>>> line.AddPoint(20,20)
>>> print (line)
LINESTRING Z (10 10 0,20 20 0)
使用线对象的SetPoint修改点坐标，例如：
>>> print (line)
LINESTRING Z (10 10 0,20 20 0)
另外，还有一些其他有用的函数，如统计所有点数量：
>>> print (line.GetPointCount())
2
又如读取0号点的x坐标和y坐标：
>>> print (line.GetX(0))
10.0
>>> print (line.GetY(0))
10.0
4.3.5. 创建线环（Ring）¶
创建线环（Ring）与创建线(Line)类似，先创建一个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。
>>> ring.CloseRings()
>>> print(ring)
LINEARRING Z (0 0 0,100 0 0,100 100 0,0 100 0,0 0 0)
可以看出，最后点坐标与起始点坐标是一致的。
4.3.6. 创建多边形（Polygon）¶
创建Polygon与创建Point与Line的区别较大。创建多边形，首先要创建ring，然后把ring添加到多边形中。下面是创建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)
这里将最后点坐标设置与起始点相同，来闭合ring：
>>> 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)
>>>
0
>>> polygon.AddGeometry(inring)
0
最后三行命令比较重要：先建立一个polygon对象，然后再添加外层ring和内层ring。
若查看polygon有几个ring，可运行下列代码：
>>> print (polygon.GetGeometryCount())
2
对于polygon对象，无法直接输出坐标，而是先获取几何对象，也就是ring。从polygon中读取ring、index的顺序与创建polygon的顺序相同。
>>> outring2 = polygon.GetGeometryRef(0)
>>> inring2 = polygon.GetGeometryRef(1)
>>> print(outring2)
LINEARRING Z (0 0 0,100 0 0,100 100 0,0 100 0,0 0 0)
>>> print(inring2)
LINEARRING Z (25 25 0,75 25 0,75 75 0,25 75 0,25 25 0)
4.3.7. 创建复合几何形状（MultiGeometry）¶
运用AddGeometry方法将普通的几何形状如MultiPoint,
MultiLineString, MultiPolygon添至复合几何形状中：
>>> multipoint = ogr.Geometry(ogr.wkbMultiPoint)
>>> point = ogr.Geometry(ogr.wkbPoint)
>>> point.AddPoint(10,10)
>>> multipoint.AddGeometry(point)
0
>>> point.AddPoint(20,20)
>>> multipoint.AddGeometry(point)
0
读取MultiGeometry的Geometry，其方法与polygon读取ring是一样的，因此Polygon是一种内置的MultiGeometry。
注意：不要删除一个已存在Feature geometry，可能会导致python崩溃，可以删除脚本运行期间创建的Geometry，可能是手动创建，或者调用其他函数自动创建的。就算这个Geometry已经创建了其它的Feature，你还是可以删除它：
例如：Polygon.Destroy()
4.3.8. 使用OGR创建数据集的几何形状¶
对于创建Geometry来讲，wkt是比较直观的。wkt是最简单的字符串格式，而且Python提供了大量的函数来处理字符串。使用wkt创建矢量数据集与拷贝创建数据集的基本步骤是一致的。首先是创建Driver，按照顺序依次创建：dataSource，Layer，Feature。需要分别创建点、线、与多边形。使用extent变量来存储多边形的四个角点，点与线坐标根据extent变量来生成，如：
extent = [400, 1100, 300, 600]
创建点数据集¶
一般情况不会生成单独的点，需要使用For循环语句来对数组参数进行遍历。
创建线数据集¶
注意坐标的格式：不同点坐标之间用“,”隔断，坐标的不同维度用空格隔断。
创建多边形数据集¶
首先使用wkt来定义一个矩形。矩形的范围由extent变量定义，矩形四个角点坐标根据extend生成，存储到wkt变量中。通过CreateGeometryFromWkt()创建一个矩形形状，并放入Feature中，最后调用Destroy()，将数据写入磁盘，这样就构成了一个矩形的数据集。
4.3.9. 使用OGR拷贝方法创建新的Shapefile¶
使用OGR拷贝创建新的矢量数据集的方法其实就是将数据读取的顺序再重复一遍。总体来说这个过程就是构建数据源->构建层->构建要素->构建几何形状->关闭数据源。
在OGR的矢量数据模型中，矢量数据分为datasource，layer，feature三个层次，在这三个层次上，OGR均有不同的复制数据的方法。
在datasource层次创建数据¶
以世界各国边界的数据作示例，运用下列语句来创建与之内容完全相同的另一套数据：
>>> from osgeo import ogr
>>> import os,math
>>> inshp = '/gdata/GSHHS_c.shp'
>>> ds = ogr.Open(inshp)
>>> driver = ogr.GetDriverByName("ESRI Shapefile")
>>> outputfile = '/tmp/xx_world_borders_copy.shp'
>>> if os.access( outputfile, os.F_OK ):
>>> driver.DeleteDataSource(outputfile)
>>> pt_cp = driver.CopyDataSource(ds, outputfile)
>>> pt_cp.Release()
>>> dir(pt_cp)[:6] + ['... ...'] + dir(pt_cp)[-3:]
Warning 1: Value 50654050.6944999993 of field area of feature 0 not successfully written. Possibly due to too larger number with respect to field width
Warning 1: Value 50654050.6944999993 of field area of feature 1 not successfully written. Possibly due to too larger number with respect to field width
Warning 1: Value 29220969.7270000018 of field area of feature 2 not successfully written. Possibly due to too larger number with respect to field width
Warning 1: Value 20154740.0899999999 of field area of feature 3 not successfully written. Possibly due to too larger number with respect to field width
Warning 1: Value 17534164.0667999983 of field area of feature 4 not successfully written. Possibly due to too larger number with respect to field width
Warning 1: Value 13919140.5654000007 of field area of feature 1736 not successfully written. Possibly due to too larger number with respect to field width
['AbortSQL',
'CommitTransaction',
'CopyLayer',
'CreateLayer',
'DeleteLayer',
'Dereference',
'... ...',
'name',
'this',
'thisown']
在这里，我们使用了CopyDataSource()的方法来对数据进行复制。这个函数的第1个参数是要进行复制的数据源，第2个参数是要生成的数据路径，并且返回一个指针。为了将数据写入磁盘，必须使用Release()方法释放此数据。
在layer层次拷贝数据¶
下面我们在layer层次来完成同样的工作。
>>> from osgeo import ogr
>>> import os,math
>>> inshp = '/gdata/GSHHS_c.shp'
>>> ds = ogr.Open(inshp)
>>> driver = ogr.GetDriverByName("ESRI Shapefile")
>>> outputfile ='/tmp/xx_world_borders_copy2.shp'
>>> if os.access( outputfile, os.F_OK ):
>>> driver.DeleteDataSource( outputfile )
>>> layer = ds.GetLayer()
>>> newds = driver.CreateDataSource(outputfile)
>>> pt_layer = newds.CopyLayer(layer, 'abcd')
>>> newds.Destroy()
>>> dir(pt_layer)[:6] + ['... ...'] + dir(pt_layer)[-3:]
Warning 1: Value 50654050.6944999993 of field area of feature 0 not successfully written. Possibly due to too larger number with respect to field width
Warning 1: Value 50654050.6944999993 of field area of feature 1 not successfully written. Possibly due to too larger number with respect to field width
Warning 1: Value 29220969.7270000018 of field area of feature 2 not successfully written. Possibly due to too larger number with respect to field width
Warning 1: Value 20154740.0899999999 of field area of feature 3 not successfully written. Possibly due to too larger number with respect to field width
Warning 1: Value 17534164.0667999983 of field area of feature 4 not successfully written. Possibly due to too larger number with respect to field width
Warning 1: Value 13919140.5654000007 of field area of feature 1736 not successfully written. Possibly due to too larger number with respect to field width
['AlterFieldDefn',
'AlterGeomFieldDefn',
'Clip',
'CommitTransaction',
'CreateFeature',
'CreateField',
'... ...',
'schema',
'this',
'thisown']
若想对layer进行拷贝，首先得有数据源。使用newds的CopyLayer()方法对图层进行拷贝，以创建newds数据源。值得注意的是，当图层拷贝完之后，需要对newds进行Destroy()操作，将数据写入磁盘。
CopyLayer()的第1个参数是OGR的Layer对象，第2个参数是要生成图层的名称。对于Shapefile来说，这个名称是没有用的，但必须给这个字符串赋变量值。
在feature层次拷贝数据¶
下列语句是在feature层次拷贝数据的方法，其实也只是多了一步操作：
>>> from osgeo import ogr
>>> import os,math
>>> inshp = '/gdata/GSHHS_c.shp'
>>> ds = ogr.Open(inshp)
>>> driver = ogr.GetDriverByName("ESRI Shapefile")
>>> outputfile ='/tmp/xx_world_borders_copy3.shp'
>>> if os.access( outputfile, os.F_OK ):
>>> driver.DeleteDataSource( outputfile )
>>> newds = driver.CreateDataSource(outputfile)
>>> layernew = newds.CreateLayer('worldcopy',None,ogr.wkbLineString)
>>> layer = ds.GetLayer()
>>> extent = layer.GetExtent()
>>> extent
(-180.0, 180.00000000000023, -89.99999999999994, 83.53036100000006)
>>> from osgeo import ogr
>>> import os,math
>>> inshp = '/gdata/GSHHS_c.shp'
>>> ds = ogr.Open(inshp)
>>> driver = ogr.GetDriverByName("ESRI Shapefile")
>>> outputfile ='/tmp/xx_world_borders_copy3.shp'
>>> if os.access( outputfile, os.F_OK ):
>>> driver.DeleteDataSource( outputfile )
>>> newds = driver.CreateDataSource(outputfile)
>>> layernew = newds.CreateLayer('worldcopy',None,ogr.wkbLineString)
>>> layer = ds.GetLayer()
>>> feature = layer.GetNextFeature()
>>> while feature is not None:
>>> layernew.CreateFeature(feature)
>>> feature = layer.GetNextFeature()
>>> newds.Destroy()
>>> dir(feature)[:6] + ['... ...'] + dir(feature)[-3:]
['__bool__',
'__class__',
'__delattr__',
'__dir__',
'__doc__',
'__eq__',
'... ...',
'__sizeof__',
'__str__',
'__subclasshook__']
特别注意：Destroy()不能省略，Destroy()除了销毁数据还有数据flush到磁盘的作用。如果没有此命令，那么刚才创建的一系列要素就不会被写入磁盘的文件中。
4.3.10. OGR属性字段的定义与使用¶
GIS数据除了图形要素之外，更重要的就是属性数据。其实GIS数据的属性可以理解成关系型数据库。只不过在这个关系数据库中，有一列特定的索引值，且数据库中的记录与图形对象也是一一对应的。
在OGR中定义属性字段与赋值¶
OGR定义字段与数据库类型时，首先需查看有哪些字段类型（表[tab:ogrfieldtype]）。
表3.3 OGR中定义的字段类型
枚举值 | 类型 | 位数 | 编码值 |
|---|---|---|---|
OFTInteger | 整型 | 32位整型 | 0 |
OFTIntegerList | 整型列表 | 32位整型 | 1 |
OFTReal | 实数型 | 双精度浮点型 | 2 |
OFTRealList | 实数型列表 | 双精度浮点型 | 3 |
OFTString | 字符串 | ASCII字符 | 4 |
OFTStringList | 字符串列表 | ASCII字符串数组 | 5 |
OFTWideString | 宽字符串 | 不再使用 | 6 |
OFTWideStringList | 宽字符串列表 | 不再使用 | 7 |
OFTBinary | 二进制型 | 原始的二进制数据 | 8 |
OFTDate | 日期型 | 9 | |
OFTTime | 时间型 | 10 | |
OFTDateTime | 日期时间型 | 11 |
我们还可以给这个矩形添加属性数据，首先定义表头，也就是定义字段，需强调的是定义一个字段需要定义其字段描述（包括字段的名称与类型），然后将这个字段描述创建到layer中。有了表头，就可以使用SetField方法在输入Feature时添加属性表内容，这时打开QGIS
ShapeFile时就可以看到结果了。若使用中文的字段值，请注意字符串字段的宽度，若超出宽度，显示内容就会出现问题。
定义字段的位置¶
当OGRFeatureDefn中无OGRFeature对象时，此方法才可以调用。此方法的优点是保证在不影响调用者响应的同时，输入的OGRFieldDefn也会被拷贝备份。
>>> import os
>>> from osgeo import ogr
>>> driver = ogr.GetDriverByName("ESRI Shapefile")
>>> extfile = '/tmp/rect3.shp'
>>> if os.access( extfile, os.F_OK ):
>>> driver.DeleteDataSource( extfile )
>>> newds = driver.CreateDataSource(extfile)
>>> layernew = newds.CreateLayer('rect3',None,ogr.wkbPolygon)
>>> fieldcnstr = ogr.FieldDefn("关于",ogr.OFTString)
>>> fieldcnstr.SetWidth(36)
>>> fieldf = ogr.FieldDefn("f",ogr.OFTReal)
>>> laydef = layernew.GetLayerDefn()
>>> laydef.AddFieldDefn(fieldcnstr)
>>> laydef.AddFieldDefn(fieldf)
>>> laydef.GetFieldCount()
2
>>> newds.Destroy()
可以看出，输出的图层属性表中并无内容。可能是这个函数只在内存中虚拟地创建了一个字段，与磁盘并无关系，或者说这个函数只是用来创建一个属性字段的框架，仅仅用来定义，在建立新的数据时可能起到参考与引用的作用。在软件处理过程中，由于软件的问题可能会导致图形对象与属性记录并非一一对应，严重时矢量数据可能会损坏。
创建MapInfo文件时可能会遇到的问题¶
以TAB格式创建的文件，在写入第一个要素前，需要在相应的驱动上为每个投影设置合理的范围，不过目前尚未有运用OGRDataSource为新文件设置缺省框架的相应机制，因此在创建新的Layer时，先暂时在mapinfo驱动上用下面的默认语句来设置框架范围：
一个经纬坐标系的框架范围BOUNDS(-180,-90)(180,90)
其他的用BOUNDS(-30000000,-15000000)(30000000,15000000)
如果在创建Layer时没有提供坐标系，可以使用以上的投影实例。但如果是地理坐标，那么此例子提供的也只是低精度的地理坐标系。你可以添加-a\_srs wgs84到ogr2ogr中，强制地转换成地理模型。
Mapinfo有很多属性限制：
只可以创建整数（Integer），实数（Real），字符串（String）字段，其他的类型均不能创建。
通常用字符串字段的宽来设置dat文件中的存储大小，因此比字段宽数长的字符串就会被限制。
若字符串字段没有分配宽时，其默认值是254。
但是BOUNDS (-30000000,-15000000)
(30000000,15000000)这个框架的限制其实并不合理。相对于北京坐标系而言，中国西安的投影坐标系得东移500公里，再加之带号，会有超过30000000的情况，又如福建在3度分带，加之北京坐标系，加之带号差不多在39-40位，加上东移的500000m，最后的坐标值变成39XXXXXX，超出了(-30000000,30000000)的范围。那么用OGR创建的Mapinfo的福建部分就是错的。但是mif却可以正常创建，其缺点就是其ACSII字符文件比较大。
4.3.11. 添加投影（Projection）¶
上述CreateLayer的第二个参数是None，其实它是一个空间参考（Spatial Reference），实则是给这些数字附有实际含义。添加了空间参考，才能依据这些数字在地图上定位。同一个Layer需要有相同的空间参考。
获取投影¶
投影使用的是Spatial Reference对象。目前Projections较为多样化，GDAL支持WKT, PROJ.4, ESPG, USGS, ESRI.prj等。我们可以从layer和Geometry中读取Projections：
>>> from osgeo import ogr
>>> ds = ogr.Open('/gdata/GSHHS_c.shp')
>>> layer = ds.GetLayer(0)
>>> spatialRef = layer.GetSpatialRef()
>>> print(spatialRef)
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AXIS["Latitude",NORTH],
AXIS["Longitude",EAST],
AUTHORITY["EPSG","4326"]]
Shapefile的投影信息一般存储在.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()
根据已有的投影创建新投影¶
创建空间参考最简单的办法是用ImportFromWkt导入Wkt中。Wkt获得最简便的方式是安装一个PostgreSQL，再安装PostGIS。在Windows系统中，有安装PostGIS的选项，而Linux系统需要自己编译。PostGIS的spatial_reference表中包含全套的空间参考，几乎所有的投影坐标信息的Wkt表达全在里面。按需拷贝，开箱即用。
假设已有Wkt，创建一个空间参考：
>>> from osgeo import osr
>>> wkt = spatialRef.ExportToWkt()
>>> spatial = osr.SpatialReference()
>>> spatial.ImportFromWkt(wkt)
0
最后用spatial替代CreateLayer的第二个参数None，完成建立空间参考的矢量数据投影。
对矢量数据进行投影转换¶
可对矢量数据要素的几何特征进行投影转换，首先先将两个Projection初始化，再创建一个CoordinateTransformation对象：
>>> targetSR = osr.SpatialReference()
>>> targetSR.ImportFromEPSG(4326) #Geo WGS84
>>> coordTrans = osr.CoordinateTransformation(spatial, targetSR)
>>> feature = layer.GetFeature(0)
>>> geom = feature.GetGeometryRef()
>>> geom.ExportToWkt()
>>> geom.Transform(coordTrans)
>>> geom.ExportToWkt()[:120] + ' ... ...'
'POLYGON ((107.544833 76.914028,106.533222 76.499556,111.093417 76.7653890000001,113.934917 75.835,113.583278 75.53119400 ... ...'
另外还需注意：要在适当的时候编辑Geometry，投影变换后最好不要再修改Geometry。一个数据源（DataSource）里面的所有Geometry都得作投影变换，可以用循环语句。
将投影写入文件¶
若单独处理投影信息，将投影信息单独写入文件时，需将你的投影写入.prj文件，用MorphToESRI()转成字符串，然后新建文本文件，写入文件即可：
>>> targetSR.MorphToESRI()
>>> file = open('/tmp/test.prj', 'w')
>>> file.write(targetSR.ExportToWkt())
>>> file.close()