Improving raster reprojection from EPSG:4326 to UTM/WGS84

[Even Rouault]

Hi,

I've enclosed a patch (ticket http://trac.osgeo.org/mapserver/ticket/2324) 
that dramatically speeds up WMS/WCS requests made with a UTM/WGS84 CRS 
(EPSG:32601 to EPSG:32660; EPSG:32701 to EPSG:32760), on raster layers whose 
data is georeferenced in LatLong WGS84 (EPSG:4326). The speed up is a factor 
typically ranging from 5 to 7.
For example, a WCS request on a DTED1 layer to UTM 31N that took 7.9s takes 
just 1.2s.

How does it work ?
I've profiled mapserver and observed that 99% of the CPU time is spent at 
doing pj_transform called from msProjTransformer. In that context, we 
reproject the coordinates from pixel centers of a whole line, expressed in 
the destination CRS (UTM) to the source CRS (WGS84). This reprojection is 
particularly slow (from a quick look to proj4, the order of complexity seems 
to be a hundred of floating point operations, involving a lot of sin/cos, 
etc...).
The idea of the optimization is that we can use a function interpolator to 
compute most of the reprojected coordinates as we have functions f and g, 
such that latitude = f(easting) and longitude = g(easting).
I've used a Chebychev polynomial interpolator as it is quite easy to code (I'm 
using a similar API to the GSL library. But I've written it from scratch, so 
it's licence is X/MIT) and gives very good results. Thanks to a test program, 
I've computed the maximum errors between the exact proj4 reprojections and 
the interpolated reprojections, for different orders of the Chebichev series. 
So I can adjust the order of the Chebichev series by looking at the spacing 
of the pixels (for example, if the spacing of the pixel is 100m, and error of 
1m is OK in this context of use.)

The conditions that must be met for the optimization to be used are :
* dest CRS is EPSG:32601 to EPSG:32660; EPSG:32701 to EPSG:32760
* source CRS is EPSG:4326
* at least 50 points to reproject at the same time (50 is quite arbitrary. It 
must be big enough to overcome the initial costs of computing the Chebychev 
interpolation coefficients)
* the set of coordinates to reproject have the same northing
* the northing is < 8.000.000 in the northern hemisphere (or > 2.000.000 in 
the souther hemisphre)
* all the eastings are between 0 and 1.000.000

In these conditions,
* with order 4 chebyshev interpolation :
   the maximum error in longitude is 2.37 arc second ~ 71m
   the maximum error in latitude is 0.058 arc second
* with order 5 chebyshev interpolation :
   the maximum error in longitude is 0.056 arc second ~ 1.5m
   the maximum error in latitude is 0.030 arc second
* with order 7 chebyshev interpolation :
   the maximum error in longitude is 3e-11 arc second
   the maximum error in latitude is 3e-4 arc second ~ 0.009 m

(The maximum errors have been computed by interating over eastings from 0 to 
1.000.000 with a step of 1.000, and over northings from 0 to 8.000.000 with a 
step of 1.000)


Of course, this could be extended to other combinations of source and 
destination CRS, but for each combination, one must study carefully the 
behaviour of the interpolator and its maximum error.


I've also thought to put the optimization directly in proj4 by adding a new 
function like :
pj_transform_with_error_control(
   projPJ src, projPJ dst, long point_count,
   int point_offset, double *x, double *y, double *z, double errorThreshold);

But I've the feeling it's maybe a bit too specific to be put in proj4.

Any comments on all this ?

Best regards,
Even