Экспорт и импорт X Y Z

[zed]

Уже описывал:

В помощь автору, по рельефу NASA.
URL с примером запроса, я уже приводил, теперь про то, что лежит в архиве:
Имя файла: координаты центра левого нижнего пикселя/элемента (уровень зума определяется из URL запроса).
Формат BIL, в нашем случае, представляет собой матрицу 150х150 элементов по 2 байта каждый. Никаких заголовков и проч. нет, блоки (по 2 байта) следуют один за другим. Размер bil файла всегда равен 45 000 байт.
Значение элемента - высота в метрах для данной точки. Принимает значения -32767..32767 (тип SmallInt в Delphi). Если значение элемента -32768, то данных по высоте для данной точки нет.
Тайлы в проекции долгота/широта, датум WGS84.
Вроде так.

[svp]

Не понял

Формат BIL, в нашем случае, представляет собой матрицу 150х150 элементов по 2 байта каждый.

Не понял, то есть это монохромный битмап в raw-формате 150x150x16bit? А почему не 256*256? Как это соотносится с тайлом? Приведите пару примеров, я что-то не понял..

[zed]

Формат bil следует рассматривать как некий бинарный файл, каждые 2 байта которого, определяют высоту над уровнем моря, для точки с определёнными координатами. Понятие "тайл" введено для определения координат внутри самого bil, а разрешение "тайла" (256*256) - для определения границ региона, описываемого данным bil файлом.
Скажем, zoom=3, один тайл (256*256) описывает площадь 90 град. по долготе и широте, и при этом значения высоты доступны с шагом 0,6 град (90/150). Т.е. если перейти к bil:
1-2-й байты - высота для точки 0,-180;
3-4-й байты - точка 0, -179.4 и т.д.
300-301-й байты точка 0,-90;
302-303-й - точка 0.6,-180...
(начальные координаты для bil (из имени файла) даны для его левого нижнего элемента, а координаты следующего элемента определяются как смещение на n градусов)
При увеличении зума, шаг уменьшается, и минимально возможный шаг = 0,0008326 град, с более точным шагом, в свободном доступе данные не распространяются (хотя, у GE, данные возможно точнее).

[zed]

У NASA, в архивах и только по Евразии ~ 5 ГБ данных

Кстати, может кто знает - а где бы нарыть столь же подробные (ну или хотя бы сравнимые с оным SRTM) данные о глубинах (bathymetry)?

На том же FTP есть данные по рекам, океанам:

Дополнительно с данными в виде отдельных слоев распространяются также данные по площадным объектам гидросети (SRTM Water Body Data - SWBD), представленные в формате данных 3-D Shapefile (шейп-формат с информацией о третьей - высотной координате).

Этот набор данных является побочным продуктом, полученным в процесс редактирования данных, осуществленным National Geospatial-Intelligence Agency США, для получения окончательного набора данных SRTM (DTED® 2). В слое объектов гидросети представлены объекты прошедшие при редактировании критерий минимальной площади, океаны, озера, водотоки. Высота н.у.м. озер постоянна. Высота н.у.м. океанов равна 0. Lake elevations were set to a constant value. Высота н.у.м. рек монотонно понижается, чтобы сохранить правильный сток. Редактирование данных производилось на базе матриц разрешения 30-м. Система координат данных также WGS 84. Горизонтальная точность - 20м (90%), вертикальная - 16м (90%). Выделенные объекты имеют следующую информацию в атрибутивной таблице:
Океаны - BA040
Озера - BH080
Реки - BH140

[svp]

Тогда так:
Пусть [x, y] -- координаты пикселя внутри тайла карты: x=0..255, y=0..255.
bil[0..149,0..149] -- матрица высот для текущего тайла (под курсором).
h[x, y] -- высота, соответствующая пикселю [x, y] тайла.
[bx, by] = [x/256*150, y/256*150] -- вещественные координаты точки внутри матрицы высот. Для этой точки нужно интерполировать высоту.
dx = bx - floor(bx); // координаты точки внутри текселя матрицы высот
dy = by - floor(bx); // 0 ≤ dx < 1; 0 ≤ dy < 1

Код: Выделить всё

h[x, y] =
  bil[ ceil(bx),  ceil(by)] * dx * dy +
  bil[floor(bx),  ceil(by)] * (1-dx) * dy +
  bil[ ceil(bx), floor(by)] * dx * (1-dy) +
  bil[floor(bx), floor(by)] * (1-dx) * (1-dy)

Эти выкладки надо бы проверить=)

P.S.
Xf =floor(X) -- округление до ближайшего целого (меньшего): Xf ≤ X.
Xc =ceil(X) -- округление до ближайшего целого (большего): Xc ≥ X.
Стандартные функции. Есть в стандартном модуле math.pas

[zed]

Пусть [x, y] -- координаты пикселя внутри тайла карты: x=0..255, y=0..255.
bil[0..149,0..149] -- матрица высот для текущего тайла (под курсором).

Проекция тайла карты и тайла bil - в общем случае не совпадают (у bil - широта/долгота, у тайлов - меркатор), поэтому их нельзя ставить в соответствие, либо должна быть "перепроецированная" матрица bil[0..n,0..m], соответствующая тайлу на конкретной широте

[svp]

Надо было сразу попрозрачнее это сказать=)
Тогда даже проще:

Пусть [x, y] -- широта и долгота; h[x, y] -- высота в этой точке.
bil[0..149,0..149] -- матрица высот.
[b0x, b0y] -- широта и долгота левого нижнего угла матрицы высот.
[bdx, bdy] -- размеры в градусах матрицы высот.
[bx, by] = [x/bdx*150, y/bdy*150] -- вещественные координаты точки внутри матрицы высот. Для этой точки нужно интерполировать высоту.

