Introduction to Look-up tables
Back in the DOS era and 80×86 CPU (I guess it’s quite old), mathematical operations for 3D Graphics were too expensive. One of the tricks to speed-up the graphics pipeline was using a pre-computed trigonometric functions implemented in a table. The benefit of using discretization of function parameter is measured by how complex the computation of function relative to accessing a table using a discretized input parameter. This method of replacing function call into looking up a table for already computed values are generally called Memoization.
Generalizing Look-up table for higher dimension
In previous example(trigonometric function) the input parameter is only one. If we think the parameter as a dimension, then the domain of the input parameter is a one dimensional object (a line). Of course the memoization technique is not only limited to one dimensional look-up table but can also be generalized to higher dimensions. In 3D Graphics (again), The memoization for higher dimension is used in texture mapping (2D) or voxel(volume pixel) representation for object modeling.
Grid shift: LUT in Coordinate transformation
Now let’s move from 3D Graphics to another area. In Geospatial domain, the same technique (look-up table) is used to transform coordinate from one Coordinate Reference System (CRS) to another CRS. Since most geographical coordinates only has two horizontal components (latitude and longitude). The look-up table is in the form of a two-dimensional grid. This grid-based transformation is also known as ‘Grid Shifts’.
- Current implementation status in GDAL/OGR (PROJ.4)
in main branch of GDAL/OGR the core coordinate transformation is done through delegation to PROJ.4 library. The grid-shift file handling and application is done internally by PROJ.4 library. Input WKT file is translated to PROJ.4 string which then executed by PROJ.4 library. There are two grid-shift that can be used: One is grid-shift in datum transform which is located in pj_datum_transform function that is called by pj_transform function. This grid-shift uses .gsb (grid shift binary, NTv1 and NTv2) file format and is used only to transform horizontal components. The other one is vertical grid shift which is inside the pj_transform function before and after the datum transform. The vertical grid-shift uses .gtx file format from NOAA and only transform vertical component (if available).
- Proposed workflow in Google Summer of Code Project (2013-2014)
Last year’s (2013) main focus of GSoC is to modify the vertical grid-shift transform to use formats that is supported as GDAL raster. The other task is also to make vertical grid-shift supports different kinds of vertical height model (projected/in-use, orthometric, ellipsoidal) via a set of grid-shifts (geoid undulation, vertical correction model) and a constant scalar parameter.
This year’s (2014) main focus is to modify grid-shift to be used not for inter-datum transform but intra-datum homogeneity correction.
Challenges in Implementation
- Analyzing the .gsb file format
During the time working for studying the grid shift. There is a question about whether the current format (NTv2) is already designed to support 3D transformation (deltaLon, deltaLat, deltaH). After examining the documentation (NTv2 Developer Guide), it is stated that the record length in NTv2 .gsb file is 16 bytes which contain 4 float32 elements (Latitude_shift, Longitude_shift, Latitude_accuracy, and Longitude accuracy). Since most of the application never uses the last two accuracy element (in meters), there is possibility using one of the last two (or both) unused record element to store Height_shift.
my conclusion regarding NTv2 is that _NTv2 grid shift_ is not yet capable of dealing with height correction grids
Removing the grid-shift transform from datum transform and adding grid-shift transform to pj_transform
Last year’s GSoC is simpler. because in order to implement the vertical grid-shift correction inside GDAL there is only required to making a modified duplicate of pj_transform function from PROJ.4 to GDAL/OGR. To maintain compatibility with the GDAL, the contributed code is moved into a new library (SpatialRef3D). The datum transform still uses the code inside PROJ.4.
However, for this year’s GSoC there are some issues. Since the PROJ.4 is used by GDAL/OGR via dynamic-linking (.dll), when pj_datum_transform is duplicated inside GDAL/OGR codebase and modified, The dependent functions is not available (not exported) when linking and some structures (LP) has different selectors between interface that is declared in the internal function and in the exported functions.
This causes more functions and macro declarations to be duplicated (including the functions that dealing directly with grid-shift files). The result is, there are more code ‘redundancies’ between the SpatialRef3D(GDAL/OGR) and the PROJ.4 library. The last commit that I pushed contains these additional codes. It’s a little bit cluttered but I’ll fixed those after mid-term evaluation.
I didn’t realized it’s already mid-term already. 😀 Let’s do better in the next term!