# Getting tile number of sinusoidal MODIS product from lat/long?

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.

• I'm looking at sn_bound_10deg and it definitely feels like some of the tiles overlap. Not sure exactly what is going on – badgley Feb 4 '18 at 0:10
• Often you can't calculate a single tile number from lat/lon bounding coordinates (try to reproject an image and you'll see why). Though the online tile calculator is open source landweb.modaps.eosdis.nasa.gov/cgi-bin/developer/tilemap.cgi – dr_times Feb 6 '18 at 11:24

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. We can code this as:

``````import math
CELLS = 2400
VERTICAL_TILES = 18
HORIZONTAL_TILES = 36
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
``````

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.