OGR坐标参考系和坐标转换教程

介绍

这个 OGRSpatialReference 和 OGRCoordinateTransformation 类分别提供表示坐标参考系（称为CRS或SRS，例如通常将地图投影与大地测量基准相关联的投影CRS）和在它们之间转换的服务。这些服务松散地基于OpenGIS坐标转换规范建模，并依赖于众所周知的文本（WKT）格式（其各种版本：OGC WKT 1、ESRI WKT、WKT2:2015和WKT2:2018）来描述坐标系。

参考文献和适用标准

• PROJ documentation: 投影方法与坐标运算

• ISO:19111 and WKT standards

• GeoTIFF Projections Transform List: 了解GeoTIFF的WKT投影公式

• EPSG Geodesy web page 也是一种有用的资源

定义地理坐标参照系

CRS封装在 OGRSpatialReference 上课。有多种方法可以将OGRSpatialReference对象初始化为有效的坐标参考系。CRS主要有两种。第一种是地理位置（位置以长/纬度测量），第二种是投影位置（如UTM-位置以米或英尺测量）。

A Geographic CRS contains information on the datum (which implies a spheroid described by a semi-major axis, and inverse flattening), prime meridian (normally Greenwich), and an angular units type which is normally degrees. The following code initializes a geographic CRS on supplying all this information along with a user visible name for the geographic CRS.

OGRSpatialReference oSRS;

oSRS.SetGeogCS( "My geographic CRS",
                "World Geodetic System 1984",
                "My WGS84 Spheroid",
                SRS_WGS84_SEMIMAJOR, SRS_WGS84_INVFLATTENING,
                "Greenwich", 0.0,
                "degree", SRS_UA_DEGREE_CONV );

备注

中的缩写CS OGRSpatialReference::SetGeogCS() 根据当前测地线术语，不合适，应理解为CRS

其中，“My geographic CRS”、“My WGS84 Spheroid”、“Greenwich”和“degree”不是键，而是用于向用户显示。但是，基准名称“World Geodetic System 1984”用作标识基准的键，应设置为EPSG注册表中的已知值，以便在坐标操作期间进行适当的基准转换。有效大地基准列表见 geodetic_datum.sql 文件。

备注

在WKT 1中，基准名称中的空格字符通常用下划线代替。WGSèu1984被用作“世界大地测量系统1984”的别名

OGRSpatialReference内置了对一些著名的CRS的支持，这些CRS包括“NAD27”、“NAD83”、“WGS72”和“WGS84”，可以在对 OGRSpatialReference::SetWellKnownGeogCS() .

oSRS.SetWellKnownGeogCS( "WGS84" );

备注

根据当前测地术语，SetWellKnownGeogCS（）中的缩写CS不合适，应理解为CRS

此外，如果EPSG数据库可用，EPSG数据库中的任何地理CRS都可以通过其GCS代码进行设置。

oSRS.SetWellKnownGeogCS( "EPSG:4326" );

为了序列化和将投影定义传输到其他包，使用了OpenGIS众所周知的坐标系文本格式。OGRSpatialReference可以从WKT初始化，也可以转换回WKT。从GDAL 3.0开始，WKT导出的默认格式仍然是OGC WKT 1。

char *pszWKT = NULL;
oSRS.SetWellKnownGeogCS( "WGS84" );
oSRS.exportToWkt( &pszWKT );
printf( "%s\n", pszWKT );
CPLFree(pszWKT);

输出：

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"]]

或以更可读的形式：

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"]]

从GDAL 3.0开始 OGRSpatialReference::exportToWkt() 方法接受选项，

char *pszWKT = nullptr;
oSRS.SetWellKnownGeogCS( "WGS84" );
const char* apszOptions[] = { "FORMAT=WKT2_2018", "MULTILINE=YES", nullptr };
oSRS.exportToWkt( &pszWKT, apszOptions );
printf( "%s\n", pszWKT );
CPLFree(pszWKT);

GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

此方法和选项在C中作为 OSRExportToWktEx() 功能。

这个 OGRSpatialReference::importFromWkt() 方法可用于从WKT CRS定义设置OGRSpatialReference。

CRS和轴顺序

前几节中省略的一个“细节”是CRS中坐标轴顺序的主题。根据ISO:19111建模，地理CRS由两个主要部分组成：大地基准面和 coordinate system . 对于2D地理CRS，坐标系轴是经度和纬度，沿这些轴的值通常用度表示（基于古法语的CRS可以使用grad）。

