[GIS] Getting tile number of sinusoidal MODIS product from lat/long

modismodis-reprojection-toolnasapython

I'm using this tool to get the tile number (horizontal/vertical) from lat/long. Is there any code,preferably python, to get the tile number ? I tried using this tool. But it gives wrong tile number. Here's my code:

import numpy as np
data = np.genfromtxt('sn_bound_10deg.txt', skip_header = 7, skip_footer = 3)
lat,lon = -25.439538,149.053418
in_tile = False
i = 0
while(not in_tile):
  #print data[i,0],data[i,1],data[i, 4],data[i, 5],data[i, 2],data[i, 3]
  in_tile = lat >= data[i, 4] and lat <= data[i, 5] and lon >= data[i, 2] and lon <= data[i, 3]
  i += 1

vert = data[i-1, 0]
horiz = data[i-1, 1]
print('Vertical Tile:', vert, 'Horizontal Tile:', horiz)

The output tile from the above code is: h30v11

Expected output and output from the online calculator is: h31v11

Here's the sn_bound_10deg.txt file.

Best Answer

MODIS assumes the Earth is a sphere with a R = 6.3781e6 m radius. The world is subdivided into 18 vertical tiles and 36 horizontal tiles, each tile has 2400 cells of 500m x 500m.

MODIS projection

We can code this as:

import math
CELLS = 2400
VERTICAL_TILES = 18
HORIZONTAL_TILES = 36
EARTH_RADIUS = 6371007.181
EARTH_WIDTH = 2 * math.pi * EARTH_RADIUS

TILE_WIDTH = EARTH_WIDTH / HORIZONTAL_TILES
TILE_HEIGHT = TILE_WIDTH
CELL_SIZE = TILE_WIDTH / CELLS

Now, what we need is some way to convert latitudes and longitudes into the MODIS grid. PyProj can do that for us.

from pyproj import Proj
MODIS_GRID = Proj(f'+proj=sinu +R={EARTH_RADIUS} +nadgrids=@null +wktext')

With that, converting from geographic coordinates can be done with

def lat_lon_to_modis(lat, lon):
    x, y = MODIS_GRID(lon, lat)
    h = (EARTH_WIDTH * .5 + x) / TILE_WIDTH
    v = -(EARTH_WIDTH * .25 + y - (VERTICAL_TILES - 0) * TILE_HEIGHT) / TILE_HEIGHT
    return int(h), int(v)

I sampled a few points with the MODLAND tile calculator and values seem to match:

if __name__ == '__main__':
    for lat, lon, h, v in [(-5, -42, 13, 9),
                           (-23.533773, -46.62529, 13, 11),
                           (0, 179, 35, 9),
                           (89, 179, 18, 0),
                           (-25.439538, 149.053418, 31, 11),]:
        print(lat, lon, h, v, lat_lon_to_modis(lat, lon))
        assert (h, v) == lat_lon_to_modis(lat, lon)

Notice we get the right output using the test coordinates: -25.439538, 149.053418

The above code and discussion is based on the information found in https://code.env.duke.edu/projects/mget/wiki/SinusoidalMODIS.

Related Question