In this post I will discuss transforming data sets from one datum to another, and focus more specifically on the change of the Australian datum from the Australian Geodetic Datum (AGD) to the Geocentric Datum of Australia (GDA).
Datum Transformations
In my previous post I mentioned that SQL Server uses an SRID to define which datum an instance of a geographic data type belongs to, and this allows SQL Server to expose spatial properties such as STLength and STArea. But something that is missing from SQL Server is that it doesn't know the differences between one datum and another, so it cannot transform data from one datum into another. This means that data from one datum cannot be used to perform spatial selections on data from another datum.
To facilitate this requirement the data needs to be in a common datum. The easiest way to achieve this is to create spatial layers that use the same datum. If your data layers use different datums, then you will need to transform the data from their current datum into the common datum.
If you have commercial GIS software such as ESRI's ArcMap, or SafeSoftware's FME products, it is a relatively simple operation to transform from one datum to another because the software knows the parameters required to transform from one datum to another, or one projection to another, but since I was just investigating spatial indexes at home, I was more interested in a cheaper option, so I started looking at ogr2ogr which is one of the utilities in FW Tools.
FW Tools is a set of open source GIS utilities and contains the ogr2ogr utility, which is part of the Geospatial Data Abstraction Library (GDAL), which is an open source library for manipulating spatial data.
ogr2ogr is a command line utility that can be used to convert data from one format to another and/or one projection or datum to another.
For instance, I could use the following command to convert a shapefile of data representing features in grid coordinates in Australian Map Grid (AMG) Zone 56, deriving from the Australian Geodetic Datum (AGD84) data to latitude longitude.
Note that, as discussed in my previous post on SQL Server Spatial Datums that the AMG projection is based on the AGD datum, so there is no datum transformation as such, just a conversion from projected coordinates to geographic coordinates.
If I wanted to convert between datums I would need to provide some transformation parameters. From what I can see ogr2ogr can accept a well known text (WKT) parameter called towgs84, which will convert the datum to the World Geocentric Datum (WGS84) but this requires knowledge of the parameters, and also means that the data can only be transformed to WGS84 using those settings. For my Australian data there is a special datum transformation I can perform that involves a transformation shift based on a distortion file, which I will discuss in further detail below.
I found an interesting "cheat sheet" post on ogr2ogr on the Boston GIS site.
To see more examples of the WKT syntax you can look at the epsg file in the proj_lib directory of the FW Tools installation area. This file will contain the definitions of a lot of the EPSG datums, so there will be examples of defining projections, and ellipsoids, as well as some of the names of ellipsoids recognised by ogr2ogr.
EPSG Datum Identifiers
As you saw above I used the EPSG switch to indicate that I wanted to transform from my projection to my datum. EPSG stands for the European Petroleum Survey Group, and the identifiers represent the same id's that SQL Server uses as the SRID to define the datum in which a geographic data feature is referenced.
You can find a listing of the EPSG datums on the spatialreference.org website, but being an Australian resident, some of the EPSG's that I am more interested in are:
AGD66 4202
If you have commercial GIS software such as ESRI's ArcMap, or SafeSoftware's FME products, it is a relatively simple operation to transform from one datum to another because the software knows the parameters required to transform from one datum to another, or one projection to another, but since I was just investigating spatial indexes at home, I was more interested in a cheaper option, so I started looking at ogr2ogr which is one of the utilities in FW Tools.
FW Tools is a set of open source GIS utilities and contains the ogr2ogr utility, which is part of the Geospatial Data Abstraction Library (GDAL), which is an open source library for manipulating spatial data.
ogr2ogr is a command line utility that can be used to convert data from one format to another and/or one projection or datum to another.
For instance, I could use the following command to convert a shapefile of data representing features in grid coordinates in Australian Map Grid (AMG) Zone 56, deriving from the Australian Geodetic Datum (AGD84) data to latitude longitude.
ogr2ogr -f "ESRI Shapefile" -s_srs EPSG:20356 -t_srs EPSG:4203 CadastreAGD84.shp CadastreAMGZone56AGD84.shp
Note that, as discussed in my previous post on SQL Server Spatial Datums that the AMG projection is based on the AGD datum, so there is no datum transformation as such, just a conversion from projected coordinates to geographic coordinates.
If I wanted to convert between datums I would need to provide some transformation parameters. From what I can see ogr2ogr can accept a well known text (WKT) parameter called towgs84, which will convert the datum to the World Geocentric Datum (WGS84) but this requires knowledge of the parameters, and also means that the data can only be transformed to WGS84 using those settings. For my Australian data there is a special datum transformation I can perform that involves a transformation shift based on a distortion file, which I will discuss in further detail below.
I found an interesting "cheat sheet" post on ogr2ogr on the Boston GIS site.
To see more examples of the WKT syntax you can look at the epsg file in the proj_lib directory of the FW Tools installation area. This file will contain the definitions of a lot of the EPSG datums, so there will be examples of defining projections, and ellipsoids, as well as some of the names of ellipsoids recognised by ogr2ogr.
EPSG Datum Identifiers
As you saw above I used the EPSG switch to indicate that I wanted to transform from my projection to my datum. EPSG stands for the European Petroleum Survey Group, and the identifiers represent the same id's that SQL Server uses as the SRID to define the datum in which a geographic data feature is referenced.
You can find a listing of the EPSG datums on the spatialreference.org website, but being an Australian resident, some of the EPSG's that I am more interested in are:
AGD66 4202
AMG Zone 56 20256
AGD84 4203
AMG Zone 56 20356
GDA94 4283
MGA Zone 56 28356
WGS84 4326
Datum Transformations in Australia
In Australia there was a very important change in datums enacted in 2000. The Inter-gevernmental Committee on Surveying and Mapping (ICSM) made a decision to change the datum used for standard mapping to use a geocentric datum. One of the main reasons for this was to provide consistency with coordinates determined from GPS, which has become much more popular since the original datums were defined (i.e. AGD66 in 1966). The explosion in the number of casual map users with smartphones is another group of users that is benefiting from the change in the datum which was unforeseen even at the time the datums were changed.
I guess one of the possible questions is that if this happened 11 years ago, why do we need to care about transforming data from one datum to another? The reason for this is that it is quite conceivable that there are data sets in general population that are still using the older datums (particularly historical data), so we need to be able to convert these data sets to our new datums to be able to overlay them correctly, or perform spatial analysis.
Migrating from AGD to GDA
AGD84 4203
AMG Zone 56 20356
GDA94 4283
MGA Zone 56 28356
WGS84 4326
Datum Transformations in Australia
In Australia there was a very important change in datums enacted in 2000. The Inter-gevernmental Committee on Surveying and Mapping (ICSM) made a decision to change the datum used for standard mapping to use a geocentric datum. One of the main reasons for this was to provide consistency with coordinates determined from GPS, which has become much more popular since the original datums were defined (i.e. AGD66 in 1966). The explosion in the number of casual map users with smartphones is another group of users that is benefiting from the change in the datum which was unforeseen even at the time the datums were changed.
I guess one of the possible questions is that if this happened 11 years ago, why do we need to care about transforming data from one datum to another? The reason for this is that it is quite conceivable that there are data sets in general population that are still using the older datums (particularly historical data), so we need to be able to convert these data sets to our new datums to be able to overlay them correctly, or perform spatial analysis.
Migrating from AGD to GDA
As I mentioned above, 11 years ago the standard datum for Australian mapping changed from the Australian Geodetic Datum (AGD66 or AGD84) to the Geocentric Datum of Australia (GDA94). The AGD66 and AGD84 datums used the same spheroid, but had a slightly different adjustment of the control points in the Australian control network. Other than the parametric differences, such as the semi-major axis and flattening, the noticeable difference between the AGD spheroid and the GDA spheroid is that the origin of AGD is not geocentric with the Earth. The GDA spheroid is the same as the WGS84 spheroid which is used by GPS sattelites, which as I said, means that measurements from GPS are usable on the GDA datum without complex transformation.
If you are interested in further reading on GDA, check out the GDA Technical Manual produced by the ICSM
Just to add a little more confusion to the AGD/GDA acronyms is the change in terminology for projections of AGD from the Australian Map Grid (AMG) to the Map Grid of Australia (MGA) for GDA. Theses standard map grids are both Universal Transverse Mercator (UTM) projections.
The UTM projection is a variation of the Mercator projection. If you use Google Maps you would be familiar with the Mercator projection - it is as though you wrapped a piece of paper around the earth as a cylinder touching the equator. Different projections have different characteristics, where a particular geometric aspect is maintained (length, angles, area) at the expense of introducing some type of distortion. You can see in the image below from Google Maps that the Mercator projection shows Greenland to be larger than Australia, but we all know that Greenland is actually smaller than Australia, so it can be seen that the Mercator projection is distorting the data particularly at the poles. In this case the area is being distorted, while angles are maintained.
On the other hand, you can see that the Mollweide projection below has maintained the relative area of Greenland and Australia, but the angles have not been maintained, showing the data as sheared, or skewed.
Standard mapping projections are defined such that they minimise the distortion while retaining their useful geometric attributes. In the case of UTM, this represents a Mercator projection, except that rather than being defined around the equator, the line touching the surface of the spheroid is vertical, along meridians of longitude, hence the "transverse" terminology. The mapping system is then broken up into small zones to minimise the amount of distortion at the edges.
The AMG and MGA mapping systems are broken into 6 degree zones - see below.
So the zones that I defined in the ESPG section above relate to the Zones I care most about, living in the South East corner of Queensland.
Transformation of AGD to GDA using a Distortion Grid
As I discussed above, I was able to use ogr2ogr to transform projected grid coordinates into the equivalent geographic coordinates of its related datum, but I was not able to transform from AGD to GDA because I didn't have any transformation parameters. If you read through the GDA Technical Manual you will find that the recommended method (and highest accuracty method) for transforming data is to use a distortion grid to translate the coordinates from one datum to another. The ICSM has released some distortion grid data that can be used to perform this translation.
The cool thing I found was that I could define switches in the ogr2ogr command line utility to use these distortion grids to transform my data.
My data was transformed from AMG Zone 56 to AGD84, in the example at the top of the post, so I could then use the following command to transform the data to GDA94.
ogr2ogr -f "ESRI Shapefile" -s_srs "+proj=longlat +ellps=aust_SA +nadgrids=C:\datum\grid84.gsb +units=m +no_defs +wktext" -t_srs EPSG:4283 CadastreGDA94.shp CadastreAGD84.shp
If I wanted to transform directly from grid coordinate data in AMG Zone 56 (on AGD84) then I could use the following command to transform the data to GDA94.
ogr2ogr -f "ESRI Shapefile" -s_srs "+proj=utm +zone=56 +south +ellps=aust_SA +nadgrids=C:\datum\grid84.gsb +units=m +no_defs +wktext" -t_srs EPSG:4283 CadastreGDA94.shp CadastreAMGZone56AGD84.shp
Note the use of the UTM projection (proj=utm), in Zone 56 (zone=56) using the australian spheroid (ellps=aust_SA).
Conclusion
As I already mentioned, you never know when you might need to transform your data to another datum in order to overlay it, so you need to know the origin of the data and the datum and/or projection it belongs to. One of the tricky things about projected data is that is just a grid based coordinate system, so quite often it is not clear which projection it belongs to, or even which datum, unless there is adequate metadata describing this.
I found that using the ogr2ogr application was a little difficult to get detailed information about the WKT parameters it uses, but once I got it working I found it to transform the data correctly (I was lucky enough to have some test data I could compare), and was actually very fast to transform.
One of the important things to ensure is that the data is being transformed as you would expect, especially before running a batch job that transforms lots of layers of data. So make sure you can test the output position of some known points to verify the transformation.
No comments:
Post a Comment
Note: Only a member of this blog may post a comment.