The order in which they are specified, that is latitude first, longitude second, or the reverse, is a constant matter of confusion and vary depending on conventions used by geodetic authorities, GIS user, file format and protocol specifications, etc. This is the source of various interoperability issues.

在GDAL 3.0之前 OGRSpatialReference 类不遵守由定义CRS的机构强制的轴顺序，因此当顺序是纬度第一、经度第二时，从WKT字符串中剥离轴顺序信息。使用ogrcoordinateransformation类的坐标转换还假定该类的Transform（）方法传递或返回的地理坐标使用经度、纬度顺序。

从GDAL 3.0开始，由定义CRS的权威授权的axis命令默认由OGRCoordinateTransformation类执行，并始终在WKT1中导出。因此，使用“EPSG:4326”或“WGS84”字符串创建的CRS使用纬度第一、经度第二轴顺序。

为了帮助从仍然使用经度、纬度顺序坐标的代码基迁移，可以将元数据信息附加到OGRSpatialReference实例，以指定为了进行坐标转换，有效传递或返回的值的顺序将是经度、纬度。为此，必须调用以下命令

oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);

争论传给 OGRSpatialReference::SetAxisMappingStrategy() 是数据轴到CRS轴的映射策略。

• OAMS_TRADITIONAL_GIS_ORDER 也就是说，对于lat/long顺序的地理CRS，数据仍然是long/lat顺序的。类似地，对于具有北距/东距顺序的投影CRS，数据仍将是东距/北距顺序的。

• OAMS_AUTHORITY_COMPLIANT means that the data axis will be identical to the CRS axis. This is the default value when instantiating OGRSpatialReference.

• OAMS_CUSTOM means that the data axes are customly defined with SetDataAxisToSRSAxisMapping().

本节讨论的地理CRS的特殊情况也适用于预计的CRS。虽然他们大多数使用东距第一，北距第二公约，一些定义在EPSG注册表使用相反的公约。

Another way to keep using the Traditional GIS order for some specific well known CRS is to calling to OGRSpatialReference::SetWellKnownGeogCS() with "CRS27", "CRS83" or "CRS84" instead of "NAD27", "NAD83" and "WGS84" respectively.

oSRS.SetWellKnownGeogCS( "CRS84" );

定义预计的CRS

A projected CRS (such as UTM, Lambert Conformal Conic, etc.) requires and underlying geographic CRS as well as a definition for the projection transform used to translate between linear positions (in meters or feet) and angular long/lat positions. The following code defines a UTM zone 17 projected CRS with an underlying geographic CRS (datum) of WGS84.

OGRSpatialReference oSRS;

oSRS.SetProjCS( "UTM 17 (WGS84) in northern hemisphere." );
oSRS.SetWellKnownGeogCS( "WGS84" );
oSRS.SetUTM( 17, TRUE );

调用 OGRSpatialReference::SetProjCS() 为投影的CRS设置用户名并确定系统已投影。这个 OGRSpatialReference::SetWellKnownGeogCS() 将地理坐标系与 OGRSpatialReference::SetUTM() 调用设置详细的投影转换参数。此时，为了创建有效的定义，上面的顺序很重要，但是在将来，对象将根据需要自动重新排序内部表示以保持有效。

小心

现在，请注意定义OGRSpatialReference的步骤顺序！

上面的定义会给出一个WKT版本，看起来像下面这样。请注意，UTM 17扩展为UTM区域的横向墨卡托定义的细节。

PROJCS["UTM 17 (WGS84) in northern hemisphere.",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG",7030]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG",6326]],
        PRIMEM["Greenwich",0,AUTHORITY["EPSG",8901]],
        UNIT["DMSH",0.0174532925199433,AUTHORITY["EPSG",9108]],
        AXIS["Lat",NORTH],
        AXIS["Long",EAST],
        AUTHORITY["EPSG",4326]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-81],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0]]

有许多投影方法，包括 OGRSpatialReference::SetTM() （横向墨卡托）， OGRSpatialReference::SetLCC() （Lambert保角圆锥曲线），以及 OGRSpatialReference::SetMercator() .

坐标参照系查询

