Convert HRAP coordinates to latitude/longitude

Questions
1

Hello all,

I’m trying to transform precipitation data in HRAP format to latitude/longitude coordinate system in Mathematica. Any help appreciated!

Thanks,

Patrice

2

Hello,

I think you are referring to the Polar Stereographic coordinates used in the Hydrologic Rainfall Analysis Project, which, if I’m not mistaken, correspond in Wolfram Language to this projection:

proj = {"Stereographic",
        "Centering" -> {90, -105},
        "GridOrigin" -> {401, 1601},
        "ReferenceModel" -> 6371.2 / 4.7625,
        "CentralScaleFactor" -> (1. + Sin[60 Degree])/2}

This is using a standard parallel of 60 degrees, which determines the central scale factor in the last line, and a central meridian of -105 degrees. It places the North pole (the centering of this polar stereographic projection) at the coordinate position {401, 1601}, the so-called “grid origin”. Finally, it uses a radius of the Earth of 6371.2 kilometers and “cells” of 4.7625 kilometers.

Then for example we can project any {lat, lon} coordinate pair:

In[]:= GeoGridPosition[GeoPosition[{40, -90}], proj]
Out[]= GeoGridPosition[{702.282, 476.602}, {"Stereographic", ...}]

and you can get the projected values out with

In[]:= %["Data"]
Out[]= {702.282, 476.602}

You can project many points at the same time (for speed) putting them all at once inside a single GeoPosition head, like RandomGeoPosition does automatically:

pos = RandomGeoPosition[Entity["Country", "UnitedStates"], 1000]
GeoGridPosition[pos, proj]

and again use …[“Data”] to extract the projected coordinates of the 1000 points.

This will produce a map of the US using that projection, showing projected-space frame ticks:

GeoGraphics[Entity["Country", "UnitedStates"],
 GeoProjection -> proj, Frame -> True, GeoGridLines -> Automatic]
3

Hi Jose,

Thanks for your explanations, gives me a lot of pointers to investigate further (I’m new to GIS).

How would I transform the other way around though, from a grid of Polar Stereographic coordinates to Lat/Lon?

4

Hi again,

Once you have the result in the form

In[]:= gridpos = GeoGridPosition[GeoPosition[{40, -90}], proj]
Out[]= GeoGridPosition[{702.282, 476.602}, {"Stereographic", ...}]

you can recover the latlon form using

In[]:= GeoPosition[gridpos]    
Out[]= GeoPosition[{40., -90.}]

The property “Data” can also be used to extract the coordinate numbers from a GeoPosition object.

5

Works like a charm, thanks a lot!