Добавьте повернутое спутниковое изображение на карту, используя OpenLayers 5

Question

Я пытаюсь добавить спутниковое изображение на карту, используя OpenLayers 5.

Проблема в том, что я не могу этого сделать, потому что я только что нашел возможность добавить изображение на карту, передавая экстент изображения (xmin, ymin, xmax, ymax), а не ограничивающий прямоугольник. Изображение должно помещаться внутри ограничительной рамки. По этой причине изображение было искажено.

Изображение находится в файле JPG (атрибут feature.properties.icon). Пример: http://exampleserver.com/220/063/353LGN00/353LGN00_thumb_large.jpg

Результат, который я хотел бы получить, выглядит примерно так:

В результате я получил:

Мой код, который добавляет это изображение на карту, следующий:

import ImageLayer from 'ol/layer/Image'
import Static from 'ol/source/ImageStatic'

...

this.olmap = new Map({
    target: 'map',
    layers: [
      baseLayerGroup, rasterLayerGroup, vectorLayer
    ],
    view: new View({
        projection: 'EPSG:4326',
        center: [ -45.8392, -3.65286 ],
        zoom: 8
    })
})

...

this.rasterLayerGroup.getLayers().push(
    new ImageLayer({
        source: new Static({
            url: feature.properties.icon,
            projection: 'EPSG:4326',
            imageExtent: [
              feature.properties.bl_longitude, feature.properties.bl_latitude,
              feature.properties.tr_longitude, feature.properties.tr_latitude
            ]  
        })
    })
)

Кто-нибудь знает, как передать ограничивающий прямоугольник изображения, а не только экстент изображения?

Заранее спасибо.

РЕДАКТИРОВАТЬ 1: решение Майка

Благодаря решению Майка мне удалось исправить ошибку, которая есть у некоторых изображений (рядом с линией экватора). По этой причине его ответ решил мою проблему, и он вставил изображение в лучшее положение, которое я ожидал в момент создания вопроса.

Однако, это решение работало со мной с изображениями рядом с линией экватора. Изображения рядом с полюсами остаются искаженными ( Изменить 2 ).

Я отправляю ниже картинку, иллюстрирующую конечный результат:

РЕДАКТИРОВАТЬ 2: Новая проблема?

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

Изображение должно помещаться внутри коробки, как показано на рисунке ниже [PS 1] :

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

Я пытался следовать объяснению Raster Reprojection [1.] На сайте Openlayers, однако мне не удалось воспроизвести его, потому что, как я сказал, обе проекции (вид и изображение) одинаковы (или они должны быть ).