一旦建立了OGRSpatialReference，就可以查询关于它的各种信息。如果是预测或地理CRS，可以使用 OGRSpatialReference::IsProjected() 和 OGRSpatialReference::IsGeographic() 方法。这个 OGRSpatialReference::GetSemiMajor() ， OGRSpatialReference::GetSemiMinor() 和 OGRSpatialReference::GetInvFlattening() 方法可用于获取有关球体的信息。这个 OGRSpatialReference::GetAttrValue() 方法可用于获取PROJCS、GEOGCS、DATUM、SPHEROID和PROJECTION名称字符串。这个 OGRSpatialReference::GetProjParm() 方法可用于获取投影参数。这个 OGRSpatialReference::GetLinearUnits() 方法可用于获取线性单位类型，并转换为米。

请注意，投影方法和参数的名称是WKT 1中的一个。

下面的代码演示了 OGRSpatialReference::GetAttrValue() 得到投影图 OGRSpatialReference::GetProjParm() 以获取投影参数。GetAttrValue（）方法在WKT文本表示中搜索与命名项关联的第一个“value”节点。使用GetProjParm（）获取投影参数时，应使用投影参数（例如SRS_PP_CENTRAL_MERIDIAN）的“define”常量。可以参考ogrspatialreference.cpp中各种投影的Set方法的代码来查找哪些参数适用于哪些投影。

const char *pszProjection = poSRS->GetAttrValue("PROJECTION");

if( pszProjection == NULL )
{
    if( poSRS->IsGeographic() )
        sprintf( szProj4+strlen(szProj4), "+proj=longlat " );
    else
        sprintf( szProj4+strlen(szProj4), "unknown " );
}
else if( EQUAL(pszProjection,SRS_PT_CYLINDRICAL_EQUAL_AREA) )
{
    sprintf( szProj4+strlen(szProj4),
    "+proj=cea +lon_0=%.9f +lat_ts=%.9f +x_0=%.3f +y_0=%.3f ",
            poSRS->GetProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
            poSRS->GetProjParm(SRS_PP_STANDARD_PARALLEL_1,0.0),
            poSRS->GetProjParm(SRS_PP_FALSE_EASTING,0.0),
            poSRS->GetProjParm(SRS_PP_FALSE_NORTHING,0.0) );
}
...

坐标变换

这个 OGRCoordinateTransformation 类用于在不同的CRS之间转换位置。新的转换对象是使用 OGRCreateCoordinateTransformation() 然后 OGRCoordinateTransformation::Transform() 方法可用于在CRS之间转换点。

OGRSpatialReference oSourceSRS, oTargetSRS;
OGRCoordinateTransformation *poCT;
double x, y;

oSourceSRS.importFromEPSG( atoi(papszArgv[i+1]) );
oTargetSRS.importFromEPSG( atoi(papszArgv[i+2]) );

poCT = OGRCreateCoordinateTransformation( &oSourceSRS,
                                          &oTargetSRS );
x = atof( papszArgv[i+3] );
y = atof( papszArgv[i+4] );

if( poCT == NULL || !poCT->Transform( 1, &x, &y ) )
    printf( "Transformation failed.\n" );
else
{
    printf( "(%f,%f) -> (%f,%f)\n",
            atof( papszArgv[i+3] ),
            atof( papszArgv[i+4] ),
            x, y );
}

There are a couple of points at which transformations can fail. First, OGRCreateCoordinateTransformation() may fail, generally because the internals recognize that no transformation between the indicated systems can be established, and will return a NULL pointer.

OGRCoordinateTransformation：：Transform（）方法本身也可能失败。这可能是由于上述问题之一造成的延迟，或者是由于一个或多个传入点的运算在数值上未定义。函数的作用是：转换成功后返回TRUE，如果转换失败则返回FALSE。出错时点数组处于不确定状态。

尽管上面没有显示，坐标转换服务可以获取三维点，并将根据球体和基准面的高程差调整高程。地理或投影CRS上给出的高程假定为椭球高度。当使用由水平CRS（地理或投影）和垂直CRS组成的复合CRS时，高程将与垂直基准（平均海平面、基于重力的等）相关。

Starting with GDAL 3.0, a time value (generally as a value in decimal years) can also be specified for time-dependent coordinate operations.

下面的示例演示如何使用与投影坐标系相同的地理坐标系方便地创建long/lat坐标系，并使用该坐标系在投影坐标和long/lat坐标系之间进行转换。返回的坐标将以经度表示，由于调用SetAxisMappingStrategy（OAMSU TRADITIONALU GISU order）而导致的纬度顺序

OGRSpatialReference    oUTM, *poLongLat;
OGRCoordinateTransformation *poTransform;

oUTM.SetProjCS("UTM 17 / WGS84");
oUTM.SetWellKnownGeogCS( "WGS84" );
oUTM.SetUTM( 17 );

