Polar reprojection (again)

[Frank Warmerdam]

Christopher Condit wrote:
> Hi Frank-
> I'm attaching a zip file with two simple mapfiles illustrating the
> problem.  There's a README.txt in zipfile that describes everything.
> Thanks for your help!

Chris,

Gee, I hope that big attachment didn't actually go to the list!

I dug into the problem and it turned out to be unusual.

The raster resampler calls PROJ.4's pj_transform() function for one
whole row of points to transform at once.  The idea is that pj_transform()
returns a success flag for each point.  But some kinds of errors are
considered to be "non-transient", and it is assumed if they are encountered
for one point they will apply to all points.

It turns out that getting errno=33 (EDOM) or errno=34 (ERANGE) from
math library functions (like atan) were considered "non-transient" by
pj_transform(), even though they are likely to only apply to some points
in a group.  Because your initial range +-10000000 was invalid to reproject,
and because the error was improperly considered non-transient the code
was assuming all the points would fail and giving up.

Using a smaller initial EXTENT (ie +-5000000) was enough to get a reason
result.  But better than that, I was able to correct the PROJ.4
pj_transform() function.

The upshot is that if you can grab the latest PROJ.4 source from CVS and
build it, you can get correct behavior.  Since I wouldn't be surprised
if that is messy for you, I have also put the modified pj_transform.c
on my web site at:

   http://www.gdal.org/tmp/pj_transform.c

If you grab that and dump over the PROJ 4.4.9 src/pj_transform.c file,
build and install it, you should be ok.

Thanks for providing a concise demonstration of the problem so I was
able to track it down fairly easily.

Best regards,
-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Programmer for Rent