Трансформация карт

[Vladimir]

Доброго дня! Есть задача похожая на http://www.prog.org.ru/topic_24274_0.html, но мне нужно трансформировать существующие тайлы гугловских карт на 2D плоскости! Т.е. есть такая проблема: если собрать тайлы в одну картинку QImage для z8 масштаба скажем и по координатам гугла выставить какую-либо видимую отметку (скажем пересечение дорог), то она не совпадет с пересечением дорог на карте, отметка будет смещена. Таким образом мне нужно трансформировать карту, чтобы координаты отметок точно совпадали с картой. В выше приведенной статье был предложен алгоритм, если его модифицировать подойдет ли он мне:
- каждый тайл представить в виде сетки точек I x J.
   - для каждой точки
     1. выполнить преобразование из проекции Меркатора на геоид WGS84
     2. выполнить преобразование из геоида WGS84 в 2D декартовые координаты
Если такой алгоритм правильный для моей задачи, то может кто подскажет формулы пересчета и как это реализовывать, т.е. необходимо попиксельно пересчитывать картинку из одних координат в другие?

[Kurles]

Вроде у гугла нумерация координат и уровней зума совпадает с OSM. Посмотри в сторону расчёта координат для Slippy map. По ссылке даны формулы для координат целых тайлов, но их не сложно привести к конкретным пикселям в нужном тайле.

[Vladimir]

Вроде у гугла нумерация координат и уровней зума совпадает с OSM. Посмотри в сторону расчёта координат для Slippy map. По ссылке даны формулы для координат целых тайлов, но их не сложно привести к конкретным пикселям в нужном тайле.

Так у меня по формулам из целых тайлов картинка и собирается, просто рассчитывается место каждого тайла на экране и он достается из нужной папки "/z%1/%2/x%3/%4/y%5.png" без всяких модификаций! Просто получается, что чем дальше от экватора в северном (южном) направлении находится точка стояния, тем больше погрешность получается (картинка вроде как должна сжиматься сверху и расширятся снизу, для северного направления), а здесь она остается просто прямоугольником, собранным из отдельных тайлов. Карты гугла вроде как в проекции Меркатора, а данная проекция искажает площади, поэтому видимо для правильного отображения необходимо делать пересчеты в другую систему координат.. вопрос в какую и как?

[Kurles]

Вроде у гугла нумерация координат и уровней зума совпадает с OSM. Посмотри в сторону расчёта координат для Slippy map. По ссылке даны формулы для координат целых тайлов, но их не сложно привести к конкретным пикселям в нужном тайле.

Так у меня по формулам из целых тайлов картинка и собирается, просто рассчитывается место каждого тайла на экране и он достается из нужной папки "/z%1/%2/x%3/%4/y%5.png" без всяких модификаций! Просто получается, что чем дальше от экватора в северном (южном) направлении находится точка стояния, тем больше погрешность получается (картинка вроде как должна сжиматься сверху и расширятся снизу, для северного направления), а здесь она остается просто прямоугольником, собранным из отдельных тайлов. Карты гугла вроде как в проекции Меркатора, а данная проекция искажает площади, поэтому видимо для правильного отображения необходимо делать пересчеты в другую систему координат.. вопрос в какую и как?

slippymaps тоже в проекции Меркатора, судя по документации. Приложи тестовый пример с проблемой, быстрей помочь получится.

[Vladimir]

минимальный проект, будет сложновато сделать..проект большой, чтобы оттуда выдрать все нюансы! вот функция загрузки карты. по сути в ней вся работа с тайлами карты!

Код:

