A Mapping and Gridding Primer: Points, Pixels, Grids, and Cells
Satellite remote sensing has made possible the collection of data over large areas of the Earth. These data are often stored in grids. Grids are an efficient means of storing data because the location of a value within the grid is implicitit is not explicitly stored in the grid. The location is also constant, which makes it easy to compare data from different sensors or different time periods.
The purpose of this document is to give a brief introduction to mapping, gridding, and the associated terminology used in the Mapx library at NSIDC. Along with comments in the source code, it also describes the format of the map and grid parameter files.
Concepts
Maps
A map projection is a procedure or mathematical formula to transform a curved surface onto a plane. The curved surface is usually the surface of the Earth and the plane is what we call a map. The derivation of a map transformation requires a mathematical representation of the surface of the Earth. The surface of constant geopotential referred to as mean sea level or the geoid is quite complicated and its mathematical form is correspondingly complex. When defining map projections then, it is common to approximate the surface of the earth with a simpler surface such as an ellipsoid or sphere.
The ellipsoid is flattened at the poles, relative to the sphere, by about 1 part in 300. For maps of large areas (continents or bigger) using the sphere introduces no significant error. When mapping smaller areas it is not uncommon to choose an ellipsoid specifically for each area; each area is mapped on an ellipsoid with a different radius and eccentricity.
Coordinates for points on the surface of the Earth are given in latitude and longitude. Geodetic latitude is the angle between a vertical ray at the point and the plane of the equator. Geocentric latitude is the angle between a ray from the center of the Earth to the point and the plane of the equator. On a sphere, the geodetic latitude is equivalent to the geocentric latitude. It's more common to use geodetic latitude. There is only one definition of longitudethe angle between a ray from the center of the Earth to the point and the plane of the Greenwich Meridian (an arbitrary reference point). Coordinates in the map plane (x,y) are defined in the usual manner with x horizontal and increasing to the right and y vertical and increasing upward.
Projections
Map transformations are referred to as projections because they are commonly visualized as geometric constructs. That is, we think of the map as a piece of paper wrapped around a globe with lines emanating from points on the globe and "projected" onto the paper. It's important to note that not all map projections can be constructed in this manner.
Map projections are specified as a set of equations giving x and y in terms of latitude and longitude. The manner in which we wrap the paper around the globe determines the form of the equations. Cylindrical projections wrap the paper around in a tube, conical projections use a cone, and azimuthal projections use a flat surface. Several parameters are used in these equations to specify exactly where the surface intersects the globe, the scale or size of the globe, and the translation and rotation of the map on the plane surface.
There is no one best map projection. Each projection has different properties and different best uses. Two of the most important characteristics of maps are whether they are conformal or equalarea. No map projection is both, and some are neither. If a map is conformal then angles within a small area are reproduced accurately. This means that shapes are preserved. A small circle on the globe will look like a small circle on the map. In formal terms, at any point on the map the scale h along a meridian of longitude is equal to the scale k along a parallel of latitude. If a map is equalarea then a small circle placed anywhere on the map will always cover the same amount of area on the globe and the product of h and k at any point is one.
For maps that are not equalarea, hk  1 gives a measure of the areal distortion. For maps that are not conformal, k/h gives the aspect ratio, which is a measure of shape distortion. For example, on the Polar Stereographic map true at 70N (a conformal map) the areal distortion varies from 6% at the pole to +29% at 45N and goes up to +276% at the equator while, by definition, the aspect ratio remains 1:1 everywhere. The Azimuthal EqualArea map over the pole has, of course, no areal distortion, while the aspect ratio varies from 1:1 at the pole to 1.17:1 at 45N and goes up to 2:1 at the equator.
Tissot's indicatrix combines the two previous measures. Tissot's indicatrix is a small ellipse which shows both the scale distortion and the maximum angular deformation at a particular point on the map. See Maling or Snyder for a complete discussion.
Grids
A grid is a rectangular array of points. Grids record regularly spaced samples over an area. When sampling over the surface of the earth a grid is determined by a map projection, a sampling interval, an origin, and the number of rows and columns. So called lat/lon or equal angle grids are equivalent to grids sampled on a Cylindrical Equidistant map projection.
A grid coordinate system is defined in the map plane with axes parallel to the rows and columns of the grid and units equal to the sampling interval. The grid sample locations are at the whole integer grid coordinate points. To conform to mathematical array conventions, the grid coordinates (r,s) start at (0,0) in the upper left corner with r increasing to the right, and s increasing downward (this also conforms to digital image processing conventions). A grid cell (i,j) is defined as the area between grid coordinates i.5 and i+.5, and j.5 and j+.5. To conform to rounding conventions, the lower bound is included in the cell while the upper bound is not (round up at .5). Note that the r coordinate corresponds to the grid column number j, and the s coordinate corresponds to the grid row number i.
Think of the grid as a map drawn on graph paper. To find a grid cell location, count graph lines over from the left and down from the top, then read the latitude and longitude off the map.
Binning
The process of sampling data on a grid is sometimes referred to as binning, especially when the sampling method is to average all data that falls into the grid cell. Other methods used to sample data are: take the closest value to the grid point (nearest neighbor) or interpolate between surrounding grid points. Also, binning doesn't need to be the average of multiple data points; it could be the minimum, the maximum, the median, or the latest value in the grid cell. The characteristics of the data should determine the sampling method.
Pictures
Pixels are analogous to grid cells. To display gridded data as a picture, define a transformation from grid coordinates to pixel coordinates and then sample the gridded data at the whole integer pixel coordinate points. The most common technique is to map grid cells onto pixels one for one. More sophisticated techniques, often referred to as "resampling", allow for scrolling, zooming, and rotating.
Parameters
Many programs we use at NSIDC depend on the mapx module to handle map projections. These include: PSQ, regrid, gridomatic, and the EASEGrid processing software. Most of these tools provide access to predefined maps, but, it is relatively easy to define your own map with a map projection parameters (.mpp) file.
Map Projection Parameters (.mpp) File Format
Each parameter is described by a keyword: value pair on a single line. A colon must delimit the end of a keyword without any intervening spaces, tabs, or newlines. Keywords are not case sensitive. Parameters can be followed by comments, which begin with a semicolon (;) or a pound sign (#): all text following a semicolon or a pound sign on a particular line is ignored. Parameters can appear in any order. The default value is used when the corresponding keyword is not found. Unrecognized keywords are silently ignored. This means that if you misspell a keyword you'll get the default value without warning. Table 1 lists the parameters that define the map projection for all projections except Integerized Sinusoidal and Universal Transverse Mercator (UTM).
Keyword  Format  Default  Description 

Map Projection  string  required field  projection name (see list below) 
Map Reference Latitude  decimal degrees  required field  reference latitude for map projection 
Map Reference Longitude  decimal degrees  required field  reference longitude for map projection 
Map Second Reference Latitude  decimal degrees  none  used by some projections (see explanation below) 
Map Rotation  decimal degrees  0.0  rotation counterclockwise 
Map Scale  float  1.0  map scale factor (radius units per map unit) 
Map Origin Latitude  decimal degrees  reference latitude  translated map origin latitude 
Map Origin Longitude  decimal degrees  reference longitude  translated map origin longitude 
Map Origin X  float  x coordinate of origin latitudelongitude transformed to map coordinates  translated map x origin in map units 
Map Origin Y  float  y coordinate of origin latitudelongitude transformed to map coordinates  translated map y origin in map units 
Map False Easting  float  0.0  offset map x origin in map units 
Map False Northing  float  0.0  offset map y origin in map units 
Map Eccentricity  float  Clark 1866 = 0.082271673  eccentricity of ellipsoid 
Map Eccentricity Squared  float  Clark 1866 = 0.00676862817822  square of eccentricity of ellipsoid 
Map Equatorial Radius  float  Ellipsoidal projections: Clark 1866 ellipsoid = 6378.2064 km  equatorial radius of ellipsoid 
Spherical projections: radius of authalic sphere to 1924 international ellipsoid = 6371.228 km  radius of sphere  
Map Polar Radius  float  Clark 1866 ellipsoid = 6356.5838 km  polar radius of ellipsoid 
Map Center Scale  float  1.0  scale factor at central meridian for Transverse Mercator 
Map Maximum Error  float  0.0  maximum allowed error in map units (0 disables map error checking) 
Table 2 lists the parameters used only by programs that actually draw maps.
Keyword  Format  Default  Description 

Map Southern Bound  decimal degrees  90.00S  bottom of map, starting point for latitude graticule 
Map Northern Bound  decimal degrees  90.00N  top of map 
Map Western Bound  decimal degrees  180.00W  left side of map, starting point for longitude graticule 
Map Eastern Bound  decimal degrees  180.00E  right side of map 
Map Graticule Latitude Interval  decimal degrees  30.  graticule spacing 
Map Graticule Longitude Interval  decimal degrees  30.  graticule spacing 
Map Graticule Label Latitude  decimal degrees  0.00N  where to label meridians 
Map Graticule Label Longitude  decimal degrees  0.00E  where to label parallels 
Map CIL Detail Level  integer  1  level of detail for coastlines, islands, and lakes 
Map BDY Detail Level  integer  0  level of detail for political boundaries 
Map RIV Detail Level  integer  0 
level of detail for rivers

NOTE: For Integerized Sinusoidal, set scale to size of each row.
Map Projection Names
 Albers Conic EqualArea
 Azimuthal EqualArea
 Azimuthal EqualArea (ellipsoid)
 Cylindrical EqualArea
 Cylindrical EqualArea (ellipsoid)
 Cylindrical Equidistant
 Integerized Sinusoidal
 Interrupted Homolosine EqualArea
 Lambert Conic Conformal (ellipsoid)
 Mercator
 Mollweide
 Orthographic
 Polar Stereographic
 Polar Stereographic (ellipsoid)
 Sinusoidal
 Transverse Mercator
 Transverse Mercator (ellipsoid)
 Universal Transverse Mercator
Keywords must be entered verbatim. Capitalization doesn't matter but white space does. Unrecognized (or mispelled) keywords are silently ignored. Map names are not case sensitive, white space and dashes are ignored, and most word permutations are acceptable. For example, EqualArea Cylindrical is equivalent to CYLINDRICAL EQUALAREA.
The reference latitude and longitude specify the original location and orientation of the map projection. The origin of the rectangular map coordinate system can be rotated and translated with the Rotation and Origin keywords. Note that a positive rotation value will rotate the coordinate system counterclockwise around the translated map origin. This will have the effect of rotating the resulting grid of pixel values clockwise with respect to the grid that would be obtained if the rotation value were 0.
For the Cylindrical EqualArea, and Polar Stereographic projections, the second reference latitude specifies the latitude of true scale. Conceptually, this is where the projection plane intersects the surface. For the Albers Conic EqualArea and Lambert Conic Conformal projections, the conic plane is secant to the surface between the reference latitude and the second reference latitude.
The equatorial radius and eccentricity specify the surface to be mapped (usually the Earth's surface). Alternatively, the square of the eccentricity rather than the eccentricity may be specified; or both the equatorial radius and the polar radius may be specified rather than either the eccentricity or the square of the eccentricity; or the polar radius and either the eccentricity or the square of the eccentricity map be specified rather than the equatorial radius. The scale and the equatorial radius determine the map units. For example if the radius is in kilometers (as is the case with the default) and the scale is 1.0, then the map units will be kilometers. For spherical projections, the eccentricity should not be specified, and the default equatorial radius is 6371.228, i.e. the authalic sphere, equivalent in surface area to the 1924 international ellipsoid, in km. For elliptical projections, the default eccentricity is 0.082271673 and the default equatorial radius is 6378.2064, i.e. the Clark 1866 ellipsoid in km. For the WGS84 ellipsoid in meters, use 0.081819190843 for the eccentricity and 6378137.0 for the equatorial radius. PostScript units are inches, so for PostScript output you would specify the scale in kilometers per inch. For a grid you can define the scale in terms of kilometers per grid cell. The gridding module grids also has parameters which specify grid cells per map unit or, alternatively, map units per grid cell. This allows the same .mpp file to be used for multiple nested grids.
Map Origin Latitude and Longitude specify the translated map origin. By default the map origin is the same as the reference latitude and longitude. Alternatively, the translated map origin can be specified by Map Origin X and Y, which take precedence over any specification of the Map Origin Latitude or Longitude. Both Map Origin X and Y must be specified if either is specified. If Map Rotation is specified, then the rotation is centered on the translated map origin. Furthermore, the specification of Grid Map Column and Row values in the Grid Parameter Definition File will locate the translated map origin at the specified column and row location within the grid.
False easting and false northing offsets are applied to the x and y map coordinates, respectively, with respect to the original map origin specified by the reference latitude and longitude. Note that Map Origin X and Y specify values after false easting and false northing offsets have been applied.
The center scale value is used only by the Transverse Mercator projection (and for the UTM projection as well). It specifies the scale factor along the central meridian.
When the maximum error value is nonzero, then error checking, is enabled. For a forward transformation (latlon to xy), the computed xy pair is inverse transformed back into a second latlon pair. For an inverse transformation (xy to latlon), the computed latlon pair is forward transformed into a second xy pair. In both cases, the distance in map units between the first and the second pairs is computed. If it is greater than the maximum error value, then both values of the returned pair (xy for forward and latlon for inverse) are set to notanumber (nan). If the maximum error value is 0 (the default), then no such error checking is performed. Note that setting maximum error to a nonzero value will slow down forward and inverse transformations by about a factor of 2.
The remaining parameters are used by programs that actually draw maps (e.g. psmap, mapenum). The southern, northern, western, and eastern bounds of the map are the starting point for the graticule and are used to speed up the search in the map outline database. They cannot be counted on to clip the map accurately. All longitudes should be in the range 180 to +360. West to east should not span more than 360 degrees. West specifies the left side of the map and east the right, not necessarily the minimum and maximum longitudes.
The latitude and longitude graticule intervals specify the spacing between graticule lines. The graticule is the overlay of latitude and longitude lines on the map. By default, a parallel is drawn every 30 degrees starting at the bottom of the map (90.00S) and a meridian is drawn every 30 degrees starting at the left side of the map (180.00W). Label latitude and longitude specify the parallel and meridian along which to draw the graticule line labels.
The last three numbers are used by the database search routine to specify the level of detail for map lines. The level refers to the number and size of features which will be selected, as opposed to resolution of each feature. The higher the number, the more detail will be included. For example, level 1 will include only the largest features, such as whole continents, while level 2 would also include some smaller features like large lakes. In either case the continental outline will appear with the same amount of detail (the same jaggedness)
Map Projection Parameters (.mpp) File Format for Integerized Sinusoidal
Table 3 lists the parameters for the Integerized Sinusoidal (ISin) projection.
Keyword  Format  Default  Description 

Map Projection  string  required  Integerized Sinusoidal 
Map Reference Longitude  decimal degrees  0.0  Longitude of the central meridian 
Map Rotation  decimal degrees  0.0  rotation counterclockwise 
Map Scale  float  1.0  map scale factor (meters per map unit) 
Map ISin NZone  int  86400  Number of equally spaced latitudinal zones 
Map ISin Justify  int  1.0=left  justify flag for ISin map 
Map Origin Latitude  decimal degrees  0.0 (= effective reference latitude)  translated map origin latitude 
Map Origin Longitude  decimal degrees  effective reference longitude  translated map origin longitude 
Map Origin X  float  x coordinate of origin latitudelongitude transformed to meters  translated map x origin in meters 
Map Origin Y  float  y coordinate of origin latitudelongitude transformed to meters  translated map y origin in meters 
Map False Easting  float  0.0  offset map x origin in meters 
Map False Northing  float  0.0  offset map y origin in meters 
Map Equatorial Radius  float  6371007.181  radius of sphere in meters 
For the ISin projection all map projection parameters have the same definitions as described in Table 3 with the following exceptions.
For the ISin projection, the Map Reference Latitude parameter value is forced to 0.0; any Map Reference Latitude value supplied by the user is ignored. The Map Reference Longitude parameter is not required for the ISin projection; the default value is 0.0.
The Map ISin NZone parameter is used only by the ISin projection; it specifies the number of equally spaced latitudinal zones (i.e. "rows"), and must be 2 or larger and even. Note that this value is a characteristic of the projection only and is independent of the actual number of rows in any grid that utilizes the ISin projection. Also note that currently only the default value of 86400 is used by all ECS MODIS products (including all 1000 m, 500 m, and 250 m products) based on the ISin projection.
The Map ISin Justify parameter is used only by the ISin projection; it is used to indicate what to do with zones with an odd number of columns. If it has a value of 0 or 1, it indicates the extra column is on the right (zero) or left (one) of the projection yaxis. If the flag is set to 2, the number of columns is calculated so there is always an even number of columns in each zone. Note that currently only the default value of 1 is used by all ECS MODIS products based on the ISin projection.
By default, setting Map Projection to Integerized Sinusoidal specifies a spherical projection using a radius of 6371007.181 meters. Note that this value is currently used for all ECS MODIS products based on the ISin projection. Only a spherical projection is supported for the ISin projection.
Map Projection Parameters (.mpp) File Format for Universal Transverse Mercator
Table 4 lists the parameters defined by the Universal Transverse Mercator (UTM) projection.
Keyword  Format  Default  Description 

Map Projection  string  required  Universal Transverse Mercator 
Map UTM Zone  integer  0  specifies value for UTM zone 
Map Reference Latitude  decimal degrees  determined by UTM zone  determines sign of UTM zone 
Map Reference Longitude  decimal degrees  determined by UTM zone  determines absolute value of UTM zone 
Map Rotation  decimal degrees  0.0  rotation counterclockwise 
Map Scale  float  1.0  map scale factor (meters per map unit) 
Map Origin Latitude  decimal degrees  0 (= effective reference latitude)  translated map origin latitude 
Map Origin Longitude  decimal degrees  effective reference longitude  translated map origin longitude 
Map Origin X  float  x coordinate of origin latitudelongitude transformed to meters  translated map x origin in meters 
Map Origin Y  float  y coordinate of origin latitudelongitude transformed to meters  translated map y origin in meters 
Map False Easting  float  500000  offset map x origin in meters 
Map False Northing  float  0 for northern hemisphere; 10000000 for southern hemisphere 
offset map y origin in meters 
Map Eccentricity  float  WGS84 = 0.081819190843  eccentricity of ellipsoid 
Map Eccentricity Squared  float  WGS84 = 0.0066943799902  square of eccentricity of ellipsoid 
Map Equatorial Radius  float  WGS84 = 6378137.0 meters  equatorial radius of ellipsoid (or radius of sphere if eccentricity = 0) in meters 
Map Polar Radius  float  WGS84 = 6356752.314245 meters  polar radius of ellipsoid in meters 
Map Center Scale  float  0.9996  scale factor at central meridian 
Map Maximum Error  float  100.0 meters  maximum allowed error in meters (0 disables map error checking) 
For the UTM projection all map projection parameters have the same definitions as described in Table 5 with the following exceptions.
The UTM zone value is used only by the UTM projection. Zones are in the range 1 to 60 for the northern hemisphere, and 1 to 60 for the southern hemisphere. If the UTM zone is 0, then the hemisphere and absolute value of the zone are determined by the reference latitude and longitude, respectively. If the UTM zone is not 0, then the reference latitude and longitude are not required and are ignored if present. The effective reference latitude value for every UTM zone is 0 (the equator); the effective reference longitude value is the central meridian of the zone starting at 177W for zones 1 and 1, 171W for zones 2 and 2, and so on until 177E for zones 60 and 60. Note that either Map UTM Zone must be nonzero or both Map Reference Latitude and Map Reference Longitude must be specified.
The default Map Origin Latitude and Longitude values are determined by the effective reference latitude value, which is always 0, and the effective reference longitude value, which is determined from the UTM zone, either as supplied in the Map UTM Zone value or as derived from the Map Reference Longitude value.
By default, setting Map Projection to Universal Transverse Mercator specifies an ellipsoidal projection using the WGS84 ellipsoid. To specify a spherical UTM projection, set Map Eccentricity to 0, and set Map Equatorial Radius to the radius of the sphere in meters.
For the UTM projection, the default false easting is 500000 (five hundred thousand) meters. For northern hemisphere UTM zones, the default false northing is 0 meters, and for southern hemisphere UTM zones, the default false northing is 10000000 (10 million) meters. Note that these values assume that Map Equatorial Radius is in meters and that Map Scale is 1.0. Note also that Map Origin X and Y specify values after false easting and false northing offsets have been applied. For example, for the UTM projection for a southern hemisphere UTM zone, specifying a Map Origin X value of 175000 will set the translated map origin to 325000 (500000  175000) meters west of the central meridian for the UTM zone, and specifying a Map Origin Y value of 5325000 will set the translated map origin to 4675000 (10000000  5325000) meters south of the equator.
The center scale value specifies the scale factor along the central meridian; the default for the UTM projection is 0.9996.
When the maximum error value is nonzero, then error checking, is enabled; the default for the UTM projection is 100.0 meters. For a forward transformation (latlon to xy), the computed xy pair is inverse transformed back into a second latlon pair. For an inverse transformation (xy to latlon), the computed latlon pair is forward transformed into a second xy pair. In both cases, the distance in map units between the first and the second pairs is computed. If it is greater than the maximum error value, then both values of the returned pair (xy for forward and latlon for inverse) are set to notanumber (nan). Note that when error checking is enabled, it will slow down forward and inverse transformations by about a factor of 2. To disable error checking, set the maximum error value to 0. Error checking is turned on by default for the UTM projection because the equations used for the projection can produce errors of more than 100 meters for points that are at a latitude of + 45 degrees and more than about 17.5 degrees of longitude beyond the center of the specified 6 degreewide UTM Zone; and in some cases, forward transformations of a latlon pair that is located far outside of the specified UTM Zone can produce an xy pair that resides within the specified UTM Zone.
Grid Parameter Definition (.gpd) File Format
Table 5 lists the grids that are defined in grid parameter definition (.gpd) files.
Keyword  Format  Default  Description 

Grid Width  integer  required field  number of columns in grid 
Grid Height  integer  required field  number of rows in grid 
Grid Map Origin Column  float  0.0  column coordinate of map origin 
Grid Map Origin Row  float  0.0  row coordinate of map origin 
Grid Cells per Map Unit  float  1.0  sets both column and row scale 
Grid Columns per Map Unit  float  1.0  sets column scale if not already set 
Grid Rows per Map Unit  float  1.0  sets row scale if not already set 
Grid Map Units per Cell  float  1.0  alternative method to set both column and row scale 
Grid Map Units per Column  float  1.0  sets column scale if not already set 
Grid Map Units per Row  float  1.0  sets row scale if not already set 
Grid MPP File  string  none  name of .mpp file, map projection parameters can also be specified in the same file with the grid definition parameters 
Limitations
The mapx library is very flexible but it does impose certain restrictions on the kinds of grids that can be defined. All grids must be rectangular arrays. Rectangular arrays are easy to use with existing scientific visualization tools, like IDL, but they may be inefficient in terms of storage space for nonrectangular data sets. If some flag value (for example zero) is stored in the unused cells, then this wasted space can be compacted by any good data compression scheme. This same invalid data flag can be used within the grid as well, leading to further storage savings. Another restriction is that each grid must be based on a map projection defined in mapx. Many of the most useful projections have already been implemented and fortunately it is easy, given a mathematical formulation of the forward and inverse transformations, to add a new projection to mapx (see source code for instructions). The mathematical formulations of most common map projections can be found in the references below.
References
American Society of Photogrammetry. 1983. Manual of Remote Sensing. Robert N. Colwell, editor. Second edition. Falls Church, VA.
Gonzalez, Rafael C. and Paul Wintz. 1987.Digital Image Processing.Second edition. Reading, MA: AddisonWesley.
Maling, D. H. 1992. Coordinate Systems and Map Projections. Second edition. Elmsford, N.Y.: Pergamon Press.
Snyder, John P. 1987. Map Projections, A Working Manual. U. S. Geological Survey Professional Paper 1395. Department of the Interior. Washington, D. C.
Last Updated December 2017