[PATCH] Improving raster reprojection from EPSG:4326 to UTM/WGS84

[rouault]

I'm enclosing a patch 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.

[warmerdam]

Even,

I'm very surprised you see 99% of time in pj_transform(). mapresample.c already does uses a linear approximating transformer, so for typical scanlines only three points are actually transformed.

Before we replace the existing mechanisms, could you perhaps isolate an example where you see this dramatic speed issue? Something I could do a bit of testing with?

Steve - if you wouldn't mind I'd be happy to take this report.

[rouault]

Frank, I see now that my patch is unnecessary with mapserver 5.0.0...

I'm sticking with mapserver 4.10.2 for daily use and have written the patch for it. I applied it blindly on mapserver trunk and it applied cleanly (just line shift), so I thought that nothing had changed in the middle time... and didn't bother to check if it still improved things really... But now I see that there was bug #1944 reported and fixed in r5805 and it's exactly the problem I've been hitting with mapserver 4.10.2. So my patch is essentially useless now. Furthermore your fix is much more simple and general than mine.

However, in ticket #1944, it's said that the fix has been applied in 4.10 branch, but I can't see it. Are you sure that it has been really applied in it ?

Anyway, I'm closing the ticket. Sorry for that. I should have really tried mapserver 5 betas >[]