void CThreadLoadMap::load()
{  
    QImage tmpImage;
    QImage tmpIm;
    QString fullPathFile;

    double lat_rad =  latitude * M_PI/180;
    double lon_rad =  longitude * M_PI/180;

    double scrRes; // разрешение экрана
    double scaleK; // масштабный коэффициент
    int startX;
    int startY;
    int tileNewSize;
    int numTilesX;
    int numTilesY;

    int Z;
    int BM0;
    int Xp;
    int Yp;

    int numX;
    int numY;

    int dx;
    int dy;

    int R0;
    int WW; // длина окружности на данной широте

    double XYRes;

    int w = szWinDraw.width(),
        h = szWinDraw.height();

    int sizeN = w*h*4;
    char *buf2N = new char[sizeN*8];

    bool fileExist = false;

    int index = 0;

    //10-20=4;30-50=3;60-100=2;110-200=1;210-400=0
    int curScale = currentScale;
    if(curScale > 200) index = 0;else
    if(curScale > 100) index = 1;else
    if(curScale >  50) index = 2;else
    if(curScale >  20) index = 3;else
                       index = 4;

    Z = 8 + index;
    BM0 = pow(2,(Z-1+8))/2; // центр мира в пикселях для данного уровня
    R0  = 6378137 *         // радиус большой полуоси эллипсоида
          cos(lat_rad);
    WW  = 2*M_PI*R0;        // длина окружности планеты на lat широте

    XYRes = WW/(pow(2,(Z-1))*256);

    scrRes = currentScale*1000.*2/(szWinDraw.height() - 2*PaddingWinDraw);
    scaleK = (XYRes / scrRes);

    // координаты на растровом массиве данного уровня в пикселах, от верхнего левого угла массива
    Xp = BM0*(1+lon_rad/M_PI);
    Yp = BM0*(1 - 0.5*log((1+sin(lat_rad))/(1-sin(lat_rad)))/M_PI);
    // номера блоков в тайловом массиве
    numX = (Xp/256);
    numY = (Yp/256);
    // смещение внутри тайла в пикселях
    dx = (Xp - numX*256);
    dy = (Yp - numY*256);

    tileNewSize = floor(256*scaleK);

    numTilesX   =  2*round(0.5+0.5*w/tileNewSize);
    numTilesY   =  2*round(0.5+0.5*h/tileNewSize);

    startX = numX - (numTilesX/2);
    startY = numY - (numTilesY/2);

    dx = dx*scaleK + abs(w - numTilesX*tileNewSize)/2;
    dy = dy*scaleK + abs(h - numTilesY*tileNewSize)/2;

    int* dest = (int*)buf2N;

    static int numb = 0;
    for(int x = startX; x < (startX+numTilesX+1); x++)
    {
      for(int y = startY; y < (startY+numTilesY+1); y++)
      {
          fullPathFile = "";

          fullPathFile = patchToMap + QString("/z%1/%2/x%3/%4/y%5.jpg")
              .arg(Z )
              .arg(x/1024)
              .arg(x)
              .arg(y/1024)
              .arg(y);

          fileExist = tmpImage.load(fullPathFile);

          if(!fileExist)
          {
              QString nameFileMap = QString("/z%1/%2/x%3/%4/y%5.png")
                      .arg(Z)
                      .arg(x/1024)
                      .arg(x)
                      .arg(y/1024)
                      .arg(y);

              fullPathFile = patchToMap + nameFileMap;

              fileExist = tmpImage.load(fullPathFile);
          }

          int width = 2*w;

          if(fileExist)
          {
              tmpIm = tmpImage.scaled(tileNewSize,tileNewSize, Qt::KeepAspectRatio,
                                                               Qt::SmoothTransformation);
              numb++;

              for(int yy = 0; yy < tmpIm.height(); yy++)
                for(int xx = 0; xx <  tmpIm.width(); xx++)
                {
                    dest[
                        ((yy)*width+xx)
                        +((x-startX)*tileNewSize+(y-startY)*tileNewSize*width)
                        ] = (int)tmpIm.pixel(xx,yy);
                }
          }
       } // end for(;;) - y
    } // end for(;;) - x

    int* src = (int*)buf2N;

    for(int y = 0; y < h; y++)
     for(int x = 0; x < w; x++)
      if(layerMap->valid(x,y))
        layerMap->setPixel(x,y,src[((x+dx) + (y+dy)*w*2)]);

    delete [] buf2N;
}

