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 REAL*8 (DOUBLE PRECISION) rather than in REAL*4, 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

*C=1+INT((X-X0)/DX), R=1+INT((Y-Y0)/DY)*.

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.