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