Developers Corner: Improving GeoTools/GeoServer raster reprojection performance

Ciao a tutti,
we hope everybody had nice seasonal holidays and is back fresh and ready to move on with another year of activity.

In this blog we’d like to present some of the work we have been doing in recent time to improve GeoTools and GeoServer raster reprojection abilities.

Raster reprojection is a quite heavy process in which every single pixel of the original image has to be mapped into a new position in the target raster. Here is a visual representation of a small set of pixels, before and after the reprojection (in this case, from WGS84 to polar stereographic):

A point by point reprojection is obviously taking quite a toll, a screen sized image has easily well over one million pixels to reproject, much more than the number of vertices one can find in the common vector map. To make things worse, the transformation to be used often involves trigonometric functions and other complex mathematical tools that visibly slow down the calculation.

The reprojection algorithm used by GeoServer and GeoTools until now worked exactly like that, it used a very accurate but also very slow approach.

Now, if you think about it, the transformation applied to one pixel and the one applied to the next one are basically the same, the difference is going to be minimal. If we only had a small bunch of pixels we could compute a simple translation and apply it to them all.
Having more pixels we could try to approximate the overall transformation with a simple linear one, which in two dimensional space is known as an affine transformation.

If the area gets larger a single affine transformation would not be good enough anymore, but we could use many of them to fit different areas of the raster.

Have a look at the following one dimensional representation to get the gist of the idea (this is known as piecewise linear function):

The trick to apply this approach if finding out where and how much to subdivide the original area into discrete pieces that can use a single affine transform: to do so we devised a simple error estimation method that proved to work well with all common map projections.

So we start by evaluating the linear transformation error over the entire raster area, if it’s larger than a certain target we split the area and evaluate again, and recurse this way until the approximation error has been put under control.

The optimized method proved to be quite effective, drastically reducing the time it takes to reproject a raster. We made some tests with GeoServer using the same data and the same reprojection test used in the FOSS4G 2010 WMS shootout measuring a six times speedup:

We also compared the images generated by the original and optimized method finding there was hardly any human noticeable difference, and very little difference that could be machine detected at all.

For all of you that want to dig into the nitty-gritty details, we have a full tech report that you can inspect.

There is still plenty that can be done to optimize further raster data reading and in memory transformations to reach even higher performance and scalability. Interested? Let us know!

The GeoSolutions team

 
  • Right during holidays I was wondering if the “approximate transformer” approach in Mapserver [1] would ever be implemented in Geotools/Geoserver… Thanks!!!

    [1] Mapserver/mapresample.c

  • He he, yeah, it’s a similar reasoning, though they do reprojection line by line and optimize that one line at a time, whilst we build up front an optimal reprojection grid instead and then pass it onto JAI.

    Besides that, the two approach are indeed similar.

  • I thought if you project source image pixels into the new coordinate system you can end up with all sorts of holes and artefacts. The right way to do it is to apply the reverse projection to the pixels in the bounding box of the reverse-projected area. Hence each destination pixel is going to have a value, and no holes. Interpolation can also be used.

  • Barry, yes, that’s exactly what’s happening, the image warping is generated by back mapping pixels from target back onto the source image (single pixel or a nxn group, depending on the interpolation chosen). So no holes.

  • Anonymous

    Would using CUDA or OpenCL be an additional way of speeding up reprojections?

  • About CUDA, we are investigating using OpenCL rather than CUDA for superior portability. Anyway at the moment the one concern we have about that approach is the applicability in a server side environment which is core for us and where we are afraid benefits could be limited (scalability VS pure speed).

  • Simone, I am rather interested speeding up things using OpenCL. Are there any efforts going on and is there a place to look at?

  • Ciao Herbert,
    I did some sparse testing in my spare time but for the moment we did not find anyone interested in in funding this development.

    Send me an email I would be happy to discuss this topic a bit more.