OK. let’s take it from the top.
The I/O API coordinate system treatment is based upon the USGS Geographic Coordinate Transform Package (GCTP), which at the time we started the project was the easiest high-quality publicly-available geographic package for us to use. Its original manual is at https://cjcoats.github.io/ioapi/GCTP.pdf and many of the I/O API coordinate concepts and parameters come from it.
The Earth is approximately an ellipsoid, in particular, because of the force of the Earth’s rotation it is a “flattened” spheroid – the rotation causes the diameter to “stretch” in the plane of the Equator (so that the diameter across the Equator is greater than the diameter pole-to-pole. A spheroid is such a description of the Earth; there have been many such descriptions proposed and used; the most common, together with their GCTP code-numbers are
- Clarke 1866: GCTP uses 0 as its index-number
- GRS 1980: 8
- WGS 84: 12
- perfect sphere: 19 (in which case we have to specify the radius elsewhere)
Generally the more recent GRS 80 or WGS 84 ellipsoids are what GIS people use; meteorologists tend to use perfect spheres since doing proper equations of motion on a more-general ellipsoid is quite “hairy” mathematics. I have seen some older data given in terms of Clarke 1866…
A map projection is a way of putting Cartesian X-Y coordinates on ellipsoids (and so needs a ellipsoid-specification to complete its definition). A Lambert Conformal Conic, for example (I/O API projection code 2, GCTP projection code 4), has the geometric algorithm:
- Form a cone that intersects the ellipsoid at two latitudes alpha and beta (and has central axis the same as the rotational N-S axis of the Earth).
- Project a geometric ray from the center of the Earth through both its surface and the surface of the cone. Map the Earth-surface intersection-point of that ray to the cone’s intersection point. If you do the arithmetic for this you’ll get into some rather messy trigonometry
- Slice the cone diametrically opposite “central longitude” gamma and unroll it flat, with the longitude-gamma’s image on the cone in the N-S, or Y direction.
- Find the point on the unrolled-cone where “origin-point” latitude-longitude (LATORIG,LONORIG) winds up, and start measuring X and Y coordinates from there – Y parallel to the longitude-gamma’s line. Usually one chooses LONORIG=gamma, and LATORIG somewhere between alpha and beta (but not always;-( )
Because of how this works, choosing different ellipsoids to represent the Earth will give somewhat-different X-Y coordinate systems; if you use Lambert projections with alpha=30, *beta=*60, and gamma=90, LONORIG=-90 (note the minus-sign, for longitudes west of the Prime Meridian LON=0), and LATORIG=45, then the Cartesian origin X=Y=0 somewhere in Illinois no matter which ellipsoid you use, but the coordinates for Portland Maine will be different by about 20 KM between choosing a spherical ellipsoid and the WGS84 ellipsoid.
Note that because of round-off problems in all the map-projection arithmetic, it is best to do everything in REAL8 (DOUBLE PRECISION) rather than in REAL4, which doesn’t have enough significant digits to be sure of giving the “right” answer…
And the I/O API’s coordinate transform routines generally do this work for you, given descriptions you don’t need to worry much about the messy mathematics.
The worst problem you’re likely to encounter is transforming back and forth between grid-cell columns and rows and Cartesian X-Y coordinates:
Suppose you have a grid with lower-left corner (X0,Y0) and cell-size (DX,DY), and a point (X,Y) with X≥X0, Y≥Y0. Then the column-
and row-numbers (C,R) for the point are
Beware that computer-languages “convert to integer” operators tend to “do the wrong thing” from a mathematical point of view: INT(-0.5) and INT(0.5) are both 0 ;-( so you have to do explicit checks for X≥X0, Y≥Y0 and say “the point is off the grid” if either one fails.