Перепроектировать плитки Web Mercator в произвольную проекцию с D3? - PullRequest
9 голосов
/ 09 июня 2019

Прошло несколько лет с тех пор, как Джейсон Дэвис поразил нас Перепроизведенными растровыми плитками - эта карта перестала работать, потому что Mapbox блокирует его сайт, но Mollweide Акварель и Прерванный Goode Raster остается отличным демо.

Теперь в Observable HQ я вижу документы для самых последних d3-гео-проекций и d3-тайлов , но нет современных примеров того, как делать то, что делал Джейсон: репроектировать стандарт Наборы плитки Mercator.

Как я могу заставить d3-плитку перейти в новую проекцию?

1 Ответ

2 голосов
/ 18 июня 2019

Этот ответ основан на:

Эти три ресурса также основаны друг на друге. Понимание этих трех примеров поможет понять, что происходит в моем примере ниже.

В ответе в качестве основы также используется моя медленная, продолжающаяся попытка построить библиотеку плиток .

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

Плитка веб-меркатора

Карта Меркатора, растягивающаяся на 360 градусов долготы и ~ 170 градусов широты (+/- 85 градусов), заполнит квадрат (превышение 85 градусов по широте приводит к искажению, выходящему из-под контроля, и включению полюсов не рекомендуется, так как полюсы находятся на +/- бесконечности на проецируемой плоскости).

Этот квадрат большей части мира для веб-картографического сервиса (с плитками Меркатора) имеет уровень масштабирования 0. Карта имеет площадь 2 ^ 0 и высоту 2 ^ 0.

Если мы разделим этот квадрат на сетку из двух квадратов на два квадрата, у нас будет уровень масштабирования 1. Карта имеет размер 2 ^ 1 на 2 ^ 1 квадратов.

Следовательно, уровень масштабирования определяет количество квадратов и высоту карты: 2 ^ zoomLevel. Если каждый квадрат имеет одинаковый размер в пикселях, то каждое увеличение уровня масштабирования на единицу увеличивает ширину пикселя в мире в 2 раза.

К счастью для нас, к северу от ~ 85 градусов нет земли, и мы не часто показываем Антарктиду, поэтому этот квадрат подходит для большинства приложений веб-картографирования. Тем не менее, это означает, что если мы перепроектируем плитки веб-меркатора на что-либо, что показано выше этих широт, у нас будет лысая точка:

world map with bald spot

Wᴇʙ-Mᴇʀᴄᴀᴛᴏʀ ʀᴇᴘʀᴏᴊᴇᴄᴛᴇᴅs ʀᴇᴘʀᴏᴊᴇᴄᴛᴇᴅ ғᴏʀ ᴀ Mᴇʀᴄᴀᴛᴏʀ ᴘʀᴏᴊᴇᴄᴛɪᴏɴ ᴛʜᴀᴛ ʜᴀs ʙᴇᴇɴ ʀᴏᴛᴀᴛᴇᴅ ᴛᴏ sʜᴏᴡ ᴛʜᴇ Nᴏʀᴛʜ Pᴏʟᴇ.

Наконец, листы веб-меркатора отображаются в проецируемом пространстве, которое является предсказуемым и регулярным по отношению к листам. Если мы перепроектируем плитки, мы можем увеличивать или уменьшать проецируемое пространство для каждой плитки, мы должны помнить об этом. На изображении выше, плитки вокруг Северного полюса перепроектированы так же намного меньше, чем те, что дальше на юг. Плитки не обязательно будут иметь одинаковый размер после проецирования.

Пластина репроекции и повторной выборки Carree

Самая большая проблема для перепроектирования плиток веб-сервиса - это время, а не только время, потраченное на понимание прогнозов и чтение таких ответов.

Функции проецирования - это сложные трудоемкие операции, которые должны выполняться для каждого визуализируемого пикселя. Все примеры d3, которые я видел, используют процесс (или близкий вариант), как показано здесь для фактического перепроецирования и повторной выборки. Этот пример будет работать, только если исходное изображение проецируется с Plate Carree . Процесс выглядит следующим образом:

  1. Создать пустое новое изображение.
  2. Для каждого пикселя в новом изображении взять его местоположение в пикселях и инвертировать его (используя желаемую проекцию), чтобы получить широту и долготу.
  3. Определите, какой пиксель в исходном изображении перекрывает эту долготу и широту.
  4. Возьмите информацию из этого пикселя в исходном изображении и назначьте ее соответствующему пикселю в новом изображении (пиксель на шаге 2)

Когда исходное изображение использует проекцию Plate Carree, нам не нужна d3-geoProjection, соотношение междупрогнозируемые и непроецированные координаты. Например: если изображение имеет высоту 180 пикселей, каждый пиксель представляет 1 градус широты. Это означает, что шаг 3 не занимает много времени по сравнению с шагом 2 и projection.invert (). Вот функция Майка для шага 3:

