PS_Corner_Points

{
 "metadata": {
  "name": "Verify Center Points-South"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "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."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "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."
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Contents of S3B.gpd and Sps.mpp"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "    Sps.mpp\t\t/* map projection parameters\t# SSM/I Brightness Temp Grids */\n",
      "    316 332\t\t/* columns rows\t\t\t# South Polar Region */\n",
      "    4 4\t\t/* grid cells per map unit\t\t# 25km */\n",
      "    157.5 173.5\t/* map origin column, row\t\t# */\n",
      "\n",
      "    PolarStereographicEllipsoid\n",
      "    -90.0 0.0 -70.0\t\t/* lat0 lon0 lat1 */\n",
      "    0.0\t\t\t/* rotation */\n",
      "    100.0\t   \t\t/* scale (km/map unit) */\n",
      "    -90.00\t0.00\t\t/* center lat lon */\n",
      "    -90.00\t-20.00\t\t/* lat min max */\n",
      "    -180.00\t180.00\t\t/* lon min max */\n",
      "    10.00 15.00\t\t/* grid */\n",
      "    0.00\t0.00\t\t/* label lat lon */\n",
      "    1 0 0\t\t\t/* cil bdy riv */\n",
      "    6378.273\t\t/* Earth equatorial radius (km) */\n",
      "    0.081816153\t\t/* eccentricity */"
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "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."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def translate_grid_to_mapcoords(r, s):\n",
      "   # From S3B.gpd\n",
      "    map_origin_col = 157.5 # grid \n",
      "    map_origin_row = 173.5\n",
      "    cols_per_map_unit = 4\n",
      "    rows_per_map_unit = 4\n",
      "    \n",
      "    # From Sps.mpp\n",
      "    scale = 100.0 #km/map unit\n",
      "    \n",
      "    # From mapx/mapx.c/inverse_mapx\n",
      "    # u,v -> rotated and translated map coordinates in map units\n",
      "    u = (r - map_origin_col) / cols_per_map_unit\n",
      "    v = -(s - map_origin_row) / rows_per_map_unit\n",
      "    u *= scale\n",
      "    v *= scale\n",
      "    print \"u (x) = {u}  v (y) = {v}\".format(u= u, v= v)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Convert the origin as defined in the gpd to show we are using the correct calculations."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print \"Origin (157.5, 173.5):\"\n",
      "translate_grid_to_mapcoords(157.5, 173.5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Origin (157.5, 173.5):\n",
        "u (x) = 0.0  v (y) = -0.0\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "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."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print \"Upper left grid point: (0, 0):\"\n",
      "translate_grid_to_mapcoords(0, 0)\n",
      "print \"Upper right grid point: (315.0, 0):\"\n",
      "translate_grid_to_mapcoords(315.0, 0)\n",
      "print \"Lower right grid point: (315.0, 331.0):\"\n",
      "translate_grid_to_mapcoords(315.0, 331.0)\n",
      "print \"Lower left rid point: (0, 331.0):\"\n",
      "translate_grid_to_mapcoords(0, 331.0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Upper left grid point: (0, 0):\n",
        "u (x) = -3937.5  v (y) = 4337.5\n",
        "Upper right grid point: (315.0, 0):\n",
        "u (x) = 3937.5  v (y) = 4337.5\n",
        "Lower right grid point: (315.0, 331.0):\n",
        "u (x) = 3937.5  v (y) = -3937.5\n",
        "Lower left rid point: (0, 331.0):\n",
        "u (x) = -3937.5  v (y) = -3937.5\n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Notice the above conversions using the corners of the grid as whole numbered rows and columns yields map coordinates that differ from those listed on the PS page by 12.5 which is 1/2 the size of a map unit (25km).  Therefore the grid points in the data itself are in the CENTER of the grid cells."
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "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. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print \"Upper left grid point: (-0.5, -0.5):\"\n",
      "translate_grid_to_mapcoords(-0.5, -0.5)\n",
      "print \"Upper right grid point: (315.5, -0.5):\"\n",
      "translate_grid_to_mapcoords(315.5, -0.5)\n",
      "print \"Lower right grid point: (315.5, 331.5):\"\n",
      "translate_grid_to_mapcoords(315.5, 331.5)\n",
      "print \"Lower left grid point: (-0.5, 331.5):\"\n",
      "translate_grid_to_mapcoords(-0.5, 331.5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Upper left grid point: (-0.5, -0.5):\n",
        "u (x) = -3950.0  v (y) = 4350.0\n",
        "Upper right grid point: (315.5, -0.5):\n",
        "u (x) = 3950.0  v (y) = 4350.0\n",
        "Lower right grid point: (315.5, 331.5):\n",
        "u (x) = 3950.0  v (y) = -3950.0\n",
        "Lower left grid point: (-0.5, 331.5):\n",
        "u (x) = -3950.0  v (y) = -3950.0\n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Another way to see this is to convert geographic coordinates to grid rows and columns using gtest mapx utility."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "    hwilcox@hwilcox:~/nsidc0080/mapx> gtest\n",
      "\n",
      "    enter .gpd file name: S3B.gpd   \n",
      "    > assuming old style fixed format file\n",
      "\n",
      "    gpd: S3B.gpd\n",
      "    mpp:Sps.mpp\n",
      "\n",
      "    forward_grid:\n",
      "    enter lat lon: -39.23 317.76\n",
      "    col,row = -0.500440 -0.505930    status = 0\n",
      "    lat,lon = -39.230000 -42.240000    status = 1\n",
      "    enter lat lon: -39.23 42.24\n",
      "    col,row = 315.500440 -0.505930    status = 0\n",
      "    lat,lon = -39.230000 42.240000    status = 1\n",
      "    enter lat lon: -41.45 135.00\n",
      "    col,row = 315.488839 331.488839    status = 1\n",
      "    lat,lon = -41.450000 135.000000    status = 1\n",
      "    enter lat lon: -41.45 225.00\n",
      "    col,row = -0.488839 331.488839    status = 1\n",
      "    lat,lon = -41.450000 -135.000000    status = 1\n",
      "    enter lat lon: \n",
      "\n",
      "    inverse_grid:\n",
      "    enter col row: 0 0 \n",
      "    lat,lon = -39.364869 -42.232570    status = 1\n",
      "    col,row = 0.000000 0.000000    status = 1\n",
      "    enter col row: 157.5 173.5\n",
      "    lat,lon = -90.000000 0.000000    status = 1\n",
      "    col,row = 157.500000 173.500000    status = 1"
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "You can see from these conversions that the geographic coordinates given in the table result in the column and row grid values associated with the outer edges of the grid cells."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}