связка географических координат с сеткой MODIS

[xxxr1]

Здравствуйте. Возник вопрос: как связать координаты географические с координатами сетки MODIS(http://gis-lab.info/qa/modis-tiles.html)? Ну вот например мне необходимы географические координаты ячейки h21v04... Получается, если я мыслю правильно, мне нужно знать координаты 4-х углов ячейки... Думаю, есть формулы для перевода пикселя от обозначений сетки MODIS к географическим координатам. Поделитесь ссылками или мыслями. Спасибо.

[Александр Мурый]

А чем не устраивает эта разграфка?

[Максим Дубинин]

да, эта сетка в географической системе координат.
данные MODIS которые распространяются тайлами - в синусоидальной (SIN)

сетка специально не делалась в SIN, так как любая ГИС умеет переводить из широты-долготы в SIN.

таким образом, если вам нужна своя реализация, то просто нужно уметь переходить от ГСК к синусоидальной проекции и обратно.

[xxxr1]

Максим, мне действительно нужна своя реализация, но я понятия не имею как перейти от SIN к ГСК. Может быть можете помочь какой-нибудь ссылкой по моему вопросу?

[Максим Дубинин]

синусоидальная проекция достаточно простая, просто почитайте про нее в интернете
начните отсюда: http://en.wikipedia.org/wiki/Sinusoidal_projection

[xxxr1]

Спасибо. Я несколько продвинулся в этом вопросе. А не умеет ли делать это MODIS Reprojection Tool ???

[Максим Дубинин]

Смотря что "это". Просто перепроецировать растр из одной СК в другую - нет, для этого лучше использовать gdal. Перепроецировать с импортом из HDF - да.

http://gis-lab.info/qa/modisimport.html

[xxxr1]

Спасибо. Я пошел несколько иным путем. Для перевода координаты из SIN в ГСК я использовал SIN. Но образовалась следующая проблема. А именно - погрешность при вычислениях.
Итак: изначально я написал простенькую программку, которая выводила бы "проекционную" координату каждого пикселя ячейки MODIS.
Формулы:

k = 926.625433055556 !размер километрового пикселя
n = 1200 !пикселей в ячейке MODIS
x = (k*n*(h-18)) !"будущая" долгота
y = (k*n * (9-v)) !"будущая" широта

h - номер ячейки по горизонтали
v - номер ячейки по вертикали

Для проверки были взяты 4 точки. Крайние координаты ячейки h23v03
№ h v
1 23 4
2 23 3
3 24 4
4 24 3

после выполнения алгоритма, который был описан выше получились результаты (расположены соответственно):
5559752.5983333355 5559752.5983333355
5559752.5983333355 6671703.118000003
6671703.118000003 5559752.5983333355
6671703.118000003 6671703.118000003

Далее я преобразовал получившиеся значения с помощью cs2cs:

cs2cs -f %.7f \
+proj=sinu +lat_0=0.0000000000 +lon_0=0.0000000000 \
+a=6371007.181 +b=6371007.181 +no_defs +to_meter=1.0 \
+to +proj=latlon

Получились следующие значения:
77.7861913 50.0000000 0.0000000
100.0000000 60.0000000 0.0000000
93.3434296 50.0000000 0.0000000
120.0000000 60.0000000 0.0000000
Первый столбец соответственно догота, второй широта.

Для сравнения я использовал данные XML файла продукта MOD09GA ячейки h23v03 (http://reverb.echo.nasa.gov/reverb/):
<GPolygon><Boundary>
<Point><PointLongitude>77.4577</PointLongitude><PointLatitude>49.7712</PointLatitude></Point>
<Point><PointLongitude>100.1725</PointLongitude><PointLatitude>60.1087</PointLatitude></Point>
<Point><PointLongitude>121.7250</PointLongitude><PointLatitude>59.8479</PointLatitude></Point>
<Point><PointLongitude>93.9503</PointLongitude><PointLatitude>49.6337</PointLatitude></Point>
</Boundary></GPolygon>

Данные схожы, но видна определенная погрешность. Возможно ошибка в параметрах +a=6371007.181 +b=6371007.181, но насколько я знаю используются именно такие значения (подставлял другие - результат все равно расходится).

Жду ваших идей. Спасибо.