var q = ((90 - φ) / 180 * dy | 0) * dx + ((180 + λ) / 360 * dx | 0) << 2;

Время, необходимое для шага 2, связано с проекцией, используемой для перепроецированного изображения. Все примеры, которые я видел, используют d3.geoProjection.invert() для второго шага в приведенном выше списке - определение местоположения пикселя в новом изображении и определение его широты и долготы. Не все прогнозы рождаются равными. Цилиндрические проекции обычно превосходят конические, а конические проекции, как правило, не выполняют азимутальные проекции. Я также видел некоторые странные различия в скорости между d3v4 и d3v5 с точки зрения projection.invert () time:

relative time for projection.invert

Tɪᴍᴇ ᴛᴏ ᴀ ᴀ ᴘᴏɪɴᴛ ᴏɴ ᴀ ᴍᴀᴘ ᴡɪᴛʜ D3 (ᴄᴏɴᴠᴇʀᴛ ғʀᴏᴍ ᴘɪxᴇʟs ᴛᴏ ʟᴀᴛ / ʟᴏɴɢ). Iᴛ ɪs ᴜɴᴄʟᴇᴀʀ ᴡʜʏ D3ᴠ4 ᴡᴏᴜʟᴅ ʙᴇ ғᴀsᴛᴇʀ.

А для полноты, вот более широкий диапазон проекций, найденных в d3-геопроекциях:

enter image description here

Cᴏᴍᴘᴀʀɪɴɢ ᴛɪᴍᴇ ғᴏʀ ᴘʀᴏᴊᴇᴄᴛɪᴏɴ.ɪɴᴠᴇʀᴛ ғᴏʀ ᴀᴅᴅɪᴛɪᴏɴᴀʟ ᴘʀᴏᴊᴇᴄᴛɪᴏɴs

Примечания

  • Этот подход может столкнуться с проблемами, описанными здесь , эти проблемы могут быть решены с помощью ответов там - но различные решения потребуют дополнительного времени для обработки.

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

  • Если базовое изображение содержит текст, текст можно поворачивать или иным образом манипулировать, чтобы сделать его менее читаемым.

  • Пример Майка для отдельного изображения, для плиток, процесс в некоторой степени изменен, мы сейчас создаем несколько изображений, и это требует знания границ каждой исходной плитки и границ каждой перепроектированной плитки в градусах, а также в единицах плитки для первого и пикселей для второго - второстепенные детали.

Перепроецирование и ресэмплинг Web Mercator

Когда я начал смотреть на этот вопрос, я посмотрел на решение / пример Алана Маккончи в качестве ссылки. На то, чтобы заметить, потребовалось некоторое время, но шаг 3 в этом примере (и я полагаю, что работа Джейсона Дэвиса тоже) не учитывает плитку веб-Меркатора при повторной выборке - только при определении границ плитки. Но соотношение между пикселями на оси y больше не является линейным, как это было в Plate Carree.

Это означает, что плитки размещены в нужных местах, но выборка рассматривает ось Y как линейную в каждой плитке. Это искажение было наиболее заметно при отображении всего мира с низким уровнем масштабирования плитки (на полпути вниз / вверх по плитке), и может быть тем, о чем говорит Алан, когда упоминает о нечетном сжатии.

Решение состоит в том, чтобы правильно спроецировать широту для каждой пары широта / долгота в шаге 3 выше. Это добавляет время, всегда больше времени - функция включает в себя Math.atan и Math.exp, но разница не должна быть слишком большой. В работе Алана и Джейсона это делается с помощью простой формулы (но используется только для границ тайлов, а не для каждого пикселя):

Math.atan(Math.exp(-y * Math.PI / 180)) * 360 / Math.PI - 90;

В моем примере ниже я только что использовал d3.geoMercator(), чтобы сделать более четкий коэффициент масштабирования, используя проекцию, включающую одну дополнительную операцию для преобразования координаты x.

В противном случае 4-х шаговый процесс остается прежним.

Поиск правильных плиток

Я видел только один чистый подход к поиску плиток для показа, d3.quadTile Джейсона Дэвиса, видел здесь . Я полагаю, что Алан МакКончи использует незавершенную версию, которая может быть изменена иным способом . Существует также этот репозиторий github для другой версии d3.quadTiles, которая очень похожа.

Для McConchie / Davies d3.quadTile, с учетом проекции с экстентом клипа (не углом скрепления) и глубиной тайла, вытянет все тайлы, которые пересекают экстент вида.

