Display Layers provided in Mercator projection

Jakub Balhar
5 min readApr 7, 2018
TimeScan over Prague in the EPSG:3857

Recently we needed to display layers provided only in the Mercator projection but using the standard WMS parameter schema. The final URL looks like this:

http://utep.it4i.cz/geoserver/gwc/service/wms?SERVICE=WMS&&request=GetMap&version=1.1.1&transparent=TRUE&layers=ESA_UTEP:GlobalTimeScan2015_30_RGB_01_11_08&styles=&format=image/png&width=256&height=256&srs=EPSG:3857&bbox=10018754.171394622,-1787187.9289878213,20037508.342789244,299486.13645630283

The prerequisities for this example is the proj4 and the Web World Wind in the version 0.9.0. In the example I am using bundled version of proj4 exposing proj4 into the global space. It isn’t best practice, but it will get the job done. You can include the proj4js from cdnjs: https://cdnjs.cloudflare.com/ajax/libs/proj4js/2.4.4/proj4.js

The Web World Wind already contains support for the layers projected in the EPSG:3857 projection. The relevant class is MercatorTiledImageLayer. This Layer handles the reprojection of the received imagery and correct size of the sector to ask for in the WGS84 space. To support layers with the schema shown above, the only missing thing is to provide mapSizeForLevel() function in the extended Layer and the URL Builder that reprojects the sector coordinates to the EPSG:3857 ones.

If you want to understand more about how this is done. There are two key parts that differs from standard layers. The createTile function which transforms the Sector to the size which represents the EPSG:3857 tiles in the WGS84 projection.

// Overridden from TiledImageLayer. Computes a tile's sector and creates the tile.
// Unlike typical tiles, Tiles at the same level do not have the same sector size.
MercatorTiledImageLayer.prototype.createTile = function (sector, level, row, column) {
var mapSize = this.mapSizeForLevel(level.levelNumber),
swX = WWMath.clamp(column * this.imageSize, 0, mapSize),
neY = WWMath.clamp(row * this.imageSize, 0, mapSize),
neX = WWMath.clamp(swX + (this.imageSize), 0, mapSize),
swY = WWMath.clamp(neY + (this.imageSize), 0, mapSize),
x, y, swLat, swLon, neLat, neLon;

x = (swX / mapSize) - 0.5;
y = 0.5 - (swY / mapSize);
swLat = 90 - 360 * Math.atan(Math.exp(-y * 2 * Math.PI)) / Math.PI;
swLon = 360 * x;

x = (neX / mapSize) - 0.5;
y = 0.5 - (neY / mapSize);
neLat = 90 - 360 * Math.atan(Math.exp(-y * 2 * Math.PI)) / Math.PI;
neLon = 360 * x;

sector = new Sector(swLat, neLat, swLon, neLon);

return TiledImageLayer.prototype.createTile.call(this, sector, level, row, column);
};

The second part is reprojecting the received imagery from the EPSG:3857 to WGS84

// Unproject the retrieved image.
for (var n = 0; n < 1; n++) {
for (var y = 0; y < image.height; y++) {
sy = 1 - y / (image.height - 1);
lat = sy * sector.deltaLatitude() + sector.minLatitude;
g = WWMath.gudermannianInverse(lat);
dy = 1 - (g - tMin) / (tMax - tMin);
dy = WWMath.clamp(dy, 0, 1);
srcRow = Math.floor(dy * (image.height - 1));
for (var x = 0; x < image.width; x++) {
kSrc = 4 * (x + srcRow * image.width);
kDest = 4 * (x + y * image.width);

destImageData.data[kDest] = srcImageData.data[kSrc];
destImageData.data[kDest + 1] = srcImageData.data[kSrc + 1];
destImageData.data[kDest + 2] = srcImageData.data[kSrc + 2];
destImageData.data[kDest + 3] = srcImageData.data[kSrc + 3];
}
}
}

MercatorUrlBuilder

Most of the code in the URL Builder is the same as in WmsUrlBuilder. It creates the URL containing standard WMS parameters such as Service, Method and all the others that are supported. If you need custom params update the URL Builder.

EOX Sentinel-2 Layer over Prague

The code relevant to this use case with Mercator is the transform function. It receives the coordinates of the sector to display in the WGS84 projection and using proj4 libray reprojects them to the EPSG:3857 projection. What follows is full implementation of a builder that can be directly used.