На картинке пример смещения точек. Черная - это то место, где хотелось бы видеть отметку, а красная - это точка вбитая по координатам гугла (то место где должно быть пересечение дорог).

[Kurles]

Посмотрел мельком код - слишком у тебя лишнего вроде нагорожено, для получения тайловых координат из широты и долготы достаточно следующей процедуры:

Код

C++ (Qt)
QPointF tileForCoordinate(qreal lat, qreal lng, int zoom)
{
    qreal zn = static_cast<qreal>(1 << zoom);
    qreal tx = (lng + 180.0) / 360.0;
    qreal ty = (1.0 - log(tan(lat * pi / 180.0) +
                          1.0 / cos(lat * pi / 180.0)) / pi) / 2.0;
    return QPointF(tx * zn, ty * zn);
}
 

Что бы узнать координаты пикселя внутри тайла - достаточно дробную часть умножить на 256.

Ну и до кучи процедуры для обратного преобразования:

Код

C++ (Qt)
qreal longitudeFromTile(qreal tx, int zoom)
{
    qreal zn = static_cast<qreal>(1 << zoom);
    qreal lat = tx / zn * 360.0 - 180.0;
    return lat;
}
 
qreal latitudeFromTile(qreal ty, int zoom)
{
    qreal zn = static_cast<qreal>(1 << zoom);
    qreal n = pi - 2 * pi * ty / zn;
    qreal lng = 180.0 / pi * atan(0.5 * (exp(n) - exp(-n)));
    return lng;
}
 

Процедуры проверенные, по крайней мере для OSM карт работают на ура.

[Vladimir]

Функция написана давно и не мной, черт ногу сломит, что там делается! но в целом она работает, за исключением того что вылез касяк со смещением.. и не понятно то ли это тайлы не совсем правильно размещаются на экране, то ли нужно в целом пересчитывать картинку по другому, ибо сами карты того требуют (Меркатор -> WGS84).. Kurles а что по вашему там лишнее?

[Kurles]

Функция написана давно и не мной, черт ногу сломит, что там делается! но в целом она работает, за исключением того что вылез касяк со смещением.. и не понятно то ли это тайлы не совсем правильно размещаются на экране, то ли нужно в целом пересчитывать картинку по другому, ибо сами карты того требуют (Меркатор -> WGS84).. Kurles а что по вашему там лишнее?

Посмотрел немного внимательнее - всё равно не очень понятно, что там делается. Эта процедура пытается всю карту планеты загрузить в картинку? Может проще самому эту часть переписать? И вполне возможно, что действительно тайлы не правильно размещаются на экране, так как смещение не только до широте, но и по долготе.

[Vladimir]

Посмотрел немного внимательнее - всё равно не очень понятно, что там делается. Эта процедура пытается всю карту планеты загрузить в картинку? Может проще самому эту часть переписать? И вполне возможно, что действительно тайлы не правильно размещаются на экране, так как смещение не только до широте, но и по долготе.

Нет не всю.. с начала создается временный буфер для склейки выбранных тайлов с диска

Код:

int sizeN = w*h*4;
char *buf2N = new char[sizeN*8];

затем выбирается масштаб карт (соотносится масштаб карт гугла Z с масштабом в программе)

Код:

//10-20=4;30-50=3;60-100=2;110-200=1;210-400=0
int curScale = currentScale;
if(curScale > 200) index = 0;else
if(curScale > 100) index = 1;else
if(curScale >  50)  index = 2;else
if(curScale >  20)  index = 3;else
                          index = 4;

здесь видимо считается масштабный коэффициент, соотношение точки на мировой карте для данного масштаба и широты к точке на экране..

Код:

    Zoom = 8 + index;
    BM0  = pow(2,(Zoom-1+8))/2;     // центр мира в пикселях для данного уровня
    R0   = 6378137 * cos(lat_rad);  // радиус большой полуоси эллипсоида
    WW   = 2*M_PI*R0;               // длина окружности планеты на lat широте

    XYRes = WW/(pow(2,(Zoom-1))*256);

    scrRes = currentScale*1000.*2/(szWinDraw.height() - 2*PaddingWinDraw);
    scaleK = (XYRes / scrRes);

