I have spent 3 days on this issue, and it appears more difficult than I initially though. The problem is that Bursa-Wolf parameters for datum changes are encoded into a WGS84ConversionInfo object. Unfortunately, this object (part of OGC 2001-09 specification) express only datum change from an arbitrary CS to WGS84, not a change from an arbitrary CS to an other arbitrary CS. At a first look, it doesn't looks like a problem since we can turn any transformation like 4817 --> 4273 --> WGS84 to a more direct transformation 4817 --> WGS84. All we need to do is to express the 4817 --> 4273 transformation as a matrix, express the 4273 --> WGS84 transformation as a matrix too and multiply them together. Expressing a WGS84ConversionInfo as a matrix is easy:

http://modules.geotools.org/cts-coordtrans/apidocs/org/geotools/cs/WGS84ConversionInfo.html#getAffineTransform()
Unfortunately, if we take two matrix like that and multiply them together, the result is more complex. For example the diagonal in the 3x3 upper left part does not contains anymore identical values (the "S" constant in the above-cited link), and is not anymore symmetric with opposite sign, as the WGS84ConversionInfo object assumes. Consequently, the result of the matrix multiplication does not fit into a WGS84ConversionInfo object.

Even if the result of concatenating two ConversionInfo object does not fit perfectly in a new WGS84ConversionInfo object, it is possible to make it fit approximatively. We can probably fit a WGS84ConversionInfo in the least-square sense (like a line fitting in a scatter plot). I have been unable to resolve fully the equations (an analytic software like Maple or Mathematica would help me; I just do too much mistakes when the equations are longer than a few terms). But for the part I have been able to solve, its look like a feasible approach. The reason is that the matrix multiplication produces expressions of the kind x + x² where x is an angle in radians, and all "x" terms fit perfectly in WGS84ConversionInfo. The only terms that don't fit are "x²" (or "xy", or "yz", etc.). Those terms can be safely ignored if |x| <<< 1. In Bursa-Wolf parameters, angles ("x") are typically less than 1E-4 radians. Not sure if it is small enough, but it may be.

So the choices are:

1) Go ahead with a naïve approximation in WGS84ConversionInfo as described above. The transformations will be approximative. The errors will be smaller if the rotation angles are small, and null if there is no rotation.

2) Replace the WGS84ConversionInfo object by something more general, capable to target arbitrary datum rather than exclusively WGS84. Then update CoordinateTransformationFactory to use this information. This is the cleanest and most accurate way (since no approximation is performed after the matrix multiplication), but would requires more work than I can afford right now. It would be better handled during the process of refactoring cts-coordtrans to the latest ISO 19107 specifications.

Even 1 as a temporary hack will probably need at least one day to implement.

http://modules.geotools.org/cts-coordtrans/apidocs/org/geotools/cs/WGS84ConversionInfo.html#getAffineTransform()

Unfortunately, if we take two matrix like that and multiply them together, the result is more complex. For example the diagonal in the 3x3 upper left part does not contains anymore identical values (the "S" constant in the above-cited link), and is not anymore symmetric with opposite sign, as the WGS84ConversionInfo object assumes. Consequently, the result of the matrix multiplication does not fit into a WGS84ConversionInfo object.

Even if the result of concatenating two ConversionInfo object does not fit perfectly in a new WGS84ConversionInfo object, it is possible to make it fit approximatively. We can probably fit a WGS84ConversionInfo in the least-square sense (like a line fitting in a scatter plot). I have been unable to resolve fully the equations (an analytic software like Maple or Mathematica would help me; I just do too much mistakes when the equations are longer than a few terms). But for the part I have been able to solve, its look like a feasible approach. The reason is that the matrix multiplication produces expressions of the kind x + x² where x is an angle in radians, and all "x" terms fit perfectly in WGS84ConversionInfo. The only terms that don't fit are "x²" (or "xy", or "yz", etc.). Those terms can be safely ignored if |x| <<< 1. In Bursa-Wolf parameters, angles ("x") are typically less than 1E-4 radians. Not sure if it is small enough, but it may be.

So the choices are:

1) Go ahead with a naïve approximation in WGS84ConversionInfo as described above. The transformations will be approximative. The errors will be smaller if the rotation angles are small, and null if there is no rotation.

2) Replace the WGS84ConversionInfo object by something more general, capable to target arbitrary datum rather than exclusively WGS84. Then update CoordinateTransformationFactory to use this information. This is the cleanest and most accurate way (since no approximation is performed after the matrix multiplication), but would requires more work than I can afford right now. It would be better handled during the process of refactoring cts-coordtrans to the latest ISO 19107 specifications.

Even 1 as a temporary hack will probably need at least one day to implement.