MODIS geolocation

Geolocation vs Data fields grid

250m vs 1km

For example in the MYD02QKM products, the geolocation is given at 1km resolution and the reflective bands at 250m. For a given 1km pixel, assuming that :
    i_tk_1km is the index of a 1km pixel along track
    i_sc_1km is the index of a 1km pixel along scan ( ie cross track ) 
The indexes of this pixel in the 250m grid is given by the relation :
i_tk_250m = 1.5 + 4 * i_tk_1km
i_sc_250m =       4 * i_sc_1km
A graphical representation of this relation is :
Modis 250m to 1km grid-scan Modis 250m to 1km grid-scan
Modis 250 to 1km

500m vs 1km

In the MYD02HKM, the geolocation is given at 1km resolution and the reflective bands at 500m. The indexes of a 1km pixel in the 250m grid is given by the relation :
i_tk_500m = 0.5 + 2 * i_tk_1km
i_sc_500m =       2 * i_sc_1km
A graphical representation of this relation is :
Modis 1 km to 500mModis 1km to 500m
Modis 1km to 500m legend

1km vs 5km

In the MYD021KM, MYD05 and MYD06 products, the geolocation is given at 5km resolution while some data fields are given at 1km. The indexes of a 5km pixel in the 1km grid is given by the simple relation :
i_tk_1km = 2 + 5 * i_tk_5km
i_sc_1km = 2 + 5 * i_sc_5km

Be careful, the cross-track dimension is composed of 270 5km pixels for the MYD021KM product but only of 269 pixels for the MYD05 and MYD06 products.
A graphical representation of this relation is :

MYD021KM
Modis 5km to 1km Modis 5km to 1km
MYD05 - MYD06
Modis 5km to 1km legend

Subpixel interpolation strategy

To build the fine resolution grid, the basic question is to obtain the ( Latitude, Longitude ) of a [ i_tk_fine, i_sc_fine ] pixel. You will find at this end of this chapter two python scripts that do the job for the 5km -> 1km and 1km -> 250m grids. Here is how it is implemented in those scripts ( rem : the same method is used in CALTRACK products )
In the following, we will take the 5km -> 1km case as an example, and i_tk and i_sc stand respectively for along-track and along-scan indexes.

Bounding coarse resolution pixels

First we search the 4 coarse resolution pixels P1 to P4 that bound the fine resolution one P.
REM: in border of scan, the fine pixel can be outside of the coarse pixel footprint. In this case, its position will be extrapolated to the nearest coarse pixels.

Coarse to fine grid interpolation strategy

Interpolation along-track

Here, P' and P" are the 1km pixels on the along-track line than the 5km pixels and on the same along-scan line than P. Their ( lat, long ) positions are interpolated using the positions of P1/P4 and P2/P3.
One reasonable approximation done here is that 2 successive scan lines are parallels.

Coarse to fine grid interpolation strategy 2

Source Code

Here is 2 python scripts that manage the process described above :
  • one for computing a 250m pixel geolocation based on 1km datasets, as in MYD02QKM products
  • one for 1km based on 5km, as in MYD05, MYD06 and MYD021KM
Both are printed below, but are also contained in the downloadable package : Modis_coarse_to_fine_geolocation.v1.0.0.tar.gz

