How the test works
Project every lon/lat sample with a candidate CRS and fit x and y as straight-line functions
of the column and row (an affine transform, six numbers). What the fit misses, taken back into
index space, is the residual in cells. A grid is regular in that CRS when the residual is down
at the precision of the numbers it shipped with, float32 lon/lat being good to about a metre.
The affine soaks up an offset, a scale, a rotation and a shear, so only parameters that bend
the grid need searching. Polar stereographic has none: lat_ts is a scale and lon_0 a rotation.
Lambert conformal conic has one, the cone constant. Geostationary has two, the sub-satellite
longitude and the height. A rotated pole has two, the pole. These are found with Nelder-Mead;
then lon_0 is turned to make the grid north-up, and the cell size is checked against the radii
and lat_ts values that make it a round number, which is how the Hughes ellipsoid and
lat_ts = 70 come back out of the NSIDC arrays.
The result is a CRS and a geotransform, which is all a reader needs (GDAL's
-a_srs and -a_ullr, or a grid spec);
the two 2D arrays can be dropped. The idea is Michael Sumner's "regular in what?" question,
asked of many datasets over the years.