4

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
3

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.