[Proj] Natural Earth projection

Thu Mar 8 18:46:20 EST 2012

In Reply to:
Oscar van Vlijmen ovv at hetnet.nl
Mar 4, 2012, at 8:33 AM

> Could be my error, but as closer I get to a pole, the more inaccurate the 
> inverse gets.
> The error amounts to several degrees. Useless!
> 
> Examples:
> 
> R = 6371008.7714 m

…

> lat, lon = 89.0, 120.0
> -> x,y = 6604995.340672151, 9042363.459174621
> inverse, back to lat, lon: 86.95668, 113.64953

I'm unable to reproduce the problem. Below is the code I use (JavaScript from http://www.shadedrelief.com/NE_proj/code.html  embedded in a HTML file). I apologize for posting JavaScript and HTML on this list.

Bernie

<html>
<head>

<script type="text/javascript">
var MAX_Y = 0.8707 * 0.52 * Math.PI;

function forward(lon, lat, xy) {
	var lat2 = lat * lat;
	var lat4 = lat2 * lat2;
	xy[0] = lon * (0.8707 - 0.131979 * lat2 + lat4 * (-0.013791 + lat4 * (0.003971 * lat2 - 0.001529 * lat4)));
	xy[1] = lat * (1.007226 + lat2 * (0.015085 + lat4 * (-0.044475 + 0.028874 * lat2 - 0.005916 * lat4)));
};

function inverse(x, y, lonlat) {
	var yc, tol, y2, y4, f, fder;

	if(y > MAX_Y) {
		y = MAX_Y;
	} else if(y < -MAX_Y) {
		y = -MAX_Y;
	}

	// Newton's method
	yc = y;
	for(; ; ) {
		y2 = yc * yc;
		y4 = y2 * y2;
		f = (yc * (1.007226 + y2 * (0.015085 + y4 * (-0.044475 + 0.028874 * y2 - 0.005916 * y4)))) - y;
		fder = 1.007226 + y2 * (0.015085 * 3 + y4 * (-0.044475 * 7 + 0.028874 * 9 * y2 - 0.005916 * 11 * y4));
		yc -= tol = f / fder;
		if(Math.abs(tol) < 0.0000000001) {
			break;
		}
	}
	lonlat[1] = yc;
	y2 = yc * yc;
	lonlat[0] = x / (0.8707 + y2 * (-0.131979 + y2 * (-0.013791 + y2 * y2 * y2 * (0.003971 - 0.001529 * y2))));
};
</script>
</head>

<body>

<script type="text/javascript">
	var R = 6371008.7714;
	var lat = 89 / 180 * Math.PI, lon = 120 / 180 * Math.PI;
	var xy = [];

	forward (lon, lat, xy);
	document.write((xy[0] * R).toFixed(10) + " / " + (xy[1] * R).toFixed(10) + "<br>");
	inverse (xy[0], xy[1], xy);
	document.write((xy[0] / Math.PI * 180).toFixed(10) + " / " + (xy[1] / Math.PI * 180).toFixed(10));
</script>

</body>