poLongLat = oUTM.CloneGeogCS();
poLongLat->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);

poTransform = OGRCreateCoordinateTransformation( &oUTM, poLongLat );
if( poTransform == NULL )
{
    ...
}

...

if( !poTransform->Transform( nPoints, x, y, z ) )
...

高级坐标变换

OGRCreateCoordinateTransformation() under-the-hood may determine several candidate coordinate operations transforming from the source CRS to the target CRS. Those candidate coordinate operations each have their own area of use. When Transform() is invoked, it will determine the most appropriate coordinate operation based on the coordinates of the point to transform and area of use. For example, there are several dozens of possible coordinate operations for the NAD27 to WGS84 transformation.

如果已知要转换的坐标所在的感兴趣区域的边界框，则可以将其指定为限制候选坐标操作以考虑：

OGRCoordinateTransformationOptions options;
options.SetAreaOfInterest(-100,40,-99,41);
poTransform = OGRCreateCoordinateTransformation( &oNAD27, &oWGS84, options );

对于必须使用特定坐标操作的情况，可以将其指定为PROJ字符串（以+PROJ=pipeline开头的单步操作或多步字符串）、描述协调操作的WKT2字符串或urn:ogc:def:coordinate operation:EPSG：：XXXX urn

OGRCoordinateTransformationOptions options;

// EPSG:8599, NAD27 to WGS 84 (46), 1.15 m, USA - Indiana
options.SetCoordinateOperation(
    "+proj=pipeline +step +proj=axisswap +order=2,1 "
    "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
    "+step +proj=hgridshift +grids=conus "
    "+step +proj=hgridshift +grids=inhpgn.gsb "
    "+step +proj=unitconvert +xy_in=rad +xy_out=deg +step "
    "+proj=axisswap +order=2,1", false );

// or
// options.SetCoordinateOperation(
//      "urn:ogc:def:coordinateOperation:EPSG::8599", false);

poTransform = OGRCreateCoordinateTransformation( &oNAD27, &oWGS84, options );

备用接口

在ogr_srs_api.h中定义了坐标系服务的C接口，并且可以通过osr.py模块使用Python绑定。方法是C++方法的类似的类似方法，但是对于C++方法来说，C和Python绑定缺少。

C绑定

typedef void *OGRSpatialReferenceH;
typedef void *OGRCoordinateTransformationH;

OGRSpatialReferenceH OSRNewSpatialReference( const char * );
void    OSRDestroySpatialReference( OGRSpatialReferenceH );

int     OSRReference( OGRSpatialReferenceH );
int     OSRDereference( OGRSpatialReferenceH );

void OSRSetAxisMappingStrategy( OGRSpatialReferenceH,
                                OSRAxisMappingStrategy );

OGRErr  OSRImportFromEPSG( OGRSpatialReferenceH, int );
OGRErr  OSRImportFromWkt( OGRSpatialReferenceH, char ** );
OGRErr  OSRExportToWkt( OGRSpatialReferenceH, char ** );
OGRErr  OSRExportToWktEx( OGRSpatialReferenceH, char **,
                        const char* const* papszOptions );

OGRErr  OSRSetAttrValue( OGRSpatialReferenceH hSRS, const char * pszNodePath,
                        const char * pszNewNodeValue );
const char *OSRGetAttrValue( OGRSpatialReferenceH hSRS,
                            const char * pszName, int iChild);

OGRErr  OSRSetLinearUnits( OGRSpatialReferenceH, const char *, double );
double  OSRGetLinearUnits( OGRSpatialReferenceH, char ** );

int     OSRIsGeographic( OGRSpatialReferenceH );
int     OSRIsProjected( OGRSpatialReferenceH );
int     OSRIsSameGeogCS( OGRSpatialReferenceH, OGRSpatialReferenceH );
int     OSRIsSame( OGRSpatialReferenceH, OGRSpatialReferenceH );

OGRErr  OSRSetProjCS( OGRSpatialReferenceH hSRS, const char * pszName );
OGRErr  OSRSetWellKnownGeogCS( OGRSpatialReferenceH hSRS,
                            const char * pszName );

OGRErr  OSRSetGeogCS( OGRSpatialReferenceH hSRS,
                    const char * pszGeogName,
                    const char * pszDatumName,
                    const char * pszEllipsoidName,
                    double dfSemiMajor, double dfInvFlattening,
                    const char * pszPMName ,
                    double dfPMOffset ,
                    const char * pszUnits,
                    double dfConvertToRadians );