Остальное останется неизменным:

dx = bx - floor(bx); // координаты точки внутри текселя матрицы высот
dy = by - floor(bx); // 0 ≤ dx < 1; 0 ≤ dy < 1

h[x, y] =
bil[ ceil(bx), ceil(by)] * dx * dy +
bil[floor(bx), ceil(by)] * (1-dx) * dy +
bil[ ceil(bx), floor(by)] * dx * (1-dy) +
bil[floor(bx), floor(by)] * (1-dx) * (1-dy)

Поподробнее с несколькими примерами опишите как именуются файлы с матрицами высот.

[Parasite]

На том же FTP есть данные по рекам, океанам

Ммда? Ок, благодарю. Ознакомлюсь.
Ранее как-то руки не доходили укачать оттуда всё...

[zed]

Надо было сразу попрозрачнее это сказать=)

Ну уж не знаю, куда прозрачней

Уже описывал:

В помощь автору, по рельефу NASA.
...
Тайлы в проекции долгота/широта, датум WGS84.
Вроде так.

Поподробнее с несколькими примерами опишите как именуются файлы с матрицами высот.

Имя файла Y_X.bil ( 0114_0165.bil, 0114_0166.bil). Опять же, подводный камень: X,Y - для левого нижнего пикселя. Если обычный тайл описывает квадрат X,Y (верх-лево); X+1,Y+1 (низ-право); то bil: X,Y-1 (верх-лево); X+1,Y (низ-право) вроде не перепутал куда там +/-1?

Изменяющаяся часть запроса: L=5&X=328&Y=230
L - уровень зума, L=0 соответствует zoom=8 (для меньшего зума данных нет, хотя для L=-1 файл скачивается, но содержит только нули) максимально-доступный L=11 (zoom=19).
X,Y - номер строки и столбца - стандартные имена тайлов для NASA (для всех типов данных). Нумерация с верхнего левого угла. Т.е. если в расчётах, гугловские имена qrst мы переводили в X,Y, то здесь уже никаких переводов не требуется, и X,Y можно получить прямо из имени. Только не знаю, нужно учитывать проекцию или нет, но скорее всего нет и X,Y - те же, что и для гугла.

[svp]

Кратко приведу диалог из лички, ибо полезно будет, если кто-то ещё выскажется.

Имя файла Y_X.bil ( 0114_0165.bil, 0114_0166.bil).

Что здесь X и Y? Как получаются 0114_0165?
Это широта и долгота? или что?

Изменяющаяся часть запроса: L=5&X=328&Y=230
L - уровень зума, L=0 соответствует zoom=8 (для меньшего зума данных нет, хотя для L=-1 файл скачивается, но содержит только нули) максимально-доступный L=11 (zoom=19).
X,Y - номер строки и столбца - стандартные имена тайлов для NASA (для всех типов данных). Нумерация с верхнего левого угла. Т.е. если в расчётах, гугловские имена qrst мы переводили в X,Y, то здесь уже никаких переводов не требуется, и X,Y можно получить прямо из имени. Только не знаю, нужно учитывать проекцию или нет, но скорее всего нет и X,Y - те же, что и для гугла.

Нам надо по координатам и зуму получать ссылку и имя файла, то есть надо получать X и Y. Если я парвильно понял, то ширина высотного тайла в градусах наразных широтах будет разной. Как её вычислять?
Не понятно, как из широты, долготы и зума получить X и Y.

Почему?
X,Y считаем стандартной функцией из SAS:

Код: Выделить всё

function GLonLat2Pos(Ll:TExtendedPoint;Azoom:byte;MT:PMapType):Tpoint;
var z,c:real;
begin
result.x:=round(zoom[Azoom]/2+ll.x*(zoom[Azoom]/360));
case MT.projection of
  1: begin
      z:=sin(Ll.y*deg);
      c:=(zoom[Azoom]/(2*Pi));
      result.y:=round(zoom[Azoom]/2-0.5*ln((1+z)/(1-z))*c);
     end;
  2: begin
      z:=sin(Ll.y*deg);
      c:=(zoom[Azoom]/(2*Pi));
      result.y:=round(zoom[Azoom]/2-c*(ArcTanh(sin(ll.y*deg))-MT.exct*ArcTanh(MT.exct*sin(ll.y*deg))) )
     end;
3: result.y:=round(zoom[Azoom]/2-ll.y*(zoom[Azoom]/360));
end;
end;

проекция 3 - (широта/долгота).
Впрочем, с проекцией bil следует уточнить эксперементально (bil что раздаются с ftp - 100% в проекции широта/долгота, а на счёт "наших" bil я начал сомневаться: на запрос l=0 x=0 y=0 - выдали файл с высотами, хотя в проекции широта/долгота такого файла быть не должно, потому что это 180,-180 )

Дело в том, что используемый в этом в этой функции массив Zoom[] зависит от размера тайла:

Код: Выделить всё

  zoom:array [1..24] of longint = (256,512,1024,2048,4096,8192,16384,32768,65536,
                                   131072,262144,524288,1048576,2097152,4194304,
                                   8388608,16777216,33554432,67108864,134217728,
                                   268435456,536870912,1073741824,2147483647);


Достаточно ли просто сделать поправку на размер тайла и всё? То есть можно ли использовать ту же формулу только с учётом, что:

Код: Выделить всё

  zoom:array [1..24] of longint = (150,300,600,1200,2400,4800,9600,19200,38400,
                                   76800,153600,307200,614400,1228800,2456700,
                                   4915200,9830400,19660800,39321600,78643200,
                                   157286400,314572800,629145600,1258291200);