import WorldWind from '@nasaworldwind/worldwind';/** 
* @param serviceAddress {String} The full URL of the service providing the imagery. An example could be http://utep.it4i.cz/geoserver/gwc/service/wms
* @param layerNames {String[]} Array of the names of the layers to ask for. Example could be ["ESA_UTEP:GlobalTimeScan2015_30_RGB_01_11_08"]
* @param styleNames {String[]} Array of the names of the styles used for the layers. The length of layerNames should be equal to the length of styleNames
* @param wmsVersion {String} Version of the WMS service to use for retrieval. Examples 1.1.0, 1.1.1, 1.3.0
*/
class
MercatorUrlBuilder extends WorldWind.WmsUrlBuilder {
constructor(serviceAddress, layerNames, styleNames, wmsVersion) {
super(serviceAddress, layerNames.join(','), styleNames.join(','), wmsVersion);
}
urlForTile(tile, imageFormat) {
let sector = tile.sector;
let sb = WorldWind.WmsUrlBuilder.fixGetMapString(this.serviceAddress);
if (sb.search(/service=wms/i) < 0) {
sb = sb + "service=WMS";
}
sb = sb + "&request=GetMap";
sb = sb + "&version=" + this.wmsVersion;
sb = sb + "&transparent=" + (this.transparent ? "TRUE" : "FALSE");
sb = sb + "&layers=" + this.layerNames;
sb = sb + "&styles=" + this.styleNames;
sb = sb + "&format=" + imageFormat;
sb = sb + "&width=" + tile.tileWidth;
sb = sb + "&height=" + tile.tileHeight;

if(this.isWms130OrGreater) {
sb = sb + "&crs=EPSG:3857";
} else {
sb = sb + "&srs=EPSG:3857";
}
sb = this.transform(sb, sector);
sb = sb.replace(" ", "%20");

return sb;
}
transform(sb, sector) {
let source = new proj4.Proj('EPSG:4326');
let dest = new proj4.Proj('EPSG:3857');

let southWestOld = new proj4.Point( sector.minLongitude, sector.minLatitude );
let northEastOld = new proj4.Point( sector.maxLongitude, sector.maxLatitude );

let southWestNew = proj4.transform(source, dest, southWestOld);
let northEastNew = proj4.transform(source, dest, northEastOld);

sb = sb + "&bbox=";
sb = sb + southWestNew.x + "," + southWestNew.y + ",";
sb = sb + northEastNew.x+ "," + northEastNew.y;

return sb;
}
}

export default MercatorUrlBuilder;

MercatorLayer

This layer is more of a convenience class which bundle together the implementation of mapSizeForLevel and creation of the MercatorUrlBuilder. It is a working implementation in the ES6 though, so you can directly use it in your project.

import WorldWind from '@nasaworldwind/worldwind'; 
import MercatorUrlBuilder from './MercatorUrlBuilder';

/**
* @param displayName {String} Name of this layer. It is mostly an internal property and as such could be anything.
* @param options {Object}
* @param options.imageSize {Number} Dimension of the tile. This layer expects the tiles to be square and as such expects only size and not width and height. For the layer I display the correct size is 256. Probably you will need 256 or 512, other are quite rare.
* @param options.service {String} The full URL of the service providing the imagery. An example could be http://utep.it4i.cz/geoserver/gwc/service/wms
* @param options.layerNames {String[]} Array of the names of the layers to ask for. Example could be ["ESA_UTEP:GlobalTimeScan2015_30_RGB_01_11_08"]
* @param options.styleNames {String[]} Array of the names of the styles used for the layers. The length of layerNames should be equal to the length of styleNames
* @param options.version {String} Version of the WMS service to use for retrieval. Examples 1.1.0, 1.1.1, 1.3.0
*/
class MercatorLayer extends WorldWind.MercatorTiledImageLayer {
constructor(displayName, options) {
super(new WorldWind.Sector(-85.05, 85.05, -180, 180), new WorldWind.Location(85.05, 180), 19, "image/png", displayName, options.imageSize, options.imageSize);
this.urlBuilder = new MercatorUrlBuilder(options.service, options.layerNames, options.styleNames, options.version);

// This is used inside of the mechanism to generate correctly visible sector.
this.imageSize = options.imageSize;

// Create a canvas we can use when unprojecting retrieved images.
this.destCanvas = document.createElement("canvas");
this.destContext = this.destCanvas.getContext("2d");
}
mapSizeForLevel(levelNumber){
return 256 << (levelNumber + 1);
}
}
export default MercatorLayer;

This way you should be able to use such layers with even limited understanding of the underlying projections and math, if that is not what you are interested in.

--

--

Jakub Balhar

I always try to find answers to the complex questions. Now exploring the world of Mainframes in the Broadcom inc.