Gdal source code analysis (11) OGR projection description

I. Introduction

For more information, see http://www.gdal.org/ogr/osr_tutorial.html.

OgrspatialreferenceClass andOgrcoordinatetransformationClass is mainly used to provide the definition coordinate system (projection and level) and the conversion coordinate. Both classes are based on the coordinate conversion instructions of OpenGIS and use well
Known text format to express the coordinate system.

Some information about the OpenGIS coordinate system and spatial reference coordinate abstract model files can be found on the OGC (Open Geospatial Consortium) website. Geotiff projection conversion list (geotiff
Projections transform list) web pages help you better understand WKT rules, and epsg websites are also useful materials.

Ii. define the geographic coordinate system

Coordinate System usageOgrspatialreferenceClass. Several types of initialization are provided hereOgrspatialreferenceClass. There are two main types of Coordinate Systems, the first is the geographical coordinate system (the location information is expressed by latitude and longitude), and the second is the projection coordinate system (such as UTM-universal horizontal axis mocato projection, location Information is expressed in meters or feet ).

The information contained in a geographic coordinate system has a (which contains a torball body represented by the reciprocal of the long half axis and flat rate), a central Meridian (usually the primary meridian, that is, 0-degree longitude line Greenwich), there is also an angle measurement unit, using degrees rather than radians. If such information is contained, a valid geographic coordinate system can be constructed.

OGRSpatialReference oSRS;oSRS.SetGeogCS( "Mygeographic coordinate system",               "WGS_1984",               "My WGS84 Spheroid",               SRS_WGS84_SEMIMAJOR, SRS_WGS84_INVFLATTENING,               "Greenwich", 0.0,                "degree", SRS_UA_DEGREE_CONV );

In the above Code, the names "My geographic coordinate system", "My WGS84 spheroid", "Greenwich" and "degree" are not keywords, these are mainly used to explain to users. However, "wgs_1984" is a key word defining the earth benchmark. Note: The Earth benchmark here must be an effective Earth benchmark! (The preceding strings are randomly specified and used for display. The strings at the position wgs_1984 must be valid and cannot be named randomly, ).

OgrspatialreferenceYou can use common strings to create a common coordinate system. These strings include "nad27", "nad83", "wgs72", and "WGS84, the setwellknowngeogcs () function is used. The usage is shown below.

oSRS.SetWellKnownGeogCS( "WGS84" );

In addition, you can use the GCS Code defined by the epsg database to define a geographic coordinate system, such:

oSRS.SetWellKnownGeogCS( "EPSG:4326" );

To facilitate association with other libraries, the WKT format of OpenGIS can be converted here. SimilarlyOgrspatialreferenceYou can use a WKT for initialization or export the information in the WKT format.

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

The preceding statement outputs the following content:

GEOGCS["WGS84",DATUM["WGS_1984",SPHEROID["WGS84",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]]

Adjust the above string format to look better:

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

FunctionOgrspatialreference: importfromwkt ()A wkt-defined coordinate system can be used to constructOgrspatialreferenceClass Object.

3. Define the projection Coordinate System

A projection coordinate system (such as UTM and lanbert equi-angle cone projection) must be built on a geographic coordinate system. In the projection coordinate system, coordinate points are expressed by length units such as meters or feet, and can also be expressed by angle coordinates of longitude and latitude. The following defines a projection Coordinate System of the UTM 17th band based on the WGS84 Earth Reference elliptical body.

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

Call the setprojcs () function to set the name of the projection coordinate system, use the setwellknowngeogcs () function to specify the geographic coordinate system, and call the setutm () function to set the projection conversion parameter information. After completing these tasks, a valid projection coordinate system is defined.Note the order of defining ogrspatialreference here!

The projection defined above is represented in the following format using WKT. Note that utm17 will be represented by the cross-axis mokto projection parameter.

PROJCS["UTM 17 (WGS84) in northernhemisphere.",   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]]

There are many ways to set the projection coordinates, including settm () (horizontal mocato projection), setkp () (lanbert equi-angle cone projection), and setmercator ().

Iv. Obtain Coordinate System Information

OnceOgrspatialreferenceObject creation, you can obtain various information in it, for example, you can useOgrspatialreference: isprojected ()AndOgrspatialreference: isgeographic ()Method to Determine whether the coordinate system is a geographic coordinate system or a projection coordinate system. Use FunctionsOgrspatialreference: getsemimajor (),Ogrspatialreference: getsemiminor ()AndOgrspatialreference: getinvflattening ()The method can be used to obtain information about the elliptical body, including the long half axis, the short half axis, and the reciprocal of the flat rate. UseOgrspatialreference: getattrvalue ()The method can be used to obtain the name strings of projcs, geogcs, datum, spheroid, and projection. UseOgrspatialreference: getprojparm ()Method To obtain projection parameter information. UseOgrspatialreference: getlinearunits ()You can obtain the length unit type and convert it to meters.

