Mapping latlng values to a point on a HRAP grid

I have a copy of some NWS rainfall data loaded into a PostGIS database. Given a lat lng point, I want to find the rainfall for that given day.

The original data are in XMRG format and projected in the Hydrologic Rainfall Analysis Project (HRAP) grid coordinate system, a polar stereographic projection true at 60°N / 105°W.

Based on a description of the HRAP system I found from the NWS, the calculations for locating a point on the HRAP grid based on lat lng value are as follows:

RE = (EARTHR*(1+SIN(60)))/xmesh
R = (RE*COS(xlat))/(1+SIN(xlat))
WLONG = xlon+75
X = R*SIN(WLONG)+401
Y = R*COS(WLONG)+1601

EARTHR is the radius of the earth (6371.2 KM)
xmesh is the mesh length at 60° latitude (4.7625 KM)
xlat is the latitude of point to be converted (decimal degrees)
xlon is the longitude of point to be converted (decimal degrees)

I wrote up a quick javascript implementation to test it:

var deg2rad = function (deg) { return deg * (Math.PI / 180); };

var hrap = function (xlon, xlat) {
var rxlon, rxlat, re, r, wlong, rwlong, x, y;
var earthr = 6371.2;
var xmesh = 4.7625;

with (Math) {
re = (earthr * (1 + sin(PI / 3))) / xmesh;
r = (re * cos(rxlat)) / (1 + sin(rxlat));

wlong = xlon + 75;
x = round(r * sin(rwlong) + 401);
y = round(r * cos(rwlong) + 1601);
}

return [x, y];
};

I ran a couple data points from the 2013-05-17 data set through to test it, with the following results:

Lng        Lat      NWS X  NWS Y  Calc X  Calc Y
-76.5633  39.7478    958    571     369    2771
-124.6132  42.7528     34    572    -431    2308

Those numbers are nowhere close to correct, and I don't know why. Is the document wrong? Did I mess up the implementation? Is the test data bad? Am I missing a simpler solution? I'd appreciate it if someone could point out where this is going wrong.

• Javascript sin and cos functions take radians not degrees. You should start with calculating zmesh for 60 north, you'll see it's not right in JS, also it's missing brackets (wolframalpha.com/input/…) so both JS and the document may be wrong. – Jakub Kania May 26 '13 at 16:02
• I fixed the deg/rad problem. Those missing parentheses concerned me too, but is that relevant? The only thing zmesh is used for in the mapping is to set the value of xmesh. – Brad Koch May 28 '13 at 14:24
• Well, if there is one thing missing there's bigger chance something else is missing if you ask me, feel free to compare it with the code I posted below. SR org also gives wrong results so I guess HRAP is just weird: spatialreference.org/ref/sr-org/4695 – Jakub Kania May 28 '13 at 19:41

I made a port of part of this fortran code:

function hrap(lon,lat){
var deg2rad = function (deg) { return deg * (Math.PI / 180); };
var stlat = 60.0;
var clon = 15.0;
var rad = 6371.2;
with (Math){