1/*! \page osr_tutorial OGR Projections Tutorial 2 3\section osr_tutorial_intro Introduction 4 5The OGRSpatialReference, and OGRCoordinateTransformation classes provide 6services to represent coordinate systems (projections and datums) and to 7transform between them. These services are loosely modelled on the 8OpenGIS Coordinate Transformations specification, and use the same 9Well Known Text format for describing coordinate systems. 10 11Some background on OpenGIS coordinate systems and services can be found 12in the Simple Features for COM, and Spatial Reference Systems Abstract Model 13documents available from the <a href="http://www.opengeospatial.org">Open Geospatial Consortium</a>. 14The <a href="http://www.remotesensing.org/geotiff/proj_list">GeoTIFF Projections Transform List</a> 15may also be of assistance in 16understanding formulations of projections in WKT. The 17<a href="http://www.epsg.org">EPSG</a> 18Geodesy web page is also a useful resource. 19 20\section osr_tutorial_cs Defining a Geographic Coordinate System 21 22Coordinate systems are encapsulated in the OGRSpatialReference class. There 23are a number of ways of initializing an OGRSpatialReference object to a 24valid coordinate system. There are two primary kinds of coordinate systems. 25The first is geographic (positions are measured in long/lat) and the second 26is projected (such as UTM - positions are measured in meters or feet). 27 28A Geographic coordinate system contains information on the datum (which implies 29an spheroid described by a semi-major axis, and inverse flattening), prime 30meridian (normally Greenwich), and an angular units type which is normally 31degrees. The following code initializes a geographic coordinate system 32on supplying all this information along with a user visible name for the 33geographic coordinate system. 34 35\code 36 OGRSpatialReference oSRS; 37 38 oSRS.SetGeogCS( "My geographic coordinate system", 39 "WGS_1984", 40 "My WGS84 Spheroid", 41 SRS_WGS84_SEMIMAJOR, SRS_WGS84_INVFLATTENING, 42 "Greenwich", 0.0, 43 "degree", SRS_UA_DEGREE_CONV ); 44\endcode 45 46Of these values, the names "My geographic coordinate system", "My WGS84 47Spheroid", "Greenwich" and "degree" are not keys, but are used for display 48to the user. However, the datum name "WGS_1984" is used as a key to identify 49the datum, and there are rules on what values can be used. NOTE: Prepare 50writeup somewhere on valid datums! 51 52The OGRSpatialReference has built in support for a few well known coordinate 53systems, which include "NAD27", "NAD83", "WGS72" and "WGS84" which can be 54defined in a single call to SetWellKnownGeogCS(). 55 56\code 57 oSRS.SetWellKnownGeogCS( "WGS84" ); 58\endcode 59 60Furthermore, any geographic coordinate system in the EPSG database can 61be set by it's GCS code number if the EPSG database is available. 62 63\code 64 oSRS.SetWellKnownGeogCS( "EPSG:4326" ); 65\endcode 66 67For serialization, and transmission of projection definitions to other 68packages, the OpenGIS Well Known Text format for coordinate systems is 69used. An OGRSpatialReference can be initialized from well known text, or 70converted back into well known text. 71 72\code 73 char *pszWKT = NULL; 74 75 oSRS.SetWellKnownGeogCS( "WGS84" ); 76 oSRS.exportToWkt( &pszWKT ); 77 printf( "%s\n", pszWKT ); 78\endcode 79 80gives something like: 81 82<pre> 83GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563, 84AUTHORITY["EPSG",7030]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG",6326]], 85PRIMEM["Greenwich",0,AUTHORITY["EPSG",8901]],UNIT["DMSH",0.0174532925199433, 86AUTHORITY["EPSG",9108]],AXIS["Lat",NORTH],AXIS["Long",EAST],AUTHORITY["EPSG", 874326]] 88</pre> 89 90or in more readable form: 91 92<pre> 93GEOGCS["WGS 84", 94 DATUM["WGS_1984", 95 SPHEROID["WGS 84",6378137,298.257223563, 96 AUTHORITY["EPSG",7030]], 97 TOWGS84[0,0,0,0,0,0,0], 98 AUTHORITY["EPSG",6326]], 99 PRIMEM["Greenwich",0,AUTHORITY["EPSG",8901]], 100 UNIT["DMSH",0.0174532925199433,AUTHORITY["EPSG",9108]], 101 AXIS["Lat",NORTH], 102 AXIS["Long",EAST], 103 AUTHORITY["EPSG",4326]] 104</pre> 105 106The OGRSpatialReference::importFromWkt() method can be used to set an 107OGRSpatialReference from a WKT coordinate system definition. 108 109\section osr_tutorial_proj Defining a Projected Coordinate System 110 111A projected coordinate system (such as UTM, Lambert Conformal Conic, etc) 112requires and underlying geographic coordinate system as well as a definition 113for the projection transform used to translate between linear positions 114(in meters or feet) and angular long/lat positions. The following code 115defines a UTM zone 17 projected coordinate system with an underlying 116geographic coordinate system (datum) of WGS84. 117 118\code 119 OGRSpatialReference oSRS; 120 121 oSRS.SetProjCS( "UTM 17 (WGS84) in northern hemisphere." ); 122 oSRS.SetWellKnownGeogCS( "WGS84" ); 123 oSRS.SetUTM( 17, TRUE ); 124\endcode 125 126Calling SetProjCS() sets a user 127name for the projected coordinate system and establishes that the system 128is projected. The SetWellKnownGeogCS() associates a geographic coordinate 129system, and the SetUTM() call sets detailed projection transformation 130parameters. At this time the above order is important in order to 131create a valid definition, but in the future the object will automatically 132reorder the internal representation as needed to remain valid. For now 133<b>be careful of the order of steps defining an OGRSpatialReference!</b> 134 135The above definition would give a WKT version that looks something like 136the following. Note that the UTM 17 was expanded into the details 137transverse mercator definition of the UTM zone. 138 139<pre> 140PROJCS["UTM 17 (WGS84) in northern hemisphere.", 141 GEOGCS["WGS 84", 142 DATUM["WGS_1984", 143 SPHEROID["WGS 84",6378137,298.257223563, 144 AUTHORITY["EPSG",7030]], 145 TOWGS84[0,0,0,0,0,0,0], 146 AUTHORITY["EPSG",6326]], 147 PRIMEM["Greenwich",0,AUTHORITY["EPSG",8901]], 148 UNIT["DMSH",0.0174532925199433,AUTHORITY["EPSG",9108]], 149 AXIS["Lat",NORTH], 150 AXIS["Long",EAST], 151 AUTHORITY["EPSG",4326]], 152 PROJECTION["Transverse_Mercator"], 153 PARAMETER["latitude_of_origin",0], 154 PARAMETER["central_meridian",-81], 155 PARAMETER["scale_factor",0.9996], 156 PARAMETER["false_easting",500000], 157 PARAMETER["false_northing",0]] 158</pre> 159 160There are methods for many projection methods including SetTM() (Transverse 161Mercator), SetLCC() (Lambert Conformal Conic), and SetMercator(). 162 163\section osr_tutorial_query Querying Coordinate System 164 165Once an OGRSpatialReference has been established, various information about 166it can be queried. It can be established if it is a projected or 167geographic coordinate system using the OGRSpatialReference::IsProjected() and 168OGRSpatialReference::IsGeographic() methods. The 169OGRSpatialReference::GetSemiMajor(), OGRSpatialReference::GetSemiMinor() and 170OGRSpatialReference::GetInvFlattening() methods can be used to get 171information about the spheroid. The OGRSpatialReference::GetAttrValue() 172method can be used to get the PROJCS, GEOGCS, DATUM, SPHEROID, and PROJECTION 173names strings. The OGRSpatialReference::GetProjParm() method can be used to 174get the projection parameters. The OGRSpatialReference::GetLinearUnits() 175method can be used to fetch the linear units type, and translation to meters. 176 177The following code (from ogr_srs_proj4.cpp) demonstrates use 178of GetAttrValue() to get the projection, and GetProjParm() to get projection 179parameters. The GetAttrValue() method searches for the first "value" 180node associated with the named entry in the WKT text representation. 181The #define'ed constants for projection parameters (such as 182SRS_PP_CENTRAL_MERIDIAN) should be used when fetching projection parameter 183with GetProjParm(). The code for the Set methods of the various projections 184in ogrspatialreference.cpp can be consulted to find which parameters apply to 185which projections. 186 187\code 188 const char *pszProjection = poSRS->GetAttrValue("PROJECTION"); 189 190 if( pszProjection == NULL ) 191 { 192 if( poSRS->IsGeographic() ) 193 sprintf( szProj4+strlen(szProj4), "+proj=longlat " ); 194 else 195 sprintf( szProj4+strlen(szProj4), "unknown " ); 196 } 197 else if( EQUAL(pszProjection,SRS_PT_CYLINDRICAL_EQUAL_AREA) ) 198 { 199 sprintf( szProj4+strlen(szProj4), 200 "+proj=cea +lon_0=%.9f +lat_ts=%.9f +x_0=%.3f +y_0=%.3f ", 201 poSRS->GetProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0), 202 poSRS->GetProjParm(SRS_PP_STANDARD_PARALLEL_1,0.0), 203 poSRS->GetProjParm(SRS_PP_FALSE_EASTING,0.0), 204 poSRS->GetProjParm(SRS_PP_FALSE_NORTHING,0.0) ); 205 } 206 ... 207\endcode 208 209\section osr_tutorial_transform Coordinate Transformation 210 211The OGRCoordinateTransformation class is used for translating positions 212between different coordinate systems. New transformation objects are 213created using OGRCreateCoordinateTransformation(), and then the 214OGRCoordinateTransformation::Transform() method can be used to convert 215points between coordinate systems. 216 217\code 218 OGRSpatialReference oSourceSRS, oTargetSRS; 219 OGRCoordinateTransformation *poCT; 220 double x, y; 221 222 oSourceSRS.importFromEPSG( atoi(papszArgv[i+1]) ); 223 oTargetSRS.importFromEPSG( atoi(papszArgv[i+2]) ); 224 225 poCT = OGRCreateCoordinateTransformation( &oSourceSRS, 226 &oTargetSRS ); 227 x = atof( papszArgv[i+3] ); 228 y = atof( papszArgv[i+4] ); 229 230 if( poCT == NULL || !poCT->Transform( 1, &x, &y ) ) 231 printf( "Transformation failed.\n" ); 232 else 233 printf( "(%f,%f) -> (%f,%f)\n", 234 atof( papszArgv[i+3] ), 235 atof( papszArgv[i+4] ), 236 x, y ); 237\endcode 238 239There are a couple of points at which transformations can 240fail. First, OGRCreateCoordinateTransformation() may fail, 241generally because the internals recognise that no transformation 242between the indicated systems can be established. This might 243be due to use of a projection not supported by the internal 244PROJ.4 library, differing datums for which no relationship 245is known, or one of the coordinate systems being inadequately 246defined. If OGRCreateCoordinateTransformation() fails it will 247return a NULL. 248 249The OGRCoordinateTransformation::Transform() method itself can 250also fail. This may be as a delayed result of one of the above 251problems, or as a result of an operation being numerically 252undefined for one or more of the passed in points. The 253Transform() function will return TRUE on success, or FALSE 254if any of the points fail to transform. The point array is 255left in an indeterminate state on error. 256 257Though not shown above, the coordinate transformation service can 258take 3D points, and will adjust elevations for elevation differents 259in spheroids, and datums. At some point in the future shifts 260between different vertical datums may also be applied. If no Z is 261passed, it is assume that the point is on the geoide. 262 263The following example shows how to conveniently create a lat/long coordinate 264system using the same geographic coordinate system as a projected coordinate 265system, and using that to transform between projected coordinates and 266lat/long. 267 268\code 269 OGRSpatialReference oUTM, *poLatLong; 270 OGRCoordinateTransformation *poTransform; 271 272 oUTM.SetProjCS("UTM 17 / WGS84"); 273 oUTM.SetWellKnownGeogCS( "WGS84" ); 274 oUTM.SetUTM( 17 ); 275 276 poLatLong = oUTM.CloneGeogCS(); 277 278 poTransform = OGRCreateCoordinateTransformation( &oUTM, poLatLong ); 279 if( poTransform == NULL ) 280 { 281 ... 282 } 283 284 ... 285 286 if( !poTransform->Transform( nPoints, x, y, z ) ) 287 ... 288\endcode 289 290\section osr_tutorial_apis Alternate Interfaces 291 292A C interface to the coordinate system services is defined in 293ogr_srs_api.h, and Python bindings are available via the osr.py module. 294Methods are close analogs of the C++ methods but C and Python bindings 295are missing for some C++ methods. 296 297<h3>C Bindings</h3> 298 299\code 300typedef void *OGRSpatialReferenceH; 301typedef void *OGRCoordinateTransformationH; 302 303OGRSpatialReferenceH OSRNewSpatialReference( const char * ); 304void OSRDestroySpatialReference( OGRSpatialReferenceH ); 305 306int OSRReference( OGRSpatialReferenceH ); 307int OSRDereference( OGRSpatialReferenceH ); 308 309OGRErr OSRImportFromEPSG( OGRSpatialReferenceH, int ); 310OGRErr OSRImportFromWkt( OGRSpatialReferenceH, char ** ); 311OGRErr OSRExportToWkt( OGRSpatialReferenceH, char ** ); 312 313OGRErr OSRSetAttrValue( OGRSpatialReferenceH hSRS, const char * pszNodePath, 314 const char * pszNewNodeValue ); 315const char *OSRGetAttrValue( OGRSpatialReferenceH hSRS, 316 const char * pszName, int iChild); 317 318OGRErr OSRSetLinearUnits( OGRSpatialReferenceH, const char *, double ); 319double OSRGetLinearUnits( OGRSpatialReferenceH, char ** ); 320 321int OSRIsGeographic( OGRSpatialReferenceH ); 322int OSRIsProjected( OGRSpatialReferenceH ); 323int OSRIsSameGeogCS( OGRSpatialReferenceH, OGRSpatialReferenceH ); 324int OSRIsSame( OGRSpatialReferenceH, OGRSpatialReferenceH ); 325 326OGRErr OSRSetProjCS( OGRSpatialReferenceH hSRS, const char * pszName ); 327OGRErr OSRSetWellKnownGeogCS( OGRSpatialReferenceH hSRS, 328 const char * pszName ); 329 330OGRErr OSRSetGeogCS( OGRSpatialReferenceH hSRS, 331 const char * pszGeogName, 332 const char * pszDatumName, 333 const char * pszEllipsoidName, 334 double dfSemiMajor, double dfInvFlattening, 335 const char * pszPMName , 336 double dfPMOffset , 337 const char * pszUnits, 338 double dfConvertToRadians ); 339 340double OSRGetSemiMajor( OGRSpatialReferenceH, OGRErr * ); 341double OSRGetSemiMinor( OGRSpatialReferenceH, OGRErr * ); 342double OSRGetInvFlattening( OGRSpatialReferenceH, OGRErr * ); 343 344OGRErr OSRSetAuthority( OGRSpatialReferenceH hSRS, 345 const char * pszTargetKey, 346 const char * pszAuthority, 347 int nCode ); 348OGRErr OSRSetProjParm( OGRSpatialReferenceH, const char *, double ); 349double OSRGetProjParm( OGRSpatialReferenceH hSRS, 350 const char * pszParmName, 351 double dfDefault, 352 OGRErr * ); 353 354OGRErr OSRSetUTM( OGRSpatialReferenceH hSRS, int nZone, int bNorth ); 355int OSRGetUTMZone( OGRSpatialReferenceH hSRS, int *pbNorth ); 356 357OGRCoordinateTransformationH 358OCTNewCoordinateTransformation( OGRSpatialReferenceH hSourceSRS, 359 OGRSpatialReferenceH hTargetSRS ); 360void OCTDestroyCoordinateTransformation( OGRCoordinateTransformationH ); 361 362int OCTTransform( OGRCoordinateTransformationH hCT, 363 int nCount, double *x, double *y, double *z ); 364\endcode 365 366<h3>Python Bindings</h3> 367 368\code 369class osr.SpatialReference 370 def __init__(self,obj=None): 371 def ImportFromWkt( self, wkt ): 372 def ExportToWkt(self): 373 def ImportFromEPSG(self,code): 374 def IsGeographic(self): 375 def IsProjected(self): 376 def GetAttrValue(self, name, child = 0): 377 def SetAttrValue(self, name, value): 378 def SetWellKnownGeogCS(self, name): 379 def SetProjCS(self, name = "unnamed" ): 380 def IsSameGeogCS(self, other): 381 def IsSame(self, other): 382 def SetLinearUnits(self, units_name, to_meters ): 383 def SetUTM(self, zone, is_north = 1): 384 385class CoordinateTransformation: 386 def __init__(self,source,target): 387 def TransformPoint(self, x, y, z = 0): 388 def TransformPoints(self, points): 389\endcode 390 391\section osr_tutorial_impl Internal Implementation 392 393The OGRCoordinateTransformation service is implemented on top of the 394<a href="http://www.remotesensing.org/proj">PROJ.4</a> library originally 395written by Gerald Evenden of the USGS. 396 397*/ 398