The following code is an example of obtaining the coordinate system information, which is from the ogr_srs_proj4.cpp file. Use getattrvalue () to obtain the projection name, and use getprojparm () to obtain the projection parameter information. The first value node of the getattrvalue () function is the description of the WKT string. Macro definitions of projection parameters (such as srs_pp_central_meridian) and functions getprojparm () can be used together to obtain projection parameters. For more usage instructions, refer to the relevant code in the ogrspatialreference. cpp file.

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));}...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));}...

V. Coordinate Transformation

OgrcoordinatetransformationClass can be used for Coordinate Transformation in different coordinate systems. Functions availableOgrcreatecoordinatetransformation ()Create a new coordinate conversion object, and then useOgrcoordinatetransformation: Transform ()Method.

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 );

The reasons for point-to-point conversion failure are as follows,Ogrcreatecoordinatetransformation ()Function execution fails because the Conversion Relationship between the specified two projections is unknown. This may be a trial of projection not supported by the proj.4 library. The conversion relationships between different elliptical bodies are not defined, or one of the coordinate systems is not completely defined. If the FunctionOgrcreatecoordinatetransformation ()If the execution fails, the return value is null.

Ogrcoordinatetransformation: Transform ()The method page may fail to be executed. The possible cause is also the above problem, or there is a number not defined in the conversion point. If the function transform () is successfully executed, true is returned. If one point fails to be converted, false is returned.

Coordinate conversion can also support coordinate transformation of three-dimensional points. The elevation values are automatically adjusted based on the datum of different elliptical regions. This can be used for Coordinate Transformation in different vertical datum planes. If there is no Z value, the system will think that all vertices are on the level.

The following example shows how to convert a point in a projection coordinate system to a point in the geographic coordinate system of the projection to a coordinate indicated by longitude and latitude.

OGRSpatialReference   oUTM, *poLatLong;OGRCoordinateTransformation *poTransform;oUTM.SetProjCS("UTM 17 /WGS84");oUTM.SetWellKnownGeogCS("WGS84" );oUTM.SetUTM( 17 );poLatLong = oUTM.CloneGeogCS(); poTransform = OGRCreateCoordinateTransformation( &oUTM,poLatLong );if( poTransform == NULL ){    ...} ...if( !poTransform->Transform( nPoints, x,y, z ) )...

Vi. interfaces in other languages

The OGR space reference provides a C language interface, which is defined in the ogr_srs_api.h file. The Python interface is defined in the OSR. py file. The names of all interfaces are similar to those of C ++, but some methods in C and Python are not provided.

1. C Language Interface

typedef void *OGRSpatialReferenceH;                              typedef void *OGRCoordinateTransformationH;OGRSpatialReferenceH OSRNewSpatialReference( const char *);void   OSRDestroySpatialReference( OGRSpatialReferenceH );int    OSRReference( OGRSpatialReferenceH );int    OSRDereference( OGRSpatialReferenceH );OGRErr OSRImportFromEPSG( OGRSpatialReferenceH, int );OGRErr OSRImportFromWkt( OGRSpatialReferenceH, char ** );OGRErr OSRExportToWkt( OGRSpatialReferenceH, char ** );OGRErr OSRSetAttrValue( OGRSpatialReferenceH hSRS, const char * pszNodePath,                         const char *pszNewNodeValue );const char *OSRGetAttrValue( OGRSpatialReferenceH hSRS,                             const char * pszName, intiChild);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 *pszParmName,                        double dfDefault,                        OGRErr * );OGRErr OSRSetUTM( OGRSpatialReferenceH hSRS, int nZone, int bNorth );int    OSRGetUTMZone( OGRSpatialReferenceH hSRS, int *pbNorth );OGRCoordinateTransformationHOCTNewCoordinateTransformation( OGRSpatialReferenceH hSourceSRS,                                OGRSpatialReferenceHhTargetSRS );void OCTDestroyCoordinateTransformation( OGRCoordinateTransformationH );int OCTTransform(OGRCoordinateTransformationH hCT,                  int nCount, double *x, double*y, double *z );

2. Python Interface

class osr.SpatialReference   def __init__(self,obj=None):   def ImportFromWkt( self, wkt ):   def ExportToWkt(self):   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):

7. Internal implementation

Ogrcoordinatetransformation depends on the proj.4 library. Therefore, to use the coordinate conversion content, gdal must be bound to proj4 during compilation for coordinate conversion. If you want to use the gdal Coordinate Transformation and re-projection related algorithms, you must have support for the proj4 Library; otherwise, the conversion will fail. For details about how to compile proj4 and bind gdal to proj4, refer to the previous blog.