Я отправляю ниже GeoJSON, который содержит информацию, относящуюся к изображению выше. Изображение можно найти в «properties.icon» (http://www.dpi.inpe.br/newcatalog/tmp/MOD13Q1/2018/MOD13Q1.A2018017.h13v14.jpg).. Координаты bbox можно найти в «geometry.coordinates» или в «properties.bl_latitude», «properties.bl_longitude», «properties.br_latitude» и скоро. «bl» означает «нижний левый», «br» означает «нижний правый», «tl» означает «верхний левый», а «tr» означает «верхний правый». Эти координаты внутри «свойств» одинаковы внутри «geometry.coordinates».

{
  "geometry": {
    "coordinates": [
      [
        [
          -77.7862,
          -50
        ],
        [
          -100,
          -60
        ],
        [
          -80,
          -60
        ],
        [
          -62.229,
          -50
        ],
        [
          -77.7862,
          -50
        ]
      ]
    ],
    "type": "Polygon"
  },
  "properties": {
    "alternate": "http://www.dpi.inpe.br/opensearch/v2/granule.json?uid=MOD13Q1.A2018017.h13v14",
    "auxpath": null,
    "bitslips": null,
    "bl_latitude": -60,
    "bl_longitude": -100,
    "br_latitude": -60,
    "br_longitude": -80,
    "centerlatitude": -55,
    "centerlongitude": -80.0038,
    "centertime": null,
    "cloud": 0,
    "cloudcovermethod": "M",
    "dataset": "MOD13Q1",
    "date": "2018-01-17T00:00:00",
    "enclosure": [
      {
        "band": "evi",
        "radiometric_processing": "SR",
        "type": "MOSAIC",
        "url": "http://www.dpi.inpe.br/newcatalog/tmp/MOD13Q1/2018/MOD13Q1.A2018017.h13v14.006.2018033223827.hdf"
      },
      {
        "band": "ndvi",
        "radiometric_processing": "SR",
        "type": "MOSAIC",
        "url": "http://www.dpi.inpe.br/newcatalog/tmp/MOD13Q1/2018/MOD13Q1.A2018017.h13v14.006.2018033223827.hdf"
      },
      ...
    ],
    "icon": "http://www.dpi.inpe.br/newcatalog/tmp/MOD13Q1/2018/MOD13Q1.A2018017.h13v14.jpg",
    "id": "http://www.dpi.inpe.br/opensearch/v2/granule.json?uid=MOD13Q1.A2018017.h13v14",
    "orbit": 0,
    "path": 14,
    "provider": "OP_CBERS1",
    "row": 13,
    "satellite": "T1",
    "sensor": "MODIS",
    "title": "MOD13Q1.A2018017.h13v14",
    "tl_latitude": -50,
    "tl_longitude": -77.7862,
    "tr_latitude": -50,
    "tr_longitude": -62.229,
    "type": "IMAGES",
    "updated": "2018-03-01T18:51:56",
    "via": "http://www.dpi.inpe.br/opensearch/v2/metadata/MOD13Q1.A2018017.h13v14"
  },
  "type": "Feature"
}

Будет ли у кого-нибудь новая идея?

[PS 1]: Исходный код, который соответствует изображению, помещается в bbox, это код Leaflet [2.], и я отправляю его ниже:

var map = L.map('map').setView([-15.22, -53.23], 5)

...

var anchor = [
    [feature.properties.tl_latitude, feature.properties.tl_longitude],
    [feature.properties.tr_latitude, feature.properties.tr_longitude],
    [feature.properties.br_latitude, feature.properties.br_longitude],
    [feature.properties.bl_latitude, feature.properties.bl_longitude]
]

layer._quicklook = L.imageTransform(feature.properties.icon, anchor).addTo(map)

[1.] https://openlayers.org/en/latest/doc/tutorials/raster-reprojection.html

[2.] https://github.com/ScanEx/Leaflet.imageTransform

Answer 1

Если координаты совпадают с координатами фотографии, а jpg, который содержит повернутую фотографию, находится в EPSG: 4326 (т.е. выровнен по меридианам и параллелям), тогда вам нужна ограничительная рамка, содержащая все углы фотографии

import {boundingExtent} from 'ol/extent';

....

this.rasterLayerGroup.getLayers().push(
    new ImageLayer({
        source: new Static({
            url: feature.properties.icon,
            projection: 'EPSG:4326',
            imageExtent: boundingExtent([
              [feature.properties.bl_longitude, feature.properties.bl_latitude],
              [feature.properties.br_longitude, feature.properties.br_latitude],
              [feature.properties.tl_longitude, feature.properties.tl_latitude],
              [feature.properties.tr_longitude, feature.properties.tr_latitude]
            ])  
        })
    })
)

Однако ваш верхний скриншот имеет сам поворот jpg. Если это необходимо, проекция не EPSG: 4326, и вам необходимо определить собственную проекцию для обработки поворота.

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

var properties = {
    "bl_latitude": -60,
    "bl_longitude": -100,
    "br_latitude": -60,
    "br_longitude": -80,
    "centerlatitude": -55,
    "centerlongitude": -80.0038,
    "icon": "https://www.mikenunn.net/demo/MOD13Q1.A2018017.h13v14.jpg",
    "tl_latitude": -50,
    "tl_longitude": -77.7862,
    "tr_latitude": -50,
    "tr_longitude": -62.229,
};

function overlaySource ( properties ) {

    var projection = ol.proj.get('EPSG:3857');  // leaflet projection

    var extentSize = [0, 0, 4096, 4096];  // arbitary extent for the projection transforms
    var size0 = extentSize[2];
    var size1 = extentSize[3];

    var url = properties.icon;

    var bl = ol.proj.transform([properties.bl_longitude, properties.bl_latitude], 'EPSG:4326', projection);
    var tl = ol.proj.transform([properties.tl_longitude, properties.tl_latitude], 'EPSG:4326', projection);
    var br = ol.proj.transform([properties.br_longitude, properties.br_latitude], 'EPSG:4326', projection);
    var tr = ol.proj.transform([properties.tr_longitude, properties.tr_latitude], 'EPSG:4326', projection);

    function normalTransform(coordinates, output, dimensions) {
        var dims = dimensions || 2;
        for (var i=0; i<coordinates.length; i+=dims) {

            var left = bl[0] + (tl[0]-bl[0]) * coordinates[i+1]/size1;
            var right = br[0] + (tr[0]-br[0]) * coordinates[i+1]/size1;

            var top = tl[1] + (tr[1]-tl[1]) * coordinates[i]/size0;
            var bottom = bl[1] + (br[1]-bl[1]) * coordinates[i]/size0;

            var newCoordinates0 = left + (right-left) * coordinates[i]/size0;
            var newCoordinates1 = bottom + (top-bottom) * coordinates[i+1]/size1;

            c = ol.proj.transform([newCoordinates0, newCoordinates1], projection, 'EPSG:3857');

            //console.log(coordinates[i] + ' ' + coordinates[i+1] + ' ' + c[0] + ' ' + c[1]);

            coordinates[i] = c[0];
            coordinates[i+1] = c[1];

        }
        return coordinates;
    }

    function rotateTransform(coordinates, output, dimensions) {
        var dims = dimensions || 2;
        for (var i=0; i<coordinates.length; i+=dims) {

            c = ol.proj.transform([coordinates[i], coordinates[i+1]], 'EPSG:3857', projection);

            var left = bl[0] + (tl[0]-bl[0]) * (c[1]-bl[1]) /(tl[1]-bl[1]);
            var right = br[0] + (tr[0]-br[0]) * (c[1]-br[1]) /(tr[1]-br[1]);

            var top = tl[1] + (tr[1]-tl[1]) * (c[0]-tl[0])/(tr[0]-tl[0]);
            var bottom = bl[1] + (br[1]-bl[1]) * (c[0]-bl[0])/(br[0]-bl[0]);

            var newCoordinates0 = (c[0]-left)*size0/(right-left);
            var newCoordinates1 = (c[1]-bottom)*size1/(top-bottom);

            //console.log(coordinates[i] + ' ' + coordinates[i+1] + ' ' + newCoordinates0 + ' ' + newCoordinates1);

            coordinates[i] = newCoordinates0;
            coordinates[i+1] = newCoordinates1;

        }
        return coordinates;
    }

    var rotatedProjection = new ol.proj.Projection({
        code: 'EPSG:' + url + 'rotated',
        units: 'm',
        extent: extentSize
    });
    ol.proj.addProjection(rotatedProjection);

    ol.proj.addCoordinateTransforms('EPSG:3857', rotatedProjection,
        function(coordinate) {
            return rotateTransform(coordinate);
        },
        function(coordinate) {
            return normalTransform(coordinate);
        }
    );

    ol.proj.addCoordinateTransforms('EPSG:4326', rotatedProjection,
        function(coordinate) {
            return rotateTransform(ol.proj.transform(coordinate, "EPSG:4326", "EPSG:3857"));
        },
        function(coordinate) {
            return ol.proj.transform(normalTransform(coordinate), "EPSG:3857", "EPSG:4326");
        }
    );

    return new ol.source.ImageStatic({
        projection: rotatedProjection,
        imageExtent: extentSize,
        url: url
    });

}

var tileLayer = new ol.layer.Tile({
    source: new ol.source.XYZ({
        attributions: [
            'Powered by Esri',
            'Source: Esri, DigitalGlobe, GeoEye, Earthstar Geographics, CNES/Airbus DS, USDA, USGS, AeroGRID, IGN, and the GIS User Community'
        ],
        //attributionsCollapsible: false,
        url: 'https://services.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}',
        maxZoom: 23
    })
});

var imageLayer = new ol.layer.Image({
    source:  overlaySource( properties ),
    opacity: 0.7
})

var map = new ol.Map({
    layers: [tileLayer, imageLayer],
    target: 'map',
    logo: false,
    view: new ol.View()
});

var imageProj = imageLayer.getSource().getProjection();
map.getView().fit(ol.proj.transformExtent(imageProj.getExtent(), imageProj, map.getView().getProjection()), {constrainResolution: false});

html, body, .map {
    margin: 0;
    padding: 0;
    width: 100%;
    height: 100%;
}

<link href="https://cdn.rawgit.com/openlayers/openlayers.github.io/master/en/v5.3.0/css/ol.css" rel="stylesheet" />
<script src="https://cdn.rawgit.com/openlayers/openlayers.github.io/master/en/v5.3.0/build/ol.js"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/proj4js/2.5.0/proj4.js"></script>
<div id="map" class="map"></div>

Это мой метод KML, но это простое вращение прямоугольника на заданный угол, а не деформация его в четырехугольник, где только две из сторон параллельны.

function kmlOverlaySource ( kmlExtent, // KMLs specify the extent the unrotated image would occupy
                            url,
                            rotation,
                            imageSize,
) {

    // calculate latitude of true scale of equidistant cylindrical projection based on pixels per degree on each axis

    proj4.defs('EPSG:' + url, '+proj=eqc +lat_ts=' +
                              (Math.acos((ol.extent.getHeight(kmlExtent)/imageSize[1])
                                        /(ol.extent.getWidth(kmlExtent)/imageSize[0]))*180/Math.PI) +
                              ' +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs');

    if (ol.proj.proj4 && ol.proj.proj4.register) { ol.proj.proj4.register(proj4); } // if OL5 register proj4

    // convert the extents to source projection coordinates

    var projection = ol.proj.get('EPSG:' + url);
    var projExtent = ol.proj.transformExtent(kmlExtent, 'EPSG:4326', projection);

    var angle = -rotation * Math.PI/180;

    function rotateTransform(coordinates, output, dimensions) {
        var dims = dimensions || 2;
        for (var i=0; i<coordinates.length; i+=dims) {
            var point = new ol.geom.Point([coordinates[i],coordinates[i+1]]);
            point.rotate(angle, ol.extent.getCenter(projExtent));
            var newCoordinates = point.getCoordinates();
            coordinates[i] = newCoordinates[0];
            coordinates[i+1] = newCoordinates[1];
        }
        return coordinates;
    }

    function normalTransform(coordinates, output, dimensions) {
        var dims = dimensions || 2;
        for (var i=0; i<coordinates.length; i+=dims) {
            var point = new ol.geom.Point([coordinates[i],coordinates[i+1]]);
            point.rotate(-angle, ol.extent.getCenter(projExtent));
            var newCoordinates = point.getCoordinates();
            coordinates[i] = newCoordinates[0];
            coordinates[i+1] = newCoordinates[1];
        }
        return coordinates;
    }

    var rotatedProjection = new ol.proj.Projection({
        code: 'EPSG:' + url + 'rotated',
        units: 'm',
        extent: projExtent
    });
    ol.proj.addProjection(rotatedProjection);

    ol.proj.addCoordinateTransforms('EPSG:4326', rotatedProjection,
        function(coordinate) {
            return rotateTransform(ol.proj.transform(coordinate, 'EPSG:4326', projection));
        },
        function(coordinate) {
            return ol.proj.transform(normalTransform(coordinate), projection, 'EPSG:4326');
        }
    );

    ol.proj.addCoordinateTransforms('EPSG:3857', rotatedProjection,
        function(coordinate) {
            return rotateTransform(ol.proj.transform(coordinate, 'EPSG:3857', projection));
        },
        function(coordinate) {
            return ol.proj.transform(normalTransform(coordinate), projection, 'EPSG:3857');
        }
    );

    return new ol.source.ImageStatic({
        projection: rotatedProjection,
        url: url,
        imageExtent: projExtent
    });

}

var tileLayer = new ol.layer.Tile({
    source: new ol.source.XYZ({
        attributions: [
            'Powered by Esri',
            'Source: Esri, DigitalGlobe, GeoEye, Earthstar Geographics, CNES/Airbus DS, USDA, USGS, AeroGRID, IGN, and the GIS User Community'
        ],
        //attributionsCollapsible: false,
        url: 'https://services.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}',
        maxZoom: 23
    })
});

// these would normally be parsed from a KML file
var kmlExtent = [8.433995415151397, 46.65804355828784, 9.144871415151389, 46.77980155828784];
var url = 'https://raw.githubusercontent.com/ReneNyffenegger/about-GoogleEarth/master/kml/the_png_says.png'
var rotation = 30;
var imageSize = [];

var img = document.createElement('img');
img.onload = imageLoaded;
img.src = url;

function imageLoaded() {

    imageSize[0] = img.width;
    imageSize[1] = img.height;

    var imageLayer = new ol.layer.Image({
        source: kmlOverlaySource(kmlExtent, url, rotation, imageSize),
    });

    var map = new ol.Map({
        layers: [tileLayer, imageLayer],
        target: 'map',
        logo: false,
        view: new ol.View()
    });

    var imageProj = imageLayer.getSource().getProjection();
    map.getView().fit(ol.proj.transformExtent(imageProj.getExtent(), imageProj, map.getView().getProjection()), {constrainResolution: false});

}

html, body, .map {
    margin: 0;
    padding: 0;
    width: 100%;
    height: 100%;
}

<link href="https://cdn.rawgit.com/openlayers/openlayers.github.io/master/en/v5.3.0/css/ol.css" rel="stylesheet" />
<script src="https://cdn.rawgit.com/openlayers/openlayers.github.io/master/en/v5.3.0/build/ol.js"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/proj4js/2.5.0/proj4.js"></script>
<div id="map" class="map"></div>