здесь идет расчет координат тайлов и их количества в временном буфере

Код:

// координаты на растровом массиве данного уровня в пикселах, от верхнего левого угла массива
    Xp = BM0*(1+lon_rad/M_PI);
    Yp = BM0*(1 - 0.5*log((1+sin(lat_rad))/(1-sin(lat_rad)))/M_PI);

    // номера блоков в тайловом массиве
    numX = (Xp/256);
    numY = (Yp/256);
    // смещение внутри тайла в пикселях
    dx = (Xp - numX*256);
    dy = (Yp - numY*256);

    tileNewSize = floor(256*scaleK);

    numTilesX = 2*round(0.5+0.5*w/tileNewSize);
    numTilesY = 2*round(0.5+0.5*h/tileNewSize);

    startX = numX - (numTilesX/2);
    startY = numY - (numTilesY/2);

    dx = dx*scaleK + abs(w - numTilesX*tileNewSize)/2;
    dy = dy*scaleK + abs(h - numTilesY*tileNewSize)/2;

Здесь по рассчитанным координатам достаем тайлы с диска и помещаем в буфер

Код:

int* dest = (int*)buf2N;

    static int numb = 0;
    for(int x = startX; x < (startX+numTilesX+1); x++)
    {
      for(int y = startY; y < (startY+numTilesY+1); y++)
      {
          fullPathFile = "";

          fullPathFile = patchToMap + QString("/z%1/%2/x%3/%4/y%5.jpg")
              .arg(Zoom )
              .arg(x/1024)
              .arg(x)
              .arg(y/1024)
              .arg(y);

          fileExist = tmpImage.load(fullPathFile);

          if(!fileExist)
          {
              QString nameFileMap = QString("/z%1/%2/x%3/%4/y%5.png")
                      .arg(Zoom)
                      .arg(x/1024)
                      .arg(x)
                      .arg(y/1024)
                      .arg(y);

              fullPathFile = patchToMap + nameFileMap;

              fileExist = tmpImage.load(fullPathFile);
          }

          int width = 2*w;

          if(fileExist)
          {
              tmpIm = tmpImage.scaled(tileNewSize,tileNewSize, Qt::KeepAspectRatio,
                                                               Qt::SmoothTransformation);
              numb++;

              for(int yy = 0; yy < tmpIm.height(); yy++)
                for(int xx = 0; xx <  tmpIm.width(); xx++)
                {
                    dest[
                        ((yy)*width+xx)
                        +((x-startX)*tileNewSize+(y-startY)*tileNewSize*width)
                        ] = (int)tmpIm.pixel(xx,yy);
                }
          }
       } // end for(;;) - y
    } // end for(;;) - x

последнее преобразование это, копирование из буфера в картинку, которая рисуется на экране

Код:

 for(int y = 0; y < h; y++)
     for(int x = 0; x < w; x++)
      if(layerMap->valid(x,y))
        layerMap->setPixel(x,y,src[((x+dx) + (y+dy)*w*2)]);

Пока писал сам немного понял но все равно не понятно, где и в чем может быть косячина..

[Kurles]

Посмотрел немного внимательнее - всё равно не очень понятно, что там делается. Эта процедура пытается всю карту планеты загрузить в картинку? Может проще самому эту часть переписать? И вполне возможно, что действительно тайлы не правильно размещаются на экране, так как смещение не только до широте, но и по долготе.

...
Пока писал сам немного понял но все равно не понятно, где и в чем может быть косячина..

Что то где-то у тебя не правильно считается. Не поленился и тестовое приложение по-быстрому сварганил, всё ок.

[Vladimir]

Что то где-то у тебя не правильно считается. Не поленился и тестовое приложение по-быстрому сварганил, всё ок.

Спасибо БОЛЬШОЕ! Буду разбираться..