## Verify that the corner map coordinates in the polar stereographic grid page (http://nsidc.org/data/polar_stereo/ps_grids.html) represent the outer edges of grid cells.¶

### The table on the PS page shows map coordinates (x, y) and geographic (latitude/longitude). To convert from geographic to map coordinates it's a more complicated set of equations which are a little long to deal with here. If there is a mistake in the implementations of that set of equations, it would take digging into the mapx code and the USGS standard equations that were implemented in that code.¶

#### Contents of S3B.gpd and Sps.mpp¶

Sps.mpp     /* map projection parameters    # SSM/I Brightness Temp Grids */
316 332     /* columns rows         # South Polar Region */
4 4     /* grid cells per map unit      # 25km */
157.5 173.5 /* map origin column, row       # */

PolarStereographicEllipsoid
-90.0 0.0 -70.0     /* lat0 lon0 lat1 */
0.0         /* rotation */
100.0           /* scale (km/map unit) */
-90.00  0.00        /* center lat lon */
-90.00  -20.00      /* lat min max */
-180.00 180.00      /* lon min max */
10.00 15.00     /* grid */
0.00    0.00        /* label lat lon */
1 0 0           /* cil bdy riv */
6378.273        /* Earth equatorial radius (km) */
0.081816153     /* eccentricity */

#### For me it is easiest to see the center vs edge values considering not geographic coordinates (latitude/longitude) to map coordinates, but using grid coordinates and map coordinates so below is the conversion between grid columns and rows to map coordinate x (u) and y (v) values.¶

In [10]:
def translate_grid_to_mapcoords(r, s):
# From S3B.gpd
map_origin_col = 157.5 # grid
map_origin_row = 173.5
cols_per_map_unit = 4
rows_per_map_unit = 4

# From Sps.mpp
scale = 100.0 #km/map unit

# From mapx/mapx.c/inverse_mapx
# u,v -> rotated and translated map coordinates in map units
u = (r - map_origin_col) / cols_per_map_unit
v = -(s - map_origin_row) / rows_per_map_unit
u *= scale
v *= scale
print "u (x) = {u}  v (y) = {v}".format(u= u, v= v)


#### Convert the origin as defined in the gpd to show we are using the correct calculations.¶

In [11]:
print "Origin (157.5, 173.5):"
translate_grid_to_mapcoords(157.5, 173.5)

Origin (157.5, 173.5):
u (x) = 0.0  v (y) = -0.0


#### Convert the corners of the grid (width = 316 cols, height = 332 rows) to map coordinates. Map coordinates are the X (km) / Y (km) values shown in the table on the PS page.¶

In [12]:
print "Upper left grid point: (0, 0):"
translate_grid_to_mapcoords(0, 0)
print "Upper right grid point: (315.0, 0):"
translate_grid_to_mapcoords(315.0, 0)
print "Lower right grid point: (315.0, 331.0):"
translate_grid_to_mapcoords(315.0, 331.0)
print "Lower left rid point: (0, 331.0):"
translate_grid_to_mapcoords(0, 331.0)

Upper left grid point: (0, 0):
u (x) = -3937.5  v (y) = 4337.5
Upper right grid point: (315.0, 0):
u (x) = 3937.5  v (y) = 4337.5
Lower right grid point: (315.0, 331.0):
u (x) = 3937.5  v (y) = -3937.5
Lower left rid point: (0, 331.0):
u (x) = -3937.5  v (y) = -3937.5


#### If we adjust the grid points to reflect a shift of 1/2 to move the points away from the center to the outer edge, we get the map coordinate values listed in the PS table.¶

In [13]:
print "Upper left grid point: (-0.5, -0.5):"
translate_grid_to_mapcoords(-0.5, -0.5)
print "Upper right grid point: (315.5, -0.5):"
translate_grid_to_mapcoords(315.5, -0.5)
print "Lower right grid point: (315.5, 331.5):"
translate_grid_to_mapcoords(315.5, 331.5)
print "Lower left grid point: (-0.5, 331.5):"
translate_grid_to_mapcoords(-0.5, 331.5)

Upper left grid point: (-0.5, -0.5):
u (x) = -3950.0  v (y) = 4350.0
Upper right grid point: (315.5, -0.5):
u (x) = 3950.0  v (y) = 4350.0
Lower right grid point: (315.5, 331.5):
u (x) = 3950.0  v (y) = -3950.0
Lower left grid point: (-0.5, 331.5):
u (x) = -3950.0  v (y) = -3950.0


#### Another way to see this is to convert geographic coordinates to grid rows and columns using gtest mapx utility.¶

[email protected]:~/nsidc0080/mapx> gtest

enter .gpd file name: S3B.gpd
> assuming old style fixed format file

gpd: S3B.gpd
mpp:Sps.mpp

forward_grid:
enter lat lon: -39.23 317.76
col,row = -0.500440 -0.505930    status = 0
lat,lon = -39.230000 -42.240000    status = 1
enter lat lon: -39.23 42.24
col,row = 315.500440 -0.505930    status = 0
lat,lon = -39.230000 42.240000    status = 1
enter lat lon: -41.45 135.00
col,row = 315.488839 331.488839    status = 1
lat,lon = -41.450000 135.000000    status = 1
enter lat lon: -41.45 225.00
col,row = -0.488839 331.488839    status = 1
lat,lon = -41.450000 -135.000000    status = 1
enter lat lon:

inverse_grid:
enter col row: 0 0
lat,lon = -39.364869 -42.232570    status = 1
col,row = 0.000000 0.000000    status = 1
enter col row: 157.5 173.5
lat,lon = -90.000000 0.000000    status = 1
col,row = 157.500000 173.500000    status = 1