double  OSRGetSemiMajor( OGRSpatialReferenceH, OGRErr * );
double  OSRGetSemiMinor( OGRSpatialReferenceH, OGRErr * );
double  OSRGetInvFlattening( OGRSpatialReferenceH, OGRErr * );

OGRErr  OSRSetAuthority( OGRSpatialReferenceH hSRS,
                        const char * pszTargetKey,
                        const char * pszAuthority,
                        int nCode );
OGRErr  OSRSetProjParm( OGRSpatialReferenceH, const char *, double );
double  OSRGetProjParm( OGRSpatialReferenceH hSRS,
                        const char * pszParamName,
                        double dfDefault,
                        OGRErr * );

OGRErr  OSRSetUTM( OGRSpatialReferenceH hSRS, int nZone, int bNorth );
int     OSRGetUTMZone( OGRSpatialReferenceH hSRS, int *pbNorth );

OGRCoordinateTransformationH
OCTNewCoordinateTransformation( OGRSpatialReferenceH hSourceSRS,
                                OGRSpatialReferenceH hTargetSRS );

void OCTDestroyCoordinateTransformation( OGRCoordinateTransformationH );

int OCTTransform( OGRCoordinateTransformationH hCT,
                int nCount, double *x, double *y, double *z );

OGRCoordinateTransformationOptionsH OCTNewCoordinateTransformationOptions(;

int OCTCoordinateTransformationOptionsSetOperation(
    OGRCoordinateTransformationOptionsH hOptions,
    const char* pszCO, int bReverseCO);

int OCTCoordinateTransformationOptionsSetAreaOfInterest(
    OGRCoordinateTransformationOptionsH hOptions,
    double dfWestLongitudeDeg,
    double dfSouthLatitudeDeg,
    double dfEastLongitudeDeg,
    double dfNorthLatitudeDeg);

void OCTDestroyCoordinateTransformationOptions(OGRCoordinateTransformationOptionsH);

OGRCoordinateTransformationH
OCTNewCoordinateTransformationEx( OGRSpatialReferenceH hSourceSRS,
                                OGRSpatialReferenceH hTargetSRS,
                                OGRCoordinateTransformationOptionsH hOptions );

Python包

class osr.SpatialReference
    def __init__(self,obj=None):
    def SetAxisMappingStrategy( self, strategy ):
    def ImportFromWkt( self, wkt ):
    def ExportToWkt(self, options = None):
    def ImportFromEPSG(self,code):
    def IsGeographic(self):
    def IsProjected(self):
    def GetAttrValue(self, name, child = 0):
    def SetAttrValue(self, name, value):
    def SetWellKnownGeogCS(self, name):
    def SetProjCS(self, name = "unnamed" ):
    def IsSameGeogCS(self, other):
    def IsSame(self, other):
    def SetLinearUnits(self, units_name, to_meters ):
    def SetUTM(self, zone, is_north = 1):

class CoordinateTransformation:
    def __init__(self,source,target):
    def TransformPoint(self, x, y, z = 0):
    def TransformPoints(self, points):

历史和实施注意事项

Before GDAL 3.0, the OGRSpatialReference class was strongly tied to OGC WKT (WKT 1) format specified by Coordinate Transformation Services (CT) specification (01-009), and the way it was interpreted by GDAL, which various caveats detailed in the OGC WKT坐标系问题 page. The class mostly contained an in-memory tree-like representation of WKT 1 strings. The class used to directly implement import and export to OGC WKT 1, WKT-ESRI and PROJ.4 formats. Reprojection services were only available if GDAL had been build against the PROJ library.

从GDAL 3.0开始 PROJ >=6.0库已成为GDAL的必需依赖项。PROJ 6内置了对OGC WKT 1、ESRI WKT、OGC WKT 2:2015和OGC WKT 2:2018表示的支持。PROJ 6还实现了ISO-19111/OGC抽象主题2的C++对象类层次结构“参照坐标”标准。因此，OGRSpatialReference类被修改为主要充当PROJ PJ*CRS对象之上的包装器，并尝试尽可能地从OGC WKT 1表示中抽象出来。但是，为了向后兼容，有些方法仍然需要OGC WKT 1特定的参数或返回值。OGRSpatialReference类的设计仍然是单片的。想要直接和细粒度访问CRS表示的用户可能希望直接使用PROJ 6 C或C++ API。