In addition to the two-dimensional coordinate transformations described in lesson 7, the CoordTransform object is also capable of performing three-dimensional coordinate transformations. Such a transformation can be performed when both the source and target coordinate systems have a valid VerticalReference.
In order to compile the code samples provided in this lesson, it will be necessary to include the "WDataSource.h" header file.
Since all of the definitions we will be using reside in the DataSource, we will start by creating and loading a DataSource. For more information about the DataSource object, see lesson 4.
DataSource * dataSource = 0;
dataSource = DataSource::CreateDataSource();
dataSource->LoadFile(L"c:\\bluemarble\\common\\data\\geocalc.xml", false, L"c:\\bmg\\geocalcpbw\\data\\custom.xml");
Performing a three-dimensional coordinate transformation is very similar to performing a two-dimensional coordinate transformation. As before, we must start by constructing a CoordTransform object, which will facilitate this transformation. In order to construct a CoordTransform, we will need to retrieve some CoordSys definitions and some VerticalReferences from the DataSource. Since we want this to be a non-trivial lesson, let's pick some interesting objects. First, we want two different types of CoordSys objects, a ProjectedCoordSys and a GeodeticCoordSys. Second, we want to select CoordSys objects that use different HorizontalDatums so that we will need to use some DatumShifts. Finally, we should select two different VerticalReferences so that we need to perform a vertical transformation. These objects will suffice for this lesson:
GeodeticCoordSys * source = 0;
source = dataSource->GetGeodeticCoordSys(L"BMG", L"SIRGAS_GEODETIC");
VerticalReference * sourceVR = 0;
sourceVR = dataSource->GetVerticalReference(L"BMG", L"CVD_GEOCOL04");
source->set_VerticalReference(sourceVR);
ProjectedCoordSys * target = 0;
target = dataSource->GetProjectedCoordSys(L"BMG", L"UTM-18N");
VerticalReference * targetVR = 0;
targetVR = dataSource->GetVerticalReference(L"BMG", L"ELLIPSOID_HEIGHT");
target->get_InnerGeodetic().set_VerticalReference(targetVR);
Here we have selected the Sirgas GeodeticCoordSys with the Colombia Geoid Model of 2004 VerticalReference as our source coordinate system, and the UTM 18N ProjectedCoordSys with the ellipsoid height VerticalReference as our target coordinate system. So we are now ready to construct our CoordTransform:
CoordTransform * ct = 0;
ct = CoordTransform::CreateCoordTransform(*source, *target);
Since our source and target coordinate systems use different HorizontalDatums, we will need to perform one or more DatumShifts. There are multiple ways to find the DatumShifts needed for a CoordTransform, and Lesson 6 describes these methods. For the purpose of this lesson, we will just retrieve a single DatumShift from the DataSource and put it into our CoordTransform::DatumShifts collection:
DatumShift * ds = 0;
ds = dataSource->GetDatumShift(L"BMG", L"SIRGAS_to_WGS84");
ct->get_DatumShifts().AddBack(*ds);
It is also possible to intentionally not use a DatumShift by not populating the DatumShifts property. Users should be aware that some coordinate shift may occur when shifting between ProjectedCoordSys' even though no DatumShift is being used. This can occur if the source and target HorizontalDatums use different Ellipsoids.
Since our source and target coordinate systems use different VerticalReferences, we will need to perform a vertical transformation. This requires that we call the CoordTransform::InitializeVerticalTransform method before we try to transform any points with this CoordTransform.
Most of the HeightModels and VerticalDatums supported by GeoCalc require additional data files that were not installed with GeoCalc. All of these files can be obtained free of charge from our website. A description of the files needed by each HeightModel and VerticalDatum, as well as a link through which they can be downloaded, can be found on the description page for each HeightModel and each VerticalDatum. For the transformation described in this lesson, we will require the files for the Colombia Geoid Model of 2004. Without these files, the InitializeVerticalTransform method will throw a GeoCalcException.
Now that we have the necessary files, we can call the InitializeVerticalTransform method:
try
{
ct->InitializeVerticalTransform();
}
catch(GeoCalcException & ex)
{
if(ex.get_ErrorCode() == GeoCalcException::Code::FileNotFound)
{
// Then we do not have all of the files needed by our VerticalReferences
return;
}
else if(ex.get_ErrorCode() == GeoCalcException::Code::IllegalVerticalTransform)
{
// Then we were unable to initialize the vertical transform.
return;
}
}
The InitializeVerticalTransform method may throw a GeoCalcException with an ErrorCode of IllegalVerticalTransform. This indicates that GeoCalc is not capable of performing the requested vertical transformation. There could be several reasons for this, but the most common cause of this exception is that one of our VerticalReferences uses a HeightModel that uses a HorizontalDatum that is not used by either the source or target CoordSys. However, if you have set things up in the same way as is done in this lesson, you should not receive this exception (nor should you receive any GeoCalcException for that matter).
Our CoordTransform is now set up and ready to perform a three dimensional coordinate transformation. As with the two-dimensional transformation as described in lesson 7, we have the option of transforming either a single point or a list of points. Here we will only transform a single point. If you would like to transform many points at once, please see the section 'Transforming a Collection of Points' in lesson 7. The method of transforming multiple three-dimensional points is the same as the method for transforming multiple two-dimensional points.
We must take care when selecting a point to transform, because the Colombia Geoid Model of 2004 is only defined from 5 to 15 north latitude and from 66 to 80 west longitude. If we try to transform a point outside of that area, we will receive a GeoCalcException. Here is the GeodeticPoint we will transform, along with the expected value of the resulting ProjectedPoint:
|
|
Sirgas Geodetic |
|
|
UTM 18N |
|
Latitude |
10 Degrees |
|
North |
1105578.58919244 Meters |
|
Longitude |
-76 Degrees |
|
East |
390399.22748554 Meters |
|
Height |
123 Meters |
|
Height |
120.75025305 Meters |
Here is how we perform the transformation:
CoordPoint * inPt = 0;
inPt = ct->get_SourceCoordSys().get_PointStyle().CloneCoordPoint();
CoordPoint * outPt = 0;
outPt = ct->get_TargetCoordSys().get_PointStyle().CloneCoordPoint();
inPt->set_InUnits(-76, 10, 123);
try
{
if(! ct->ConvertAhead(*inPt, *outPt))
{
// ConvertAhead failed
return;
}
}
catch(GeoCalcException & ex)
{
if(ex.get_ErrorCode() == GeoCalcException::Code::InPointWrongType)
{
// Then inPt has the wrong ClassType
return;
}
else if(ex.get_ErrorCode() == GeoCalcException::Code::OutPointWrongType)
{
// Then outPt has the wrong ClassType
return;
}
else if(ex.get_ErrorCode() == GeoCalcException::Code::InPointOutOfBounds)
{
// Then inPt is not a valid point within the SourceCoordSys
return;
}
else if(ex.get_ErrorCode() == GeoCalcException::Code::OutPointOutOfBounds)
{
// Then outPt is not a valid point within the TargetCoordSsy
return;
}
}
double east, north, height;
outPt->get_InUnits(east, north, height);
The values of east, north, and height should now match the corresponding values shown above.
Finally, we must clean up after ourselves and free the memory that we have allocated in this lesson using the Disposal::Dispose method:
if(ct) Disposal::Dispose(ct);
if(source) Disposal::Dispose(source);
if(sourceVR) Disposal::Dispose(sourceVR);
if(target) Disposal::Dispose(target);
if(targetVR) Disposal::Dispose(targetVR);
if(ds) Disposal::Dispose(ds);
if(inPt) Disposal::Dispose(inPt);
if(outPt) Disposal::Dispose(outPt);
if(dataSource) Disposal::Dispose(dataSource);