В решении / примере Алана Маккончи уровень масштабирования зависит от масштаба проекции, но это не обязательно самый мудрый: у каждой проекции есть разные коэффициенты масштабирования, масштаб 100 в одной шкале будет показать другую степень, чем шкала 100 на другом. Кроме того, взаимосвязь между значениями масштаба и размером карты в цилиндрической проекции может быть линейной, в то время как нецилиндрические проекции могут иметь нелинейные зависимости между размером карты и масштабом.

Я немного изменил этот подход - я использую масштабный коэффициент, чтобы определить начальную глубину плитки, а затем уменьшаю эту глубину плитки, если число плиток, возвращаемое d3.quadTile, превышает определенное число:

geoTile.tileDepth = function(z) {
    // rough starting value, needs improvement:
    var a = [w/2-1,h/2]; // points in pixels
    var b = [w/2+1,h/2];
    var dx = d3.geoDistance(p.invert(a), p.invert(b)) ; // distance between in radians      
    var scale = 2/dx*tk;
    var z = Math.max(Math.log(scale) / Math.LN2 - 8, 2);
    z = Math.min(z,15) | 0;

    // Refine:
    var maxTiles = w*h/256/128;
    var e = p.clipExtent();
    p.clipExtent([[0,0],[w,h]])
    while(d3.quadTiles(p, z).length > maxTiles) {
        z--;
    }
    p.clipExtent(e);

    return z;
}

Затем, используя d3.quadTile, я вытаскиваю соответствующие плитки:

geoTile.tiles = function() {
    // Use Jason Davies' quad tree method to find out what tiles intercept the viewport:
    var z = geoTile.tileDepth();
    var e = p.clipExtent(); // store and put back after.

    p.clipExtent([[-1,-1],[w+1,h+1]]) // screen + 1 pixel margin on outside.
    var set = d3.quadTiles(p, Math.max(z0,Math.min(z,z1))); // Get array detailing tiles
    p.clipExtent(e);

    return set;
}

Сначала я подумал, что вытягивание плиток с разной глубины масштабирования (чтобы учесть расхождения в размерах перепроектированных плиток) было бы идеальным: но это столкнется с проблемами, такими как толщина линий в растре, а также прерывистыми аннотаций.

Принятие работы Джейсона и Алана

Я беру набор плиток, сгенерированный выше, с помощью geoTile.tiles() и запускаю его через цикл ввода / обновления / выхода, используя координату плитки (в координатах плитки, строке, столбце, глубине масштабирования) в качестве ключа, добавляя элементы image родителю g или svg. При загрузке изображений, как только изображение загружается, мы вызываем функцию загрузки, чтобы выполнить реальное воспроизведение. Это почти не изменилось по сравнению с Джейсоном и Аланом, я рассмотрел следующие проблемы, которые я видел в этом коде:

  • Ресэмплинг не учитывал веб-меркатор (упомянутый выше)
  • Глубина плитки выбрана не очень хорошо (упомянуто выше)
  • Плитки были перепроектированы как холсты, помещенные в div, а не в SVG - создав два родительских контейнера, по одному для каждого типа объектов: плитки или вектора.