1km -> 250m

 # -*- coding: utf-8 -*-
 
 """
 Description :
     Give the geolocation of a MODIS 250m pixel specified by its indices [i_250m_along_track, i_250m_cross_track] for products where the 
     (latitude, longitude) datasets are given at 1km resolution
     For the interpolation of the 250m positions, it takes care of borders of scan. In this cases, the position will be extrapolated.
     It also deals with the change date meridian crossing
 
 Usage :
     python modis_250m_subpixel_geolocation.py <infile> <i_250m_along_track> <i_250m_cross_track>
 
     It will print out the ( lat, lon ) geolocation of the given pixel
 
     Optionally, this script can also display a figure of the pixels around the requested one.
 
 Prerequisites :
     [REQUIRED]
         - python >= 2.5
         - numpy
         - pyhdf
 
     [OPTIONAL]
         When enabling the plot of the grid ( by setting __DEBUG__ to True ), you need to have installed the matplotlib library 
         with a fonctionnal graphical toolkit support
 
 Author :
     CGTD-ICARE/UDEV Nicolas PASCAL ( nicolas.pascal-at-icare.univ-lille1.fr  )
 
 License :
     This file must be used under the terms of the CeCILL.
     This source file is licensed as described in the file COPYING, which
     you should have received as part of this distribution.  The terms
     are also available at
     http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 
 History :
 --------
 v0.1.0 : 2009/11/06
  - creation
 """
 
 import math
 import sys
 
 import warnings
 warnings.filterwarnings("ignore")
 
 import numpy as np
 from pyhdf import SD
 
 __DEBUG__ = False
 if __DEBUG__ is True :
     import matplotlib.pyplot as plt
     import matplotlib.ticker as ticker
 
 # --- 1km -> 250m ---
 sz_sc_1km  = 1354
 sz_sc_250m = 5416
 
 def itk_1km_to_250m ( i_tk_1km ) :
     """
     return the 250m grid index along track of a 1km pixel
     """
     return 1.5 + 4. * i_tk_1km
 
 def itk_250m_to_1km ( i_tk_250m ) :
     """
     return the 1km grid index along track of a 250m pixel
     """
     return ( i_tk_250m - 1.5 ) / 4.
 
 def isc_1km_to_250m ( i_sc_1km ) :
     """
     return the 250m grid index cross track of a 1km pixel
     """
     return 4. * i_sc_1km
 
 def isc_250m_to_1km ( i_sc_250m ) :
     """
     return the 1km grid index cross track of a 250m pixel
     """
     return i_sc_250m / 4.
 
 def get_250m_pix_pos ( itk_250m, isc_250m, lat_1km, lon_1km ) :
     """
     return the (lat,lon) of a 250m pixel specified with its indexes, when the geolocation datasets are given at 1km resolution
     @param itk_250m grid index of the 250m pixel along-track
     @param isc_250m grid index of the 250m pixel cross-track
     @param lat_1km latitudes dataset at 1km resolution
     @param lon_1km longitudes dataset at 1km resolution
     @return the ( lat, lon ) of the 250m pixel  [ itk_250m, isc_250m ]
     """
     # change date meridian case sentinel
     is_crossing_change_date = False
 
     # check 250m indexes validity
     sz_tk_1km = lat_1km.shape [0]
     sz_tk_250m = 4 * sz_tk_1km
     if ( isc_250m < 0 ) or ( isc_250m > sz_sc_250m - 1 ) :
         raise ValueError ( "Invalid scan index %d. Must be in range [%d,%d]"%(isc_250m,0,sz_sc_250m-1) )
     if ( itk_250m < 0 ) or ( itk_250m > sz_tk_250m - 1 ) :
         raise ValueError ( "Invalid track index %d. Must be in range [%d,%d]"%(itk_250m,0,sz_tk_250m-1) )
 
     # --- set the bounding 1km pixels to take for interpolation
     # set the (track,scan) indexes of the 1km pixel in the 250m grid
     itk_1km = itk_250m_to_1km ( itk_250m )
     isc_1km = isc_250m_to_1km ( isc_250m )
 
     #print "i_250m=[%.2f, %.2f] -> i_1km=[%.2f, %.2f]"%(itk_250m, isc_250m, itk_1km,isc_1km)
 
     # the width of one scan, in number of pixels, at 250m resolution
     w_scan_250m = 40
     # - extra/interpolation 1km pixels along track
     if (   itk_250m % w_scan_250m ) <= 1 : # extrapolate start of scan
         itk_top_1km    = math.ceil ( itk_1km )
         itk_bottom_1km = itk_top_1km + 1
     elif ( itk_250m % w_scan_250m ) >= 38 : # extrapolate end of scan
         itk_top_1km    = math.floor ( itk_1km )
         itk_bottom_1km = itk_top_1km - 1
     else : # general case : middle of scan
         itk_top_1km    = math.floor ( itk_1km )
         itk_bottom_1km = itk_top_1km + 1
     # - extra/interpolation 1km pixels along track
     if ( isc_1km >= 1353. ) : # extrapolate end of scan line
         isc_left_1km  = math.floor ( isc_1km ) - 1
         isc_right_1km = math.floor ( isc_1km )
     else : # general case : interpolation
         isc_left_1km  = math.floor ( isc_1km )
         isc_right_1km = isc_left_1km + 1
 
     #print "itk_top_1km=%d itk_bottom_1km=%d isc_left_1km=%d isc_right_1km=%d"%(itk_top_1km, itk_bottom_1km, isc_left_1km, isc_right_1km)
 
     # --- set the 1km track lines position ; left border ---
     lat_left_250m, lon_left_250m   = get_y_pos_1km_to_250m (
                 isc_left_1km, itk_top_1km, itk_bottom_1km,
                 lat_1km, lon_1km,
                 itk_250m )
     # --- set the 1km track lines position ; right border ---
     lat_right_250m, lon_right_250m = get_y_pos_1km_to_250m (
                 isc_right_1km, itk_top_1km, itk_bottom_1km,
                 lat_1km, lon_1km,
                 itk_250m )
 
     #print "left_250m=[%f,%f] right_250m=[%f,%f]"%(lat_left_250m, lon_left_250m, lat_right_250m, lon_right_250m)
 
     # check for change date meridian case
     if  abs ( lon_right_250m  - lon_left_250m  ) > 180. :
         is_crossing_change_date = True
         # all negative longitudes will be incremented of 360 before interpolation
         if lon_left_250m < 0. :
             lon_left_250m += 360.
         if lon_right_250m < 0. :
             lon_right_250m += 360.
 
     # for each track line position, interpolate along scan to retrieve the 250m geolocation
     lat, lon = get_x_pos_1km_to_250m ( lat_left_250m,  lon_left_250m,  isc_left_1km,
                                      lat_right_250m, lon_right_250m, isc_right_1km,
                                      isc_250m )
     #print "geolocation = [%f, %f]"%(lat,lon)
 
     # in case of change date crossing, turn values > 180. to negative
     if lon > 180.:
         lon -= 360.
     elif lon < -180.:
         lon += 360.
     return lat, lon
 
 def get_x_pos_1km_to_250m ( lat_left_250m,  lon_left_250m,  isc_left_1km,
                             lat_right_250m, lon_right_250m, isc_right_1km,
                             isc_250m
                         ) :
     """
     retrieve the position of a 250m pixel set by its cross-track index and the position of 2 along-track lines set by 2 sucessive 1km
     pixels
     PRECONDITION : change date meridian case must have already been treated, and so lon_left_250m or lon_right_250m can have values > 180.
     @param lat_left_250m latitude of the 250m along track pixel that is the intersection between cross-track line at index isc_250m and
     the along-track line that joins the 2 sucessive bounding 1km pixels on the left of the point to interpolate
     @param lon_left_250m longitude of the 250m along track pixel that is the intersection between cross-track line at index isc_250m and
     the along-track line that joins the 2 sucessive bounding 1km pixels on the left of the point to interpolate
     @param isc_left_1km cross-track index at 1km of the left border
     @param lat_right_250m latitude of the 250m along track pixel that is the intersection between cross-track line at index isc_250m and
     the along-track line that joins the 2 sucessive bounding 1km pixels on the right of the point to interpolate
     @param lon_right_250m longitude of the 250m along track pixel that is the intersection between cross-track line at index isc_250m and
     the along-track line that joins the 2 sucessive bounding 1km pixels on the right of the point to interpolate
     @param isc_right_1km cross-track index at 1km of the right border
     @return (lat,lon) of the pixel at [itk_250m, isc_1km], interpolated between the 2 successive 1km pixels p1 and p2
     """
     # make sure P1 and P1 are 2 successive pixels along-track
     if abs ( isc_left_1km - isc_right_1km ) != 1 :
         raise ValueError ( "The 2 borders are on the same along-track line and must be successive" )
 
     # coordinates of left and right border in the 250m grid
     isc_left_250m  = isc_1km_to_250m ( isc_left_1km  )
     isc_right_250m = isc_1km_to_250m ( isc_right_1km )
 
     #print "isc_left_1km=%d isc_right_1km=%d"%(isc_left_1km, isc_right_1km)
     #print "isc_left_250m=%d isc_right_250m=%d"%(isc_left_250m, isc_right_250m)
     #print "isc_250m_min=%d isc_250m_max=%d"%(isc_250m_min,isc_250m_max)
 
     # linear interpolation on the position
     alpha_lon = ( lon_right_250m - lon_left_250m ) / ( isc_right_250m - isc_left_250m )
     alpha_lat = ( lat_right_250m - lat_left_250m ) / ( isc_right_250m - isc_left_250m )
 
     lat = lat_left_250m + alpha_lat * ( isc_250m - isc_left_250m )
     lon = lon_left_250m + alpha_lon * ( isc_250m - isc_left_250m )
 
     return lat, lon
 
 def get_y_pos_1km_to_250m ( isc_1km, itk_p1_1km, itk_p2_1km,
                             lat_1km, lon_1km,
                             itk_250m ) :
     """
     return the position of 250m pixel defined by its along-track index between 2 successive 1km pixels on a same cross-track line
     @warning If the 2 pixels are crossing the changing date meridian, the negative longitude will be returned
     with an increment of 360 to manage correctly the interpolation
     @param isc_1km index of the cross-track line in the 1km grid
     @param itk_p1_1km index of the first along-track pixel in the 1km grid
     @param itk_p2_1km index of the second along-track pixel in the 1km grid
     @param lat_1km latitudes dataset at 1km resolution
     @param lon_1km longitudes dataset at 1km resolution
     @param itk_250m index of the along-track 250m along-track pixel to interpolate
     @return (lat,lon) of the pixel at [itk_250m, isc_1km], interpolated between the 2 successive 1km pixels p1 and p2
     """
     # make sure P1 and P1 are 2 successive pixels along-track
     if abs ( itk_p1_1km - itk_p2_1km ) != 1 :
         raise ValueError ( "P1 and P2 are on the same cross-track line and must be successive" )
 
     # lat, lon of the 1km bounding pixels
     p1_1km_lat = lat_1km [ itk_p1_1km, isc_1km ]
     p1_1km_lon = lon_1km [ itk_p1_1km, isc_1km ]
     p2_1km_lat = lat_1km [ itk_p2_1km, isc_1km ]
     p2_1km_lon = lon_1km [ itk_p2_1km, isc_1km ]
 
     # check for change date meridian particular case
 
     # change date meridian case sentinel
     if abs ( p1_1km_lon - p2_1km_lon ) > 180. :
         if p1_1km_lon < 0. :
             p1_1km_lon += 360.
         elif p2_1km_lon < 0. :
             p2_1km_lon += 360.
 
     # coordinates of p1, p2 in the 250m grid
     itk_p1_250m = itk_1km_to_250m ( itk_p1_1km )
     itk_p2_250m = itk_1km_to_250m ( itk_p2_1km )
 
     #print "itk_p1_250m=%f itk_p2_250m=%f"%( itk_p1_250m, itk_p2_250m )
 
     # linear interpolation on the position
     alpha_lon = ( p2_1km_lon - p1_1km_lon ) / ( itk_p2_250m - itk_p1_250m )
     alpha_lat = ( p2_1km_lat - p1_1km_lat ) / ( itk_p2_250m - itk_p1_250m )
     lon = p1_1km_lon + alpha_lon * ( itk_250m - itk_p1_250m )
     lat = p1_1km_lat + alpha_lat * ( itk_250m - itk_p1_250m )
 
     # in case of change date crossing, turn values > 180. to negative
     if lon > 180. :
         lon -= 360.
     elif lon < -180. :
         lon += 360.
 
     return lat, lon
 
 def get_1km_pix_to_plot ( itk_250m, isc_250m, lat_1km, lon_1km ) :
     """
     return the arrays of (lat,lon) and of [i,j] labels of the 1km pixelsto draw arounf the (itk_250m, isc_250m) pixel
     """
     v_itk_1km = [ ]
     v_isc_1km = [ ]
     sc_250m_width = 40
     itk_1km = round ( itk_250m_to_1km ( itk_250m ) )
     isc_1km = max ( 0, min ( sz_sc_1km - 1, round ( isc_250m_to_1km ( isc_250m ) ) ) )
 
     if ( itk_250m % sc_250m_width ) <= 2 : # start of scan
         #itk_1km = round ( itk_1km )
         if isc_250m <= 2 :
             v_itk_1km = [   itk_1km,     itk_1km,     itk_1km,
                             itk_1km + 1, itk_1km + 1, itk_1km + 1,
                             itk_1km + 2, itk_1km + 2, itk_1km + 2 ]
             v_isc_1km = [  isc_1km,     isc_1km + 1, isc_1km + 2,
                             isc_1km,     isc_1km + 1, isc_1km + 2,
                             isc_1km,     isc_1km + 1, isc_1km + 2 ]
         elif isc_250m >= 5410 :
             v_itk_1km = [  itk_1km,     itk_1km    , itk_1km,
                             itk_1km + 1, itk_1km + 1, itk_1km + 1,
                             itk_1km + 2, itk_1km + 2, itk_1km + 2 ]
             v_isc_1km = [  isc_1km,     isc_1km - 1, isc_1km - 2,
                             isc_1km,     isc_1km - 1, isc_1km - 2,
                             isc_1km,     isc_1km - 1, isc_1km - 2 ]
         else :
             v_itk_1km = [  itk_1km,     itk_1km    , itk_1km,
                             itk_1km + 1, itk_1km + 1, itk_1km + 1,
                             itk_1km + 2, itk_1km + 2, itk_1km + 2 ]
             v_isc_1km = [  isc_1km,     isc_1km - 1, isc_1km + 1,
                             isc_1km,     isc_1km - 1, isc_1km + 1,
                             isc_1km,     isc_1km - 1, isc_1km + 1 ]
     elif ( itk_250m % sc_250m_width ) >= 37 : # end of scan
         #itk_1km = round ( itk_1km )
         if isc_250m <= 2 :
             v_itk_1km = [  itk_1km,     itk_1km    , itk_1km,
                            itk_1km - 1, itk_1km - 1, itk_1km - 1,
                            itk_1km - 2, itk_1km - 2, itk_1km - 2 ]
             v_isc_1km = [  isc_1km,     isc_1km + 1, isc_1km + 2,
                            isc_1km,     isc_1km + 1, isc_1km + 2,
                            isc_1km,     isc_1km + 1, isc_1km + 2 ]
         elif isc_250m >= 5410 :
             v_itk_1km = [  itk_1km,     itk_1km    , itk_1km,
                             itk_1km - 1, itk_1km - 1, itk_1km - 1,
                             itk_1km - 2, itk_1km - 2, itk_1km - 2 ]
             v_isc_1km = [  isc_1km,     isc_1km - 1, isc_1km - 2,
                            isc_1km,     isc_1km - 1, isc_1km - 2,
                            isc_1km,     isc_1km - 1, isc_1km - 2 ]
         else :
             v_itk_1km = [  itk_1km,     itk_1km    , itk_1km,
                             itk_1km - 1, itk_1km - 1, itk_1km - 1,
                             itk_1km - 2, itk_1km - 2, itk_1km - 2 ]
             v_isc_1km = [  isc_1km,     isc_1km - 1, isc_1km + 1,
                             isc_1km,     isc_1km - 1, isc_1km + 1,
                             isc_1km,     isc_1km - 1, isc_1km + 1 ]
     else : # middle of scan
         #itk_1km = round ( itk_1km )
         if isc_1km == 0 :
             v_itk_1km = [  itk_1km,     itk_1km    , itk_1km,
                             itk_1km - 1, itk_1km - 1, itk_1km - 1,
                             itk_1km + 1, itk_1km + 1, itk_1km + 1 ]
             v_isc_1km = [  isc_1km,     isc_1km + 1, isc_1km + 2,
                             isc_1km,     isc_1km + 1, isc_1km + 2,
                             isc_1km,     isc_1km + 1, isc_1km + 2 ]
         elif isc_1km == 1353 :
             v_itk_1km = [  itk_1km,     itk_1km    , itk_1km,
                             itk_1km - 1, itk_1km - 1, itk_1km - 1,
                             itk_1km + 1, itk_1km + 1, itk_1km + 1 ]
             v_isc_1km = [  isc_1km,     isc_1km - 1, isc_1km - 2,
                             isc_1km,     isc_1km - 1, isc_1km - 2,
                             isc_1km,     isc_1km - 1, isc_1km - 2 ]
         else :
             v_itk_1km = [  itk_1km,     itk_1km    , itk_1km,
                             itk_1km - 1, itk_1km - 1, itk_1km - 1,
                             itk_1km + 1, itk_1km + 1, itk_1km + 1 ]
             v_isc_1km = [  isc_1km,     isc_1km - 1, isc_1km + 1,
                             isc_1km,     isc_1km - 1, isc_1km + 1,
                             isc_1km,     isc_1km - 1, isc_1km + 1 ]
     npix = len(v_itk_1km)
     # build the labels array
     v_label = [ "[%d,%d]"%(v_itk_1km[i],v_isc_1km[i]) for i in xrange ( npix ) ]
 
     return v_itk_1km, v_isc_1km, v_label
 
 def get_bounds_1km_to_250m( itk_1km, isc_1km ) :
     """
     return the 250m pixel indexes limits in the 1km pixel [ itk_1km, isc_1km ] footprint
     """
 
     # set the (track,scan) indexes of the 1km pixel in the 1km grid
     itk_250m = itk_1km_to_250m ( itk_1km )
     isc_250m = isc_1km_to_250m ( isc_1km )
 
     # set the 250m indexes of pixels to interpolate along track
     itk_250m_min    =  int ( itk_250m - 1.5 )
     itk_250m_max    =  int ( itk_250m + 1.5 )
 
     # general case : 2 interpolations done along scan : [isc-1, isc] then [isc, isc+1]
     isc_250m_min  = isc_250m - 2
     isc_250m_max  = isc_250m + 2
     if   ( isc_1km == 0 ) :
         isc_250m_min = 0
     elif ( isc_1km >= sz_sc_1km - 1 ) :
         isc_250m_max  = isc_250m + 3
 
     #print itk_1km, itk_250m_min, itk_250m_max
     #print isc_1km, isc_250m_min, isc_250m_max
 
     return itk_250m_min, itk_250m_max, isc_250m_min, isc_250m_max
 
 ##############  MAIN  ###############
 def main():
     """
     Main script
     """
     if len ( sys.argv ) != 4 :
         raise ValueError ( "Invalid number of arguments.\nUsage : python modis_1km_to_250m_geolocation.py <infile> <itrack_250m> <iscan_250m>" )
 
     infile     =       sys.argv [ 1 ]
     itk_250m   = int ( sys.argv [ 2 ] )
     isc_250m   = int ( sys.argv [ 3 ] )
     #print itk_250m, isc_250m
     #print itk_250m_to_1km ( itk_250m )
     #itk_1km = max ( 0, math.floor ( itk_250m_to_1km ( itk_250m ) ) + 1 )
     #isc_1km = math.floor ( isc_250m_to_1km ( isc_250m ) )
     #print itk_1km, isc_1km
 
     # -----------------------------------
     # --- Load the 1km (lat,lon) data ---
     # -----------------------------------
     hdf  = SD.SD ( infile )
 
     # read the latitudes of the points
     sds  = hdf.select ( "Latitude" )
     lat_1km = sds.get ( )
     sds.endaccess()
 
     # read the longitudes of the points
     sds  = hdf.select ( "Longitude" )
     lon_1km = sds.get ( )
     sds.endaccess()
 
     hdf.end()
 
     # ---------------------------------------------
     # --- compute the (lat,lon) of a 250m pixel ---
     # ---------------------------------------------
     lat, lon = get_250m_pix_pos ( itk_250m, isc_250m, lat_1km, lon_1km )
     print "%f\t%f"%(lat,lon)
 
     if __DEBUG__ :
         # ----------------------------------
         # --- set the 1km pixels to plot ---
         # ----------------------------------
         sz_tk_1km, sz_sc_1km = lat_1km.shape
         #if itk_1km >= sz_tk_1km :
             #itk_1km = sz_tk_1km - 1
         #if isc_1km >= sz_sc_1km :
             #isc_1km = sz_sc_1km - 1
         v_itk_1km, v_isc_1km, v_label = get_1km_pix_to_plot ( itk_250m, isc_250m, lat_1km, lon_1km )
 
         # set the lat,lon of those points
         v_lat_1km = lat_1km [ v_itk_1km, v_isc_1km ]
         v_lon_1km = lon_1km [ v_itk_1km, v_isc_1km ]
 
         # ----------------------------------
         # --- plot 1km pixels
         # ----------------------------------
         # - init the plot where will be drawn the pixels
         fig = plt.figure ()
         ax  = plt.subplot (111)
 
         plt.xlabel ('Longitude' )
         plt.ylabel ('Latitude' )
         plt.axis   ( 'equal' )
         ax.get_xaxis().set_major_formatter( ticker.FormatStrFormatter('%.2f') )
         ax.get_yaxis().set_major_formatter( ticker.FormatStrFormatter('%.2f') )
 
         plt.axis   ( [ v_lon_1km.min(), v_lon_1km.max(), v_lat_1km.min(), v_lat_1km.max() ] )
 
         #print "v_lon_1km "+str(v_lon_1km)
         #print "v_lat_1km "+str(v_lat_1km)
 
         # - plot the 250m pixel
         plt.scatter ( [ lon ], [ lat ], color='m', marker='x', linewidth = 1 )
         ax.text     ( lon, lat, "(%d,%d)"%( itk_250m, isc_250m ), color = 'm' )
 
         # - plot the 1km pixels
         npix = len ( v_lat_1km )
         for i in xrange ( npix ) :
             _lat   = v_lat_1km [ i ]
             _lon   = v_lon_1km [ i ]
             label = v_label   [ i ]
             sz_marker = 1
             if i == 0 :
                 # set a biggest marker for pixel to interpolate
                 sz_marker = 5
             plt.scatter ( [ _lon ], [ _lat ], color='r', marker='o', linewidth = sz_marker, label = "_" )
             ax.text     ( _lon, _lat, label )
 
         # ------------------------------------------
         # --- Compute and plot inner 250m pixels ---
         # ------------------------------------------
         itk_1km = round ( itk_250m_to_1km ( itk_250m ) )
         isc_1km = min ( sz_sc_1km - 1, round ( isc_250m_to_1km ( isc_250m ) ) )
 
         itk_250m_min, itk_250m_max, isc_250m_min, isc_250m_max = get_bounds_1km_to_250m( itk_1km, isc_1km )
         v_lat_250m = []
         v_lon_250m = []
         v_itk_250m = []
         v_isc_250m = []
 
         for itk_250m in xrange ( itk_250m_min, itk_250m_max + 1 ) :
             for isc_250m in xrange ( isc_250m_min, isc_250m_max + 1 ) :
                 v_itk_250m.append ( itk_250m )
                 v_isc_250m.append ( isc_250m )
                 # intert/extrapolate the lat, lon of this 250m pixel using the 1km grid
                 lat, lon = get_250m_pix_pos ( itk_250m, isc_250m, lat_1km, lon_1km )
                 v_lat_250m.append ( lat )
                 v_lon_250m.append ( lon )
         v_itk_250m = np.array ( v_itk_250m )
         v_isc_250m = np.array ( v_isc_250m )
         v_lat_250m = np.array ( v_lat_250m )
         v_lon_250m = np.array ( v_lon_250m )
 
         plt.scatter ( v_lon_250m [ v_lon_250m != -9999. ].flat, v_lat_250m [ v_lat_250m != -9999. ].flat, \
                     color='b', marker='+', linewidth = 1, label = "_" )
 
         plt.show ()
 
 if __name__ == "__main__":
   main()
 
