4

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;

    rxlat = deg2rad(xlat);

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

        wlong = xlon + 75;
        rwlong = deg2rad(wlong);
        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
5
+50

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){
var sfactor = (1+sin(deg2rad(stlat)))/(1+sin(deg2rad(lat)));
var R = rad*cos(deg2rad(lat))*sfactor;
var x = R*cos(deg2rad(lon+clon));
var y = R*sin(deg2rad(lon+clon));
var hrapx = x/4.7625 + 401;
var hrapy = y/4.7625 + 1601;
return [hrapx,hrapy];
}}

Here's a working example.

  • This is phenomenal Jakub! Digging into a random fortran reference implementation? Well deserved. – Brad Koch May 28 '13 at 20:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.