Add a rotated satellite image on the map using OpenLayers 5

Question

I'm trying to add a satellite image on my map using OpenLayers 5.

The problem is that I'm not able to do this, because I've just found an option to add an image on the map passing the image extent (xmin, ymin, xmax, ymax) and not the bounding box. The image should fit inside the bounding box. For that reason, the image was distorted.

The image is in JPG file (attribute feature.properties.icon). Example: http://exampleserver.com/220/063/353LGN00/353LGN00_thumb_large.jpg

The result that I would like is something like this:

The result that I've got was that:

My code that adds this image on the map is the following:

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
            ]  
        })
    })
)

Would someone know how to pass the image bounding box instead of just the image extent?

Thank you in advance.

EDIT 1: Mike's solution

Through Mike's solution I was able to fix a bug that some images have (near to the equator line). For that reason, his answer solved my problem and it inserted the image in a better position that I was expecting in the moment that I created the question.

However, this solution worked to me with images near to the equator line. Images next to the poles stay distorted (Edit 2).

I send below a picture illustrating the final result:

EDIT 2: New problem?

I was testing some images and I have discovered a new bug. Now I have discovered that the image should fit inside the bounding box. If the image does not fit inside the bbox, it stays distorted, such as the print that I send below illustrating.

The image should fit inside the bbox like in the image below [PS 1]:

I believe that it can be a problem of reprojection, but I don't know, because both view projection and image projection is EPSG:4326.

I tried to follow the explanation about Raster Reprojection[1.] on Openlayers site, however I was not able to reproduce it, because, as I said, both projections (view and image) are the same (or they should be).

I send below the GeoJSON that contains the information related to the image above. The image can be found in "properties.icon" (http://www.dpi.inpe.br/newcatalog/tmp/MOD13Q1/2018/MOD13Q1.A2018017.h13v14.jpg). The bbox coordinates can be found in "geometry.coordinates" or in "properties.bl_latitude", "properties.bl_longitude", "properties.br_latitude" and so on. "bl" means "bottom left", "br" means "bottom right", "tl" means "top left" and "tr" means "top right". These coordinates inside "properties" are the same inside "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"
}

Would someone have a new idea?

[PS 1]: The original code that does the image fits inside the bbox is a Leaflet code [2.] and I send it below:

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

Do you have lon/lat values for tl and br or a rotation value? — Mike, Mar 26 '19 at 20:35
Yes, I do lon/lat values for the four points of the image bounding box. They are "bl_latitude, bl_longitude, br_latitude, br_longitude, tl_latitude, tl_longitude, tr_latitude, tr_longitude". A rotation value I do not have. Should I calculate it by the bbox values? — rmmariano, Mar 27 '19 at 11:38
"Ask a question", on OpenLayers site, redirects to StackOverflow instead of GIS SE. — rmmariano, Mar 29 '19 at 17:25

Mike · Answer 1 · 2019-03-29T16:31:21.500

If the coordinates are those of the photo and the jpg which contains the rotated photo is in EPSG:4326 (i.e. aligned to meridians and parallels) then you need a bounding box containing all of the corners of the photo

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]
            ])  
        })
    })
)

However your top screenshot has the jpg itself rotated. If that is desired the projection isn't EPSG:4326 and you would need to define a custom projection to handle the rotation.

I've managed to get something close, but simply stretching the image to fit the polygon doesn't give the exact alignment at the side that the leaflet method does

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>

This is my KML method, but that is a simple rotation of a rectangle by a specified angle, not warping it into a quadrilateral where only two of the sides are parallel.

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>

Your solution worked to my images that are next to equator lines, but the ones that are next to the poles stay distorted, unfortunately. Now I have discovered that the images should fit inside the bounding box passed (bl_longitude, bl_latitude, ...). Would you have a new tip that you would like to share? Thank you so much for your reply. — rmmariano, Mar 27 '19 at 18:30
Can you post the original image for that part of South America and all the coordinates that go with it. Clearly it's going to need a rotated projection of some kind. I've managed those for KML ground overlays designed for Google Earth https://renenyffenegger.ch/notes/tools/Google-Earth/kml/index to work with OpenLayers https://i.stack.imgur.com/68bYZ.jpg so something similar may be possible with yours. — Mike, Mar 27 '19 at 20:45
I updated the question including what you would like. Really, both cases seem similar. How did you do to rotate the PNG image? In the question I included a Leaflet code that plots that South America image on the map. Thank you so much. — rmmariano, Mar 28 '19 at 12:32
I've adding working snippets for an attempt at reproducing the leaflet result and for the KML overlay.. — Mike, Mar 29 '19 at 16:34
I found this old question. https://stackoverflow.com/questions/52904587/displaying-georeferenced-images-using-openlayers-5 Maybe the affline.js library can give a better transform — Mike, Apr 05 '19 at 09:38

Add a rotated satellite image on the map using OpenLayers 5

1 Answers1

Linked