5km -> 1km
 
 # -*- coding: utf-8 -*-
 
 """
 Description :
     Give the geolocation of a MODIS 1km pixel specified by its indices [i_1km_along_track, i_1km_cross_track] for products where the 
    (latitude, longitude) datasets are given at 5km resolution
     For the interpolation of the 1km positions, it takes care of borders of scan. In this cases, the position will be extrapolated.
     It also deals with the change date meridian crossing.
 
 Limitation :
     This script works well with MYD05 and MYD06 products but needs a minor adaptations for fully running on MY021KM products.
     Due to the difference of along-scan dimension size ( 271 pixels for MYD021KM and 270 for others ), the pixels with a 1km scan index
     in range [1349,1352] are always extrapolated using the 5km scan pixels [268,269].whereas they can be interpolated between [269,270] in
     the MYD021KM case.
 
 Usage :
     python modis_5km_to_1km_geolocation.py <infile> <i_1km_along_track> <i_1km_cross_track>
 
     It will print out the ( lat, lon ) geolocation of the given pixel.
 
     Optionally, this script can also display a figure of the pixels around the requested one.
 
 Prerequisites :
     [REQUIRED]
         - python >= 2.5
         - numpy
         - pyhdf
 
     [OPTIONAL]
         When enabling the plot of the grid ( by setting __DEBUG__ to True ), you need to have installed the matplotlib library with 
         a fonctionnal graphical toolkit support
 
 Author :
     CGTD-ICARE/UDEV Nicolas PASCAL ( nicolas.pascal-at-icare.univ-lille1.fr  )
 
 License :
     This file must be used under the terms of the CeCILL.
     This source file is licensed as described in the file COPYING, which
     you should have received as part of this distribution.  The terms
     are also available at
     http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 
 History :
 --------
 v0.1.0 : 2009/11/25
  - creation
 """
 
 import math
 import sys
 
 import warnings
 warnings.filterwarnings("ignore")
 
 import numpy as np
 from pyhdf import SD
 
 __DEBUG__ = False
 if __DEBUG__ is True :
     import matplotlib.pyplot as plt
     import matplotlib.ticker as ticker
 
 # --- 5km -> 1km ---
 sz_sc_5km = 270 # number of 5km pixels along scan
 
 def itk_5km_to_1km ( i_tk_5km ) :
     """
     return the 1km grid index along track of a 5km pixel
     """
     return 2 + 5 * i_tk_5km
 
 def itk_1km_to_5km ( i_tk_1km ) :
     """
     return the 5km grid index along track of a 1km pixel
     """
     return ( i_tk_1km - 2. ) / 5.
 
 def isc_5km_to_1km ( i_sc_5km ) :
     """
     return the 1km grid index cross track of a 5km pixel
     """
     return 2. + 5. * i_sc_5km
 
 def isc_1km_to_5km ( i_sc_1km ) :
     """
     return the 5km grid index cross track of a 1km pixel
     """
     return  ( i_sc_1km - 2. ) / 5.
 
 def get_1km_pix_pos ( itk_1km, isc_1km, lat_5km, lon_5km ) :
     """
     return the (lat,lon) of a 1km pixel specified with its indexes, when the geolocation datasets are given at 5km resolution
     @param itk_1km grid index of the 1km pixel along-track
     @param isc_1km grid index of the 1km pixel cross-track
     @param lat_5km latitudes dataset at 5km resolution
     @param lon_5km longitudes dataset at 5km resolution
     @return the ( lat, lon ) of the 1km pixel  [ itk_1km, isc_1km ]
     """
     # check 1km indexes validity
     sz_tk_5km = lat_5km.shape[0] + 1
     sz_tk_1km = 5 * sz_tk_5km
     sz_sc_1km = ( 5 * sz_sc_5km ) + 6
     if ( isc_1km < 0 ) or ( isc_1km > sz_sc_1km - 1 ) :
         raise ValueError ( "Invalid scan index %d. Must be in range [%d,%d]"%(isc_1km,0,sz_sc_1km-1) )
     if ( itk_1km < 0 ) or ( itk_1km > sz_tk_1km - 1 ) :
         raise ValueError ( "Invalid track index %d. Must be in range [%d,%d]"%(itk_1km,0,sz_tk_1km-1) )
 
     # --- set the bounding 5km pixels to take for interpolation
     # set the (track,scan) indexes of the 5km pixel in the 1km grid
     itk_5km = itk_1km_to_5km ( itk_1km )
     isc_5km = isc_1km_to_5km ( isc_1km )
 
     #print "i_1km=[%.2f, %.2f] -> i_5km=[%.2f, %.2f]"%(itk_1km, isc_1km, itk_5km, isc_5km)
 
     # the width of one scan, in number of pixels, at 1km resolution
     w_scan_1km = 10.
     # - extra/interpolation 5km pixels along track
     if ( itk_1km % w_scan_1km ) <= 2 : # extrapolate start of scan
         itk_top_5km    = math.ceil ( itk_5km )
         itk_bottom_5km = itk_top_5km + 1
     elif ( itk_1km % w_scan_1km ) >= 7 : # extrapolate end of scan
         itk_top_5km    = math.floor ( itk_5km )
         itk_bottom_5km = itk_top_5km - 1
     else : # general case : middle of scan
         itk_top_5km    = math.floor ( itk_5km )
         itk_bottom_5km = itk_top_5km + 1
     # - extra/interpolation 5km pixels cross track
     if ( isc_1km <= 2 ) : # extrapolate start of scan line
         isc_left_5km  = 0
         isc_right_5km = 1
     elif ( isc_5km >= ( sz_sc_5km - 1 )  ) : # extrapolate end of scan line
         isc_left_5km  = sz_sc_5km - 2
         isc_right_5km = sz_sc_5km - 1
     else : # general case : interpolation
         isc_left_5km  = math.floor ( isc_5km )
         isc_right_5km = isc_left_5km + 1
 
     #print "itk_top_5km=%d itk_bottom_5km=%d isc_left_5km=%d isc_right_5km=%d"%(itk_top_5km, itk_bottom_5km, isc_left_5km, isc_right_5km)
 
     # --- set the 5km track lines position ; left border ---
     lat_left_1km, lon_left_1km   = get_y_pos_5km_to_1km (
                 isc_left_5km, itk_top_5km, itk_bottom_5km,
                 lat_5km, lon_5km,
                 itk_1km )
     # --- set the 5km track lines position ; right border ---
     lat_right_1km, lon_right_1km = get_y_pos_5km_to_1km (
                 isc_right_5km, itk_top_5km, itk_bottom_5km,
                 lat_5km, lon_5km,
                 itk_1km )
 
     #print "left_1km=[%f,%f] right_1km=[%f,%f]"%(lat_left_1km, lon_left_1km, lat_right_1km, lon_right_1km)
 
     # check for change date meridian case
     if  abs ( lon_right_1km  - lon_left_1km  ) > 180. :
         # all negative longitudes will be incremented of 360 before interpolation
         if lon_left_1km < 0. :
             lon_left_1km += 360.
         if lon_right_1km < 0. :
             lon_right_1km += 360.
 
     # for each track line position, interpolate along scan to retrieve the 1km geolocation
     lat, lon = get_x_pos_5km_to_1km ( lat_left_1km,  lon_left_1km,  isc_left_5km,
                                      lat_right_1km, lon_right_1km, isc_right_5km,
                                      isc_1km )
     #print "geolocation = [%f, %f]"%(lat,lon)
 
     # in case of change date crossing, turn values > 180. to negative < -180 to positive
     if lon > 180. :
         lon = lon - 360.
     elif lon < -180. :
         lon = lon + 360.
     return lat, lon
 
 def get_x_pos_5km_to_1km ( lat_left_1km,  lon_left_1km,  isc_left_5km,
                            lat_right_1km, lon_right_1km, isc_right_5km,
                            isc_1km
                         ) :
     """
     retrieve the position of a 1km pixel set by its cross-track index and the position of 2 along-track lines set by 2 sucessive 5km
     pixels
     PRECONDITION : change date meridian case must have already been treated, and so lon_left_1km or lon_right_1km can have values > 180.
     @param lat_left_1km latitude of the 1km along track pixel that is the intersection between cross-track line at index isc_1km and
     the along-track line that joins the 2 sucessive bounding 5km pixels on the left of the point to interpolate
     @param lon_left_1km longitude of the 1km along track pixel that is the intersection between cross-track line at index isc_1km and
     the along-track line that joins the 2 sucessive bounding 5km pixels on the left of the point to interpolate
     @param isc_left_5km cross-track index at 5km of the left border
     @param lat_right_1km latitude of the 1km along track pixel that is the intersection between cross-track line at index isc_1km and
     the along-track line that joins the 2 sucessive bounding 5km pixels on the right of the point to interpolate
     @param lon_right_1km longitude of the 1km along track pixel that is the intersection between cross-track line at index isc_1km and
     the along-track line that joins the 2 sucessive bounding 5km pixels on the right of the point to interpolate
     @param isc_right_5km cross-track index at 5km of the right border
     @return (lat,lon) of the pixel at [itk_1km, isc_5km], interpolated between the 2 successive 5km pixels p1 and p2
     """
     # make sure P1 and P1 are 2 successive pixels along-track
     if abs ( isc_left_5km - isc_right_5km ) != 1 :
         raise ValueError ( "The 2 borders are on the same along-track line and must be successive" )
 
     # coordinates of left and right border in the 1km grid
     isc_left_1km  = isc_5km_to_1km ( isc_left_5km )
     isc_right_1km = isc_5km_to_1km ( isc_right_5km )
 
     #print "isc_left_5km=%d isc_right_5km=%d"%(isc_left_5km, isc_right_5km)
     #print "isc_left_1km=%d isc_right_1km=%d"%(isc_left_1km, isc_right_1km)
     #print "isc_1km_min=%d isc_1km_max=%d"%(isc_1km_min,isc_1km_max)
 
     # linear interpolation on the position
     alpha_lon = ( lon_right_1km - lon_left_1km ) / ( isc_right_1km - isc_left_1km )
     alpha_lat = ( lat_right_1km - lat_left_1km ) / ( isc_right_1km - isc_left_1km )
 
     lat = lat_left_1km + alpha_lat * ( isc_1km - isc_left_1km )
     lon = lon_left_1km + alpha_lon * ( isc_1km - isc_left_1km )
 
     return lat, lon
 
 def get_y_pos_5km_to_1km (  isc_5km, itk_p1_5km, itk_p2_5km,
                             lat_5km, lon_5km,
                             itk_1km ) :
     """
     return the position of 1km pixel defined by its along-track index between 2 successive 5km pixels on a same cross-track line
     @warning If the 2 pixels are crossing the changing date meridian, the negative longitude will be returned
     with an increment of 360 to manage correctly the interpolation
     @param isc_5km index of the cross-track line in the 5km grid
     @param itk_p1_5km index of the first along-track pixel in the 5km grid
     @param itk_p2_5km index of the second along-track pixel in the 5km grid
     @param lat_5km latitudes dataset at 5km resolution
     @param lon_5km longitudes dataset at 5km resolution
     @param itk_1km index of the along-track 1km along-track pixel to interpolate
     @return (lat,lon) of the pixel at [itk_1km, isc_5km], interpolated between the 2 successive 5km pixels p1 and p2
     """
     # make sure P1 and P1 are 2 successive pixels along-track
     if abs ( itk_p1_5km - itk_p2_5km ) != 1 :
         raise ValueError ( "P1 and P2 are on the same cross-track line and must be successive" )
 
     # lat, lon of the 5km bounding pixels
     p1_5km_lat = lat_5km [ itk_p1_5km, isc_5km ]
     p1_5km_lon = lon_5km [ itk_p1_5km, isc_5km ]
     p2_5km_lat = lat_5km [ itk_p2_5km, isc_5km ]
     p2_5km_lon = lon_5km [ itk_p2_5km, isc_5km ]
 
     # check for change date meridian case
     if  abs ( p1_5km_lon  - p2_5km_lon  ) > 180. :
         # all negative longitudes will be incremented of 360 before interpolation
         if p1_5km_lon < 0. :
             p1_5km_lon += 360.
         if p2_5km_lon < 0. :
             p2_5km_lon += 360.
 
     # coordinates of p1, p2 in the 1km grid
     itk_p1_1km = itk_5km_to_1km ( itk_p1_5km )
     itk_p2_1km = itk_5km_to_1km ( itk_p2_5km )
 
     #print "itk_p1_1km=%f itk_p2_1km=%f"%( itk_p1_1km, itk_p2_1km )
 
     # linear interpolation on the position
     alpha_lon = ( p2_5km_lon - p1_5km_lon ) / ( itk_p2_1km - itk_p1_1km )
     alpha_lat = ( p2_5km_lat - p1_5km_lat ) / ( itk_p2_1km - itk_p1_1km )
     lon = p1_5km_lon + alpha_lon * ( itk_1km - itk_p1_1km )
     lat = p1_5km_lat + alpha_lat * ( itk_1km - itk_p1_1km )
 
     if lon > 180. :
         lon = lon - 360.
     elif lon < -180. :
         lon = lon + 360.
 
     return lat, lon
 
 def get_5km_pix_to_plot ( itk_5km, isc_5km, lat_5km, lon_5km ) :
     """
     return the arrays of (lat,lon) and of [i,j] labels of the 5km pixels that bound the one at (itk_5km, isc_5km).
     The first element of these arrays represent the (itk_5km, isc_5km) pixel
     """
     v_itk_5km = [ ]
     v_isc_5km = [ ]
     sc_width = 2 # width of the scan along track
     if ( itk_5km % sc_width ) == 0 : # start of scan
         v_itk_5km = [  itk_5km,     itk_5km    , itk_5km,
                        itk_5km + 1, itk_5km + 1, itk_5km + 1 ]
         if isc_5km == 0 :
             v_isc_5km = [  isc_5km,     isc_5km + 1, isc_5km + 2,
                            isc_5km,     isc_5km + 1, isc_5km + 2 ]
         elif isc_5km == 269 :
             v_isc_5km = [  isc_5km,     isc_5km - 1, isc_5km - 2,
                             isc_5km,     isc_5km - 1, isc_5km - 2 ]
         else :
             v_isc_5km = [   isc_5km,     isc_5km - 1, isc_5km + 1,
                             isc_5km,     isc_5km - 1, isc_5km + 1 ]
     elif ( itk_5km % sc_width ) == ( sc_width - 1 ) : # end of scan
         v_itk_5km = [  itk_5km,     itk_5km    , itk_5km,
                        itk_5km - 1, itk_5km - 1, itk_5km - 1 ]
         if isc_5km == 0 :
             v_isc_5km = [  isc_5km,     isc_5km + 1, isc_5km + 2,
                            isc_5km,     isc_5km + 1, isc_5km + 2 ]
         elif isc_5km == 269 :
             v_isc_5km = [  isc_5km,     isc_5km - 1, isc_5km - 2,
                            isc_5km,     isc_5km - 1, isc_5km - 2 ]
         else :
             v_isc_5km = [  isc_5km,     isc_5km - 1, isc_5km + 1,
                            isc_5km,     isc_5km - 1, isc_5km + 1 ]
 
     npix = len( v_itk_5km )
     # build the labels array
     v_label = [ "[%d,%d]"%(v_itk_5km[i],v_isc_5km[i]) for i in xrange ( npix ) ]
 
     return v_itk_5km, v_isc_5km, v_label
 
 def get_bounds_5km_to_1km( itk_5km, isc_5km ) :
     """
     return the 1km pixel indexes limits in the 5km pixel [ itk_5km, isc_5km ] footprint
     """
 
     # set the (track,scan) indexes of the 5km pixel in the 5km grid
     itk_1km = itk_5km_to_1km ( itk_5km )
     isc_1km = isc_5km_to_1km ( isc_5km )
 
     # set the 1km indexes of pixels to interpolate along track
     itk_1km_min    =  itk_1km - 2
     itk_1km_max    =  itk_1km + 2
 
     # general case : 2 interpolations done along scan : [isc-1, isc] then [isc, isc+1]
     isc_1km_min  = isc_1km - 2
     isc_1km_max  = isc_1km + 2
     # if last 5km pixel along scan, only 4 1km pixels in the 5km footprint in this direction
     if ( isc_5km == sz_sc_5km - 1 ) :
         isc_1km_max  = isc_1km + 6
 
     return itk_1km_min, itk_1km_max, isc_1km_min, isc_1km_max
 
 ##############  MAIN  ###############
 def main():
     """
     Main script
     """
     if len ( sys.argv ) != 4 :
         raise ValueError ( "Invalid number of arguments.\nUsage : python modis_5km_to_1km_geolocation.py <infile> <i_track_1km> <iscan_1km>" )
 
     infile  =       sys.argv [ 1 ]
     itk_1km = int ( sys.argv [ 2 ] )
     isc_1km = int ( sys.argv [ 3 ] )
 
     # -----------------------------------
     # --- Load the 5km (lat,lon) data ---
     # -----------------------------------
     hdf  = SD.SD ( infile )
 
     # read the latitudes of the points
     sds  = hdf.select ( "Latitude" )
     lat_5km = sds.get ( )
     sds.endaccess()
 
     # read the longitudes of the points
     sds  = hdf.select ( "Longitude" )
     lon_5km = sds.get ( )
     sds.endaccess()
 
     hdf.end()
 
     # ---------------------------------------------
     # --- compute the (lat,lon) of a 1km pixel ---
     # ---------------------------------------------
     lat, lon = get_1km_pix_pos ( itk_1km, isc_1km, lat_5km, lon_5km )
     print "%f\t%f"%( lat, lon )
 
     if __DEBUG__ :
         # ----------------------------------
         # --- set the 5km pixels to plot ---
         # ----------------------------------
         itk_5km = itk_1km_to_5km ( itk_1km )
         isc_5km = isc_1km_to_5km ( isc_1km )
 
         itk_5km = max ( 0, round ( itk_1km_to_5km ( itk_1km ) ) )
         isc_5km = max ( 0, min ( round ( isc_1km_to_5km ( isc_1km ) ), sz_sc_5km - 1 ) )
 
         sz_tk_5km = lat_5km[0].shape
         if itk_5km >= sz_tk_5km :
             itk_5km = sz_tk_5km - 1
         if isc_5km >= sz_sc_5km :
             isc_5km = sz_sc_5km - 1
         v_itk_5km, v_isc_5km, v_label = get_5km_pix_to_plot ( itk_5km, isc_5km, lat_5km, lon_5km )
 
         # set the lat,lon of those points
         v_lat_5km = lat_5km [ v_itk_5km, v_isc_5km ]
         v_lon_5km = lon_5km [ v_itk_5km, v_isc_5km ]
 
         # ----------------------------------
         # --- plot 5km pixels
         # ----------------------------------
         # - init the plot where will be drawn the pixels
         fig = plt.figure ()
         ax  = plt.subplot (111)
 
         plt.xlabel ( 'Longitude' )
         plt.ylabel ( 'Latitude' )
         plt.axis   ( 'equal' )
         ax.get_xaxis().set_major_formatter( ticker.FormatStrFormatter ( '%.2f' ) )
         ax.get_yaxis().set_major_formatter( ticker.FormatStrFormatter ( '%.2f' ) )
 
         plt.axis   ( [ v_lon_5km.min(), v_lon_5km.max(), v_lat_5km.min(), v_lat_5km.max() ] )
 
         #print "v_lon_5km "+str(v_lon_5km)
         #print "v_lat_5km "+str(v_lat_5km)
 
         # - plot the 1km pixel
         plt.scatter ( [ lon ], [ lat ], color='m', marker='x', linewidth = 1 )
         ax.text     ( lon, lat, "(%d,%d)"%( itk_1km, isc_1km ), color = 'm' )
 
         # - plot the 5km pixels
         npix = len ( v_lat_5km )
         for i in xrange ( npix ) :
             _lat   = v_lat_5km [ i ]
             _lon   = v_lon_5km [ i ]
             label = v_label   [ i ]
             sz_marker = 1
             if i == 0 :
                 # set a biggest marker for pixel to interpolate
                 sz_marker = 5
             plt.scatter ( [ _lon ], [ _lat ], color='r', marker='o', linewidth = sz_marker, label = "_" )
             ax.text     ( _lon, _lat, label )
 
         # --------------------------------------------------------------------
         # --- Compute and plot the in 1km pixels inside the 5km footprint  ---
         # --------------------------------------------------------------------
         itk_1km_min, itk_1km_max, isc_1km_min, isc_1km_max = get_bounds_5km_to_1km( itk_5km, isc_5km )
         v_lat_1km = []
         v_lon_1km = []
         v_itk_1km = []
         v_isc_1km = []
 
         for itk_1km in xrange ( itk_1km_min, itk_1km_max + 1 ) :
             for isc_1km in xrange ( isc_1km_min, isc_1km_max + 1 ) :
                 v_itk_1km.append ( itk_1km )
                 v_isc_1km.append ( isc_1km )
                 # intert/extrapolate the lat, lon of this 1km pixel using the 5km grid
                 lat, lon = get_1km_pix_pos ( itk_1km, isc_1km, lat_5km, lon_5km )
                 v_lat_1km.append ( lat )
                 v_lon_1km.append ( lon )
         v_itk_1km = np.array ( v_itk_1km )
         v_isc_1km = np.array ( v_isc_1km )
         v_lat_1km = np.array ( v_lat_1km )
         v_lon_1km = np.array ( v_lon_1km )
 
         #print "v_lon_1km "  + str ( v_lon_1km )
         #print "v_lat_1km "  + str ( v_lat_1km )
         #print "v_itk_1km " + str ( v_itk_1km )
         #print "v_isc_1km " + str ( v_isc_1km )
 
         plt.scatter ( v_lon_1km [ v_lon_1km != -9999. ].flat, v_lat_1km [ v_lat_1km != -9999. ].flat, \
                     color='b', marker='+', linewidth = 1, label = "_" )
 
         plt.show ()
 
 if __name__ == "__main__":
   main()