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 )
i_tk_250m = 1.5 + 4 * i_tk_1km i_sc_250m = 4 * i_sc_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
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
MYD05 – MYD06
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.
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.
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
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-lille.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-lille.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()