Я считаю, что мой пример с очень незначительными изменениями обратился к ним. Я также добавил более подробные комментарии для просмотра:

    function onload(d, that) { // d is datum, that is image element.

        // Create and fill a canvas to work with.
        var mercatorCanvas = d3.create("canvas")
          .attr("width",tileWidth)
          .attr("height",tileHeight);
        var mercatorContext = mercatorCanvas.node().getContext("2d");           
        mercatorContext.drawImage(d.image, 0, 0, tileWidth, tileHeight); // move the source tile to a canvas.

        //
        var k = d.key; // the tile address.
        var tilesAcross = 1 << k[2]; // how many tiles is the map across at a given tile's zoom depth?

        // Reference projection:
        var webMercator = d3.geoMercator()
          .scale(tilesAcross/Math.PI/2) // reference projection fill square tilesAcross units wide/high.
          .translate([0,0])
          .center([0,0])

        // Reprojected tile boundaries in pixels.           
        var reprojectedTileBounds = path.bounds(d),
        x0 = reprojectedTileBounds[0][0] | 0,
        y0 = reprojectedTileBounds[0][1] | 0,
        x1 = (reprojectedTileBounds[1][0] + 1) | 0,
        y1 = (reprojectedTileBounds[1][1] + 1) | 0;

        // Get the tile bounds:
        // Tile bounds in latitude/longitude:
        var λ0 = k[0] / tilesAcross * 360 - 180,                     // left        
        λ1 = (k[0] + 1) / tilesAcross * 360 - 180,                   // right
        φ1 = webMercator.invert([0,(k[1] - tilesAcross/2) ])[1],     // top
        φ0 = webMercator.invert([0,(k[1] + 1 - tilesAcross/2) ])[1]; // bottom.             

        // Create a new canvas to hold the what will become the reprojected tile.
        var newCanvas = d3.create("canvas").node();

        newCanvas.width = x1 - x0,      // pixel width of reprojected tile.
        newCanvas.height = y1 - y0;     // pixel height of reprojected tile.
        var newContext = newCanvas.getContext("2d");    

        if (newCanvas.width && newCanvas.height) {
            var sourceData = mercatorContext.getImageData(0, 0, tileWidth, tileHeight).data,
                target = newContext.createImageData(newCanvas.width, newCanvas.height),
                targetData = target.data;

            // For every pixel in the reprojected tile's bounding box:
            for (var y = y0, i = -1; y < y1; ++y) {
              for (var x = x0; x < x1; ++x) {
                // Invert a pixel in the new tile to find out it's lat long
                var pt = p.invert([x, y]), λ = pt[0], φ = pt[1];

                // Make sure it falls in the bounds:
                if (λ > λ1 || λ < λ0 || φ > φ1 || φ < φ0) { i += 4; targetData[i] = 0; continue; }  
                    // Find out what pixel in the source tile matches the destination tile:
                    var top = (((tilesAcross + webMercator([0,φ])[1]) * tileHeight | 0) % 256  | 0) * tileWidth;
                    var q = (((λ - λ0) / (λ1 - λ0) * tileWidth | 0) + (top)) * 4;

                    // Take the data from a pixel in the source tile and assign it to a pixel in the new tile.
                    targetData[++i] = sourceData[q];
                    targetData[++i] = sourceData[++q];
                    targetData[++i] = sourceData[++q];
                    targetData[++i] = 255;
              }
            }
            // Draw the image.
            if(target) newContext.putImageData(target, 0, 0);
        }

        // Add the data to the image in the SVG:
        d3.select(that)
          .attr("xlink:href", newCanvas.toDataURL()) // convert to a dataURL so that we can embed within the SVG.
          .attr("x", x0)
          .attr("width", newCanvas.width)
          .attr("height",newCanvas.height)
          .attr("y", y0);
    }

Помещение в большую структуру.

Обычная карта тайлов с наложенными элементами имеет несколько систем координат:

  • Единицы плитки (3D), которые отмечают столбец, строку и уровень масштабирования каждой плитки (x, y, z соответственно)
  • Географические координаты (3D), которые обозначают широту и долготу точки на трехмерной сфере.
  • Единицы увеличения (3D), которые отслеживают трансфокатор (x, y) и масштаб увеличения (k).
  • проекционные единицы (2D), пиксельные единицы, на которые проецируются широта и долгота.

Цель любой скользкой карты - объединить эти координаты в удобную для использования систему.

Когда мы перепроектируем плитки, нам нужно добавить координатное пространство:

  • (/ a) проекция набора плиток.

Я чувствовал, что примеры не были особенно ясны в том, как они связали вместе все системы координат. Итак, я поместил вышеупомянутые методы, как вы могли видеть, в объект geoTile, который взят из личного проекта для библиотеки плиток . Целью этого является более плавная координация различных подразделений. Я не пытаюсь подключить его, он все еще находится в стадии разработки (слишком занят, чтобы действительно закончить его); тем не менее, я посмотрю, даст ли мне время возможность привести пример с d3-плиткой .

Проблемы, продвигающиеся вперед

Скорость увеличения и скорость реагирования - это самая большая проблема, которую я вижу. Чтобы решить эту проблему, я установил функцию масштабирования для запуска в конце масштабирования - это наиболее заметно при событиях панорамирования, так как обычно панорамирование непрерывно запускает функцию масштабирования через панорамирование, это можно решить путем перевода существующих изображений. Самый надежный способ использовать это, однако, на статической карте. Реализация перевода для уже нарисованных изображений была бы идеальной для событий панорамирования, а не для повторной выборки, как в настоящее время.

Анимация такой карты, вероятно, невозможна.

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

Примеры * * тысяча двести двадцать-одна К сожалению, кода слишком много для фрагмента, поэтому я сделал несколько ошибок для демонстрации. Выбор набора плиток Базовый зум / панорамирование с векторным слоем Ограниченная кастрюля с размером Fit Статическая карта У них было только минимальное тестирование. Если мне удастся закончить основную библиотеку плиток, я для этого раскошелюсь, а пока достаточно в качестве примера. Суть кода находится в geoTile.tile() в файле d3-reprojectSlippy.js, который содержит цикл ввода / обновления / выхода (довольно простой) и функцию загрузки, описанную выше. Поскольку я немного работаю над плитками, я буду постоянно обновлять этот ответ. Альтернатива

Перепроектирование плиток является громоздким и отнимает много времени. Альтернативой, если это возможно, будет создание набора плиток в желаемой проекции. Это было сделано с плитками OSM , но это также громоздко и отнимает много времени - только для создателя карты, а не для браузера.

TL; DR

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

...