OGRCoordinateTransformation 类可用在不同坐标系统间转换坐标。
新的转换对象由OGRCoordinateTransformation()创建。 创建后
OGRCoordinateTransformation::Transform()
方法就可以用来转换不同坐标系的点。
导致在转换时会出错的两个原因：
第一，OGRCreateCoordinateTransformation() 可能失败。
一般是因为程序认为两个指定坐标系统不能建立转换关系，
这可能是由于使用的投影Proj内部暂时不支持，
这样两个不一样的基准面就没有已知关系。 或者一个投影定义得不适当。
如果OGRCreateCoordinateTransformation() 失败了将返回空。
OGRCoordinateTransformation::Transform() 方法本身也可能失败。
这可能是因为一个有问题积累后形成的问题。
或者是“因为一个或者多个点因数值上未定义而省略操作”起因的后果。
Transform()函数如果成功，则返回为TRUE，
如果某个点转换失败，剩下的点就会以一个错误的状态保持。
坐标转换服务其实可以处理3D点的。 并且服务会根据椭球体和基准面上的高程差异调节高程。 将来，也有可能会在不同矢量基准面上进行的点间的转换的应用。 如果没有Z值，则假设点都在大地水准面上。
下面的例子显示了如何方便创建一个地理坐标系统 到一个基于同一个地理坐标系统的投影坐标系统。并在两个坐标系统之间转换。
In [2]:
from osgeo import gdal
from osgeo import osr
# import osr
dataset = gdal.Open("/data/gdata/geotiff_file.tif")
/opt/conda/lib/python3.12/site-packages/osgeo/gdal.py:311: FutureWarning: Neither gdal.UseExceptions() nor gdal.DontUseExceptions() has been explicitly called. In GDAL 4.0, exceptions will be enabled by default. warnings.warn(
从数据集中获取空间参考并且建立一个SpatialReference对象：
In [3]:
sr = dataset.GetProjectionRef()
osrobj = osr.SpatialReference()
osrobj.ImportFromWkt(sr)
Out[3]:
0
In [4]:
osrobj.ExportToWkt()
Out[4]:
'PROJCS["Albers_Beijing54",GEOGCS["Unknown datum based upon the Krassowsky 1940 ellipsoid",DATUM["Not_specified_based_on_Krassowsky_1940_ellipsoid",SPHEROID["Krassowsky 1940",6378245,298.3,AUTHORITY["EPSG","7024"]],AUTHORITY["EPSG","6024"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["latitude_of_center",0],PARAMETER["longitude_of_center",105],PARAMETER["standard_parallel_1",25],PARAMETER["standard_parallel_2",47],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]'
In [5]:
osrobj.MorphToESRI() # 重点，转成Esri的wkt格式
Out[5]:
0
In [6]:
osrobj.ExportToWkt()
Out[6]:
'PROJCS["Albers_Beijing54",GEOGCS["GCS_Krasovsky_1940",DATUM["D_Krasovsky_1940",SPHEROID["Krasovsky_1940",6378245.0,298.3]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Albers"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",105.0],PARAMETER["Standard_Parallel_1",25.0],PARAMETER["Standard_Parallel_2",47.0],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]'
In [7]:
osrobj.IsGeographic()
Out[7]:
0
In [8]:
osrobj.IsProjected()
Out[8]:
1
再获取一个进行比较：
In [9]:
dataset2 = gdal.Open("/data/gdata/lu75c.tif")
sr2 = dataset2.GetProjectionRef()
osrobj2 = osr.SpatialReference()
osrobj2.ImportFromWkt(sr2)
osrobj2.IsSame(osrobj)
Out[9]:
0
创建一个地理坐标系，然后和已有的坐标系统进行比较：
In [10]:
osrobj3 = osr.SpatialReference()
osrobj3.SetWellKnownGeogCS("WGS84")
osrobj3.IsSame(osrobj2)
osrobj3.IsSame(osrobj)
osrobj3.ExportToWkt()
Out[10]:
'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"]]'
In [11]:
osrobj3.ExportToMICoordSys()
Out[11]:
'Earth Projection 1, 104'
In [12]:
osrobj3.ExportToPCI()
Out[12]:
['LONG/LAT D000', 'DEGREE', (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)]
In [13]:
osrobj3.ExportToUSGS()
Out[13]:
[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), 12]
In [14]:
osrobj3.ExportToProj4()
Out[14]:
'+proj=longlat +datum=WGS84 +no_defs'
In [15]:
osrobj3.IsGeographic()
Out[15]:
1
proj.4 字符串转成Esri wkt字符串的方法 April 12, 2016barry.z proj.4 字符串转成Esri wkt字符串的方法 使用gdal/ogr库进行转换
In [16]:
import os
import sys
import string
import osgeo.osr
srs = osgeo.osr.SpatialReference()
srs.ImportFromProj4(sys.argv[1])
srs.MorphToESRI() # 重点，转成Esri的wkt格式
srs.ExportToWkt()
ERROR 1: PROJ: proj_create: unrecognized format / unknown name
Out[16]:
''
In [17]:
from osgeo import ogr, osr
driver = ogr.GetDriverByName('ESRI Shapefile')
dataset = driver.Open('/data/gdata/region_popu.shp')
layer = dataset.GetLayer()
spatialRef = layer.GetSpatialRef()
spatialRef.ExportToWkt()
Out[17]:
'PROJCS["WGS 84 / Pseudo-Mercator",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"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]'
In [18]:
spatialRef.ExportToPrettyWkt()
Out[18]:
'PROJCS["WGS 84 / Pseudo-Mercator",\n GEOGCS["WGS 84",\n DATUM["WGS_1984",\n SPHEROID["WGS 84",6378137,298.257223563,\n AUTHORITY["EPSG","7030"]],\n AUTHORITY["EPSG","6326"]],\n PRIMEM["Greenwich",0,\n AUTHORITY["EPSG","8901"]],\n UNIT["degree",0.0174532925199433,\n AUTHORITY["EPSG","9122"]],\n AUTHORITY["EPSG","4326"]],\n PROJECTION["Mercator_1SP"],\n PARAMETER["central_meridian",0],\n PARAMETER["scale_factor",1],\n PARAMETER["false_easting",0],\n PARAMETER["false_northing",0],\n UNIT["metre",1,\n AUTHORITY["EPSG","9001"]],\n AXIS["Easting",EAST],\n AXIS["Northing",NORTH],\n EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"],\n AUTHORITY["EPSG","3857"]]'
In [19]:
spatialRef.ExportToPCI()
Out[19]:
['MER D894', 'METRE', (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)]
In [20]:
spatialRef.ExportToUSGS()
Out[20]:
[5, 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), 12]
In [21]:
from pprint import pprint
pprint(spatialRef.ExportToXML())
6
ERROR 6: Unhandled projection method Mercator_1SP
In [22]:
from osgeo import ogr, osr
spatialRef = osr.SpatialReference()
spatialRef.ImportFromEPSG(26912)
spatialRef.MorphToESRI()
file = open('yourshpfile.prj', 'w')
file.write(spatialRef.ExportToWkt())
file.close()
In [23]:
from osgeo import osr
from osgeo import gdal
dataset = gdal.Open("/data/gdata/geotiff_file.tif")
从数据集中获取空间参考并且建立一个SpatialReference对象
In [24]:
sr = dataset.GetProjectionRef()
osrobj = osr.SpatialReference()
osrobj3 = osr.SpatialReference()
osrobj3.SetWellKnownGeogCS("WGS84")
osrobj3.IsSame(osrobj)
osrobj3.ExportToWkt()
Out[24]:
'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"]]'
In [25]:
osrobj3.IsGeographic()
Out[25]:
1
In [26]:
ct = osr.CoordinateTransformation(osrobj,osrobj3)
# ct.TransformPoint([590000.0,4928000.0,0])
# ct.TransformPoint(609000,4928000)
# ct.TransformPoint(609000,4914000)
# ct.TransformPoint(590000,4914000)
ERROR 1: PROJ: proj_create: unrecognized format / unknown name
下面的代码，展示了从 EPSG 2927 到 EPSG 4326 的投影转换：
In [27]:
from osgeo import ogr
from osgeo import osr
source = osr.SpatialReference()
source.ImportFromEPSG(2927)
target = osr.SpatialReference()
target.ImportFromEPSG(4326)
transform = osr.CoordinateTransformation(source, target)
point = ogr.CreateGeometryFromWkt("POINT (1120351.57 741921.42)")
point.Transform(transform)
point.ExportToWkt()
Out[27]:
'POINT (47.3488070138318 -122.598149943144)'